Table Of Content

Fixation and consensus times on a network: a unified approach G. J. Baxter,1 R. A. Blythe,2 and A. J. McKane3 1School of Mathematics, Statistics and Computer Science, Victoria University of Wellington, PO Box 600, Wellington, New Zealand 2SUPA, School of Physics, University of Edinburgh, Mayfield Road, Edinburgh EH9 3JZ, UK 3Theory Group, School of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK 9 We investigate a set of stochastic models of biodiversity,population genetics, language evolution 0 andopiniondynamicsonanetworkwithinacommonframework. Eachnodehasastate,0<xi<1, 0 withinteractionsspecifiedbystrengthsmij. Foranysetofmij wederiveanapproximateexpression 2 forthemeantimetoreach fixationorconsensus(all xi=0or1). Remarkablyinacaserelevantto language change thistime is independent of thenetwork structure. n a PACSnumbers: 05.40.-a,89.75.Hc,87.23.-n,51.10.+y J 5 Mathematical models predicting biological and social ] (b) ech clsathasattnifsgetewicaayrleepabrhesychsoiamcvisiintnsggiinsneccerunelatausnirnaeglxldpyylocnosaimomnmicoosfn[ap1cl]ta.icvOei,tnyweaiatmshpotenhcget i (a(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1))(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) mij (cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) m - of this work which is not widely appreciated is that sev- (cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) j eral seemingly distinct phenomena can be described by t a very similar models: in some cases they can even be ex- t s actly mapped into each other [2]. Examples include bio- . FIG. 1: Stochastic copying dynamics between two islands i t diversity [3], population genetics [4], opinion dynamics a andj formingpartofalargernetwork. Afteradeathoni,an m [1] and language change [5]. The common thread is that objectisreplacedeitherwithacopyfrom(a)thesameisland - objects that come in different variants are copied from i or (b) a different island j (rate mij). d one place (an “island”) to another according to some n stochastic rule. If no new variants are created in the o process (e.g., by mutation) and the number of objects at random to die and to be replaced by a copy of either c [ does not grow without bound, one variant is eventually (a) another individual picked at random from island 2, guaranteedtotakeoveranentirepopulation,orgotofix- or (b) a migrant from island 1. Process (b) is assumed 2 ation in the genetics parlance. Changes in the network to happen less frequently than (a), and as mentioned v 3 structure connecting differentislands,and the stochastic above, if (b) is absent then the final state of the system 8 rules used to choose the source and target islands, lead isonewhichcontainsindividualsofonlyonespecies. The 0 toavarietyofscalinglawsrelatingthenumberofislands number of individuals on island 2 is a constant, n, and 3 andtimetoreachfixation(see,e.g.,[6,7,8]). InthisLet- new species are created on island 1, but not island 2, by . 1 ter, we present a theoretical treatment of a very general mutation-like events. This model is a neutral theory; no 0 stochasticcopyingmodel thatincludes many previously- one species is assumed to be “fitter” than another. 8 studiedcasesandunifies the diversefixationtime results The relationship between this model and simple neu- 0 obtained so far. We also discuss a mapping to a particle : tral models of population genetics is well known [3]: in- v reactionsystemfromwhichitcanbeshownthatourpre- dividuals are analogous to genes and species are types i dictionforthe fixationtime, obtainedbymakingvarious X of genes (alleles). A more general model comprises a set approximations,can in many cases be stated as a bound r of islands labeled i = 1,...,N each of which contains n a thatsimulationsshowisoftensaturated. Wealsodiscuss genes (in reality, n individuals each containing one copy a key application of our findings—to a current theory of of the gene of interest). For simplicity we assume that new-dialect formation [9]. there is no mutation, so that no new alleles may be cre- Toestablishthebasicfeaturesofthelargeclassofmod- ated. The alleles are labeled by α = 1,...,M. The els we consider, we describe a prominent special case, dynamical processes, illustrated in Fig. 1 are as before: Hubbell’s model of biodiversity and biogeography [3]. a death onislandi followedby (a) a birth onislandi, or Here there are only two islands: a metacommunity or (b)amigrantoffspringfromislandj arrivingonislandi. mainland (island 1) and a local community (island 2). There are various ways to parametrize these dynamics. Theobjectsareindividualswhichcompeteforacommon Asin[2],welettheparent(copiedobject)betakenfrom resource(e.g.treescompetingforspace,sunlightandnu- island j with probability f . In process (b) the proba- j trients [10]) and the variants are different species. At bility the offspring (copy) lands on island i is taken to regular time intervals an individual in island 2 is picked be proportional to m , which specifies a migration rate ij 2 within a standardcontinuous-time limit described in [2], 1 but whose details are unimportant here. At any time t, 0.8 α i the state ofthe systemcanbe givenin terms ofthe frac- 0.6 tion xiα(t) of genes on island i that are allele α. Then, 0.4 α σ the islands can be represented as nodes on a network, 0.2 j ii σ whichintheopiniondynamicsandlanguagechangemod- 0 ij T T elsareindividualsiwithopinionsorlanguagevariantsα 0 0 500 time (arbitrary units) 1000 expressed with frequency x . In these social contexts, iα fixation is sometimes called consensus. FIG.2: Numericalsolutionformeansαi,αj and(co)variances σij ≡βij−αiαj onafully-connectednetworkofN =20sites. Using the formalism first developed for these systems Note the αs converge at a time T0 much less than the mean in population genetics [4], we describe the evolution in fixation time T, and at time T0, σij ≈ 0 for i 6= j. These terms of a Fokker-Planck equation. Suppose, first of all, features become more pronounced as N is increased. that there is only one island and no mutation. Then the only dynamical process is random genetic drift in which the frequencies xα(t) diffuse on the interval [0,1] due to is conserved by the dynamics. Decomposing αi(t) in the random sampling in the death/birth process. For terms of its right eigenvectors we see that it approaches simplicity, we will assume that there are only two alleles a constant independent of i as t → ∞. Since all xi(t) which have frequencies x and (1 − x), so the state of tend to 0 or 1 as t → ∞, this is the probability of the the system is described by the single stochastic variable allele fixing. From (3) we see hξ(t)i also approaches this x∈[0,1]. The probabilitythatthe systemis inthe state value in this limit and since ξ is conserved, the fixation x at time t, P(x,t), satisfies the Fokker-Planckequation probability is ξ(0) [6, 11]. ∂ P(x,t)=∂2[D(x)P(x,t)] where the diffusion constant Numerical studies suggest that convergence of the en- t x is state dependent: D(x) = x(1−x)/2 [4]. Moving to semble average αi(t) to its asymptote ξ(0) occurs on a the case of N islands, with migration rate m from j to much shorter timescale than the ultimate fixation of a ij i, the Fokker-Planck equation for P(x ,t) now reads [2] variant, which in turn governs the rate of change of the i secondmoments β (t)=hx (t)x (t)i, see Fig. 2. Eq.(1) ij i j ∂P ∂ ∂ implies for the latter = m −m [(x −x )P] ij ji i j ∂t (cid:18) ∂x ∂x (cid:19) Xhiji i j dβ ij = m β + m β +δ f (α −β ) (4) 1 N ∂2 dt ik kj jℓ iℓ ij i i ii + f [x (1−x )P] , (1) Xk Xℓ 2 i∂x2 i i Xi=1 i whilst the mean time to fixation,T,is givenby the solu- where hiji means sum over distinct pairs i,j. This may tionofabackwardversionoftheFokker-Planckequation begeneralizedtoM >2andtoincludemutation[2],but (1) [12] Eq. (1) will be sufficient for our purposes. Analysis of Eq. (1) is, on the face of it, a hopeless ∂T ∂T −1 = − (x −x ) m −m i j ij ji task since it has many degrees of freedom, xi, interact- Xhiji (cid:18) ∂xi ∂xj(cid:19) ing with arbitrary strengths m . However, much of the ij macroscopic dynamics is captured by the first and sec- 1 N ∂2T + f [x (1−x )] . (5) ond moments of xi(t). The mean, αi(t)=hxi(t)i can be 2Xi=1 i i i ∂x2i found from Eq. (1) to evolve according to The assumption that the time overwhich all the α con- i N dαi verge, T0, is much less than T (see Fig. 2) leads to the = m (α −α )≡ m α , (2) ij j i ij j dt followingapproximatetreatmentofthisequation. Weas- Xj6=i Xj=1 sume that T depends only on the state of the system at where the equivalence holds if the diagonal elements time T0, and principally through ξ(0). Changing to the m are defined to be − m . This matrix has a ξ(0)variable using Eq. (3) we find thatthe firsttermon ii j6=i ij zero eigenvalue,which wePwill assume is non-degenerate. the right-hand side of Eq. (5) vanishes [6] giving The associated right eigenvector has all elements equal to one, and the left eigenvector we denote Qi, so that N 2 d2T PThieQnimweijfi=nd0f,roamndEnqo.rm(a2l)iztehdatsutchhetehnastemPbileQi(n=ois1e-. −2=Xi=1fiQixi(T0)[1−xi(T0)]dξ(0)2 . (6) history) average of the collective variable This equation still depends on the variables x at time i N T0. We canestimate these byassumingthatcorrelations betweenthenodesareabsent,i.e. β =α α =ξ(0)2i6= ξ(t)= Q x (t) (3) ij i j i i Xi=1 j, and that the rate of change of the variance of xi is 3 sufficientlyslowthatthetime derivativeinEq. (4)when predicts a fixationtime that, inour time units, increases i = j can be neglected (see again Fig. 2). Then β at quadratically with N. ii time T0 can be estimated from Eq. (4). Wealsorecoverknownresultsforlinkdynamics[7,13], By replacing x (1−x ) in Eq. (6) by α −β at time whereanopinioniscopiedinarandomly-chosendirection i i i ii T0, we find [ξ(0)(1−ξ(0))]d2T/dξ(0)2 =−2/r where alonga randomly-chosenlink in eachtimestep. Here, we have f = k /N and m =A /(2N) which implies that i i ij ij r ≈ Q2f 2 j6=imij . (7) Qi =1/N. Inthisspecialcaseithasbeennotedthatthe i i2 P m +f mean time to fixation does not depend on the network Xi j6=i ij i structure [7, 13]. Furthermore, our general result (7) P includes the various scaling forms found in [8]. The mean fixation time is obtained by integrating the The generality of our results canbe tested more strin- equation for T(ξ(0)), with the boundary conditions gently by choosing f and m from various random dis- T(0)=T(1)=0 to give i ij tributions in such a way that the numerators and de- 2 nominators in (7) are of similar magnitudes. Simulation T(ξ(0))=− [ξ(0)lnξ(0)+(1−ξ(0))ln(1−ξ(0))] , r results (not shown) on Erdo¨s-Rényi random graphs of (8) varying densities [14] are consistent with Eqs. (7) and whichintandemwith(7)isourmainresult,derivedfora (8), except when T and T0 turn out to be of a similar largeclassofstochastic-copyingprocessesonanynetwork order in N. and arbitrary migration rates. It should be emphasized Aconcreteapplicationofthemodelistoanevolution- that the reduction from a system with (M −1) degrees ary model of language change [5]. In fact it was our ex- of freedom, to a system with only one degree of freedom perienceinapplyingthis modeltothe emergenceofNew is not imposed, but instead emerges from the dynamics. Zealand English [15] that alerted us to the possibility of For orientationandasa checkofthese results,we ver- a general analysis and to some surprising behavior for a ifythatspecificchoicesoftheparametersf andm give subclassofmodels. Asindicatedearlier,inthisparticular i ij expressions which have been previously obtained. As a applicationthe index i labels speakerswho converseand firstcheck,weexaminethecaseofvoterdynamicsonhet- so affect eachothers grammar. The probability that two erogeneous graphs described in [6, 7]. The voter model speakersinteractisgivenbyamatrixG ;thisisareflec- ij proper has each node chosen uniformly to be the recipi- tion of the topology of the network. Another factor H ij ent of a new opinion, hence all fi = 1/N. This opinion accountsfortheweightthatigivestotheutterancesofj. is selected randomly from the node’s neighbors, giving Itis the product ofthe G andH thatis equalto m . ij ij ij mij = Aij/(2Nki), where Aij = 1 if nodes i,j are con- WhileGij isbydefinitionsymmetric,Hij neednotbe. In nected,ki isthedegreeofnodei,ki = j6=iAij,andthe theapplicationtotheformationofNewZealandEnglish, factorof2N isaconsequenceoftheconPtinuous-timelimit Trudgill[9]haspostulatedthatthelargequantityofdata described in [2]. To calculate r, one needs Q , the nor- that exists on the emergence on this English variant[16] i malizedzerolefteigenvectorofthematrixm . Onefinds may be explained by assuming that social factors which ij thatthiseigenvectorequalsQ =k /(Nk¯),wherek¯isthe are modeled by H are unimportant, and therefore this i i ij meandegreeandisincludedtonormalizetheQ . Putting factor may be replaced by a constant. If this is so, the i all this together we find from (7) that r = k2/2(Nk¯)2. m is symmetric; if the model is not socially-neutral, ij This corresponds with the results of [6, 7] after taking then m will not in general be symmetric. ij into account the different choice of time units made in To test Trudgill’s theory, we assume m is symmetric ij those works. and using the results (7) and (8) check if the mean time Aninvasionprocess,wherethesourcenodeofanopin- tofixationisinaccordwiththedata. Apriorionewould ion is randomly selected in each timestep, was also con- expect that the precise nature of the speaker network sidered in [7]. This corresponds to the choice m = during the formation of New Zealand English would be ij A /(2Nk ), from which Q may again be found. How- required. However, remarkably, it turns out that in this ij j i ever, the resulting expression for r via (7) does not sim- casethetimetofixationiscompletelyindependentofthe plify unless one additionally assumes (as in [6, 7]) that network structure! This is because the symmetry of m ij node degreesareuncorrelated,A =k k /(Nk¯). Taking implies that its left and right eigenvectors are propor- ij i j f = 1/N as before, one finds from (7) the result of [7] tional to one another, and so Q = 1/N. This together i i in the limit N → ∞. Our treatment however also ex- with m ∝ f [5] implies a constant r from Eq. (7), j ij i tends to networks with strong degree correlations, e.g., whichPis confirmed by simulation results shown in Fig. 3 a star network that has a central node connected to all (although see below for some caveats). To evaluate the N −1 outer nodes, which in turn are connected only to plausibilityofthissocially-neutralmodelasamechanism the central node. For this network, (7) predicts a fixa- for new dialect formation requires in addition to this re- tion time proportional to N3, confirmed by simulation sult detailed consideration of human memory lifetimes data (not shown). By contrast, the approximation of [7] that are beyond the scope of this article and discussed 4 300 T0 ≪ T corresponds to taking the longest timescale in theory the backward-time dynamics to be that associated with me 250 fruanlldyo cmo n(n2e0c%te dp onsestiwbloer kedges connected) thefinalcoalescencereaction;thatis,thatthesubsequent n ti 200 star network relaxationof the one particle state to its equilibrium oc- o ati 150 cursonamuchshortertimescale. Ifthelasttwoparticles x n fi 100 arewell-mixed,wewouldexpectthisassumptiontohold, a me since each particle would be close to the single-particle 50 equilibrium state. 0 0 10 20 30 40 50 60 70 80 90 100 110 Insummary,wehaveprovidedanapproximatetheory, N validated numerically on a range of networks, for cal- culating the fixation time within a general model that FIG.3: Fixationtimesperspeaker withinthesocially-neutral hasfeaturessharedbyphysical,biologicalandsocialsys- utteranceselectionmodel(definedinthetext)onvariousnet- works with N nodes (see legend). tems. Since in reality specific applications will contain additionalprocesses,itisofinteresttoextendthegeneral elsewhere [15]. Ultimately, we conclude that a purely approach we have described to cater for these. We have neutral model is hard to reconcile with available empiri- alsobrieflydiscussedhowaparticularlystrikingresult— cal data. thatfixationtimeisindependentofnetworkstructure—is We remark that our main result, Eq. (7), can also be ofdirectrelevancetocurrentlinguistictheory[9,15]and obtained and interpreted within a backwards-time for- hope that the trend whereby mathematical results for mulation. Consider for example two objects of the same formalmodels of socialbehaviorareappliedin empirical type; by asking which objects these were copied from contexts will continue, and increase, in the future. at some earlier time, one can reconstruct their ances- RAB holds an RCUK Academic Fellowship. AJM tral lineages which will hop between nodes as one goes wishes to thank the EPSRC (UK) for financial support back in time. Eventually, both objects will have been under grant GR/T11784. copiedfromthe same‘parent’,causingmergerofthe lineages. This processiscalledthe coalescent inpopulation genetics [17], but by viewing the lineages as particle tra- jectories, one can also recognize it as the A+A → A particle coalescence reaction of non-equilibrium statisti- [1] C. Castellano, S. Fortunato, and V. Loreto, Statisti- cal mechanics [18] on a network. cal physicsofsocial dynamics,Preprint,arXiv:0710.3256 (2007). In this picture, the quantityr inEq.(7) is the asymp- [2] R. A. Blythe and A. J. McKane, J. Stat. Mech. P07018 totic rate at which the last two unreacted particles co- (2007). alesce. The specific expression quoted above can be in [3] S.P.Hubbell,TheUnifiedNeutralTheoryofBiodiversity fact obtained as a bound by means of a variationalprin- and Biogeography (Princeton University Press, Prince- ciple. Most usefully, this approach—whose details will ton, 2001). be presented elsewhere [19]—yields physical insight into [4] J. F. Crow and M. Kimura, An Introduction to Pop- whentheapproximationswehavemadetoobtain(7)are ulation Genetics Theory (Harper and Row, New York, 1970). valid. Thekeyassumptionis thatthe lasttwounreacted [5] G. Baxter, R. A. Blythe, W. Croft and A. J. McKane, particles must eachtypically have been able to explore a Phys. Rev.E73, 046118 (2006). large number of the sites of the network before the final [6] V. Sood and S. Redner, Phys. Rev. Lett. 94, 178701 coalescencetakesplace. Then,we regardthe particlesas (2005). well-mixed: the asymptotic probability of finding a par- [7] V. Sood, T. Antal and S. Redner, arXiv:0712.4288 ticle on site i relative to finding it on site j is assumed (2007). to be independent of the location of the other particle [8] R. A.Blythe, Theor. Popul. Biol. 71, 454 (2007). [9] P. Trudgill, New-dialect formation (Edinburgh Univer- unless both particles are on the same site (since then sity Press, 2004). theymayreact). Inpractice,wehavefoundthepresence [10] R. Condit et al., Science 295, 666 (2002). of sufficient long-range connections in the network gives [11] K. Suchecki, V. M. Eguiluz, and M. San Miguel, Euro- risetoawell-mixedtwo-particlestateatlatetimesinthe phys.Lett. 69, 228 (2005). particle-reaction process. On the other hand, a bottle- [12] C. W. Gardiner, Handbook of Stochastic Methods neck that dramatically extends the time taken to reach (Springer, Berlin, 2004). one subset of sites fromanother would likely leadto this [13] C. Castellano, AIPConf. Proc. 779, 114 (2005). [14] P. Erdös, and A. Rényi, Publ. Math. Debrecen 6, 290 assumption of well-mixedness breaking down. Although (1959). (7)wouldthenceasetobevalid,thevariationalapproach [15] G. Baxter, R. A. Blythe, W. Croft and A. J. McKane, of [19] could, in principle, be refined to take such net- submitted. work structures into account, a task we leave for future [16] E. Gordon et al, New Zealand English: Its origins and work. Finally, we remark that the pivotal assumption evolution (Cambridge University Press, 2004). 5 [17] J.Wakeley,Coalescent theory: AnIntroduction (Roberts A38, R79 (2005). & Co., Englewood, 2006). [19] R. A.Blythe, in preparation. [18] U. Täuber, M. Howard, and B. Vollmayr-Lee, J. Phys.