Spontaneous formation of synchronization clusters in homogenous neuronal ensembles induced by noise and interaction delays

0.3 MB·English

by  Igor Franovic

Spontaneous formation of synchronization clusters in homogenous neuronal ensembles
induced by noise and interaction delays
Igor Franovi´c,1 Kristina Todorovi´c,2 Nebojˇsa Vasovi´c,3 and Nikola Buri´c4,∗
1Faculty of Physics, University of Belgrade, PO Box 44, 11001 Belgrade, Serbia
2Department of Physics and Mathematics,Faculty of Pharmacy,
University of Belgrade, Vojvode Stepe 450, Belgrade, Serbia
3Department of Applied Mathematics, Faculty of Mining and Geology,
University of Belgrade, PO Box 162, Belgrade, Serbia
4Scientific Computing Lab., Institute of Physics,
University of Beograd, PO Box 68, 11080 Beograd-Zemun, Serbia
(Dated: January 31, 2012)
2
Spontaneous formation of clusters of synchronized spiking in a structureless ensemble of equal
1
stochastically perturbedexcitableneuronswithdelayedcouplingisdemonstratedforthefirsttime.
0
Theeffectisaconsequenceofasubtleinterplaybetweeninteraction delays,noiseandtheexcitable
2
character of a single neuron. Dependence of the cluster properties on the time-lag, noise intensity

n and thesynaptic strength is investigated.
a
J

0 Collective behavior in large ensembles of physiological where the exact system shows an onset of cluster states.
3 and inorganicsystems canbe reducedto that of coupled Network dynamics and the tools to analyze it– We fo-

