Quenched disorder in the contact process on bipartite sublattices
QQuenched disorder in the contact process on bipartite sublattices
M. N. Gonzaga , C. E. Fiore and M. M. de Oliveira Departamento de F´ısica e Matem´atica, CAP, Universidade Federal de S˜ao Jo˜ao del Rei, Ouro Branco-MG, 36420-000 Brazil, Instituto de F´ısica, Universidade de S˜ao Paulo, S˜ao Paulo-SP, 05314-970, Brazil
We study the effects of distinct types of quenched disorder in the contact process (CP) with acompetitive dynamics on bipartite sublattices. In the model, the particle creation depends on itsfirst and second neighbors and the extinction increases according to the local density. The clean(without disorder) model exhibits three phases: inactive (absorbing), active symmetric and activeasymmetric, where the latter exhibits distinct sublattice densities. These phases are separated bycontinuous transitions; the phase diagram is reentrant. By performing mean field analysis and MonteCarlo simulations we show that symmetric disorder destroys the sublattice ordering and thereforethe active asymmetric phase is not present. On the other hand, for asymmetric disorder (eachsublattice presenting a distinct dilution rate) the phase transition occurs between the absorbingand the active asymmetric phases. The universality class of this transition is governed by the lessdisordered sublattice.
PACS numbers: 05.50.+q,05.70.Ln,05.70.Jk,02.50.EyKeywords: Contact process, disordered systems, symmetry-breaking, absorbing state, nonequilibrium phasetransitions
I. INTRODUCTION
In nonequilibrium systems, an absorbing-state phasetransition occurs when a control parameter such as a cre-ation or annihilation rate is varied, and the system un-dergoes a phase transition from a fluctuating state to afrozen state, with no fluctuations (the “absorbing” state).Absorbing-state phase transitions have attracted consid-erable interest in recent years since they are related tothe description of several phenomena such as populationdynamics, epidemic spreading, chemical reactions, andothers [1–4]. During the last decade, several experimen-tal realizations,e.g., in turbulent liquid crystals [5], drivensuspensions [6], superconducting vortices [7] and openquantum systems [8] have highlighted the importance ofthis kind of phase transition.In analogy with equilibrium phase transitions [9], itis expected that the critical phase transitions into ab-sorbing states belong to a finite number of universal-ity classes [4]. However, a complete classification ofthese nonequilibrium classes is still lacking. In general,absorbing-state transitions in models with short-range in-teractions and lacking a conserved quantity or symmetrybeyond translational invariance belong to the directedpercolation (DP) universality class [10]. On the otherhand, models presenting two absorbing states linked byparticle-hole symmetry are known fall in the voter modeluniversality class [11]. There are also models that arefree of absorbing states but cannot achieve thermal equi-librium because their transition rates violate the detailedbalance. An example is the majority vote model [12]. Inits ordered phase a Z symmetry is spontaneously bro-ken, leading to Ising-like behavior for spatial dimensions d ≥ dν ⊥ <
2, where d is the dimensionality and ν ⊥ is the correlation lengthexponent of the pure model. In these cases, quencheduncorrelated randomness induces the emergence of rareregions, typically located λ c (0) < λ < λ c , the λ c (0) and λ c being the critical point for the pure and disorder mod-els, respectively. Although globally the whole system isconstrained in the subcritical phase, local supercriticalregions emerge due to the presence of the disorder. Thelifetime of such “active rare regions” grows exponentiallywith the domain size, usually leading to a slow dynamics,characterized for nonuniversal exponents toward the ex-tinction for some interval of the control parameter below a r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec criticality. This behavior characterizes a Griffiths phase(GP), and was verified in DP models with uncorrelateddisorder irrespective to the disorder strength [23, 24]. Inaddition, it was shown this behavior corresponds to theuniversality class of the random transverse Ising model[23, 24]. However, some kinds of correlated disorder donot alter the critical behavior [25–27].In this work, we provide a step further by investigat-ing the effects of quenched disorder in the phase diagramof the contact processes on sublattices. Following theHarris’criterium, disorder should be relevant for the ab-sorbing phase transition, since ν ⊥ = 0.734(4), for d = 2in the clean system [23]. On the other hand, for the sym-metry breaking phase transition, the Harris criterion isinconclusive, because it corresponds to a marginal case( ν = 1 for the pure Ising model in d = 2) [28]. Here, westudy distinct kinds of disorder, (i) a random homoge-neous (ii) inhomogeneous deletion of sites, in which thedisorder strength is different in each sublattice. Inter-esting, we show, through mean-field analysis and MonteCarlo simulations, that each one of the above disorderprescriptions yields completely different outcomes.The remainder of this paper is organized as follows.In the next section we review the model and analyze itsmean-field theory. In Sec. III we present and discuss oursimulation results; Sec. IV is to summarize our conclu-sions. II. MODEL AND MEAN-FIELD THEORY
Consider a stochastic interacting particle system, de-fined on a square lattice of linear size L , where each sitecan be either occupied by a particle or empty. Each par-ticle creates a new particle in one of its first-neighborsites with rate λ , and in one of its second-neighbor siteswith rate λ . Note that in such bipartite sublattice, λ isthe creation rate in the opposite sublattice, whereas λ isthe rate in the same sublattice as the replicating particle.Therefore, unequal sublattice occupancies are favored if λ > λ . An occupied site is emptied at a rate of unity(independent of its neighboring sites). In addition to theintrinsic annihilation rate of unity, an “inhibition term”proportional to the local density is included in the dy-namics. As a consequence of this term, if the occupationfraction ρ A of sublattice A is much larger than that ofsublattice B , for instance, then any particles created insublattice B will die out quickly, stabilizing the unequalsublattice occupancies.In order to typify the model properties, we evaluate themacroscopic particle densities of each sublattice A and B given by ρ A and ρ B , respectively. In the absorbing (AB)and active symmetric (AS) phases, ρ A = ρ B = 0 and ρ A = ρ B (cid:54) = 0, respectively. Hence, ρ = ρ A + ρ B is areliable order parameter for absorbing phase transitions.Conversely, for the active asymmetric (AA) phase, it isconvenient to calculate the difference of sublattice occu-pation by φ = | ρ A − ρ B | , since φ distinguishes from the AS phase, where φ = 0.The disorder is introduced by means of a fraction Γ ofrandom deletion of sites. We shall consider two cases: thesymmetric and asymmetric, in which the sublattice remo-tion is equal (Γ A = Γ B ) and different (Γ A (cid:54) = Γ B ), respec-tively. In both cases, the disorder is quenched in spaceand time, i.e, its position or strength does not changeduring the evolution of the process. In order to achievea qualitative portrait of the phase diagram, we begin byemploying the one-site mean-field theory (MFT). For alattice with coordination number q , it results in followingcoupled equations dρ A dt = − (cid:2) µq ρ B (cid:3) ρ A + ( λ ρ B + λ ρ A )(1 − Γ A − ρ A )(1)and dρ B dt = − (cid:2) µq ρ A (cid:3) ρ B + ( λ ρ A + λ ρ B )(1 − Γ B − ρ B ) . (2) λ ρ , φ Γ = 0.0 Γ = 0.1 Γ = 0.2 Γ = 0.3 λ ρ , φ FIG. 1: (Color online) Mean-field densities φ (solid curves) e ρ (dahsed) for λ = 0 . µ = 2 .
0. Inset: detail of the data closeto the absorbing transition (linear scale).
Let us first consider that the disorder is homogeneouslydistributed in both sublattices, i.e., Γ A = Γ B = Γ. In thiscase, one derives explicit solutions for the densities as ρ = 12 k (cid:115)(cid:18) λ + λ (cid:19) + 4 k [( λ + λ )(1 − Γ) − − λ + λ (3)and φ = (cid:114) − [( λ − λ )(1 − Γ) − − λ ρ − kρ ] k . (4)The mean-field densities ρ and φ , are plotted as func-tion of λ , for distinct values of disorder and fixed λ =0 . λ = λ ABS ,c (Γ), with ρ = 0, to the activesymmetric phase, where ρ > φ = 0. Note thephase transition moves for larger values of λ when in-creasing the disorder, as expected. A further increase of λ gives rise to the AS-AA and AA-AS phase transitionsat λ = λ I ,c (Γ) and λ = λ II ,c (Γ), respectively. λ λ Γ = 0.00 Γ = 0.10 Γ = 0.20 Γ = 0.30 AAABS AS
FIG. 2:
Phase diagram in the λ − λ plane for µ = 2, showingabsorbing (ABS), active-symmetric (AS) and active asymmetric(AA) phase, for distinct values of disorder Γ. The phase diagram for distinct values of Γ is shown inFig. 2. The increase of Γ enlarges the absorbing phaseand then the AB-AS phase transition moves for largervalues of λ ABS ,c . In the active phase, the phase diagramis reentrant, and we note the size of the AA phase isreduced when the disorder increases.The asymmetric disorder case, Γ A (cid:54) = Γ B , is shown inFig. 3. MFT indicates the suppression of the AS phase,signed by a smooth change of both of ρ and φ at the same λ ,c . The AS phase is not stable for any value of λ , sothat φ does not vanish for finite values of λ . So, thephase diagram only presents two phases: the absorbingphase and the active symmetric phase separated by acontinuous phase transition.In the next section, we compare the results from MFTwith those evaluated from numerical simulations. III. RESULTS AND DISCUSSIONA. Methods
We performed extensive Monte Carlo simulations ofthe model on square lattices with periodic boundaries.The simulation scheme is as follows. First, a site is se-lected at random. If the site is occupied, it creates a λ ρ , φ Γ A =0.00, Γ B =0.00 Γ A =0.10, Γ B =0.00 Γ A =0.20, Γ B =0.00 Γ A =0.30, Γ B =0.00 λ ρ , φ FIG. 3:
Mean-field densities φ (solid curves) e ρ (dahsed) for λ =0 . µ = 2 . A (cid:54) = Γ B . Inset: detailof the data close to the absorbing transition (linear scale). particle at one of its first neighbors with a probability p = λ /W or at one of its second neighbors with a prob-ability p = λ /W . Here, W = (1 + λ + λ + µn ) isthe sum of the rates of all possible events. With a com-plementary probability 1 − ( p + p ), the site is vacated.To improve the efficiency, we choose the sites from a listcontaining the currently N occ occupied sites; accordingly,the time is incremented by ∆ t = 1 /N occ after each event.For simulations in the subcritical and critical absorbingregime, we sample the quasi-stationary (QS) regime us-ing the simulation method detailed in [29]. In order todraw a comparison with the results from MFT, in allcases we take µ = 2 . B. Symmetric disorder
We begin by analysing the symmetric case, whereΓ A = Γ B = Γ. In Fig. 4, we show typical configura-tions observed on the bipartite lattice for λ = 0 . λ = 0 .
1, for clean Γ = 0 and disordered system Γ = 0 . ρ and φ on a squarelattice with linear system size L = 80. As expected, theabsorbing phase transition occurs at a higher value of λ when disorder is introduced. We also observe thatincreasing the fraction of disorder reduces the possibilityof an AA phase, since the values of φ reduces when weincrease the fraction of disorder.In order to clarify the effects of disorder on the stabilityof the AA phase, we analyze the dependence of the orderparameter φ with the linear system size L . In Fig. 6,we compare the clean model with the disordered system,with Γ = 0 .
1. Panels (a) and (b) show the order param-eter ρ and φ versus λ for distinct systems sizes. Whilefor the clean version the the AA phase and the reentrantphase transition are observed for all system sizes (see FIG. 4: (Color online) Typical configurations observed on the bi-partite lattice for λ = 0 . λ = 0 .
1, for clean Γ = 0 . .
2. Green: inerte sites; Blue: occupied sitesand Black: empty sites. Linear system size L = 80 and µ = 2 . λ ρ , φ Γ=0.0 Γ=0.1 Γ=0.2 Γ=0.3
FIG. 5:
Order-parameters ρ and φ for distinct (symmetric) Γ’s vs. λ , for µ = 2 and λ = 0 .
1. Linear system size L = 80. panel (a)), the disordered system shows a remarkablydifferent picture, in which the sublattice occupations be-come equivalent for large L . Therefore, we conclude that φ vanishes when L → ∞ and the AA phase is suppressedby the disorder. Such behavior is reinforced by analyz-ing the order parameter variance χ = L ( (cid:104) φ (cid:105) − (cid:104) φ (cid:105) ),as shown in panels ( c ) (clean) and ( d ) (disordered). Incontrast to the clean system, in which χ diverges nearbythe transitions AB-AA and AA-AB as L → ∞ , no diver-gence is verified in such case. So, in contrast to the MFTpredictions, we observe that symmetric disorder forbidsthe stability of AA phase and therefore the disorderedsystem does not show symmetry breaking.Now, let us characterize the absorbing phase transi-tion in the disordered system. For locating the criticalcreation rate λ ,c , we study the time evolution of thenumber of active particles n ( t ) starting from an initialconfiguration the simulation with one pair of neighboringactive particles. The critical value λ ,c can be estimatedas the threshold λ separating asymptotic growth fromthe decay towards the absorbing phase. ρ , φ φ, L=20 φ, L=40 φ, L=80 φ, L=160 ρ, L=160 λ λ χ (a) (b)(c) (d) FIG. 6:
Quasi-stationary densities φ and ρ vs. λ , for µ = 2 and λ = 0 .
1, for the clean system (a), and for Γ = 0 . χ of the order parameter φ vs. λ for µ = 2 and λ = 0 . . t n ( t ) λ = 2.483λ = 2.486λ = 2.489λ = 2.490λ = 2.491λ = 2.493λ = 2.494 ln(ln(t)) -40 l n ( P s ) FIG. 7:
Number of active particles n for Γ = 0 . µ = 2 and λ = 0 .
1. Inset: survival probability for λ = 2 . .
490 and2.491.
Figure 7 exemplifies the results for Γ = 0 . λ ,c = 2 . n ( t ) ∼ ln( t/t ) θ (5)and with the survival probability P s ∼ ln( t/t ) − δ . (6)From the data in Fig.7, we estimate ln( t ) = 3 . θ = 0 . θ = 0 . -3 -2 -1 P s n λ = 2.486 λ = 2.489 λ = 2.490 λ = 2.491 λ = 2.493 ln L l n ( l n ( τ )) FIG. 8: . Log-log plot of the number of ocuppied sites n as afunction of the survival probability P s obtained from spreadingsimulations. Inset: Finite-size scaling of the lifetime of the QSstate τ . for Γ = 0 . µ = 2 and λ = 0 . t -5 -4 -3 -2 -1 ρ λ = 2.30 λ = 2.42 λ = 2.47 λ = 2.48 λ = 2.485 λ = 2.488 λ = 2.490 λ = 2.492 λ = 2.50 λ = 2.52 λ z ’ FIG. 9:
Initial decay simulations: particle density ρ vs. t for µ = 2 and λ = 0 . . L = 2000). Inset: Dynamical exponent z (cid:48) as afunction of λ . Fig. 7, we find δ = 2 . δ =1 . n ( t ) and P s ( t ), eqs. (5) and (6), we find that n ∼ P − δ/θs (7)at criticality. This relation does not depend on t , and isuseful to check our estimate of the critical point. In Fig.8,we plot n as function of P s , and observe a power law thatat λ = 2 . λ = 2 .
491 ( λ =2 . δ/θ = 0 . δ and θ obtained from the data in Fig.7.This value is also in agreement with the value δ/θ =0 . τ ∼ L ψ , (8)where ψ is a universal exponent. From the data in theinset of Fig 8, we found ψ = 0 . ψ = 0 . ρ to-wards the extinction is algebraic (with non-universal ex-ponents), rather than exponential. In Fig.9, we presentresults from initial decay simulations, where the systemstarts its dynamics from a fully occupied lattice. Weobserve the existence of a range of values of λ in thesubcritical regime where the long-time behavior of thedensity decays as a power law, ρ ( t ) ∼ t − /z (cid:48) (9)with the non-universal dynamical exponents z (cid:48) following z (cid:48) ∼ | λ − λ ,c | − ψν ⊥ , (10)where λ ,c is the critical value of the control parameter λ [24]. In the inset of Fig. 9 we observe that z (cid:48) divergeswhen λ approaches λ ,c = 2 . ψν ⊥ = 0 . ψν ⊥ ≈ .
61 reported in [24] for the diluted CP.In resume, our results for the effects of symmetric dis-order are in agreement with the predictions from the Har-ris criterium, since ν ⊥ = 0.734(4), for d = 2 in the cleansystem, and therefore the disorder is relevant for the ab-sorbing phase transition [23]. Similar trends have beenobserved for lower values of Γ, although larger crossovertimes toward the infinite-random behavior are expectedin such cases [23, 24]. C. Asymmetric disorder
Now, we investigate the effects of asymmetric disorderin which disorder strengths is different in each sublattice.This is exemplified in Fig. 10 for Γ = 0 and distinctvalues of Γ . We note that even a very small asymmetricdisorder (Γ = 0 . ρ ∼ L − β/ν ⊥ and the lifetime τ ∼ L − z .For example, for Γ = 0 .
2, shown in Fig 11, we find β/ν ⊥ = 0 . z = 1 . β/ν ⊥ = 0 . z = 1 . m = (cid:104) ρ (cid:105) / (cid:104) ρ (cid:105) goes to a universal value m = 1 . m =1 . λ ρ , φ Γ A = Γ B = 0.00 Γ A = 0.00, Γ B = 0.05 Γ A = 0.00, Γ B = 0.10 Γ A = 0.00, Γ B = 0.20 FIG. 10:
For L = 80 and assymmetric disorder: Order-parameters φ and ρ for µ = 2 and λ = 0 . A and Γ B . ln(L) -4-202468 l n ( ρ ) , l n ( τ ) λ m FIG. 11:
Assymetric disorder, with Γ = 0 . = 0 . ρ (bottom) and scaledlifetime of the QS state ln τ (top), versus ln L . Inset: Moment ratio m versus λ . on a square lattice. Parameters: µ = 2 and λ = 0 . D. Phase diagrams
Our simulation results for the effects of disorder areresumed in the phase diagrams shown in Fig.12. In (a),the clean system exhibits, besides the absorbing phase,a reentrant active asymmetric phase (with φ >
0) inside the active symmetric phase. In [13], it was shown thatthe active-absorbing phase belongs to the DP universalityclass, while the AA-AS symmetry-breaking phase transi-tion is Ising-like.The effects of asymmetric disorder in the phase dia-gram are shown in Fig.12 (b) and (c). We note that theAS phase vanishes, and there are only two phases, theabsorbing and the AA phase. The universality class ofthe transition between these phases is governed by theless disordered sublattice. So, while in the case (b), thephase transition belongs to the DP universality class, inthe case (c) it belongs to the universality class of thediluted CP.Finally, if the disorder is symmetric, as in case (d),we note that the AA phase vanishes, and the absorbing-AS phase transition belongs to the class of the dilutedCP. This system exhibits activated scaling and Griffithsphases as shown in detail in section III.B. λ λ λ λ ABS
ASABS AA AA ABS (b) (c) (d) (a)
ABS AS AA FIG. 12:
Phase diagrams for (a) clean system, (b) Γ A = 0 andΓ B = 0 .
3, (c) Γ A = 0 . B = 0 .
3, and (d) Γ A = Γ B = 0 . µ = 2 in all cases). IV. CONCLUSIONS
In this work, we have investigated the effects ofquenched disorder in the phase diagram of the competi-tive contact processes on sublattices. Through mean-fieldanalysis and Monte Carlo simulations, we have studieddistinct types of disorder, (i) a random homogeneous dele-tion of sites (ii) an asymmetrical disorder, in which thedisorder strength is different in each sublattice. Interest-ing, each one of the above disorder prescriptions yieldscompletely different outcomes.We observe that in the case (i), the disorder de-stroys the asymmetric active phase, and therefore, thesymmetry-breaking phase transition. So, there are onlytwo phases, the absorbing and the active (symmetric)phases. The absorbing-active phase transition exhibitedbelongs to the universality class of the disordered con-tact process. Effects related to an infinite randomnessfixed point are observed, such as activated dynamics andGriffiths phases in the subcritical regime.A distinct behavior is observed if each sublattice hasa different disorder strength. In such case, the symmet-rical active phase is not stable, and we observe a phasetransition directly from the active asymmetric phase tothe absorbing phase. The critical behavior in this caseis governed by the sublattice with less disorder, for in-stance, if one of the sublattices has no disorder, then thephase transition will fall in the DP class.A natural extension of the present work would be thestudy of the effects of temporal disorder [31–34] on therobustness of the symmetry-breaking phase transition, and in its the critical behavior. It is important to men-tion that, besides the theoretical interest in the field ofnonequilibrium phase transitions, suppression of activityat the nearest neighbors of active sites resembles bio-logical lateral inhibition, known to be important in thevisual system of many animals [35]. Also, our work canbe useful to understand the effects of heterogeneities inextended systems showing checkerboard pattern distribu-tions such as mutually exclusive species co-occurrences[36].
Acknowledgments
This work was supported by CNPq and FAPEMIG,Brazil. [1] J. Marro and R. Dickman,
Nonequilibrium Phase Tran-sitions in Lattice Models (Cambridge University Press,Cambridge, 1999).[2] H. Hinrichsen, Adv. Phys. , 815 (2000).[3] M. Henkel, H. Hinrichsen and S. Lubeck, Non-Equilibrium Phase Transitions Volume I: AbsorbingPhase Transitions (Springer-Verlag, The Netherlands,2008).[4] G. ´Odor, Rev. Mod. Phys , 663 (2004).[5] K. A. Takeuchi, M. Kuroda, H. Chat´e, and M. Sano,Phys. Rev. Lett. , 234503 (2007).[6] L. Cort´e, P. M. Chaikin, J. P. Gollub, and D. J. Pine,Nature Physics , 420 (2008).[7] S. Okuma, Y. Tsugawa, and A. Motohashi, Phys. Rev.B , 012503 (2011).[8] R. Guti´errez, C. Simonelli, M. Archimi, F. Castellucci,E. Arimondo, D. Ciampini, M. Marcuzzi, I. Lesanovsky,and O. Morsch Phys. Rev. A , 041602(R) (2017).[9] N. Goldenfeld, Lectures on phase transitions and therenormalization group (Addison-Wesley, 1992).[10] H. K. Janssen, Z. Phys. B , 151 (1981);P. Grassberger, Z. Phys. B , 365 (1982).[11] I. Dornic, H. Chat´e, J. Chave and H. Hinrichsen, Phys.Rev. Lett. , 045701 (2001).[12] M. J. de Oliveira, J. Stat. Phys. , 273 (1992).[13] M. M. de Oliveira and R. Dickman, Phys. Rev. E ,011125 (2011).[14] S. Pianegonda and C. E. Fiore, J. Stat. Mech. 2014,P05008 (2014).[15] M. M. de Oliveira and C. E. Fiore, J. Stat. Mech. 2017,P053211 (2017).[16] H. Hinrichsen, Braz. J. Phys. , 69 (2000).[17] A. J. Noest, Phys. Rev. Lett. , 90 (1986).[18] Moreira A G and Dickman R, Phys. Rev. E R3090(1996);Dickman R and Moreira A G, Phys. Rev. E , 035701(2006).[20] M. Bramson, R Durrett, and R. Schonmann, Ann. Prob. , 960 (1991).[21] M. S. Faria, D. J. Ribeiro and S. A. Salinas, J. Stat.Mech, P01022 (2008).[22] A. B. Harris, J. Phys. C , 1671 (1974).[23] M. M. de Oliveira and S. C. Ferreira, J. Stat. MechP11001 (2008).[24] T. Vojta, A. Farquhar and M. Mast, Phys. Rev. E ,011111 (2009).[25] M. M. de Oliveira, S.G. Alves, S.C. Ferreira and R. Dick-man, Phys. Rev. E , 031133 (2008).[26] H. Barghathi and T. Vojta Phys. Rev. Lett. , 120602(2014).[27] M. M. de Oliveira, S.G. Alves and S.C. Ferreira, Phys.Rev. E , 012110 (2016).[28] P.H.L. Martins and J.A. Plascak, Phys. Rev. E ,012102 (2007).[29] M. M. de Oliveira and R. Dickman, Phys. Rev. E ,016129 (2005); M. M. de Oliveira and R. Dickman, Braz.J. Phys. 36, 685 (2006).[30] A. O Hada and M. J. de Oliveira, J. Stat. Mech. ,043209 (2017).[31] F. Vazquez, J. A. Bonachela, C. L´opez and M. A. Mu˜noz,Phys. Rev. Lett. , 235702 (2011).[32] R. Mart´ınez-Garc´ıa, F. Vazquez, C. L´opez, and M. A.Mu˜noz, Phys. Rev. E , 051125 (2012).[33] M.M. de Oliveira and C.E. Fiore, Phys. Rev. E ,052138 (2016).[34] C.E Fiore, M.M. de Oliveira and J.A. Hoyos, Phys. Rev.E , 032129 (2018).[35] W. Lytton, From Computer to Brain (Springer-Verlag,New York, 2002).[36] E.F. Connor, M.D. Collins and D. Simberloff, Ecology94