Excitable Greenberg-Hastings cellular automaton model on scale-free networks

0.15 MB·English

by  An-Cai Wu

ExcitableGreenberg-Hastings cellularautomaton model onscale-free networks
An-Cai Wu,1 Xin-Jian Xu,2 and Ying-Hai Wang1
1Instituteof Theoretical Physics, Lanzhou University, Lanzhou Gansu 730000, China
2Departamento de F´ısica da Universidade de Aveiro, 3810-193 Aveiro, Portugal
(Dated:February6,2008)
7
0 WestudytheexcitableGreenberg-Hastingscellularautomatonmodelonscale-freenetworks. Weobtained
0 analyticalexpressionsfornoexternalstimulusandtheuncoupledcase. Itisfoundthatthecurves,theaverage
2 activityF versustheexternalstimulusrater,canbefittedbyaHillfunction,butnotexactly,andthereexists
n arelationF ∝ rα for thelow-stimulusresponse, whereStevens-Hillexponent αranges fromα = 1inthe
a subcritical regime to α = 0.5 at criticality. At the critical point, the range reaches the maximal. We also
J calculatetheaverageactivityFk(r)andthedynamicrange∆k(p)fornodeswithgivenconnectivityk. Itis
1 interestingthatnodeswithlargerconnectivityhavelargeroptimalrange,whichcouldbeappliedinbiological
1 experimentstorevealthenetworktopology.

] PACSnumbers:87.10.+e,87.18.Sn,05.45.-a
h
c
e Many systems in the real world, either naturally evolved [10]. Analytical results have recently been obtained for the
m
or artificially designed, are organized in a network fashion one-dimensionalcellularautomatonmodelunderthetwo-site
- [1]. Themoststudiednetworksaretheexponentialnetworks mean-field approximation [11]. In ref. [12], Kinouchi and
t
a in which the nodes’ connectivity distribution (the probabil- Copelli studied the GHCA model on Erdo¨s-Re´nyi random
t
s ity P(k) that a node if connected to other k nodes) is expo- graph with stochastic activity propagation, and they found a
. nentially bounded. A typical model of this network is the new and important result that dynamic range maximized at
t
a Watts-Strogatz (WS) graph [2] which exhibits the “small- the critical point, whichis perhapsthe first clear exampleof
m
world”phenomenon that was observed in realistic networks signalprocessingoptimizationataphasetransition.
- [3]. On the contrary, it has been observed recently that the
d In the n-state GHCA model [9] for excitable systems,
nodes’degreedistributionofmanynetworkshavepower-law
n the instantaneous membrane potential of the i-th cell (i =
tails(scale-freeproperty,duetotheabsenceofacharacteristic
o 1,...,N) at discrete time t is represented by x (t) ∈
c valueforthedegrees)P(k)∼k−γ with2≤γ ≤3[4]. {0,1,...,n−1},n ≥ 3. Thestatex (t) = 0denotesianeu-
i
[
Scale-free (SF) networks have been widely studied in the ron at its resting (polarized) potential, xi(t) = 1 represents

1 pastfewyears,mainlybecauseoftheirconnectiontoalotof a spiking (depolarizing) neuron and xi(t) = 2,...,n − 1
v real-world structures, including networks in nature, such as account for the afterspike refractory period (hyperpolariza-
8 the cell, metabolicnetworksandthe foodweb, artificialnet- tion). There are two ways for the i-th element to go from
4
workssuchastheInternetandtheWWW,orevensocialnet- state x (t) = 0 to x (t+1) = 1: a) due to an externalsig-
2 i i
works, such assexualpartnershipnetworks[1]. The SFnet- nal, modelled here by a Poisson process with rate r (which
1
workischaracterizedwiththepower-lawdegreedistribution. impliesatransitionwithprobabilityλ=1−exp(−r∆t)per
0
7 Thustherealways existsa smallnumberof nodeswhich are time step); b) with probability p, due to a neighbor j being
0 connected to a large number of other nodes in SF networks. in the excited state in the previous time step. If x (t) ≥ 1,
i
/ This heterogeneity leads to intriguing properties of SF net- thenx (t+1) = (x (t)+1)modn,regardlessofthestimu-
t i i
a works. A lot of work has been devoted in the literature to lus. In other words, the rulesstate that a neurononlyspikes
m
thestudyofstaticpropertiesofthenetworks,whilemanyin- if stimulated,after whichitundergoesanabsoluterefractory
- terestsaregrowingfordynamicalpropertiesonthesekindof periodbeforereturningtorest. Timeis discrete. We assume
d
networks. Inthecontextofpercolation,SFnetworksaresta- ∆t=1mswhichcorrespondstotheapproximatedurationof
n
ble against random removal of nodes while they are fragile a spikeand isthe time scale adoptedforthe timestep of the
o
c underintentionalattackstargetingonnodeswithhighdegree model.Thenumberofstatesnthereforecontrolstheduration
: [5,6]. Also, thereareabsencesofepidemicthresholdsinSF oftherefractoryperiod(whichcorrespondston−2,inms).
v
networks [7] and of the kinetic effects in reaction-diffusion Inthebiologicalcontext,rcouldberelatedforexamplewith
i
X processestakingplaceonSFnetworks[8]. theconcentrationofagivenodorantpresentedtoanolfactory
r In this paper, we will study the excitable Greenberg- epithelium[13],orthelightintensitystimulatingaretina[14].
a
Hastings cellular automaton (GHCA) model [9] on the SF Weshallrefertorasthestimulusrateorintensity.
network, especially we focus here on the Baraba´si and Al- BAgraphisakindofSFnetworksandcanbeconstructed
bert (BA) graph [4]. Due to experimental data which sug- according to ref. [4]. Starting from a small number m0
gest that some classes of spiking neurons in the first layers of nodes, every time step a new vertex is added, with m
ofsensorysystemsareelectricallycoupledvia gapjunctions links that are connected to an old node i with probability
or ephaptic interactions, the GHCA model is employed to Π =k / k ,wherek istheconnectivityoftheithnode.
i i j j i
model the response of the sensory network to external stim- After iteraPting this scheme a sufficient number of times, we
uliinsomerecentworks. Two-dimensionaldeterministiccel- obtainanetworkcomposedbyN nodeswithconnectivitydis-
lular automatonmodelwas studied by computersimulations tribution P(k) ∼ k−3 and average connectivity hki = 2m.
2
Forthiskindnetwork,theabsenceofacharacteristicscalefor
the connectivity makes highly connected nodes statistically