oscillators engaged in various synchronization phenom- cus on an N-size population of all-to-all diffusively cou-
]
O ena. In terms of macroscopic coherent rhythms, it may pledFitzhugh-Nagumoneurons,whose localdynamicsis
A either be the casewhereallthe units arerecruitedinto a set by
giant component or the case of cluster states character-
.
n izedbyexactorin-phaseintra-subsetandlaginter-subset c N
nli synchronization. The spontaneousonsetofclusterstates ǫdxi =(xi−x3i/3−y+I)dt+ N X[xj(t−τ)−xi(t)]dt,
is of particular interest to neuroscience [1] for the con- j=1
[
jectured role in information encoding, as well as for par- dy =(x+b)dt+√2DdW , (1)
i i
1 ticipating in motor coordination or accompanying some
v
neurological disorders. The approach to clustering has where the activator variables x embody the membrane
9 i
mostly relied on modeling neurons as autonomous oscil- potentials, while the recovery variables y mimic the ac-
9 i
1 lators, treating separately the question of whether the tion of the K+ membrane gating channels. c denotes
6 proposedmechanisms maybe robustunder noise [2] and the synaptic strength and τ stands for the coupling de-
1. transmission delays [3]. We explore a new mechanism lay, both parameters for simplicity assumed homoge-
0 which rests on the excitable character of neuronal dy- neous across the ensemble. The √2DdWi terms rep-
2 namics and mutual adjustment between noise and time resent stochastic increments of the independent Wiener
1 delay to yield the self-organization into functional mod- processes, i.e. the white noise. For the external stimula-
:
v ules within an otherwise unstructured network. tion holds I = 0, whereas the small parameter ǫ = 0.01
Xi For the instantaneous couplings, the research on pop- warrants a clear separation between the fast and slow
ulations of excitable neurons has covered pattern forma- time scales. Selecting b = 1.05, the neurons are poised
r
a tiondueto localinhomogeneity[4],orhasinvokedasce- near the Hopf bifurcation threshold b = 1, which places
nario where noise enacts a control parameter tuning the them in an excitable regime where each possesses a sin-
dynamicsofensembleaveragesbetweenthethreegeneric gle equilibrium. An adequate stimulation, be it by the
globalregimes[5]. Distinctfromthelayoutwithcomplex noise or the interaction term, may evoke a large excur-
connection topologies, here it is demonstrated how cou- sion of membrane potential, passing through the spiking
pling delays do alter the latter landscape in a significant and refractory states before it loops back to rest.
fashion, giving rise to an effect one may dub the clus- To characterize the degree of correlation be-
ter forming time-delay-induced coherence resonance. In tween the firing events, we use primarily the
part,thestrategytoanalyzeglobaldynamicsrestsonde- interneuron spike train coherence [6] κij =
m m
riving the mean-field (MF) approximation for the exact X (k)X (k)/ X (k)X (k). This requires
Pk=1 i j pPk=1 i j
system. ThelikelygainfromtheMFtreatmentisatleast one to split the simulation period T into bins k of
twofold: except for allowing one to extrapolate what oc- length ∆ = T/m, awarding each neuron a variable
cursinthethermodynamiclimitN ,itmayserveas X (k) = 1(0), contingent on whether a spike was
i
→∞
anauxiliarymeanstodiscriminatebetweentheeffectsof triggered or not within the given time bin, respectively.
noise and time delay. Unexpectedly, the MF model un- As with all the quantities below, we have been careful
dergoes a global bifurcation at certain parameter values to exclude fromcalculations the transientbehavior. The
2
spike threshold and the time bin are set to X = 1
0 c = 0.08 c = 0.08
and ∆ = 0.008, verifying that no change of the results 0.5 (a) cc == 0 0.1.12 1.0 (b) cc == 0 0.1.12
occurred if X or ∆ were reduced. The distribution
0
0.8
of the κ values may serve to distinguish between the 0.4
ij
homogeneous and clustered network states. Another 0.6
0.3
aspect we are interested in is whether the clustered
0.4
states are monostable or coexistent with the homoge- 0.2
neous ones at the given network size. To probe this, 0.2
0.1
we have monitored if the values of the global coherence 0.001 0.002 0.003 0.004 0.005 0.006 0.001 0.002 0.003 0.004 0.005 0.006
D D
κ = 1 N κ for different realizations at the 1.0 (c) cc == 0 0.0.18 1.0 (d) cc == 0 .00.81
N(N−1) P ij c = 0.12 c = 0.12
i,j=1;i6=j 0.8 0.8
fixed parameters clumped together, expecting bunching
into distinct groups as evidence of multistable behavior. 0.6 0.6
Addressing the temporal structure of the network 0.4 0.4
states, it is useful to look into the distribution of the
0.2 0.2
local neuron jitters r [7]. They represent the normal-
i
0.001 0.002 0.003 0.004 0.005 0.006 0.001 0.002 0.003 0.004 0.005 0.006
ized variations of the interspike intervals Tk extracted D D
from x (t), r = <T2 > <T >2/ < T >, with
i i p k − k k FIG.1. (color online). Profiles of theκ(D) families ofcurves
smallervaluesindicatingbetterregularity. Themodality
over the synaptic strengths c = 0.08,0.1 and 0.12 display
and the width of the ri distribution over the population strongdependenceonthedelay,increasingfromτ =2in(a),
may serve as rough indices on how the cluster dynamics τ =6 in (b), τ =10 in (c) to τ =11 in (d). The location of
is mutually adjusted. In the final part, we analyze the ”wells” may point tothe emergence of theclustered states.
N
behaviorofthe ensembleaveragesX =1/N x and
Pi=1 i in Figs. 1(b) and 1(c) may occur for just two reasons,
N
Y = 1/N y , where the former increases if a larger
Pi=1 i as κ decreases either for the incoherent or the clustered
fraction of neurons fire in synchrony. The results for the
states. The latter alternative is supported by the co-
exact system are compared to those of the approximate
herence matrices in Fig. 3, which are discussed shortly.
MF model [9]. The latter presents a two-dimensionalset
The importance of the D - τ adjustment for the cluster-
of delayed differential equations
ing effect is also witnessed by the c-dependence within
the families in Fig. 1: the stronger the interaction term,
dX(t)
ǫ =X(t) X(t)3/3 the more salient is the picture of ”irregularity” sections
dt − −
immersed into a ”regular” curve profile. Increasing the
X(t)
1 c X(t)2+ [c 1+X(t)2]2+4D delay, the cluster states first occur, apparently monos-
− 2 n − − p − o
table, around τ =2 for the small D =0.00025, whereby
Y(t)+c[X(t τ) X(t)],
the typicalphase portrait(PP)projectionshowstwisted
− − −
dY(t) orbits with two clearly discernable segments, see Fig.
=X(t)+b, (2)
dt 2(a). These reflect the two macroscopic fractions of
the population firing alternately, such that the homo-
derived within a cumulant approach by employing the
geneous network dynamically splits into clusters of mu-
Gaussian approximation.
tually synchronized neurons, with the clusters locked in
We note that the results for the exact system refer to
antiphase. Thefrequencyentrainmentisindicatedbythe
a network of N = 200 neurons, applying independently
shape of the r distribution, which peaks sharply around
i
a method from[8] to verify no qualitative changesin the
< r > = 0.01. We tested the invariance of clustering
m
clustering behavior for larger N.
with N via the asymptotical behavior of the quantity
Results– To get a sense of what may be the param- χ2 = σX2 , where σ2 =< X(t)2 > < X(t) >2
eter ranges to admit the cluster states, we plot the c– N N1 PNi=1σx2i X t − t
families of the κ curves in dependence of D for different and σ2 =< x (t)2 > < x (t) >2 holds. If the clus-
xi i t − i t
τ. Without the delay, the curves would conform to a ter states endure, there should be a residual component
stereotype profile, where one distinguishes between the χ( ) (0,1) in the large N limit [8]. For this and
∞ ∈
three”regular”segmentsforverysmall,intermediateand the remaining cases, the onset of such a regime is found
large D, showing first a reduced κ due to incoherent os- around N 200, implying that no qualitatively novel
≈
cillations, then steady high values for the coherent ones phenomena occur above this system size. An interesting
and the decaying segment at D where the stochastic dy- observation is that the cluster configuration N ,N ,
1 2
{ }
namics prevails. However, from Fig. 1 we learn how determined by the fractions’ sizes, fluctuates around the
this is upheld for some τ, say τ = 11, but is violated ratio 2 : 1 for different stochastic realizations and ap-
manifestly at the ”cluster-resonant” values τ = 2,6,10. pears to aggregatewith enhancing N. For certain τ, the
The ”wells” seen at approximately D (0.001,0.003) 2-cluster state also emerges outside the D-region delim-
∈
3
FIG.2. GlobalPPsforthe2-clusterstatesshowtwistedLCs,
whereby the two discernable segments reflect the alternate
firingoftheneuronsubsets. TheN1/N2 ratiodependsonthe FIG. 3. (color online). Rearranged coherence matrices for
interplay of D and τ, as seen from the examples τ = 2,D = τ =5,D=0.0005,c=0.1 in (a) and τ =10,D=0.0013,c=
0.00025,c=0.1 in (a) and τ =5,D=0.0005,c=0.1 in (b). 0.1in (b)implythestrongclusterseparation inthe2-cluster
states and mixing between the clusters in the 3-cluster case.
iting the incoherent and coherent global regimes. This
Darkershading reflects higher coherence.
holdsforτ =5andD (0.0004,0.0008),wherethe clus-
nature, such that the neurons once engaged in synchro-
∈
ter layout is also such that if one is active, the other re-
nized spiking are much more likely to do so again.
mainsrefractory. Ther distributionmaintainsanarrow
i Aa understanding of clustering mechanism is revealed
form, but its maximum shifts to <r > 0.19. Though
m by comparing the typical PPs of neurons participant in
≈
oneretrievesthegeneralpicturefromabove,avarianceis
thehomogeneouscoherentstateandtheclusteredstates,
that larger τ seems to favorthe partition N /N 1:1,
1 2 see Figs. 4(a) and 4(b). A striking feature in the latter
≈
see PP in Fig. 2(b). The 1:1 ratio is preferred both for
case is a kink at the refractory branch of the slow mani-
increasing N and if the delay is set to τ =6.
fold. The appearance of a kink is the key manifestation
Theclusteredstatessofarmaybecastasstationaryin of the D τ co-effect, that consists in separating the
−
the sense of stability against neurons switching between ensemble into clusters and maintaining the proper phase
the clusters. We also report on the existence of 3-cluster difference between them. The purpose of the kink is to
statesthatmaybeconsidered”dynamical”,withtheneu- keeptheneuronsfrustratedlongenoughattherefractory
ronsabletojumptoandfromclusters. Suchanoutcome branchbeforebeingallowedtoslidedowntoitsleftknee.
arisesforthestrongernoiseD 0.0013,oncethedelayis Thismaybeimaginedasaformofalock-and-releasebe-
≈
increasedto τ =10. To underline the difference between havior, where the delay primarily gives rise to the first,
the stationary and dynamical clustered states at τ = 5 and noise to the second part. If a fraction of the en-
and τ =10, we plot side-by-side the corresponding pair- semble were to move beyond the left knee and the other
wise coherence matrices κ , see Figs. 3(a) and 3(b), were to lag behind, the split should be amplified with
ij
{ }
wherethenetworknodeshavebeenrearrangedbyahier- each population cycle, eventually becoming resilient to
archicalclusteringalgorithmaccordingto aformofmet- perturbation precisely due to trapping at the refractory
ricdistancethathasthemostcoherentnodestheclosest. branch. For trapping to be successful, the kink has to
This makes it explicit how the intercluster coherence for be placed properly, approximately where the dynamics
the 2-cluster state is virtually negligible with respect to of the representativepoint is most susceptible to pertur-
the 3-cluster case. Loosely speaking, within an unstable bation along the slow manifold. Then, for a brief pe-
three-part population division, when a certain fraction riod, due to an influence from x , the evolution of y is
i i
is firing, the other is refractory and the neurons in the locally accelerated, becoming comparable to a speed of
smallest cluster are at rest (excitable). This less clear change in the direction orthogonal to the slow manifold,
separation is also apparent when comparing the nodal driven by the spiking fraction of the population. Note
degree distributions in cases τ =5 and τ =10, obtained that the trapping interval has to be adjusted so that the
if one assumes κ to provide weights for the network entirepopulationisentrainedtoasinglefrequencyoffir-
ij
{ }
whose links stand for the correlated dynamics between ing. The latter matchesthe oneindelay-freecase,which
the neurons. For τ = 5, the bimodal degree distribution warrantsstability againstperturbations. The arguments
is clearly seen without raising the connectivity thresh- above and the numerical data seem to indicate how the
old,whereasforτ =10theinitiallysmearedthree-modal delayswherethecoherenceresonanceisfeltthestrongest
distributionrefinesaftersomethresholdingisperformed. maybeapproximatedbyτ =T /2+n T ,withT being
0 0 0
∗
The rationale of dynamical clustering may best be un- the period of coherent oscillations at τ = 0. Noise-wise,
derstoodbyanalyzingther distributioninthe 3-cluster with increasing τ, D has to be adjusted to higher values
i
state. Apartfrombeingwiderthaninthe2-clusterstate, toregulatetherelaxationfromthekinktotheslowman-
it peaks at a much smaller value < r > 0.09, im- ifold while maintaining the entrainment to the proper
m
≈
plying the more regular neuron firing. For this to hold, frequency. In parallel, for stronger D, the representa-
synchrony within the clusters has to be of intermittent tion cloud of the firing fraction tends to disperse more,
4
FIG. 4. (a) and (b) show typical PPs of neurons par- FIG.5. (coloronline). BistabilityintheMFmodel: (a)shows
ticipating the homogeneous global oscillations and clustered the trajectories converging either to the FP or the LC, de-
states, respectively. The latter are distinguished by a kink pendingontheinitialconditions,whereasin(b)aredisplayed
K, which is a signature of the D − τ co-effect. The pa- thetwo basins of attraction for τ =2,D=0.00025,c=0.1.
rameter sets are τ = 6,D = 0.0005,c = 0.1 in (a) and possibleonlyasaninterplayofexcitability,noiseandin-
τ =2,D=0.00025,c=0.1 in (b).
teraction delay. Once the phenomenon is recognized as
requiring a sufficient τ for this effect to be averagedout. causedonlybythesequalitativepropertiesonecanstudy
The interplay between D and τ is further highlighted the effects of more realistic assumptions on the distribu-
by exploring the behavior of the MF model (2). Local tion of neuronal properties and connection patterns. An
bifurcation analysis shows that the MF exhibits a suc- interesting point concerns the derived MF model, which
cession of super- and subthreshold Hopf bifurcations [9], can aid in understanding the precise roles played by D
which account for the transition from the stochastically and τ. Notably, beneath the surface lies a more strati-
stable FP to the stable LC. Still, this scenario is con- fied phenomenon, where the subtle adjustment between
fined to noise higher than here: analytical and numer- the parameters affects the number of clusters, their con-
ical means corroborate the Hopf bifurcations to emerge figuration, stationary or dynamical character, as well as
about D 0.0025 at relevant τ. Now we argue that the whether the cluster states occur monostable or coexist
≈
approximate model is in qualitative terms able to cap- with the homogenous solution at the given population
ture the clustering effect occurring for small D,c and τ. size. This framework could find application within the
Focus is on the finding that MF system predicts an on- researchon neural systems and other excitable media.
This work was supported in part by the Ministry of
set of cluster states by undergoing a global bifurcation
Education and Science of the Republic of Serbia, under
for the parameter values around τ =2,D =0.00025 and
project Nos. 171017 and 171015.
c=0.08. At the givenτ and D, for c<0.08the approx-
imate model has only the equilibrium, whereas around
c 0.08a large anda small LC are born via a fold-cycle ∗ [email protected]
≃
scenario. Note how then the PP of the MF acquires the [1] P. A. Tass, Phase Resetting in Medicine and Biology:
form qualitatively similar to those of the exact system’s StochasticModellingandDataAnalysis,(Springer,Berlin
in Figs. 2(a) and2(b). The two sections of the emerging Heidelberg, 2007).
[2] N. Buri´c and D. Todorovi´c, Phys.Rev. E, 67, 066222
MF orbit mimic the action of the fractions within the
(2003);M.Dhamala,V.K.JirsaandM.Ding,Phys.Rev.
fullpopulation. ThisstructureoftheLCbecomesunsta-
Lett, 92, 074104 (2004); S. H. Strogatz, Nature 394, 316
ble under increasing c or τ, i.e. for the stronger impact
(1998); O. V. Popovych, C. Hauptmann and P. A. Tass,
of the interaction term. Another interesting aspect to Phys. Rev.Lett. 94, 164102 (2005).
the approximatesystemis thatit showsthe complexLC [3] A. Pikovsky, A. Zaikin and M. A. de la Casa, Phys. Rev.
to coexist with the FP, viz. the basins of attraction in Lett. 88 050601 (2002); B. Lindner, J. Garcia-Ojalvo, A.
Fig. 5(b), whichisafeatureapparentlyabsentintheex- Neiman and L. Schimansky-Geier, Phys. Rep. 392, 321
(2004);N.Buri´c,K.Todorovi´candN.Vasovi´c,Phys.Rev
act model. However, the FP is located very close to the
E, 75, 026209 (2007); Q. Wang, G. Chen and M. Perc,
basins’ boundary which indicates it to be stochastically
PLoS ONE6, e15851 (2011).
unstableintheexactsystemforanarbitrarysmallnoise. [4] X. Sailer et al., Phys.Rev.E 73, 056209 (2006).
We have reported on a novel phenomenon where clus- [5] P.Kalutza,C.StregeandH.Meyer-Ortmanns,Phys.Rev.
tering within the homogeneous neural population is in- E 82, 036104 (2010).
duced by an interplay of noise and time delay. This [6] M. Yi and L. Yang, Phys. Rev. E 81, 061924 (2010); X.
paradigm is distinct from most current explanations on Sun, M. Perc, Q. Lu and J. Kurths, Chaos 20, 033116
(2010).
howtheclusteredstatesmayarise,foritdoesnottreatD
[7] A. Pikovsky, J. Kurths, Phys. Rev. Lett. 78, 775 (1997);
and τ as destabilizing and detrimental, but rather as bi- Y. Gao and J. Wang, Phys.Rev. E83, 031909 (2011).
asedtowardtheformationofdynamicalstructureinnet- [8] Neuro-Informatics and Neural Modelling, edited by F.
works that are unstructured both in terms of topology Moss and S. Gielen (Elsevier, Amsterdam, 2000), p.903.
and local parameters. The analyzed model is minimal [9] N.Buri´c,K.Todorovi´candN.Vasovi´c,Phys.Rev.E,82,
yet sufficient to display an interesting type of behavior, 037201 (2010)