Table Of ContentEmergence of self-sustained patterns in small-world excitable media
Sitabhra Sinha1, Jari Sarama¨ki2, and Kimmo Kaski2
1The Institute of Mathematical Sciences, C.I.T. Campus, Taramani, Chennai - 600 113 India
2Laboratory of Computational Engineering, Helsinki University of Technology, P.O. Box 9203, FIN-02015 HUT, Finland
(Dated: February 4, 2008)
Motivated byrecent observations that long-range connections (LRC) play a role in variousbrain
7
phenomena, we haveobserved two distinct dynamical transitions in theactivity of excitable media
0
where waves propagate both between neighboring regions and through LRC. When LRC density p
0
is low, single or multiple spiral waves are seen to emerge and cover theentire system. This state is
n 2 saeltfr-asnussittaiionnintgoaansdpraotbiaullsyt hagoaminosgtenpeeorutus,rbtaemtiopnosr.alAlytppe=riopdlcicthsteastpe.iraFlisnaarlley,suapbpovreessped=apnucd, athcteirveitiys
a ceases after a brief transient.
J
6 PACSnumbers: 05.65.+b,87.18.Bb,05.45.Jn,89.75.Hc
]
n Patternformationinexcitablemediahasrecentlybeen points towards intriguing possibilities for the functional
n
a very exciting area of research, not least because it is role of these connections. In this Letter we have investi-
-
s observed in a wide variety of natural systems, ranging gatedageneric modelof excitable media with increasing
i
d from spiral waves in mammalian brain [1] to chemical densityofrandomlong-rangeconnections,andwereport
. systems such as the Belousov-Zhabotinsky reaction [2]. theexistenceoftwoqualitativelydifferentregimesofself-
t
a Thefunctionalroleofsuchpatternsinbiologicalsystems sustained pattern formation. The correspondence of the
m
makes it imperative to understand better the conditions observed patterns with those observed in nature, e.g.,
- in which they can spontaneously emerge. For example, epileptic bursts and seizures [12], as well as their depen-
d
thesepatternshavebeenimplicatedinthegenesisoflife- dence onthe topologicalstructure of connections,makes
n
o threatening arrhythmias in the heart [3], while in the these results highly relevant in our view.
c brainthey havebeenthoughtto provideaspatialframe- The model we consider here consists of a two-
[ work for cortical oscillations [1]. Until now, work has dimensionalarrayofN×N excitable cellscoupleddiffu-
1 mostly focussed on generating patterns in excitable me- sively,
v diabyusingstochasticstimulation[4]. Inthiscase,how-
D
21 ever,theactivityisnotreallyself-sustainedasthenoiseis xit,+j1 =(1−D)f(xit,j,yti,j)+ 4 X f(xit+q,j+q,yti+q,j+q).
akin to external intervention necessary for the initiation q=±1
1
and persistence of spiral waves [5, 6]. If such patterns
1 Here D is the diffusion coefficient and the dynamics of
0 are to be seen as spontaneously emerging from arbitrary individualcellsaredescribedbyapairofvariablesx ,y ,
t t
7 initial conditions, then the pattern formation should be
evolving according to a discrete-time model of generic
0 anoutcomeoftheinternalstructureofthesystemalone.
excitable media [13]:
/
t Furthermore, small variations in this structure may re-
a sult in transitions between different dynamical regimes x = f(x ,y )=x2e(yt−xt)+k,
m t+1 t t t
characterized by distinct spatiotemporal patterns. y = g(x ,y )=ay −bx +c,
- t+1 t t t t
d
Here we have consideredexcitable systems which have where we have fixed the parameters as a = 0.89,b =
n
o a regular topology with cells communicating only with 0.6,c = 0.28, and k = 0.02. When this system is ex-
c nearest neighbors, but where there are a few random citedwithasuprathresholdstimulation,thefastvariable
:
v long-range connections, linking spatially distant cells. x shows anabrupt increase. This triggerschanges in the
i Such “small-world” topologies have been observed in a slow variable y, such that the state of the cell is grad-
X
large number of real world systems [7], and have also ually brought down to that of the resting state. Once
r
a been associatedwith self-sustained activity in a chain of excited,the cellremainsimpervious to stimulationupto
modelneurons[8]aswellasperiodicepidemicpatternsin arefractoryperiod, the durationofwhichis governedby
disease spreading [9]. However, while it has been shown the parameter a. Neighboring cells communicate exci-
that such long-rangeconnections do play a role in main- tation to each other with a strength proportional to D,
taining spiral waves that already exist [10], so far there chosen D =0.2 for most of our simulations. In addition
havebeenveryfewattemptsatshowingtheemergenceof, to the diffusive coupling, we introduce long-range con-
and transitions between, different types of self-sustained nections such that each cell receives a connection from
patterns as a resultof suchtopology. The observationof a randomly chosen cell with probability p. The strength
spontaneous pattern formationin naturalsystems where of this connection is chosen to be the same as that of
sparse long-range connections coexist with fairly regu- the nearest neighbors, i.e., D/4. These random long-
larunderlyingconnectiontopology,likeinthebrain[11], range links can be either quenched (i.e., chosen initially
2
riodicity disappears as the system settles into the spiral
wave mode. For the time series preceding the onset of
the spiral wave, we typically observe slow modulations
of the envelope of the periodic oscillations, which arise
from the interaction between the waves as well as from
the shortcut-induced excitations. Power spectral densi-
tiesofsuchserieswereobservedtoshowapower-lawlike
decay, indicating the presence of 1/f-noise.
Thespiralwaveswereobservedtobeprimarilycreated
by a shortcut-induced excitation occurring in the refrac-
tory “shadow” of a circular wave front, sparking a semi-
circularwavewhosetransmissionispartiallyhinderedby
the shadow. We verified that indeed spiral waves could
be triggered in this fashion by using externally applied
FIG. 1: Emergence of spiral waves, recorded in a simulation signal to stimulate an appropriate point in this region
of a N2 = 128×128 system, with shortcut probability p = andthen observingthe resultantpatterns. Once a spiral
0.25. After initialization, the dynamics of the system can be wave is created, it will eventually take over the dynam-
characterized by circular waves. At around t ∼ 1,500 time
ics of the system, as the successive excitation wavefronts
steps,aspiralwaveisspontaneouslycreatedandsubsequently
occur with the highest frequency compared to all other
takesoverthedynamics. Panels(a)and(b)showthestateof
excitations which will then be swept away. Evidently,
the system at different times in terms of the fast variable x.
Colorsindicateexcitationlevelofthecells. (c)Timeseriesof the probability of a spiral wave creation per unit time
theaverage activity,i.e. the fraction of cells with x>0.9. increaseswiththesystemsizeN2 andtheshortcutprob-
abilityp. Thusforlongtimesandlargesystemsizes,exci-
tationviarandomshortcutswilleventuallyalwaysresult
andkeptfixedforthe durationofthe simulation),oran- inthe formationofspiralwaves. This is corroboratedby
nealed (i.e., randomly created at each time step). While Fig. 2, where the fraction of spiral wave configurations
theresultsreportedbelowareforannealedrandomlinks, in 400simulationruns is shownas function of time t, for
we could observe no qualitative difference between these varying system sizes at p = 0.05 (panel a) and for fixed
two cases. Note that we have used both periodic and system size but increasing values of p (panel b). The
absorbing boundary conditions for the system, and ob- spatial structure of the spiral waves, i.e. an excitation
served no significant difference in the results. All the front followed by a refractory shadow, makes them very
following results were obtained for absorbing boundary robust against perturbations. For example, if the state
conditions. of the system is frozen, and the state of large areas (say,
In the simulations, the observed patterns formed in upto a quarter of the total area)or every second cell are
the system for varying values of the shortcut density p set to x = 0, the spiral wave mode is quickly recovered
can,generallyspeaking,be divided into threecategories. once the simulation is restarted.
First,belowthelowercriticalprobability,i.e.,0<p<pl, However, at high enough values of p, the shortcut-
c
the stateofthe systemafter aninitialtransientperiodis induced excitations become too numerous for sustaining
characterized by self-sustaining single or multiple spiral thespiralwavedynamics. Asalmosteverypointisliable
waves covering the entire system. Second, at p=pl, the to be excited with a frequency proportionalto its refrac-
c
spiral wave mode is suppressed and the system under- tory period, the spirals become unstable and we see a
goes a transition to a regime in which a large fraction of transition to a new regime at p = pl ≈ 0.553 (Fig. 3).
c
the system gets simultaneously active, and then refrac-
tory, in a periodic manner. Third, when the value of p
a) b)
isincreasedaboveasystem-size-dependentuppercritical
1 1
probability pu, the self-sustained activity ceases and the N=100
c 0.8 N=200 0.8
systemfallsintotheabsorbingstatewherexi,j =0,∀i,j. N=300
0.6 N=400 0.6
Inthebeginningofeachsimulationrunthesystemwas FS p=0.1
0.4 0.4 p=0.2
initialized such that x = 1 for a small number of cells, p=0.3
0.2 0.2
and x = 0 for the rest. For small shortcut probabilities p=0.4
0 0
p ≪ pl, upon starting the simulation, multiple coexist- 0 5000 10000 15000 20000 0 5000 10000
c time t time t
ing circular excitation waves were seen to emerge (see
Fig.1a), to be later takenoverby spiralwaves(Fig.1b). FIG.2: Fraction FS of systemswith spiralwavesin 400runs
TheactivitytimeseriesinFig.1c,displayingthefraction as function of time t: (a) fixed shortcut probability p=0.05
of cells where x > 0.9, shows a high frequency periodic- withvaryingsystemsize,(b)fixedsystemsizeN2 =300×300
ity corresponding to the refractory period; then, the pe- with varyingshortcut probability p.
3
1 such that after a transient, almost all cells become re-
0.8
0.8 fractory and not enough susceptible cells are left to sus-
0.7
tainthe excitation. We observedthatonce a system-size
FS0.6 lpC00..65 dependent value p=puc (N)is approached,the probabil-
0.4 N=100 ity of reaching the absorbing state xi,j =0,∀i,j grows
N=200 0.4 (cid:0) (cid:1)
0.2 NN==340000 0.15 0D.20 0.25 rapidly(see Fig.5) suchthat athighenoughvalues ofp,
the system is never seen to sustain its activity. The lim-
0
-2 0 2 4 iting behavior for N → ∞ was investigated by plotting
(p-pl)N
c pu against1/N and extrapolatingits value for 1/N →0.
c
Fitting a quadratic function yielded pu →0.86, hence it
FIG.3: ThefractionofspiralconfigurationsFS asafunction c
of `p−plc´ normalized bysystem size N, for shortcut proba- appears that the regime where activity is sustained has
bilities p around plc ≈ 0.553. FS is calculated at t = 20,000 for all system sizes an upper limit for the shortcut den-
timesteps,averagedover400runs. Inset: Dependenceofthe sityp. Notethatthiscanbeviewedasanapproximation
critical valueplc ondiffusion constantD. Thecirclesaresim- only, because there is no a priori reason to assume any
ulationresultsforN =200whilethecurveshowsfittingwith particular form for the dependence of pu on N.
plc ∼D−1.14. We have also carried out simulationscwith disorder in
the parameters describing the properties of individual
cells. For instance, we have made the parameter a that
controls the refractory period a quenched random vari-
able ranging over a small interval. We find that there is
a tendency for greater fragmentation of waves with dis-
order,andcorrespondingincreaseinthenumberofcoex-
istingspiralwaves,butotherwisenoremarkablechanges.
The robustness of our results in the presence of disorder
in the individual cellular properties underlines the rel-
evance of this study to real-world systems, where cells
are unlikely to have uniform properties. In addition, we
havelookedattheeffectofthediffusionconstantD. For
higher D the wave propagates much faster, so that, in a
fixed amount of time, the system will initiate excitation
FIG. 4: The periodic regime, recorded in a simulation of a throughmanymorelong-rangeconnectionsthanthesys-
N2 = 128×128 system, with shortcut probability p = 0.6. tem with lower D. Therefore, the higher D system will
a)-d): Snapshotsofthestateof thesystemtakenat intervals
be equivalent to a lower D system with a larger number
of ∆t=10 steps. e): Time series of the activity (see Fig. 1)
of long-range connections, i.e., higher p (Fig. 3, inset).
and (f): the corresponding power spectrum. The main peak
is at f0 ≈0.022, corresponding to a wavelength of ≈43 time The above results carry potential relevance to all nat-
steps; other peaks are harmonics at integer multiples of f0. ural systems that are excitable and have sparse, long-
range connections. Experiments carried out in the ex-
citable BZ reaction system with nonlocal coupling [14]
This value of pc was found to be independent of the sys- haveexhibitedsomeofthe featuresdescribedabove. For
tem size N2. Here, the spatial pattern becomes more example, the coexistence of transient activity and sus-
homogeneous, as a large fraction of the system becomes tained synchronized oscillations that we observe close to
simultaneously active and subsequently decays to a re- pu has been reported in the above system, and can now
c
fractory state. However, some cells not participating in be understood in terms of the model introduced here.
this wave of excitation carry on the activity to the next Recently,therehasalsobeenalotofresearchactivityon
cycle, where it again spreads through almost the entire theroleofnon-trivialnetworktopologyonbrainfunction
system. This results in a remarkably periodic behavior (e.g., Ref. [15]). In particular, studies show a connec-
of the system in time, with a large fraction of cells be- tion between the existence of small-world topology and
ing recurrently active with a period close to the refrac- epilepsy [12, 16]. This is one of the areas where our re-
tory period of the cells (see Fig. 4 (e)-(f)). Often, small sults can have a possible explanatory role. In the brain,
spiral-like waves were also observed (Fig. 4 (a)-(d)), but glial cells form a matrix of regular topology with cells
these were short-lived and spatial correlations could not communicating between their nearest neighbors through
be maintained. calcium waves. Neurons are embedded on this regular
Finally,whenpisincreasedstillfurther,theverylarge structure and are capable of creating long-range links
numberofshortcutconnectionsguaranteesalmostsimul- between spatially distant regions. It is now known that
taneous spread of excitation to nearly all cells. As a re- neurons and glialcells can communicate with each other
sult, the dynamics of the system tends to “burn out” through calcium waves [17]. Therefore, the aggregate
4
1 that is driven by the system’s own internal architecture.
1
Most important of all, the system exhibits a non-trivial
0.8
0.9 N=100 transition point at which the pattern goes from the spa-
C0.6 N=150
F 0.4 N=200 tial domain, with multiple coexisting spiral waves where
N=300
pc,u0.8 0.2 NN==345000 the global activity level remains more or less uniform,
0-0.1 0 0.1 to the temporal domain, where the global activity level
(p-p )
0.7 c,u shows large oscillations as a large fraction of cells be-
comes simultaneously active and then refractory, with a
0.6 strict periodicity. The connection to biological phenom-
0 0.005 0.01
1/N ena, most importantly, to calcium waves in the brain,
and the possibility of the functional role of small-world
FIG. 5: Inset: Fraction of configurations FC where activity connections in epileptic seizures and bursts, is expected
has ceased at t = 20,000 time steps in 100 simulation runs
as function of (p−puc). Main: The upper critical shortcut to make our study of special relevance.
probability puc as function of inverse system dimension 1/N
(◦). Thesolid linedisplaysafittedquadraticfunction,where
puc (N →∞)≈0.86.
[1] X. Huang, W. C. Troy, Q. Yang, H. Ma, C. R. Laing, S.
J. Schiff and J.-Y. Wu,J. Neurosci. 24, 9897 (2004).
system of neurons and glial cells can be seen as a small-
[2] A. N. Zaikin and A. M. Zhabotinsky, Nature (London)
world network of excitable cells. Here the neurons, that 225, 535 (1970).
are outnumbered by glial cells approximately by one or-
[3] R. A. Gray, A. M. Pertsov, and J. Jalife, Nature (Lon-
der of magnitude, form the sparse, long-range connec- don) 392, 75 (1998); F. X. Witkowski, L. J. Leon, P. A.
tions. Recently, the role of neural-glial communication Penkoske, W. R. Giles, M. L. Spano, W. L. Ditto, and
in epilepsy has been investigated (e.g., Ref. [18]), but A. T. Winfree, Nature(London) 392, 78 (1998).
the observation of spiral intercellular calcium waves in [4] J. Garc´ıa-Ojalvo and L. Schimansky-Geier, Europhys.
Lett. 47, 298 (1999).
the hippocampus [11], and the similarity of the patterns
[5] P. Jung and G. Mayer-Kress, Chaos 5, 458 (1995).
seen in our model with the observedfeatures of epileptic [6] Z. Hou and H.Xin, Phys.Rev. Lett. 89, 280601 (2002).
seizuresandbursts,makesitespeciallyappealingtopos- [7] D. J. Watts and S. H. Strogatz, Nature (London), 393,
tulate the role of small-world topology in the generation 440 (1998).
of epilepsy. Experimental verification of this suggested [8] A.Roxin,H.RieckeandS.A.Solla,Phys.Rev.Lett.92
scenario can be performed through calcium imaging in (2004). See also, S. Yonker and R. Wackerbauer, Phys.
glial-neuronalco-culture systems. Rev. E73 026218 (2006).
[9] M. Kuperman and G. Abramson, Phys. Rev. Lett 86,
The relationship between the fraction of long-range
2909 (2001).
connections and normal brain function is indicated by
[10] D.He,G.Hu,M.Zhan,W.RenandZ.Gao, Phys.Rev.
experimental evidence not just for epilepsy but for var- E 65, 055204(R) (2002).
ious other mental disorders as well. A recent study has [11] M. E. Harris-White, S. A. Zanotti, S. A. Frautschy and
foundthatthebrainsofpatientswithclinicallydiagnosed A. C. Charles, J. Neurophysiol. 79, 1045 (1998).
schizophrenia, depression or bipolar disease have lower [12] T. I. Netoff, R. Clewley, S. Arno, T. Keck and J. H.
glia to neuron ratio compared to normal subjects [19]. White, J. Neurosci. 24, 8075 (2004).
[13] D. R. Chialvo, Chaos, Solitons and Fractals 5, 461
Therefore, understanding the dynamical ramifications of
(1995).
increasing shortcuts in an excitable media can motivate
[14] M.Tinsley,J.Cui,F.V.Chirila,A.Taylor,S.Zhongand
experimentsthathavethepotentialsignificanceofaiding K. Showalter, Phys. Rev. Lett. 95, 038306 (2005); A. J.
clinical breakthroughs in treating a whole class of men- Steele, M. Tinsley and K. Showalter, Chaos 16, 015110
taldisorders. This is alsoconnectedwith the questionof (2006).
theevolutionarysignificanceforincreasinggliatoneuron [15] V. M. Egu´ıluz, D. R. Chialvo, G. A. Cecchi, M. Baliki
ratio with brain size [20]. It is known that excitable me- andA.V.Apkarian,Phys.Rev.Lett.94,018102(2005).
dia of larger dimensions are more likely to exhibit spiral [16] B. Percha, R. Dzakpasu and M. Zˆochoswski, Phys. Rev.
E 72, 031909 (2005).
waves[21]. Therefore,decreasingthe fractionof neurons
[17] A. C. Charles, Dev.Neurosci. 16, 196 (1994).
(and therefore, long-range shortcuts) could be Nature’s [18] S. Nadkarni and P. Jung, Phys. Rev. Lett. 91, 268101
way of ensuring dynamical stability for neural activity.
(2003).
Toconclude,wehaveinvestigatedspontaneouspattern [19] R. A. Brauch, M. A. El-Masri, J. C. Parker and R. S.
formation in excitable media with small-world connec- El-Mallakh, J. Affect.Disord. 91, 87 (2006).
tions. Ourresultsshowthecreationofnon-trivialspatio- [20] A. Reichenbach,Glia 2, 71 (1989).
temporal patterns similar to those seen in many real-life [21] S. Sinha, A. Pande and R. Pandit, Phys. Rev. Lett. 86,
3678 (2001).
systems. These patterns are formed through dynamics
