A piecewise linear model of self-organized hierarchy formation
aa r X i v : . [ n li n . AO ] S e p A piecewise linear model of self-organized hierarchy formation
Tomoshige Miyaguchi, ∗ Takamasa Miki, and Ryota Hamada Department of Mathematics, Naruto University of Education, Tokushima 772-8502, Japan (Dated: September 15, 2020)The Bonabeau model of self-organized hierarchy formation is studied by using a piecewise linearapproximation to the sigmoid function. Simulations of the piecewise-linear agent model show thatthere exist two-level and three-level hierarchical solutions, and that each agent exhibits a transitionfrom non-ergodic to ergodic behaviors. Furthermore, by using a mean-field approximation to theagent model, it is analytically shown that there are asymmetric two-level solutions, even though themodel equation is symmetric (asymmetry is introduced only through the initial conditions), andthat linearly stable and unstable three-level solutions coexist. It is also shown that some of thesesolutions emerge through supercritical-pitchfork-like bifurcations in invariant subspaces. Existenceand stability of the linear hierarchy solution in the mean-field model are also elucidated.
I. INTRODUCTION
Hierarchy formation has been intensively studied in awide range of animal species: insects [1], fish [2–4], birds[5], and mammals [6] including even humans [7, 8]. Ithas been considered that not only differences in the priorattributes of individuals such as weight and aggressive-ness but also social interactions between individuals areimportant in the hierarchy formation [3, 9]. In fact, itis known that an individual who won an earlier contesthas a higher probability of winning later contests thanan individual who lost the earlier contest (winner-losereffects) [10]. Positive feedback generated through sucheffects might enhance the formation of hierarchies in an-imal groups.To elucidate such a feedback mechanism in hierar-chy formations, a mathematical model is proposed byBonabeau et al [11]. The Bonabeau model consists of N agents, and each agent i ( i = 1 , . . . , N ) is character-ized by a variable F i ( t ), where t is time. F i ( t ), whichis called strength or fitness in the literature, is referredto as a dominance score (DS) in this paper [5]. After acontest between two agents i and j , F i ( t ) increases if theagent i wins and decreases if i loses [the same rule is ap-plied to F j ( t )]. A greater value of F i ( t ) means a higherprobability to win a contest. In addition, the agents areassumed to perform random walks on a two-dimensionalsquare lattice L × L , and a contest occurs when two agentsmeet; thus, the density of the agents ρ = N/L is a pa-rameter, which controls the frequency of contests. Inaddition to these pairwise interactions, F i ( t ) is assumedto show a relaxation according to a differential equation dF i ( t ) /dt = − µ tanh (cid:0) F i ( t ) (cid:1) .It is found that, as the density ρ increases, theBonabeau model shows a transition from an egalitarianstate in which all F i ( t ) are equal to a hierarchical statein which F i ( t ) = F j ( t ) for some i = j . The Bonabeaumodel is one of the basic models of the hierarchy forma-tion, and it is compared with experimental observations ∗ [email protected] of hierarchies in animal groups [5, 12]. Hierarchical struc-tures can be well described by the Bonabeau model [12],but some discrepancies are also reported [5].Since the Bonabeau model is a simple model, manymodified versions have been proposed to make it morerealistic. In Refs. [13, 14], another feedback mechanismand an asymmetric rule are incorporated into the dynam-ics of F i ( t ). This generalized model (a Stauffer version)is analyzed in Ref. [15], and it was found that the egal-itarian solution is always stable, while a two-level sta-ble solution (a hierarchical solution) appears at a criticalparameter value through a saddle-node bifurcation. Inaddition, a model with a simpler relaxation dynamics dF i ( t ) /dt = − µF i ( t ) is also analyzed in Ref. [15], and itwas found that a similar transition occurs but in this casethe bifurcation is supercritical-pitchfork type. Recently,an asymmetric model is intensively studied in Ref. [16];In this asymmetric model, each agent has an intrinsic pa-rameter called a talent, which can be considered as a priorattribute of that agent. Moreover, two modified modelsare proposed in Refs.[17–19]: a timid-society model anda challenging-society model. In the timid-society model,an agent can choose a vacant site when it moves andthereby it can avoid a contest; in the challenging-societymodel, the agent chooses the strongest neighbor as anopponent.In contrast to these modifications trying to incorporaterealistic features, there are also works intending to sim-plify the Bonabeau model [20–22]. In these studies, theDSs of agents are assumed to attain only integer values,and the DS of the winner increases by one and that of theloser does not change. The dynamics can be describedby a partial differential equation in a continuum limit.This model also shows a transition from the egalitariansolution to a hierarchical solution.In spite of this diversity of models of the hierarchy for-mation, understanding of the original Bonabeau model isstill limited. For example, in the Bonabeau model withrelaxation dynamics dF i ( t ) /dt = − µF i ( t ), it is found thatthe egalitarian solution is stable at low densities (at smallvalues of ρ ). This egalitarian solution becomes unstableat ρ = ρ c , and two-level stable solutions appear through asupercritical-pitchfork bifurcation [15]. But, it seems im-possible to rigorously derive the stable range of this two-level solution (some approximation is necessary). Thisdifficulty stems mainly from nonlinearity of the sigmoidfunction employed in the Bonabeau model (See Sec. II).In this paper, we propose another simplified versionof the Bonabeau model by introducing a piecewise linearfunction in place of the sigmoid function. Piecewise linearapproximations are often used in the studies of nonlineardynamical systems. In fact, even for systems in whichrigorous approaches are difficult, more detailed analysisis possible for piecewise linear versions [23–26]. Here, wederive the stable ranges of two-level and three-level solu-tions for the piecewise-linear model. Moreover, we foundthat asymmetric two-level solutions exist even thoughthe system is symmetric (asymmetry is introduced onlythrough the initial conditions). It is also shown that var-ious stable and unstable solutions coexist.This paper is organized as follows. In Sec. II, we definetwo piecewise-linear models of the hierarchy formation:an agent model and a mean-field model. In Sec. III, linearstability analysis for steady solutions (i.e., fixed points[27]) of the mean-field model is presented. In Sec. IV, atransition from ergodic to non-ergodic behaviors in theagent model is numerically studied. Finally, Sec. V is de-voted to a discussion, in which we suggest possible gen-eralizations of the agent model. II. MODELS
In this section, we introduce two models of the self-organized hierarchy formation. The first model is re-ferred to as an agent model, and the second as a mean-field model. It is shown that the mean-field model is agood approximation of the agent model in a weak inter-action limit.
A. Agent model
Let us suppose that there are N agents, and each agent i ( i = 1 , . . . , N ) is characterized by a real number F i ( t ),which is referred to as the DS at time t [11]. F i ( t ) is ameasure of strength or fitness of the agent i [9, 15], andchanges through interactions with other agents; Here-after, the interaction between two agents is referred toas a contest. Firstly, we define the dynamics just at thecontest; Secondly, we define inter-contest dynamics byusing a Poisson process.
1. Contest dynamics
Let us define the dynamics of F i ( t ) at the contests. Atrandom time t = t n , two agents i and j contest witheach other, where i and j are randomly chosen from the gain η − η ηf ( x )0 probability(a) gain ηq ( x ) − η q ( − x ) ηf ( x )mean0probability(b) FIG. 1. Probabilities of the gain F i ( t + n ) − F i ( t − n ) (a) for themodel defined in Eq. (1), and (b) for the Bonabeau model.Here, x stands for the difference of the DSs of two contestants,e.g., x = F i ( t − n ) − F j ( t − n ). Then, q ( x ) [ q ( − x )] is the probabilitythat the agent i wins (loses) the contest with j . The meangain µ is expressed as µ = η [ q ( x ) − q ( − x )] = ηf ( x ) [Eq. (8)]. N agents. In this contest, the values of F i ( t ) and F j ( t )change as F i ( t + n ) = F i ( t − n ) + ηf (cid:0) F i ( t − n ) − F j ( t − n ) (cid:1) + ξ i ( t n ) , (1) F j ( t + n ) = F j ( t − n ) + ηf (cid:0) F j ( t − n ) − F i ( t − n ) (cid:1) + ξ j ( t n ) , (2)where t − n and t + n are the times just before and after thecontest, respectively. A gain from winning or losing thecontest is defined by the difference of the DSs beforeand after the contest, F i ( t + n ) − F i ( t − n ), which is equiv-alent to the sum of the second and the third terms inthe right side of Eq. (1). Therefore, the parameter η controls the amount of the gain, thereby characterizingthe impact of the contest result. Moreover, ξ i ( t ) is arandom variable following the normal distribution withmean 0 and variance σ , and satisfies an independenceproperty h ξ i ( t n ) ξ j ( t m ) i = δ ij δ nm σ . The gain is thus arandom variable, and its probability density is illustratedin Fig. 1(a).In Eqs. (1) and (2), f ( x ) is a non-linear function similarto the sigmoid function. In this paper, we assume it hasthe following piecewise linear form: f ( x ) = − x ≤ − F ) , x F ( − F ≤ x ≤ F ) , x ≥ F ) , (3)where F characterizes the scale of F i ( t ), and it can beremoved by rescaling [See Appendix A]. Note also that x stands for the DS difference of two contestants, e.g., x = F i ( t − n ) − F j ( t − n ) [See Eq. (1)]. This function f ( x ) isa piecewise-linear approximation to the function f b ( x ) = 11 + e − x/F −
11 + e x/F . (4)This function f b ( x ) is employed in the original Bonabeaumodel [11]. Due to the nonlinearity in f b ( x ), theoreticalanalysis of the Bonabeau model is difficult except fora few simple steady solutions. For the piecewise linearapproximation given by Eq. (3), however, more detailedanalysis of hierarchical solutions is possible due to itssimplicity.The first and the second terms on the right side ofEq. (4) have a simple probabilistic interpretation. Letus define a function q ( x ) as q ( x ) := [ f b ( x ) + 1] /
2, thenEq. (4) can be expressed as f b ( x ) = q ( x ) − q ( − x ). If x is given by x = F i ( t − n ) − F j ( t − n ), the first term q ( x )is the winning probability of i against j , and the sec-ond term q ( − x ) is the losing probability of i against j .A similar interpretation is possible also for our model[Eq. (3)]; If we rewrite f ( x ) as f ( x ) = q ( x ) − q ( − x ) with q ( x ) = [ f ( x )+1] /
2, then q ( x ) [ q ( − x )] is the probability ofwinning (losing). Thus, ηf ( x ) in the right side of Eq. (1)is the mean gain of the agent i through the contest with j . Apart from this difference in f ( x ) and f b ( x ), Eqs. (1)and (2) are still slightly different from the Bonabeaumodel, for which the dynamics is given by F i ( t + n ) = F i ( t − n ) ± η (5)with the plus sign if the agent i wins and the minus sign ifit loses (a similar equation holds for the opponent). Theprobabilities of winning and losing are given by q ( x ) and q ( − x ) defined above. Thus, the gain of the contest is arandom variable following a dichotomous distribution asshown in Fig. 1(b).The present model shown in Fig. 1(a) can be consid-ered as a coarse-grained version of the Bonabeau model.This is because a sum of several gains, each followingthe dichotomous distribution in Fig. 1(b), should followa continuous distribution similar to the one in Fig. 1(a)by virtue of the central limit theorem [28]. Therefore,our model might well be plausible for some species forwhich the same pair of individuals contest in succession[5].More precisely, if we assume that the same pair con-tests T times in succession in the Bonabeau model, thenoise terms in Eqs. (1) and (2) can be considered as small.In fact, the mean and the variance of the sum of the di-chotomous gains approximately become µ ≈ T η [ q ( x ) − q ( − x )] , (6) σ ≈ T η q ( x ) q ( − x ) . (7)Let us rescale η by replacing it with η/T , we obtain µ ≈ η [ q ( x ) − q ( − x )] , (8) σ ≈ η T q ( x ) q ( − x ) . (9)This is the situation shown in Fig. 1(a). From Eqs. (8)and (9), it is found that, if the timescale T is large, the − F i α ( t ) /F t i m e t / T r a n k α F i α ( t ) /F − . .
02 0 10 20 30initial state(a)8 agents F i α ( t ) / F rank α −
202 0 10 20 30final state (b)10 agents y F i α ( t ) / F rank α FIG. 2. (Left) Two-level hierarchy formation in the agentmodel with N = 32. Time evolution of the DS profile F i α ( t )is displayed as a function of the rank α and time t . Here, i α ( t ) is the agent index of which rank is α at time t . Theparameters η and γ are set as η = 10 − F and γη = 1 . ρ c with ρ c given by Eq. (19). (Right) The initial and final DSprofiles are shown in (a) and (b), respectively. standard deviation σ can be considered as small com-pared with the mean value µ . Therefore, in the following,we study the simplest case σ = 0 and neglect the noiseterms ξ i ( t ) and ξ j ( t ) in Eqs. (1) and (2). Note also thatthis noiseless model can be considered as a simplificationof the Bonabeau model in that the random dichotomousgains ± η in the Bonabeau model are replaced by its meanvalue µ given in Eq. (8) [See Fig. 1(b)].
2. Inter-contest dynamics
In addition to the dynamics just at the contests,we should define the inter-contest dynamics. We as-sume that the contests occur at random times t = t , · · · , t n , · · · (we set t = 0 for convenience). In theBonabeau model, the agents are assumed to perform ran-dom walks, and the times t n are determined by randomencounters of the agents [11]. However, the random walkmodel introduces non-trivial correlations in the sequenceof the intervals τ n := t n − t n − ( n = 1 , , . . . ).Here, however, we assume that these intervals τ n aremutually independent random variables, and follow theexponential distribution: w ( τ ) = γ a e − γ a τ , (10)where γ a is the interaction rate and its inverse 1 /γ a isthe mean of τ . Thus, the inter-contest dynamics is thePoisson process [28], and simplifies the model dynam-ics thanks to the independence of the intervals τ n . Inthe original Bonabeau model, the contest is consideredas a diffusion-limited reaction, while the Poisson processmight arise from a reaction-limited random walk.Note that γ a dt is the mean number of contests in thetime interval dt . Then, the mean number of contests inwhich the agent i involves is γ a dt × /N . Therefore, letus define γ := 2 γ a /N , which is an interaction rate for a − − − F i α ( t ) /F t i m e t / T r a n k α F i α ( t ) /F − . − . . .
020 10 20 30(a) initial20 agents F i α ( t ) / F rank α −
202 0 10 20 30(b) final y
18 agents F i α ( t ) / F rank α FIG. 3. (Left) Three-level hierarchy formation in the agentmodel with N = 32. Time evolution of the DS profile F i α ( t )is displayed as a function of the rank α and time t . Theparameters η and γ are set as η = 10 − F and γη = 1 . ρ c with ρ c given by Eq. (19). (Right) The initial and final DSprofiles are shown in (a) and (b), respectively. single agent. In Appendix A, we show that η and γ a (or γ ) completely characterize the agent model.Relaxation of the dominance relationship is observed inexperiments of animal groups. For example, in Ref. [3],a group of fish is assembled to form a hierarchy, theneach individual in the group is separated for long time,and finally they are assembled to form a hierarchy again.This second hierarchy is often different from the first,and thus it is considered that individual fish forgets theearlier dominance relationship.Therefore, in the meantime of the contests in ourmodel, F i ( t ) is assumed to decay. As a relaxation dy-namics, we employ the following differential equation: dF i ( t ) dt = − F i ( t ) T . (11)Here, T > ξ i ( t ) (i.e., we set σ = 0). Theinitial condition F i (0) is weakly stratified into two andthree groups as shown in Figs. 2(a) and 3(a), respec-tively. At long times, the DS profile F i ( t ) converges tostratified profiles slightly different from the initial profiles(but there are some fluctuations at the final states dueto stochastic dynamics, i.e., the random sampling of thecontestants, and the random intervals τ n ). As shown inFig. 2(b), the final state is an asymmetric two-level pro-file, whereas in Fig. 3(b), the final state is a symmetricthree-level profile. In addition, even if the parameters arethe same, there are several final profiles depending onlyon the initial conditions and realizations of the stochasticdynamics. Therefore, it is conjectured that several stableprofiles coexist at the same parameter values. B. Mean-filed model
To analyze the stable profiles in the agent model, themean-field model has been employed in previous works[11]. In contrast to the agent model, which is a stochasticmodel, the mean field model is deterministic and thusdescribed by ordinary differential equations. Here, let usapply the mean-field approximation to the agent modelintroduced in the previous subsection.If 1 /γ ≪ T , there are many contests between theagent i and the other agents in the time scale T . In ad-dition, if η ≪ F [29], then F i ( t ) does not change greatly(compared with F ) in each contest. Under these as-sumptions, changes of F i ( t ), denoted as δF i ( t ), due tocontests in the interval ( t, t + δt ) (1 /γ ≪ δt ≪ T ) canbe approximated as η γδt X k =1 f (cid:0) F i ( t ) − F j k ( t ) (cid:1) ≈ γδt ηN ′ N X j =1 j = i f (cid:0) F i ( t ) − F j ( t ) (cid:1) , (12)where N ′ is defined as N ′ := N − j k is the indexof the k -th contestant of i in the interval δt .By incorporating the relaxation term [Eq. (11)], thedynamics of F i ( t ) can be described by the ordinary dif-ferential equations dF i ( t ) dt ≈ − F i ( t ) T + γηN ′ N X j =1 j = i f (cid:0) F i ( t ) − F j ( t ) (cid:1) . (13)This equation (13) is the same form as the Bonabeau’smean-field model [11], in which the function f ( x ) is givenby Eq. (4). Here, however, we employ the piecewise linearfunction given in Eq. (3). III. LINEAR STABILITY ANALYSIS OFMEAN-FIELD MODEL
In this section, we study steady solutions of the mean-field model with T = 1: dF i ( t ) dt = − F i ( t ) + ρN ′ N X j =1 j = i f (cid:0) F i ( t ) − F j ( t ) (cid:1) , (14)where ρ ≥ ρ = γη . Moreover, f ( x ) inEq. (14) is assumed to be given by Eq. (3) with F = 1[i.e., Eq. (A4) in Appendix A]. In Appendix A, it is shownthat such simplifications do not lead to loss of generality,and that ρ is the only parameter of the mean-field model.In the figures, however, we give the units explicitly.It can be shown that the total DS defined by S ( t ) = P Ni =1 F i ( t ) follows the equation dS/dt = − S , and thus S ( t ) decays to zero as t → ∞ . Therefore, any stablesteady solution F i ( t ) ≡ F ∗ i satisfies P Ni =1 F ∗ i = 0. A. Single-level solution (egalitarian solution)
It is easy to see that F i ( t ) ≡ i = 1 , . . . , N ) is asteady solution of Eq. (14) for any values of ρ ≥
0. TheJacobian of the right side of Eq. (14) at this solution isgiven by a circulant matrix J ,N ( a ) = a b · · · b bb a . . . b b ... . . . . . . . . . ... b b . . . a bb b · · · b a , (15)where a := ρ/ − b := − ρ/ (2 N ′ ). For later use,the Jacobian is denoted as J ,N ( a ) to indicate that it is a(matrix-valued) function of a . The eigenvalues of J ,N ( a )are given by λ = a + ( N − b, a − b (16)= − , N N ′ ρ − . (17)The multiplicity of the first eigenvalue a + ( N − b is 1and that of the second a − b is N − F i ( t ) ≡ ρ satisfies ρ < N ′ N =: ρ c , (18)where we define the critical value ρ c . Note that, for thegeneral case with F = 1 and T = 1, this definitionbecomes ρ c := 2 N ′ N F T . (19)For ρ > ρ c , the single-level solution is unstable. Thisis consistent with the corresponding result in Ref. [11].Note that N − λ = ρ/ρ c − ρ = ρ c . B. Two-level solution
Two-level asymmetric solutions have been studied inprevious works [15, 16], but in these studies, asymmetryis incorporated directly into the model equations. Here,however, we show that there exist stable asymmetric so-lutions even in the symmetric model given by Eq. (14)(Asymmetry is incorporated through the initial condi-tion). Moreover, linearly stable ranges in terms of ρ arederived for these asymmetric solutions.Let us study steady two-level solutions with asymme-try: F i ( t ) ≡ ( F u ( i ≤ m ) , − F l ( i > m ) , (20) ρρ c m . − m = 12 F i ( t ) / F agent index i − m = 8 F i ( t ) / F agent index i FIG. 4. (Left) Phase diagram ( ρ vs m ) of two-level stablesolutions. The total number of the agents is N = 32. Onthe horizontal solid lines, the two-level solutions are stable.Dashed lines are theoretical prediction Eq. (25). Arrows in-dicate the parameter values used in the right figures. (Right)Examples of asymmetric two-level solutions obtained by nu-merical simulations. The density ρ is set as ρ/ρ c = 1 .
3. Thenumber of the upper-level agents m is (a) m = 12, and (b) m = 8. A weakly hierarchical state F i (0), similar to the oneshown in Fig. 2(a), is used as the initial condition, which issorted as F i (0) > F j (0) for i < j . Note that this order of F i ( t ) does not change with t in the mean-field model [i.e., i α ( t ) ≡ α ], thus the agent index i is used as the horizontalaxis. where the constants F u , F l are positive F u , F l > m = 1 , . . . , N/
2. This parameter m is the number ofthe upper-level agents [ F i ( t ) ≡ F u ]; m = N/ m > N/ m < N/ F i ( t ) → − F i ( t ). However, we omit these cases m > N/ F u and F l can be determined by settingthe right side of Eq. (14) zero. Thus, we obtain ( ρ N − mN ′ f ( F u + F l ) − F u = 0 ,ρ mN ′ f ( F u + F l ) − F l = 0 . (21)If F u + F l ≤
2, then f ( F u + F l ) = ( F u + F l ) /
2, andtherefore we have F u = F l = 0 from Eq. (21) . Thus, F u + F l > F u = ρ N − mN ′ , F l = ρ mN ′ . (22)Since F u + F l = 2 ρ/ρ c >
2, the two-level solution existsfor ρ > ρ c .Linear stability analysis can be carried out in the sameway as the previous subsection. The Jacobian at the two-level steady solutions is given by J ,m = (cid:18) J ,m ( a ) OO J ,N − m ( a ′ ) (cid:19) , (23)where J ,m ( a ) is an m × m matrix of the form of Eq. (15)but with a = ρ ( m − / (2 N ′ ) − b is the same as that inEq. (15)], J ,N − m ( a ′ ) is an ( N − m ) × ( N − m ) matrix with a ′ = ρ ( N − m − / (2 N ′ ) −
1, and O is a zero matrix. Byusing Eq. (16), it is easy to find the eigenvalues of J ,m as λ = − , ρ m N ′ − , ρ N − m N ′ − , (24)with multiplicities 2 , m −
1, and N − m −
1, respectively.Therefore, the two-level stable solution with m exists for ρ satisfying 1 < ρρ c < NN − m . (25)Thus, at ρ = ρ c , the steady solution of Eq. (14) changesabruptly from F i ( t ) ≡ f ( x ) isnot differentiable. Even for ρ/ρ c > N/ ( N − m ), thetwo-level solutions with m exist, but they are unstablebecause the third eigenvalue in Eq. (24) becomes positive.In Fig. 4, the ranges of ρ where a stable two-level so-lution exists are displayed by horizontal lines. The sym-metric solution ( m = N/
2) has the widest stable range;the stable range is shorter for stronger asymmetry (i.e.,for smaller values of m ). Asymmetric solutions shownin Figs. 4 (a) and (b) are obtained by numerical simu-lations; these solutions resemble the result for the agentmodel shown in Fig. 2(b).As shown in Fig. 4(Left), the two-level solutions ( m =1 , . . . , N −
1) appear simultaneously at ρ = ρ c throughbifurcations of the pitchfork type (though there is adiscontinuity). This can be easily checked by setting F i ( t ) = F u ( t ) for ( i ≤ m ), and F i ( t ) = − F l ( t ) for( i > m ); this form of the trajectory F i ( t ) is a solutionin an invariant two-dimensional subspace. If we define∆ F ( t ) = F u ( t ) + F l ( t ), it is easy to show that d ∆ F ( t ) dt = − ∆ F ( t ) + 2 ρρ c f (cid:0) ∆ F ( t ) (cid:1) . (26)Examining the functional form of the right side [as a func-tion of ∆ F ( t )] for ρ < ρ c and ρ > ρ c , it is found that thebifurcation at ρ = ρ c is the superciritical pitchfork type[27]. Note however that this analysis in invariant sub-spaces is insufficient for a proof of the linear stability. SeeAppendix B for a similar argument on three-level solu-tions, for which two stable solutions appear also throughthe superciritical pitchfork bifurcations, but they are un-stable in some directions perpendicular to the invariantsubspaces. C. Three-level solution
There are also many three-level steady solutions, andthus here we focus only on symmetric three-level solu-tions. In this subsection, steady solutions of the following ρρ c m
321 1 . − m = 10 F i ( t ) / F agent index i − m = 6 F i ( t ) / F agent index i FIG. 5. (Left) Phase diagram ( ρ vs m ) of three-level stablesolutions. The total number of the agents is N = 32. Onthe horizontal solid lines, the three-level solutions are sta-ble. Dashed lines are theoretical prediction Eq. (31). Ar-rows indicate the parameter values used in the right figures.(Right) Examples of three-level solutions obtained by numer-ical simulations. The density ρ is set as ρ/ρ c = 1 .
5. Half thenumber of the middle-level agents m is (a) m = 10, and (b) m = 6. Weakly hierarchical states, similar to the one shownin Fig. 3(a), are used as initial conditions. form are shown to be stable: F i ( t ) ≡ F (cid:0) ≤ i ≤ N − m (cid:1) , (cid:0) N − m < i ≤ N + m (cid:1) , − F (cid:0) N + m < i ≤ N (cid:1) . (27)Here, 2 m is the number of the middle-level agents, forwhich F i ( t ) ≡
0; therefore m should satisfy 0 < m 2. Moreover, the constant F is assumed to satisfy F > F < 2, there exist some steady solutions,but they are linearly unstable. See Appendix B). Substi-tuting Eq. (27) into the right side of Eq. (14), we foundthat F = ρρ c + ρmN ′ . (28)Since we assume F > 2, the steady solution [Eq. (28)]exists for ρ > N ′ / ( N + 2 m ).The Jacobian of these steady solutions is given by J F > ,m = J ,N/ − m ( a ) O OO J , m ( a ′ ) OO O J ,N/ − m ( a ) , (29)where J ,N/ − m ( a ) is an ( N/ − m ) × ( N/ − m ) matrix ofthe form of Eq. (15) with a = ρ ( N ′ − − m ) / (4 N ′ ) − J ,m ( a ′ ) is an 2 m × m matrix with a = ρ (2 m − / (2 N ′ ) − 1. By using Eq. (16), we obtain the eigenval-ues of J F > ,m as λ = − , ρ N − m N ′ − , ρ mN ′ − , (30)with multiplicities 3 , N − m − 2, and 2 m − 1, respec-tively. Therefore, the three-level stable solution with m [Eq. (27)] exists for ρ satisfying2 NN + 2 m < ρρ c < max (cid:18) NN − m , N m (cid:19) . (31)Even for ρ/ρ c larger than this upper bound, the three-level solutions exist, but they are unstable because thesecond or the third eigenvalues in Eq. (30) become posi-tive.In Fig. 5, the ranges of ρ where the stable three-levelsolutions exist are displayed by horizontal lines. Thewidest stable range is at m = N/ 6, at which the threelevels have the equal numbers of agents (the exampleshown in the figure is for N = 32, and thus N/ N is a multiple of 3, there is a steady solutionfor which each level has N/ D. N -level solution (linear hierarchy) Linear hierarchies are frequently observed in animalsocieties. In a linear hierarchy, if an individual A domi-nates B and B dominates C, then A dominates C [3] (i.e.,a transitive relationship). At high values of ρ , there existsa steady N -level solution, in which each agent has a dif-ferent values of F i ( t ). This completely stratified solutionis reminiscent of the linear hierarchy.Here, let us assume the following solution F i ( t ) ≡ F − FN ′ ( i − , ( i = 1 , . . . , N ) , (32)where F is a constant to be determined, and we alsoassume F > N ′ . In order that the above F i ( t ) is a steadysolution, i.e., dF i ( t ) /dt ≡ F should satisfy F = ρ. (33)In the derivation, we used the assumption F > N ′ as F i ( t ) − F j ( t ) = 2 F ( j − i ) /N ′ > F/N ′ > 2, where i < j .From F > N ′ and Eq. (33), ρ should also satisfy ρ > N ′ ,or ρρ c > N . (34)Thus, the N -level solution exists only at large ρ .The stability of the N -level solution is easy to prove.The Jacobian of this steady solution is simply given by J N = − I , where I is the N × N identity matrix. There-fore, the N -level solution [Eq. (32)] is stable. IV. ERGODICITY IN AGENT MODEL As shown in Figs. 2 and 3, the agent model behavessimilarly to the mean-field model for 1 /γ ≪ T and η ≪ F . But, if these conditions are not fulfilled, the agentmodel behaves differently from the mean-field model. Inthis section, the dependence of the agent model on theseparameters γ and η is numerically studied.As a quantity characterizing the dynamics of the agentmodel, we use the standard deviation σ ( γ, η ) of the time-averaged DS, F i , defined as µ ( γ, η ) := 1 N N X i =1 F i , (35) σ ( γ, η ) := 1 N N X i =1 [ F i − µ ( γ, η )] . (36)The time averaged DS, F i , is defined as F i := 1 T Z T F i ( t ) dt, (37)where the dynamics of F i ( t ) is given by Eq. (1) with theparameters η and γ . Thus, the standard deviation σ ( γ, η )also depends on these parameters.If the system is ergodic, a time average tends to asingle value, which is equal to the ensemble average, in along time limit ( T → ∞ ). In the agent model [Eqs. (1)and (2)], all the agents are equivalent, and therefore thelimiting value is the same for all the agents, and thusit follows that σ ( γ, η ) vanishes at T → ∞ . Accordingly, σ ( γ, η ) can be used as a parameter of ergodicity breaking[30, 31]. In Fig. 6, we set γη is constant (i.e., ρ = γη isconstant) to fix the corresponding mean-field model [SeeEq. (13)], and numerically obtain the variance σ ( ρ/η, η )as a function of η .As shown in Fig. 6 (Left), the standard deviation σ ( ρ/η, η ) is far away from zero for small values of η . Infact, the agents are separated into two groups as shownin the inset of Fig. 6(Left); these two groups correspondto the two-level solution in the mean-field model with m = N/ η , the members of thesetwo groups rarely change in the course of time evolution,as shown in Fig. 6(a), where two typical trajectories F i ( t )are displayed.For large values of η , the agents are still separated intotwo groups again [See the inset of Fig. 6(Left)], but theagents frequently move from one group to the other asshown in Figs. 6(b) and (c). Accordingly, all the time av-erages F i ( i = 1 , . . . , N ) tend to zero as T increases, andtherefore the standard deviation σ ( ρ/η, η ) also vanishesas shown in Fig. 6 (Left). The transitions of the agentsfrom one group to the other occur, because the impactof each contest becomes significant for large η [though atime average of this effect, given by γη , is the same in allthe numerical simulations in Figs. 6 (a)–(c)]. It shouldbe also noted that, even though the time average F i van-ishes at large η , a hierarchy exists in snapshots F i ( t ) asshown in the inset of Fig. 6(Left), where the agents areseparated into two groups, and thus the system is notegalitarian. . . . . . . 40 0 . 02 0 . 04 0 . − 202 0 10 20 30 F i α /F S t a nd a r dd e v i a t i o n σ ( ρ /η , η ) η/F rank α − 202 (a) − 202 (b) − F i ( t ) / F time t/T FIG. 6. (Left) Standard deviation σ ( ρ/η, η ) vs η in the agentmodel with N = 32. The value of ρ is fixed as ρ = 1 . ρ c ,for which the egalitarian solution F i ( t ) ≡ contests occur.The inset is a snapshot of F i α ( t ) vs rank α for η = 0 . F (circle), η = 0 . F (triangle), and η = 0 . F (square).(Right) Typical trajectories of F i ( t ) for (a) η = 0 . F , (b) η = 0 . F , and (c) η = 0 . F . In (a), two trajectories aredisplayed, whereas a single trajectory is displayed in (b) and(c). At small η , the ergodicity seems to be violated asshown in Fig. 6(a). However, it is probable that it justtakes too long time to observe transitions of the agentsfrom one group to the other, and thus the ergodicitymight not be violated. This is because a sequence of con-tests at large η which causes a transition of an agent canbe possible, in principle, to occur even at small η (thoughthe probability of occurrence of such sequence of contestsis quite small). Therefore, the observed violation of theergodicity might well be just apparent. V. DISCUSSION Since the appearance of the seminal paper [11], theBonabeau model has been employed to explain experi-mental data of animal hierarchy formations, and manymodified versions have been proposed [13, 14, 17–19].But, understanding of the original Bonabeau model hasnot been far from satisfactory due to difficulty in treat-ing its nonlinearity. In this paper, a piecewise linear ver-sion of the Bonabeau model was introduced. By usingthe mean-field approximation, it was shown that thereare many asymmetric solutions, and that coexistence ofthe stable solutions takes place. In addition, an appar-ent transition in ergodic behaviors is found in the agentmodel.Our model assumed that encounters of the agents arecompletely random. Namely, at each contest time t n , theagents i and j are randomly chosen from the N agents.But, it is known that if the agents i and j contest, then these agents i and j are more likely to contest in thenext contest event than other agents [5]. Remarkably, itis also shown in Ref. [5] that the persistent time duringwhich the same individuals successively contest follows apower-law distribution. Such a non-Markovian memoryeffect can be easily implemented in the agent model, byintroducing a persistent-time distribution [31, 32] w p (˜ τ ) ≃ a ˜ τ α ( τ → ∞ ) , (38)where a and α are positive constants. We choose asequences of persistent times ˜ τ , ˜ τ , · · · , each following w p (˜ τ ), and define renewal times as ˜ t n := P nk =1 ˜ τ k , atwhich the contestants change. In each interval [˜ t n − , ˜ t n ],the same agents i and j contest. This generalized modelshould be studied in future works.The linear hierarchy, frequently observed in animal so-cieties, is characterized by the transitive relationship (SeeSec. III D); however, intransitive relationships are alsoobserved by suppressing group processes [3]. Such intran-sitive relationships cannot be described by the Bonabeaumodel, because it is always transitive from its defini-tion; i.e., if F i ( t ) > F j ( t ) and F j ( t ) > F k ( t ), then F i ( t ) > F k ( t ). To describe the intransitive relation-ships, it is necessary to introduce an anti-symmetric ma-trix F ij ( t ) which describes the dominance relationshipbetween i and j . In the Bonabeau model, F ij ( t ) could bedefined by F ij ( t ) := F i ( t ) − F j ( t ), but the matrix F ij ( t )cannot be described by a single vector in general.Therefore, future work is needed to develop a gener-alized model for F ij ( t ), and to elucidate how the transi-tive relationship [i.e., if F ij ( t ) > F jk ( t ) > 0, then F ik ( t ) > 0] emerges (or self-organizes). In such a general-ized model, a bystander effect should be incorporated, inaddition to the winner/loser effects [3, 4]. The bystandereffect is a mechanism that an individual who witnessesa contest between other individuals is influenced by theresult of that contest; the witness might learn its statusvicariously by observing contests between other individ-uals [4]. Without such a bystander effect, intransitiverelationships should be frequently observed [3].Finally, we neglect the noise terms in Eqs. (1) and (2)in this paper, and thus contest dynamics is purely de-terministic except the random choice of two contestants.In real societies, however, contestants have some randomfactors such as their physical conditions. Therefore, thenoise terms might be important and should be studied infuture works. Appendix A: Rescaling In this Appendix, the agent model and the mean-field-model are transformed into simpler forms by introduc-ing rescaled variables. Let us define the rescaled (non-dimensional) variables as¯ t = tT , ¯ F i (¯ t ) = F i ( t ) F . (A1)Then, Eqs. (1) and (2) can be rewritten as (we omit thenoise terms)¯ F i (¯ t + n ) = ¯ F i (¯ t − n ) + ¯ η ¯ f ( ¯ F i (¯ t − n ) − ¯ F j (¯ t − n )) , (A2)¯ F j (¯ t + n ) = ¯ F j (¯ t − n ) + ¯ η ¯ f ( ¯ F j (¯ t − n ) − ¯ F i (¯ t − n )) , (A3)where ¯ η and ¯ f ( x ) are defined respectively as ¯ η = η/F and ¯ f ( x ) = − x < − , x ( − ≤ x ≤ , x > . (A4)The exponential distribution of the intervals τ [Eq. (10)] is also rescaled as¯ w (¯ τ ) = ¯ γ a e − ¯ γ a ¯ τ , (A5)where ¯ γ a = γ a T . The relaxation dynamics [Eq. (11)] issimply given by d ¯ F i (¯ t ) d ¯ t = − ¯ F i (¯ t ) . (A6)Therefore, the two parameters ¯ η and ¯ γ a completely char-acterize the agent model.Similarly, by using the transformations in Eq. (A1),the mean-field model in Eq. (13) becomes d ¯ F i (¯ t ) d ¯ t = ¯ ρN ′ N X j =1 j = i ¯ f (cid:0) ¯ F i (¯ t ) − ¯ F j (¯ t ) (cid:1) − ¯ F i (¯ t ) . (A7)where ¯ ρ is defined as ¯ ρ = ρT /F . Thus, ¯ ρ is the onlyparameter of the mean-field model. Note that, even if ¯ ρ is constant, corresponding parameter values in the agentmodel (¯ η and ¯ γ a ) are not uniquely determined, because¯ ρ = ¯ γ ¯ η with ¯ γ := 2¯ γ a /N . Appendix B: Unstable three-level solution In Sec. III C, we study stable three-level solutions[Eq. (27)], but there also exist unstable three-level so-lutions. In this Appendix, we show that the three-levelunstable solutions emerge simultaneously at ρ = ρ c , andthese unstable solutions become stable at some values of ρ > ρ c .First, it is easy to show that a steady three-level solu-tion of the form of Eq. (27) does not exist for 0 < F < F > < F < 2. For1 < F < 2, the three-level solution is given by F = ρ N − mN ′ − ρm , (B1) . . . . . m = 0 (two levels) m = 1 m = 2 m = 3 m = 4 m = 5 m = 6 m = 7891011121314 F / F ρ/ρ c FIG. 7. Dominance score F for two-level and three-level sym-metric solutions as functions of ρ . The number of agents areset as N = 32. The numbers in the figure are the correspond-ing values of m . The dashed curves are unstable solutionsgiven by Eq. (B1), and the solid curves are stable two-leveland three level solutions given by Eq. (22) with m = N/ m still exist, but they are unstable (SeeSecs. III B and III C). These unstable ranges are not shownfor brevity. with m = 0 , . . . , N/ − m = 0 corresponds to thesymmetric two-level solution). Accordingly, the range of ρ satisfying 1 < F < < ρρ c < NN + 2 m . (B2)Note that the upper bound is equivalent to the lowerbound of the stable three-level solution [Eq. (31)]. In fact,the stability of the three-level solutions for each value of m changes at ρ/ρ c = 2 N/ ( N + 2 m ) as shown below. Thebifurcation at this point might be subcritical-pitchforktype (however, we should note that details of this bifur-cation remain still unclear).Next, the stability of the three-level solutions[Eq. (B1)] is elucidated. In this case, the Jacobian matrixis given by J 1, and B is an ( N/ − m ) × m matrix, all the elements of whichare the same and given by − ρ/ (2 N ′ ). B t is the transposeof B .After a somewhat lengthy but elementary calculation,we obtain four eigenvalues of the Jacobian in Eq. (B3).0 − . − − . − − . . . . − − − d F ( t ) / d t F ( t ) ρ/ρ c = 0 . ρ/ρ c = 0 . ρ/ρ c = 1 . ρ/ρ c = 1 . FIG. 8. dF ( t ) /dt in Eq. (B7) vs F ( t ) for four different valuesof ρ : ρ/ρ c = 0 . . . . N and m are set as N = 32and m = 6. For ρ < ρ c , dF ( t ) /dt is monotonically decreasing,and thus the origin F ( t ) ≡ ρ >ρ c , however, the origin is unstable, and two stable fixed points,that correspond to Eq. (B1), appear. Note that these stablefixed points are stable only in the invariant subspace, andunstable in some directions perpendicular to this subspace. Two of the four eigenvalues are given by λ = ρρ c − , ρ N + 2 m N ′ − , (B4)with multiplicities 2 m and N − m + 1), respectively.The first eigenvalue in Eq. (B4) is positive because ofEq. (B2). Therefore, the three-level solutions in Eq. (B1)are unstable. Note however that the first eigenvalue doesnot exist for m = 0 (i.e., for the two-level symmetric solu-tion), because the multiplicity becomes zero, and there-fore it is not contradicting with the fact that the two-level solution is stable (See Sec. III B). Note also thatthe second eigenvalue is negative, and it does not exist,if m = N/ − 1. The remaining two eigenvalues are givenby λ = − , ρ mN ′ − . (B5)These eigenvalues are simple and negative.A phase diagram of the symmetric two- and three-levelsolutions are displayed in Fig. 7. These steady solutionsemerge at ρ = ρ c , but only two-level solution is stable,and all the three-level solutions are unstable. For F > F i ( t ) = F ( t ) (cid:0) ≤ i ≤ N − m (cid:1) , (cid:0) N − m < i ≤ N + m (cid:1) , − F ( t ) (cid:0) N + m < i ≤ N (cid:1) , (B6)where m = 0 , , . . . , N/ − 1, and F ( t ) can be eitherpositive or negative. The time evolution equation for F ( t ) is obtained by inserting Eq. (B6) into Eq. (14) as dF ( t ) dt = − F ( t ) + ρN ′ (cid:20) mf (cid:0) F ( t ) (cid:1) + N − m f (cid:0) F ( t ) (cid:1)(cid:21) = (cid:16) ρρ c − (cid:17) F ( t ) [0 < F ( t ) < , (cid:0) ρ mN ′ − (cid:1) F ( t ) + ρ N − m N ′ [1 < F ( t ) < , − F ( t ) + ρ N +2 m N ′ [2 < F ( t )] , (B7)where the equation only for F ( t ) > F ( t ) < f ( x ) given in Eq. (A4) is an odd func-tion. Note also that the slope ρm/N ′ − ρ > ρ c , corresponds to the secondeigenvalue in Eq. (B5).From the first equation in the right side of Eq. (B7),the single-level solution F ( t ) ≡ ρ < ρ c andunstable for ρ > ρ c . The two-level ( m = 0) and three-level ( m > 0) solutions emerge at ρ = ρ c simultaneously,and they are stable because of ρm/N ′ − < ρ > ρ c .Due to the symmetry, − F ( t ) is also a solution in theinvariant subspaces, and thus there are two stable fixedpoints in each invariant subspace with m .This bifurcation is readily understood by a phase dia-gram [27] shown in Fig. 8, in which dF ( t ) /dt in Eq. (B7)is displayed as a function of F ( t ). It is clear that the bi-furcation at ρ = ρ c can be considered as a superciriticalpitchfork type. Although the bifurcations are pitchforktype and thus the two emerged fixed points are stable,these fixed points except the two-level solutions ( m = 0)are stable only in the invariant subspaces; in fact, theyare unstable in some directions perpendicular to the sub-spaces, because the first eigenvalue in Eq. (B4) is positive. [1] P. S. Oliveira and B. Holldobler, Behav. Ecol. Sociobiol. , 385 (1990).[2] C. Goessmann, C. Hemelrijk, and R. Huber, Behav.Ecol. Sociobiol. , 418 (2000).[3] I. D. Chase, C. Tovey, D. Spangler-Martin, and M. Man-fredonia, Proc. Natl. Acad. Sci. U.S.A , 5744 (2002).[4] L. Grosenick, T. S. Clement, and R. D. Fernald, Nature , 429 (2007).[5] W. B. Lindquist and I. D. Chase, Bulletin of mathemat-ical biology , 556 (2009).[6] R. M. Wittig and C. Boesch, Int. J. Primatology , 847(2003).[7] R. C. Savin-Williams, J. Youth Adolescence , 75 (1980).[8] C. Garandeau, I. Lee, and C. Salmivalli, J. Youth Ado- lescence , 1123 (2014).[9] C. Castellano, S. Fortunato, and V. Loreto,Rev. Mod. Phys. , 591 (2009).[10] Y. Hsu, R. L. Earley, and L. L. Wolf, Biol. Rev. , 33(2006).[11] E. Bonabeau, G. Theraulaz, and J.-L. Deneubourg,Physica A , 373 (1995).[12] E. Bonabeau, G. Theraulaz, and J.-L. Deneubourg, Bul-letin of mathematical biology , 727 (1999).[13] D. Stauffer, Int. J. Mod. Phys. C , 237 (2003).[14] D. Stauffer and J. Sa Martins, Advances in Complex Sys-tems , 559 (2003).[15] L. Lacasa and B. Luque,Physica A: Statistical Mechanics and its Applications , 472 (2006).[16] M. P´osfai and R. M. D’Souza, Phys. Rev. E , 020302.[17] T. Odagaki and M. Tsujiguchi, Physica A , 435(2006).[18] M. Tsujiguchi and T. Odagaki, Physica A , 317(2007).[19] T. Okubo and T. Odagaki, Phys. Rev. E , 036105(2007).[20] E. Ben-Naim and S. Redner, J. Stat. Mech. , L11002(2005). [21] E. Ben-Naim, F. Vazquez, and S. Redner, Euro. Phys.J. B , 531 (2006).[22] E. Ben-Naim, S. Redner, and F. Vazquez, Europhys.Lett. , 30005 (2007).[23] R. L. Devaney, Physica D , 387 (1984).[24] S. Tasaki and P. Gaspard, J. Stat. Phys. , 803 (2002).[25] T. Miyaguchi, Prog. Theore. Phys. , 31 (2006).[26] T. Miyaguchi and Y. Aizawa, Phys. Rev. E , 066201(2007).[27] S. Strogatz, Nonlinear dynamics and chaos: with appli-cations to physics, biology, chemistry, and engineering (Westview Press, Cambridge, MA, 1994).[28] W. Feller, An Introduction to Probability Theory and itsApplications , 2nd ed., Vol. II (Wiley, New York, 1971).[29] More precisely, ηγδt ≪ F should be satisfied. If we canchoose δt satisfying this condition with 1 /γ ≪ δt ≪ T ,the approximation in Eq. (12) is valid.[30] Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev.Lett. , 058101 (2008).[31] T. Miyaguchi and T. Akimoto, Phys. Rev. E , 032130(2013).[32] T. Miyaguchi, T. Uneyama, and T. Akimoto, Phys. Rev.E100