significant, and induces strong fluctuationsin the connectiv-
0.15
itydistributionwhichcannotbeneglected.
n = 4
Letρ (s)bethedensitiesofneuronswhichareinstatesat
t 5
timet. We havethenormalizationcondition, n−1ρ (s) = 0.1 6
s=0 t
1. SincethedynamicsoftherefractorystateisPdeterministic,
theequationsfors≥2aresimply
F
0.05
ρt+1(2) = ρt(1)
ρt+1(3) = ρt(2)
.
.
.
ρt+1(n−1) = ρt(n−2). (1)
4 5 6 7 8 9 10
Byimposingthestationaritycondition,wehave
1/p
ρ (0)=1−(n−1)ρ (1), (2)
t t
Following ideas developed by Pastor-Satorras et al. [7], to FIG.1: Intheabsenceofstimulus,theaverageactivityF asafunc-
takeintoaccountstrongfluctuationsintheconnectivitydistri- tionof1/pforBAnetworksofsizeN =104andm0 =m=4.The
bution,weconsidertherelativedensityρ k(1)ofactivenodes linearbehavioronthesemi-logarithmicscaleprovesanexponential
t
withgivenconnectivitykbythefollowingequation behaviorpredictedbyEq.(10).
∂ ρ k(1)=−ρ k(1)+λρ k(0)+kp(1−λ)ρ k(0)Θ(ρ (1)),
t t t t t t
(3) whichyieldsthesolution
where Θ(ρ (1)) is the probability that any given link point
t
to an active node and is assumed as a function of the total λ λ 1 maΘ(λ,p) maΘ(λ,p)
density of exciting nodes ρ (1). In the steady state, ρ (1) is Θ(λ,p)= +( − ) ln ,
t t b b n−1 b maΘ(λ,p)+b
just a functionof λ and p. Thus, the probabilityΘ becomes (8)
animplicitfunctionofλandp. Byimposingthestationarity where a = (n − 1)(1 − λ)p and b = 1 + (n − 1)λ. We
condition,∂tρtk(1)=0,weobtain can solve the above equation to obtain the solution Θ(λ,p).
Theoretically,combingEqs. (4) , (6) and(8), we can obtain
λ+(1−λ)pkΘ(λ,p)
ρk(1)= . (4) the final result of ρ(1). In some special cases the analytical
1+(n−1)[λ+(1−λ)pkΘ(λ,p)]
expressionscanbeobtained.
This set of equations show that the higher the node connec- i)noexternalstimulus,i.e.,λ=0orr =0. Wecanobtain
tivity,thehighertheprobabilitytobeinaspikingstate. This
inhomogeneity must be taken into account in the computa- e−1/mp
tion of Θ(λ,p). Indeed, the probability that a link points to Θ(0,p)= (1−e−1/mp)−1. (9)
(n−1)mp
anodewithq linksisproportionaltoqP(q). Inotherwords,
a randomlychosen link is more likely to be connectedto an Combing Eqs. (4), (6), and (9), we find the solution of the
excitingnodewithhighconnectivity,yieldingtherelation densityofactivenodeswhenthereisnoexternalstimulus,
kP(k)ρk(1)
Θ(λ,p)= qP(q) . (5) ρ(1)∼ 1 e−1/mp. (10)
Xk q n−1
P
Since ρk(1) is on its turn a functionof Θ(λ,p), we obtaina ii)the uncoupledcase, i.e., p = 0. CombingEqs. (4)and
self-consistency equation that allows to find Θ(λ,p) and an (6),wehave
explicitformfor Eq. (4). Finally, we can evaluate the order
parameter(persistence)ρ(1)usingtherelation λ
ρ(1)= , (11)
1+(n−1)λ
ρ(1)= P(k)ρk(1). (6)
X
k which is independentof the network’stopology,so it is also
For the BA graph, the full connectivity distribution is given have been obtained in other networks, such as one dimen-
by P(k) = 2m2k−3, where m is the minimum number of sional case [11]. When the external stimulus intensity r is
connection at each node. By noticing that the average con- verysmall,ρ(1)∼r−1.
∞
nectivityishki= kP(k)dk =2m,Eq.(5)gives Wedefinetheaverageactivity
m
∞R
1 λ+(1−λ)pkΘ(λ,p)
Θ(λ,p)=m , 1 T
Zm k21+(n−1)[λ+(1−λ)pkΘ(λ,p)] F = ρt(1), (12)
(7) T Xt=1
3

21
-1
10
20
-2
) 10 )
1 B 19
- ms -3 0.12 d
10
(
F( F00.08 18
-4 0.04
10
0.00
-5 0.0 0.1 p 0.2 0.3 17
10
-4 -3 -2 -1 0 1
10 10 10 10 10 10
0.00 0.02 0.04 0.06 0.08 0.10 0.12
-1
r(ms ) p
FIG.3:Dynamicrangevs.branchingprobabilityp.Thecurveisob-
FIG.2: Responsecurves(meanfiringratevs. stimulusrate). Points tainedbycalculatingthedatafromFig.2. Inthesubcriticalregime,
representsimulationresultswithN = 104,m0 = m = 4,n = 5 thedynamicrange∆(p)increasesmonotonicallywithp,Inthesu-
statesand T = 103 ms, from p = 0 top = 0.12 (inintervals of percriticalregime,∆(p)decreaseswhenpincreases.Thereisamax-
0.02).ThesecurvesarepowerlawsF ∝rαwithα=1(subcritical) imalrangepreciselyatthecriticalpointpc =0.06.
andα=1/2(critical).Inset:spontaneousactivityF0vs.branching
probabilityp,thecriticalpointispc=0.06.
[r0.1,r0.9] is found from its corresponding response interval
[F0.1,F0.9],whereFx =F0+x(Fmax−F0).
where T is a large time window (of the order of 104 time Figure 3 depicts the dynamic range ∆ versus the branch-
steps). Inthe stationarystate, itis obviouslythatF = ρ(1). ing probability p. In the subcritical regime, the dynamic
Toconfirmthepictureextractedfromtheaboveanalytictreat- range∆(p)increasesmonotonicallywithp, Inthesupercrit-
ment,weperformnumericalsimulationsontheBAnetwork. ical regime, ∆(p) decreases when p increases. There is a
Figure 1 showsthe behaviorof the averageactivity F in the maximalrangepreciselyatthecriticalpoint. Thisresultwas
absence of stimulus. All the plots decays with an exponent alsofoundinErdo¨s-Re´nyirandomgraph[12]. Kinouchiand
formF ∼ exp(−c/mp), where c is a constant. The numer- Copelliexplainedtheresultas“Inthesubcriticalregime,sen-
icalvalueobtainedc = 0.226isin goodagreementwith the sitivityisenlargedbecauseweakstimuliareamplifieddueto
theoreticalpredictionc=1/m=0.25. activity propagationamong neighbours. As a result, the dy-
InFig. 2,weshowtheaverageactivityF versustheexter- namic range ∆(p) increases monotonically with p. In the
nalstimulusrater fordifferentbranchingprobabilityp. The supercritical regime, the spontaneous activity F0 masks the
insetshowsthespontaneousactivityF0 versusthebranching presenceofweakstimuli,therefore∆(p)decreases.Theopti-
probabilitypintheabsenceofstimulus(r =0),thereisacrit- malregimeoccurspreciselyatthecriticalpoint”[12].Wealso
icalpointp = 0.06,onlyatp > p thereisaself-sustained performsimulationsonothercomplexnetworks,suchasWS
c c
activity, F > 0. The curves F(r) could be fitted by a Hill networksand randomSF networkswith differentγ, and ob-
function, but are not exactly Hill, and there exists the rela- tainthesamebehavioroftheaverageactivityF asitisshown
tion F ∝ rα for the low-stimulus response. It is found that inFig. 2,theoptimalregimeoccuringpreciselyatthecritical
theStevens-Hillexponentαchangesfromα = 1inthesub- point. We state thatthisphenomenonistheuniversalbehav-
criticalregimetoα = 0.5atcriticality. Thisimportantpoint ior of excitable GHCA model on complex networks and the
wasalsoreportedbyref. [12]wherethenetworkisanErdo¨s- explanationprovidedbyKinouchiandCopelliisalsosuitable
Re´nyirandomgraph. Wealsonoticethatapparentexponents forothertypesofcomplexnetworks.
between0.5and1.0areobserved[11]iffinitesizeeffectsare In SF networks, there exist strong fluctuations in the con-
present,thatis,ifN issmall. nectivitydistribution. Itisworthytoinvestigatethebehavior
As a function of the stimulus intensity r, networks have of the averageactivity Fk for nodeswith given connectivity
a minimum response F0 (= 0 for the subcritical and critical k.InFig. 4weplotthequantityFkversustheexternalstimu-
cases) and a maximum response F (due to the absolute lusintensityrfordifferentbranchingprobabilityp. Figs.4(a)
max
nature of the refractory period, F = 1/n, which can be and(b)correspondto the node’sconnectivityk = 4 and36,
max
obtainedbysettingλ=1inEqs. (11)). Ref. [12]definesthe respectively.ItisfoundthatthecurvesFk(r)havethesimilar
dynamicrange∆=10log(r0.9/r0.1)asthestimulusinterval propertyofF(r)whichwasshowninFig. 2,i.e.,theStevens-
(measuredindB)wherevariationsinrcanberobustlycoded Hillexponentαchangesfromα=1inthesubcriticalregime
by variationsin F, discardingstimuliwhichare too weak to toα=0.5atcriticality.
bedistinguishedfromF0ortooclosetosaturation.Therange Finally, we calculate the dynamic range ∆k(p) of nodes
4
with givenconnectivityk, and the result is shown in Fig. 5.
Thephenomenonthattheoptimalregimeoccurspreciselyat

thecriticalpointrecurs.Itisnotablethattheoptimaldynamic
-1 -1
10 10 range for nodes with given connectivity k increases with k,
i.e., for the node with larger connectivity, its corresponding
-2 -2 largeroptimal range. One can investigateevery node’sopti-
)10 10
1 maldynamicrangeandcalculatetheirfluctuations,since the
-
ms -3 -3 fluctuationsofnode’soptimaldynamicrangereflectthefluc-
10 10
( tuations of the node’s connectivityin the network. To some
k
F extent, we can discoverythe topologyof the network by in-
-4 -4
10 10 vestigatingthedynamicsonit.
(a) (b) In summary, we have investigated the behavior of the ex-
-4 -3 -2 -1 0 1 -4 -3 -2 -1 0 1 citableGHCAmodel[9]onBAnetworks. Wefoundoutthat
10 10 10 10 10 10 10 10 10 10 10 10
the curves of the average activity F as a function of the ex-
-1
r(ms ) ternalstimulusrater canbefittedbyaHillfunction,butare
not exactly Hill, and there exists a relation F ∝ rα for the
FIG.4: AverageactivityFk fornodeswithgivenconnectivitykvs. low-stimulusresponse. TheStevens-Hillexponentαchanges
external stimulus rate r. (a) for connectivity k = 4, and (b) for fromα =1inthesubcriticalregimetoα = 0.5atcriticality.
connectivity k = 36. Simulations are performed with N = 104, There is a maximal range precisely at the critical point. We
m0 = m = 4,n = 5statesandT = 104 ms,fromp = 0top = alsoobservedtheseresultsnumericallyinotherkindofcom-
0.12(inintervalsof0.02). ThesecurvesarepowerlawsFk ∝ rα plex networks. We conclude that these phenomena are the
withα=1(subcritical)andα=1/2(critical). universal behaviors of the excitable GHCA model on com-
plexnetworks. Duetostrongfluctuationsintheconnectivity
distributionontheBAgraph,wecalculatedtheaverageactiv-
ityFk(r)andthedynamicrange∆k(p)fornodeswithgiven

connectivityk. ThetwoquantitiesFk(r)and∆k(p)havethe
32
similar behavior as that of F(r) and ∆(p), respectively. It
30 k = 4 is interesting that nodes with larger connectivity have larger
8
28 optimal range. This property could be applied in biological
15
26 26 experimentrevealingthenetworktopology.
)
B 36
d 24
( 22
20
18
16
0.00 0.02 0.04 0.06 0.08 0.10 0.12

p
FIG.5: Dynamicrange∆k(p)ofnodeswithgivenconnectivityk This work was supported by the Fundamental Research
vs. branching probabilityp. Pointsrepresent theresultscalculated Fund for Physics and Mathematics of Lanzhou Univer-
fromsimulationwithN = 104 sites,m0 = m = 4,n = 5states sity under Grant No. Lzu05008. X.-J. Xu acknowl-
andT =104 ms. Dynamicrange∆k(p)isoptimizedatthecritical edges financial support from FCT (Portugal), Grant No.
pointpc =0.06,althoughtheconnectivitykisdifferent. SFRH/BPD/30425/2006.
[1] R.AlbertandA.-L.Baraba´si,Rev.Mod.Phys.74, 47(2002); 378(2000).
S.N. Dorogovtsev and J.F.F. Mendes, Adv. Phys. 51, 1079 [6] R.Cohen,K.Erez,D.ben-Avraham,andS.Havlin,Phys.Rev.
(2002);M.E.J.Newman,SIAMRev.45,167(2003). Lett.85,4626(2000).
[2] D.J.WattsandS.H.Strogatz,Nature393,440(1998). [7] R.Pastor-SatorrasandA.Vespignani,Phys.Rev.Lett.86,3200
[3] D.J.Watts,SmallWorlds: TheDynamicsofNetworksBetween (2001); Phys. Rev. E 63, 066117 (2001); ibid 65, 036104
OrderandRandomness (PrincetonUniversityPress,NewJer- (2002).
sey,1999). [8] L.K. Gallos and P. Argyrakis, Phys. Rev. Lett. 92, 138301
[4] A.-L.Baraba´siandR.Albert, Science286, 509(1999); A.-L. (2004).
Baraba´si,R.Albert,andH.Jeong,PhysicaA272,173(1999). [9] J.M.GreenbergandS.P.Hastings,SIAMJ.Appl.Math.34,515
[5] R.Albert,H,Jeong,andA.-L.Baraba´si,Nature(London)406, 1978.
5
[10] M.CopelliandO.Kinouchi,PhysicaA349,431(2005). BioSystems58,133(2000).
[11] L.S.FurtadoandM.Copelli,Phys.Rev.E73,011907(2006). [14] M.R. Deans, B. Volgyi, D.A. Goodenough, S.A. Bloomfield,
[12] O.KinouchiandM.Copelli,Nat.Phys.2,348(2006). andD.L.Paul,Neuron36,703(2002).
[13] J.-P.Rospars,P.La´nsky´,P.Duchamp-Viret,andA.Duchamp,