Synchronization and long-time memory in neural networks with inhibitory hubs and synaptic plasticity

6.8 MB·

by  Elena Bertolotti

Synchronization and long-time memory in neural networks with inhibitory hubs and
synaptic plasticity
Elena Bertolotti,1,2,∗ Raffaella Burioni,1,2,† Matteo di Volo,3,4,5,‡ and Alessandro Vezzani1,6,§
1Dipartimento di Fisica e Scienza della Terra, Universita` di Parma, via G.P. Usberti, 7/A-43124, Parma, Italy
2INFN, Gruppo Collegato di Parma, via G.P. Usberti, 7/A-43124, Parma, Italy
3Group for Neural Theory, Depart´ement des Etudes Cognitives, Ecole Normale Sup´erieure, Paris, France
4Centro Interdipartimentale per lo Studio delle Dinamiche Complesse, via Sansone, 1-50019 Sesto Fiorentino, Italy
5Indiana University-Purdue University, 420 University Blvd., Indianapolis, Indiana 46202, USA
7 6IMEM–CNR, Parco Area delle Scienze 37/A-43124 Parma, Italy
1
Weinvestigatethedynamicalroleofinhibitoryandhighlyconnectednodes(hub)insynchronizationand
0
input processing of leaky-integrate-and-fire neural networks with short term synaptic plasticity. We take
2
advantageofaheterogeneousmean-fieldapproximationtoencodetheroleofnetworkstructureandwetune

n thefractionofinhibitoryneuronsfI andtheirconnectivityleveltoinvestigatethecooperationbetweenhub
a features and inhibition. We show that, depending on fI, highly connected inhibitory nodes strongly drive
J the synchronization properties of the overall network through dynamical transitions from synchronous to
asynchronousregimes. Furthermore,ametastableregimewithlongmemoryofexternalinputsemergesfora
8
specificfractionofhubinhibitoryneurons,underliningtheroleofinhibitionandconnectivity alsoforinput
1
processinginneuralnetworks.

]
n
n I. INTRODUCTION cooperation between inhibition and connectivity structureas
- a drivingof synchronization properties in neural networks.
s
i Neural ensembles feature a wide range of dynamical pat- In this paper we investigate theoretically and numerically
d terns, that can be suitably analyzed through collective ob- theroleandtheinterplayofthefractionofinhibitoryneurons
.
t servables, like local field potentials generated by single units andoftheirhubcharacter(theirrelativeconnectivitywithre-
a
firing activity [1–3]. These observables typically display os- spectoftherestofthenetworkunits)insynchronizationand
m
cillations on different frequency ranges, that encode specific functional properties of a neural network. In particular, we
- structuresofsynchronizationintheactivityofsingleneurons considerherealeaky-integrate-and-fire(LIF)neuralnetwork
d
[4, 5]. Synchronization among neurons has been shown tobe composed by inhibitory and excitatory neurons with a short
n
relevantfornormalneuralfunctionsandalsotobeassociated, term synaptic plasticity mechanism, developed by Tsodyks
o
whenanomalous, with neuraldisorders [6,7]. Thereforeneu- and colleagues (Refs. [18–20]). In order toinvestigate the in-
c
ral synchronization has been so far widely investigated, both terplaybetweeninhibitionandconnectivity,wetunetwomain
[
from experimental and theoretical pointsof view [8, 9]. parameters: thefractionofinhibitoryneuronsinthenetwork,

1 Severalfeaturesofneuralnetworkshavebeenshowntoin- f ,andtheirhubcharacterquantifiedbythedistancebetween
I
v fluencesynchronizationofneuralpopulations,e.g.,neuralex- the average connectivity hki of inhibitory (I) and excitatory
6 citability [10], the presence of inhibitory component [11] and (E)populations,∆=hk i−hk i. Weshowhowtheinterplay
I E
5 the structure of connections [12]. In recent years, a crucial ofthesetwoingredientsgivesrisetoawiderangeofdynamical
0 role of inhibitory component in neural ensembles has been regimes, with specific synchronization properties and differ-
5 proposed to reproduce specific patterns observed in cortical entabilitytoprocessexternalinputs. Thenetworkdynamics
0
regions of the brain. Indeed, an optimal balance between is typically organized in aquasi–synchronous regime for very
.
1 excitatory and inhibitory neurons is required to obtain pat- smallfractionsofinhibitoryneuronsandturnstoacompletely
0 terns of activity observed in different experimental setups, asynchronous state when the network is mainly inhibitory.
7 and an optimal fraction of inhibitory components has also Nevertheless,increasingthefractionofinhibitorycomponents
1 beensuggested tooptimizetheperformances ofaneuralnet- doesnotsimplydecreaseneuralsynchronization. Indeed,the
: work [13, 14]. synchronizationofthenetworkincreasesforintermediateval-
v
Also the structure of connections plays a central role in uesoff ,reachinganalmost synchronousregimeforanopti-
i I
X neurons synchronization and in their capacity to store and malbalancebetweenexcitation andinhibition. Furthermore,
process information [15]. Recently, evidence of hubs neurons closetothisbalance,wefindanintriguingmetastabledynam-
r
a in the brain, i.e., neurons with high connectivity, has been icalregime,whichturnsouttobethemosteffectiveforstoring
putforward [16]. Theseneuronsarecrucial inorganizing the information of external inputs. As a result of this analysis,
synchronization of neurons they are connected with. More- the fraction of inhibitory neurons and their hub character is
over,theactivityofasinglehubneuroncanstronglyinfluence showntocombinenontriviallyfortheemergingnetworksyn-
the firing pattern of the overall population. A recent paper chronization and inputprocessing capability.
by Bonifazi et al. [17] has put into evidence that hub neu-
The model we consider is described in details in Sec. II,
rons are typically inhibitory, suggesting an unifying view of
where we also implement a heterogeneous mean-field model
(HMF) reproducing the finite size dynamics for sufficiently
high connectivity (see, e.g., [21, 22]). Such a HMF model
turns out to be very useful to represent the role of connec-
∗ Electronicaddress: elena.bertolotti1@fis.unipr.it tivity distribution and to reduce the computational cost of
† Electronicaddress: raffaella.burioni@fis.unipr.it numerical simulations. An extended analysis of the HMF
‡ Electronicaddress: [email protected] model ability to reproduce finite size dynamics is then re-
§ Electronicaddress: alessandro.vezzani@fis.unipr.it ported in Sec. VI. In Sec. III we discuss the synchronization
2
properties of the network as a function of f and of the hub theinactiveresourcesz arerecovered,returningavailablein
I ij
characterofinhibitoryneurons∆andweanalyzethedynam- atimeτi. Actually,whenneuronj spikes,itactivatesafrac-
r
icaltransitionsfrom partially synchronous,synchronous,and tion u of the available resources x . If the postsynaptic
ij ij
asynchronous regimes. In Sec. IV we describe the dynamical neuron i is excitatory, u = U is constant in time. Oth-
ij
features of the metastable regime emerging for specific val- erwise, if i is inhibitory, u evolves according the following
ij
ues of f and ∆. In Sec. V we investigate how the network facilitation mechanism:
I
processes external inputs, showing a high capacity of storing
u (t)
information in the metastable regime. Section VI is devoted u˙ (t)=− ij + U [1−u (t)]S (t). (6)
ij τ f ij j
to the comparison between the HMF results and finite size f
networkssimulations,toshowthatthetwodynamicsbecome
Accordingly,thevariable u decaysin atime τ , while it in-
ij f
similarbyincreasingnetworkconnectivity,asexpected. How-
creasesinstantaneouslyofanamountU atthespikingevent.
f
ever, this limit is not trivial, since in the quasi-synchronous Moreover, forpostsynapticinhibitory neurons,τi istypically
and asynchronous regimes the effect of a finite connectivity r
much smaller than in the excitatory case. The equations are
is a noise superimposed to the HMF behavior, while in the thenclosedbythesynapticcurrentIsyn(t)inEq.(1),defined
balance regime the HMF dynamics is obtained considering i
as thesum overall theconnected active resources y :
ij
metastable states in finite size networks, whose lifetime di-
verges with the system connectivity. Eventually, in Sec. VII, Isyn(t) = g ǫ y (t), (7)
we discuss conclusions and perspectives. i hki ij ij
Xj6=i
where g is the strength of the synaptic coupling and ǫ is
ij
II. FROM THE FINITE SIZE MODEL TO THE theconnectivitymatrix (ǫ =0iftherearenolinksbetween
ij
HMF FORMULATION i and j, ǫ = 1 if the presynaptic neuron j is excitatory,
ij
ǫ =−1ifj isinhibitory). Thenumberofsynapsesoutgoing
ij
from neuronj isk = |ǫ |. Similarly,theinputdegree
We consider a neural network where each node represents j,out i ij
k = |ǫ | corresponds to the number of neurons from
an excitatory or inhibitory neuron and directed links model i,in j ij P
thesynapsesconnectingthecorrespondingnervouscells. Neu- which i receives synaptic inputs. Accordingly, hki is the av-
P
ral dynamics is described by a LIF oscillator (Refs. [23, 24]), erage in-degree of the graph. Equations (1)-(7) can be very
sothat theevolution ofthemembranepotentialv (t)ofneu- efficientlysimulated throughan event-drivenmap,exploiting
i
ron i is given by the integrability of the model between two subsequent firing
events(see, e.g., Refs. [25, 26]).
v˙i = a−vi+Iisyn(t), (1) In this paper we will study the case where the inhibitory
neurons are hubs of the networks, i.e., their connectivity is
where a is the leakage current and Isyn(t) is the synaptic
i much larger than theaverage. We are interested in networks
current, encoding the coupling with other neurons. Without
with hubs displaying a large numberof both in and out con-
loss of generality we have rescaled variables and parameters,
nections. Accordingly, we impose a simple correlation be-
that are now expressed in adimensional units. Wheneverthe
tweenin-andout-degree,sothatk = k =k ,foreach
membranepotentialv reachesthethresholdvth = 1,aspike i,out i,in i
i i node i. In particular we consider a “configuration model,”
is instantaneously sent toall thepostsynaptic connections of
i.e., random networks with fixed degree distribution. As we
neuron i and v is reset to vr = 0. For a > 1 the neuron
i i are interested in networks with very large nodes degree, that
emits periodically spikeseven if it does not interact with the
is with k → ∞, we can approximate the discrete degree k
i i
restofthenetworkanditissaidtoevolveinaspikingregime.
with a continuous variable k and we can call P (k) [P (k)]
E I
WedescribethecouplingamongneuronswiththeTsodyks-
thedegreep.d.f.ofexcitatory(inhibitory)neurons. Asacon-
Uziel-Markram (TUM) model for short term synaptic plas-
sequence, the fraction of excitatory (inhibitory) neurons f
E
ticity, already presented in [19]. Each synapse connecting
(f ) is
I
neuron j to neuron i is endowed with synaptic transmitters
or resources, that can be found in three different states. Ac-
f = P (k)dk,
cordingly,foreachsynapse(i,j)afractionofavailable,active, E E
andinactiveresources,x ,y ,z ,respectively,isdefinedand Z (8)
ij ij ij
their time evolution equations read fI = PI(k)dk,
Z
y (t)
y˙ij(t) = − iτj + uij(t)xij(t)Sj(t), (2) where fE +fI = 1 and the connectivity p.d.f. for a generic
in neuron is P (k)+P (k).
E I
z (t)
x˙ (t) = ij − u (t)x (t)S (t), (3) Forpurelyexcitatorynetworksnumericalsimulationsshow,
ij τri ij ij j asafunctionofg andτin,differentdynamicalphasesranging
x (t) + y (t) + z (t) = 1, (4) fromanasynchronousregimetoapartiallyandatotallysyn-
ij ij ij
chronousone[27]. Inparticular,in thepartially synchronous
where the last equation is the resources normalization and
state,allneuronsin thenetwork arrangeintotwosubgroups:
S (t) represents the spike train produced by the presynaptic
j thelockedones,withaperiodicdynamicsandacommon pe-
neuron j:
riod,andtheunlockedones,with anaperiodic dynamicsand
Sj(t) = δ[t−tj(n)], (5) different periods. The only topological property, apart from
n finitesizefluctuations,thatsurprisinglydistinguishesalocked
X
wheret (n)isthetimeatwhichneuronj emitsitsnthpulse. neuron from an unlocked one is its input degree k .
j i,in
In between two consecutive spikes by neuron j the active re- Figure 1 shows that the same property is preserved for
sourcesy exponentiallydecaywithatimeconstantτ ,while networks with both excitatory and inhibitory neurons, as
ij in
3
1.4 5000 neuron i with the average over the whole graph. We then
4000 take into account that the contribution of each neuron has
1.2
to be weighted with its out-degree, that is, in a graph with
3000
SIik 1 i12000000 ykii,obuytk=i/khik,iin.=Ackcio,rwdeinngelye,dwteocmaunltwiprliyteeachmicroscopicfield
I0.8
h0.6 0 9121 9122 912t3 9124 9125 9126 Yi(t)≈ hkkiiN1 N kjhykji(t) = hkkiiY(t), (10)
j=1
X
0.4
100 150 200 250 300 350 400 wherewehavedefinedtheglobalquantityorfieldY(t),which
k
no longer depends on the single neuron i. Such a mean-field
approximation is exact in the limit of k → ∞ and it can
i
FIG. 1. hISIik as a function of the degree k for a network of berigorously treated formassivegraphs,wherethedegreeki
N = 5000 neurons, which evolves according Eqs. (1)-(7). In the
scaleswiththenumberofneuronsN. Furthermore,theHMF
inset we show a raster plot of the corresponding dynamics, where
approximation allows toreplace thedynamicsfor each single
nodesareorderedalongtheverticalaxisaccordingtotheirdegree
nodewiththedynamicsofeachclassofnodeswithdegreek .
k; a dot means that neuron i has fired at time t. Here and in all i
otherplotstimeisdimensionless. TheinhibitoryfractionfI is0.1. In particular, in a network with both excitatory and in-
Blue (green) data refer to the excitatory (inhibitory) nodes. The hibitory neurons, for each class of degree we introduce two
excitatory average degree hkEi is 100, the inhibitory one hkIi is variables vkE and vkI describing the membrane dynamics of
350 and the standard deviation σ of both Gaussian p.d.f.s PE(k) excitatory and inhibitory neuron class respectively. Indeed,
andPI(k)is10. duetothefacilitationmechanism,thesynapticcurrentreach-
ing an inhibitory neuron is different from the current of an
excitatory neuron and this gives rise to a different evolution
also shown in [28]. Data are obtained from a network of of the membrane potential. In this framework, four different
N = 5000 neurons, f = 0.1, where both P (k) and P (k) heterogeneous mean fields need to be defined, i.e., the field
I E I
are Gaussian with standard deviation σ = 10 and average produced and received by excitatory and inhibitory neurons
values hk i = 100 and hk i = 350 respectively. In the inset, respectively [28]. If we denote with ∗,† the indexes E,I for
E I
the raster plot represents the spike event of each node as a excitation and inhibition, thefour global fields are
functionoftimeandnodesareorderedaccordingtotheirde-
greeki. Weobservethattheexcitatoryneurons(bluedatain Y (t) = P∗(k)kyk(†,∗) dk, (11)
Fig. 1) with small connectivity ki are actually characterized †,∗ hki
Z
by a periodic and phase-locked dynamics, while all theother
nodes evolve aperiodically (excitatory nodes with k > 115 wherey(†,∗) representstheactiveresources of nodeswith de-
i k
and all inhibitory nodes). The dynamical properties of the greek andoftype∗outgoingtowardsanodeoftype†. Simi-
ith neuron are synthetically represented by the time aver- larlyx(†,∗) andz(†,∗) aretheavailableandinactiveresources.
k k
age hISIi of the interspike interval (ISI is defined as the time Accordingly, the synaptic fields Y (t) and Y (t) received
E I
lapse between two consecutive firing events). Plotting hISIi respectively byexcitatory and inhibitory neuronsare
as a function of the node degree, locked or periodic excita-
tory neurons with low degree (i.e., with ki <∼ hkEi ) belong YE(t)=+ YEE(t) − YEI(t), (12)
to the initial horizontal plateau while, for all other aperiodic Y (t)=+ Y (t) − Y (t),
I IE II
neurons, hISIi depends on k . Therefore, also in this frame-
i
work, the dynamical properties of a neuron seem to depend where the signs + and − take into account the excitatory
on the network topology, apart from finite size fluctuations, and inhibitory nature of the presynaptic neuron. The final
only through thedegree of thenode. equations for thenetwork dynamics in the HMF formulation
We thus expect that a heterogeneous mean-field approach are therefore
(Refs. [29, 30]), which actually gets rid of connectivity pat- g
v˙†(t)=a−v†(t)+ kY (t),
terns and takes into account only the neurons connectivity k k hki †
k andthein-degreedistribution,couldbeeffectivelyapplied
i y(†,∗)(t)
also in oursetup. Thefield Yi(t)received byneuron icanbe y˙(†,∗)(t)=− k +u(†,∗)(t)x(†,∗)(t)S∗(t),
written as k τin k k k
y(†,∗)(t) z(†,∗)(t)
1 N k 1 ki z˙(†,∗)(t)= k − k ,
Yi(t)= hki ǫijyij(t)= hkiik yj(t), (9) k τin τr† (13)
Xj=1 i Xj=1 x(†,∗)(t)+y(†,∗)(t)+z(†,∗)(t)=1,
hi,ji k k k
u(E,∗)(t)=U,
wherewehaveputintoevidencethedegreeoftheithneuron. k
Thenotationhi,jibelowthesecondsummeansthatonlythe u(I,∗)(t)
coupleswhereiandj arenearest-neighborarenowincluded. u˙(kI,∗)(t)=− kτ +Uf[1−u(kI,∗)(t)]Sk∗(t).
f
Indeed, each row of the matrix ǫ is a vector with (N −
ij
k ) elements equal to 0 and k elements equal to 1 and here Unlessotherwisespecified,thevaluesofallparametersap-
i i
we want to consider only the k nonzero elements, that is pearing in the equations are set according to accepted phe-
i
those j indices such that ǫ = 1. Then we make a mean- nomenological values (Refs. [31, 32]), after a proper time
ij
field approximation on the right-hand side of Eq. (9), i.e., rescaling, as follows: τ =0.2, τI =17·τ , τE =133·τ ,
in r in r in
we approximate the average over the k neighboring sites of τ =33.25, a=1.3, g=30, U =U =0.5.
i f f
4
Integrating the set of Eqs. (13) through the same steps 3
0.1
usedforthedynamics(1)-(3),asimilarevent-drivenmapcan 2.5 1.5 0.15
be obtained. In order to simulate the dynamics of the HMF 0.17
model we need to perform a sampling of the k values from 2 1.4 0.214
the p.d.f P(k). We point out that numerical simulations are i 0.22
now very efficient, since, through an importance sampling of SIk1.5 1.3 80 90 100 110 120 00..2259
the degree k, say we sample M values of k, M ∼ O(102) is Ih
1
sufficient to obtain simulations wherethefluctuationsdueto
discretization of the distribution are very small. In Sec. VI
0.5
we will compare the dynamics of finite large networks with 100 150 200 250 300 350 400
the results of the HMF approach and we will show that two k
different behaviors emerge. In the quasi-periodic regime, at (a)
small f , the effect of a finite degree is simply a dynamical
I 500
noisesuperimposedtothemean-fieldsolution. Moreinterest-
ingly,weputintoevidencethatinthebalanceregionthesyn-
400
chronousmean-fieldsolutionisunstableforthefinitenetwork
and an asynchronous state is always reached asymptotically; 300
however,inthiscasethelifetimeofthesynchronousstatedi-
i
verges with the connectivity, recovering also in this case the 200
HMF dynamicsfor large enough connectivities.
100
f =0.1 f =0.22 f =0.29
I I I
0
III. SYNCHRONIZATION EFFECTS OF 1943 1945 1947 117107 117109 117111 1707 1709 1711
t
INHIBITORY HUBS
(b)
In this section we show numerical HMF results, putting
into evidence the synchronization effects of inhibitory hubs. FIG.2. (a)hISIkiasafunctionofthedegreek fornetworkswith
differentinhibitoryfractionsfI (seethelegend). ∆isfixedat250
Asthehubnatureofaneurondependsonlyonitsdegree,we and the number M of connectivity classes is 500. In the inset we
needtochoose suitable p.d.f.s PE(k) andPI(k). Forsimplic- show a zoom of the excitatory plateau. In panel (b) we display
itywefixthetwop.d.f.stobebothGaussians with standard snapshots of the raster plots for different inhibitory fractions, to
deviation σ = 10 and averages hk i and hk i for excitatory represent the three main dynamical regimes: partial synchroniza-
E I
andinhibitorynodes,whilethedistancebetweenthetwodis- tion for fI = 0.1, total synchronization for fI = 0.22, and an
tributions is ∆ = hkIi−hkEi. The larger ∆ is, the more asynchronous state for fI = 0.29. A dot in raster plots of HMF
hub the inhibitory neurons are. In particular, in our simu- simulations represents the spike event of connectivity class i at
timetandclasses areorderedalong the vertical axis accordingto
lations we fix hk i = 100 and we vary f and ∆, that tune
E I theirdegreek. Thedegreeclassesbelow(above)theblacklineare
the fraction of inhibitory nodes and how strong is their hub
excitatory (inhibitory). These plots refer to different dynamical
character, respectively. From a technical ground, the contin-
regimes, however a time unit in the abscissa axis corresponds to
uous distributions PE(k) and PI(k) are efficiently sampled, aboutoneglobaloscillationofthenetwork(afterthistimeinterval
obtainingM discreteconnectivityclasses, each oneincluding allneuronshavefiredatleastonce).
all neuronswith degree k∈(k ,k ) (i = 1, ..., M −1).
i i+1
First let us fix ∆=250 in order to observe the dynamical
configurations of the network as a function of fI. The raster work,wecalculatetheKuramotoparameterR(Refs.[27,33]),
plotsinFig.2(b)havebeenobtainedforthreedifferentvalues definedas follows:
of this parameter. When f is close to 0 the network is par-
I
tially synchronized and neurons are divided into two groups: R= 1 M e2πi[(t−tnk)/(tnk+1−tnk) , (14)
excitatory neurons with low degree are locked, that is they *(cid:12)M (cid:12)+
evolveaccordingaperiodicdynamicsandtheysharethesame (cid:12) Xk=1 (cid:12)
(cid:12) (cid:12)
period, whileall theothers(excitatorynodeswith higherde- wherethevertical(cid:12)barsrepresentthemodulus(cid:12)ofthecomplex
(cid:12) (cid:12)
gree and inhibitory hubs) are unlocked, displaying different numberandtheanglebracketsthetemporalaveragevalueof
periods and an aperiodic dynamics. The plot of hISIki as a that quantity, while the exponent is the phase we assign to
function of k, as shown in Fig. 2(a), synthetically represent the connectivity class k. The closer to 0 (1) R becomes, the
the dynamical properties of each connectivity class. Locked more asynchronous (synchronous) the dynamics is. M is the
neurons, whose hISIki does not depend on theirdegree, give total number of classes, tnk is the time of the nth impulse of
risetoaplateau,aswecanobserveinthatfigureforfI =0.1 classk,andtnk+1isthetimeofthefollowing(n+1)thimpulse.
and k<106. Foreach ∆,we varyf from 0to avaluethatcorresponds
I
ForlargervaluesoffI,thenetworksynchronizationgrows. tothetotally asynchronousnetworkandwecomputetherel-
Indeed, the plateau in Fig. 2(a) in the excitatory region ative Kuramoto parameter R, namely the parameter 1−R,
stretches, while a new inhibitory plateau appears. After the plotted in Fig. 3 as a function of f in semilogarithmic scale
I
network has reached the total synchronization for a particu- to make the transition from the asynchronous to the syn-
lar valueof fI (which for ∆=250 is 0.214), with all neurons chronous regime more evident. For all values of ∆ ≥ 0, R
locked and hISIki belonging to two identical plateaus, the firstincreases towards 1, meaningthat,at small valuesoffI,
network becomes more and more asynchronousfor larger fI. increasing the fraction of inhibitory neurons really improves
In order to quantify the synchronization level of the net- thesynchronizationofthewholenetwork. Then,thesynchro-
5
Writing Y (t) according to Eqs. (12) and (11),we obtain
†
Y (t)−Y (t)=0,
†,E †,I
P (k)ky(†,∗) P (k)ky(†,∗)
E k dk= I k dk,
hki hki
Z Z
y y (15)
P (k)kdk= P (k)kdk,
hki E hki I
Z Z
hk if =hk if ,
E E I I
hk i(1−f )=(∆+hk i)f ,
E I E I
where in the second step we have used y(†,∗) = y(∀k) and
k
we factorized them out of the integrals, since the network is
totally synchronous and all microscopic variables are equal.
FIG.3. Networksynchronizationfromthemeasureof1−R(where
Simplifyingy/hkifromthethirdequation,weobtainthatthe
RistheKuramotoparameter)asafunctionoffI fordifferentval-
uesof∆(seethelegend),insemilogarithmicscale. Inallnetworks left integral represents the average value of PE(k) multiplied
M is equal to 500. The closer to 0 the parameter 1−R is, the by fE and similarly the right integral [see Eq. (8)]. Solving
more synchronous the network. The grey-shaded area represents the last equation for fI we get the inhibitory fraction fIB
theintervalofinhibitoryfractions,whichhasbeenreallyobserved for the balance regime. We verified indeed that, by fixing
in some brain tissues. Without hub nodes (∆ = 0, black data), the inhibitory fraction to the predicted balance value, the
thenetworkreachesahighsynchronizationleveloutofthisregion, Kuramoto parameter is maximum and very close to 1: R =
differently from 160 <∼ ∆ <∼ 500. Therefore, if we want to build 0.987−0.997forthedifferentvaluesof∆,asshowninFig.3.
synchronousnetworkswithinhibitoryfractionsclosetorealvalues,
Furthermore,inthebalanceregimeallneuronsarelockedand
theinhibitoryneuronsshouldbehubs,withanaveragedegree3-6
characterizedbythesameperiod,thatcanbecomputedfrom
timesgreaterthantheexcitatoryaveragedegree.
the membrane potential equation [first Eq. in (13)] without
the field term. The period of these free oscillators is simply
hISI i=ln[a/(a−1)]=1.466(∀k).
k
nization level reaches a peak, where the network is totally Outofthisbalancephase,oneofthetwomechanics,excita-
synchronized and R is very close to 1. We notice that for toryorinhibitory,exceedstheother. Wemeasuretherelative
different valuesof ∆ thispeak correspondsto adifferent and weight of theexcitatory and theinhibitory component intro-
specific value for the inhibitory fraction. Finally, R quickly ducingthe following quantities:
decays towards 0 and the network becomes asynchronous.
Y (t)−Y (t)
If inhibitory neurons do not feature a higher connectivity WE = YEE(t)+YEI(t) ,
with respect to the excitatory ones (∆=0), we observe that (cid:28) EE EI (cid:29) (16)
Y (t)−Y (t)
theinhibitory-drivensynchronousregimeisobservedonlyfor W = IE II ,
high values of fI. Nevertheless, experimental observations I (cid:28)YIE(t)+YII(t)(cid:29)
suggestthatthefractionofinhibitoryneuronsismuchsmaller
where the angle brackets represent again the temporal aver-
(f ∼10−30%)[34,35]. Actually,weshowhowbyincreasing
I age. W refers to the field received by an excitatory neuron,
E
thehubcharacterofinhibitoryneuronsonecan observesuch
while W refers to the field received by an inhibitory neu-
I
a synchronized regime for much smaller values of f , reach-
I ron. If W ≃W ≃1, the dynamics of the whole network is
E I
ingphenomenologicalvalues(grey-shadedregioninFig.3)for
purelyexcitatoryandinhibitoryneuronsdonotsubstantially
2×102 <∼∆<∼103. Apartfrom thequantitativereliability of contributetotheglobalfields. InthelimitofW ≃W ≃−1
E I
ourresultswithrespecttotheexperimentallysuggestedfrac-
the dynamics is totally driven by the inhibitory mechanism,
tionofinhibitoryneurons,Fig.3showshowthehubcharacter
eventhoughthislimitistheoretical, sinceapurelyinhibitory
of inhibitory neurons is a fundamental ingredient to allow a
networkisnotinaspikingregime. InFig.4weplotW and
E
small fraction of inhibitory nodes to drive the synchroniza-
W forthevaluesof∆andf consideredinFig.3. Weimme-
I I
tion of the overall network, thus underlining the importance diatelyobservethatforthebalancevaluefB oftheinhibitory
of inhibitory neurons connectivity in affecting network syn- I
fraction we get W =W =0.
E I
chronization. The narrow region of f where the network is
I However,weseethattheincreaseoff doesnotnecessarily
I
highly synchronized is a consequence of a dynamical balance
lead toamonotonousincrease intheinhibitorycontribution,
betweentheexcitatoryandtheinhibitorymechanisms,which
but some unexpected oscillations appear, due to nontrivial
will be discussed in the nextsection.
equilibrium mechanisms between the excitatory and the in-
hibitory dynamics outside thebalance region.
The Kuramoto plot in Fig. 3 can be now compared with
whatemergesinFig.4. Atfirst,foreach∆,asthefractionof
IV. BALANCE REGIME inhibitoryhubsincreases andthenetworkdynamicsbecomes
more inhibitory (in Fig. 4 W and W are indeed decreas-
E I
ing), the synchronization improves, starting from a partially
In the totally synchronous phase, all neurons follow the synchronous evolution, until the totally synchronous state is
same dynamics and the evolution of v†(t) no longer depends reachedforthebalancevalueoff . Then,thefurtherincrease
k I
on k. This means that the last term in the first equation off leadstotheprevalenceoftheinhibitorydynamics,since
I
in (13) should be 0, i.e., the excitatory and the inhibitory both W and W become negative, and it gives rise to asyn-
E I
contributionsin thedefinition of Y (t) cancel each otherout. chronousregimes. Byincreasing∆,thisdynamicaltransition
†
6
1 ∆ = 0 500
∆ = 160
400
∆ = 250
WE,I0.5 ∆∆ == 5100000 i300
200
0
100
0 0.1 0.2 0.3 0.4 0.5 026978 26982 26986 26990 26994 26998
fI t
FIG.4. Relativeweightsoftheexcitatoryandtheinhibitorycom- FIG. 6. Raster plot of a network with M = 500 connectivity
ponentindrivingthenetworksdynamicswithM =500connectiv- classes, ∆ = 250 and fI = 0.1. At a given time (see the red
ityclasses,fromthemeasureofWE andWI fordifferentvaluesof arrow), we apply a synchronizing external stimulus to 30% of ex-
∆(seethelegend)andfI. ThebalancevaluefIB oftheinhibitory citatoryneurons,randomlychosen. Afterabouttenglobaloscilla-
fractionforwhichWE =WI =0is0.5for∆=0,0.28for∆=160, tions, the stimulated nodes recover the initial dynamics. The de-
0.22for∆=250,0.143for∆=500and0.083for∆=1000. greeclassesbelow(above)theblacklineareexcitatory(inhibitory).
Light(dark)bluedotsrepresentthespikeeventsofstimulated(not
stimulated)nodes.
500
400
starting from random initial conditions, an evidence of the
300 slow dynamics. For this reason, we point out that in Fig. 3
i the value of the Kuramoto parameter of networks with this
200 metastabledynamicshasbeencomputedreferringtothetotal
synchronous state, namely to the state the network assumes
100
using total synchronized initial conditions. The robustness
of the metastable dynamics is also preserved, even when we
0
290704 290706 290708 290710 581491 581493 581495 581497 modifythediscretizationforthedegreeclassesfromtheinitial
t t
continuous p.d.f. P(k), which is another relevant parameter
for the simulations.
FIG.5. Rasterplotsofametastablepartiallysynchronousstatefor
thenetworkwithM =500,∆=250,andfI =0.214,startingfrom
randominitialconditions. Theseplotsaretwosnapshots,takenat
V. RESPONSE TO AN EXTERNAL STIMULUS
different times during the same dynamical evolution. The degree
classesbelow(above)theblacklineareexcitatory (inhibitory).
Inthissection weinvestigatetheinputprocessing features
of thenetwork byanalyzing thereaction of our system to an
from a synchronous to an asynchronous configuration occurs external stimulus. In particular, we focus on the transient
in a narrower interval of fI. dynamicsatthestimulusoffsetasthespacetrajectory ofthe
When ∆ has intermediate values, i.e., 150<∆<500, the system in this transient has been supposed to be crucial for
region immediately before the balance (for example 0.2 < thestimulusdetection[36]. Inparticular,thelongerthetime
f < 0.214 for ∆ = 250) is characterized by a highly to return to the unperturbed state is, the better the system
I
metastable regime of partial synchronization, where almost is able to detect and distinguish theinput [37].
allexcitatoryneuronsarelockedwithverysimilarfrequencies We will show that the metastable states, like those in
butdifferentphasesandthereforetheysplitintodifferentsyn- Fig. 5, can be affected by the external perturbation longer
chronized groups, as we can see in the raster plots in Fig. 5. than the other dynamical regimes. Indeed, as the configu-
In this regime, inhibitory neurons are unlocked, even if they ration of the synchronous excitatory groups depends on the
spikealmostsimultaneouslywithoneoftheexcitatorylocked initialconditionsandchangesonlong-timescales,anexternal
groups, but the group changes in time. Moreover, the sys- synchronous stimulus, creating a new group, can produce an
temis metastable. Indeed,eventheconfiguration ofthesyn- alteration of the network dynamics for long-time intervals.
chronousexcitatory groups changes on verylong-time scales, We consider ∆ = 250 and vary f , so that the different
I
typically of thousands of oscillations; for example the raster dynamical regimes can be explored. We stimulate only the
plots in Fig. 5 are different snapshots taken during an evolu- excitatory nodes, since they are the most abundant in real
tion of the same network. As one can expect, the structure networks, as it has been experimentally observed, and there-
ofthesemetastablestatesdependsontheinitialconditionsof fore they are more likely to be affected by external pertur-
the microscopic variables. In particular, if one imposes syn- bations with respect to inhibitory cells. In all networks, we
chronousinitialconditions(potentialsandsynapticresources randomly choose 30% of the excitatory neurons and we syn-
are equal for all neurons), the dynamics reaches a configura- chronizethemartificially; e.g., wecan forcethesimultaneous
tion where a single group of synchronous neurons is present firing of the stimulated neurons, when the excitatory field
and still some are unlocked. This configuration seems to be Y (t) reaches its minimum value (see Fig. 6). Such forced
E
theasymptoticstablesituationsince,atvariancewithFig.5, synchronization, that may represent external visual or odor
wedonotobserveasystemmetastabilityonlong-timescales. input,couldencodetheintensityofthestimuli: thehigherthe
However,weneverreachsuchatotalsynchronousstatewhen numberof synchronousunits thehighertherelease of neuro-
7
100 100
) 0.1
−Rt(10−5 0 Rt() 00..11
1 0.1 − 0.1
10−10 0.214 (metastable) 110−10 0.214
0.23 0.214
0.214
0.28
0.214
10−2 100 t 102 104 10−2 100 102 104
t
FIG. 7. Transient dynamics at the stimulus offset: return to the 100
original synchronization level after the application of an external
input. We plot the parameter 1−R(t) of the stimulated nodes
as a function of time, for networks with M = 500, ∆ = 250, and
sdtiiffmeruelnuts visalaulewsayofsfaIpp(lsieeedtahtetle=ge0ndto).3T0%heoefxetxercnitaaltosryyncnheruornoonuss, Rt()10−5 10%, 0.1
randomlychosen. − 20%, 0.1
1 30%, 0.1
10%, 0.214
20%, 0.214
transmitter y dueto thespikingof forced neurons. After the 10−10 30%, 0.214
perturbation, we compute the reduced Kuramoto parameter 10−2 100 102 104
R(t) of the neurons which has been initially boosted to 1 by t
theinput. Thenwemeasurethetimerequiredtoreturntothe
originalvalueandthereforehowlongittakesforthenetwork FIG. 8. Response of the network to external stimuli applied at
toabsorbtheexternalstimulusandcomebacktotheoriginal differenttimesandwithdifferentintensity. Inbothfigures,weplot
dynamicalscenario. InFig.7weplottheevolutionofthepa- the parameter 1−R(t) of the stimulated nodes as a function of
rameter1−R(t)asafunctionoftime,fordifferentf . Atthe timeandweconsidernetworkswithM =500and∆=250. Green
I
beginning,1−R(t)iscloseto0,becausethestimulusroughly (red) shades data refer to fI =0.1(0.214). Upper panel: for each
improves the phase coherence of the stimulated group; then network,wechoosefourdifferenttimestoswitchonastimuluswith
the same intensity, i.e., on 30% of excitatory neurons, randomly
itincreases,inafastorslowwaydependingontheinhibitory
chosen. All data are shifted sothat the onset of the perturbation
fraction whichispresent,anditreturnstotheoriginal value,
is at t= 0. Lower panel: we compare the transient dynamics for
represented by theplateau on theright of Fig. 7.
stimuliappliedon10%,20%,and30%oftheexcitatoryneuronsat
We observe that, when the inhibitory fraction is too small t=0.
or too high, the network quickly destroys the new synchro-
nization imposed bythestimulusandafterabout tenoscilla-
tions the initial state is recovered. In the network with the timerequiredtoreturntotheoriginalsynchronizationlevelis
metastable dynamics, the duration of the perturbation is 20 at least an order of magnitude longer than in other regimes,
times longer. Therefore, in a partially synchronized configu- as in Fig. 7.
ration, the metastable dynamics of the network stores infor-
We verify the robustness of our results showing that the
mation of a synchronousperturbation for long-time scales.
response of the network to a given stimulus is independent
In Fig. 7, the initial growth of the reduced Kuramoto pa- of the specific time the perturbation is switched on. In the
rameter follows a simple power law, namely 1−R(t) ∼ t2, upperpanelofFig.8werandomlychoosefourdifferenttimes
which can beexplained byexpandingR(t) for small t: to apply the same stimulus to the same network, so now the
perturbation does not necessarily occur at the minimum of
R(t)= 1 M eitnkt+−1t−nktnk the excitatory fields YE(t). In all trials we obtain similar
(cid:12)M (cid:12) transient dynamicsand the relaxation times are consistent.
(cid:12) Xk=1 (cid:12) The same qualitative results are also obtained considering
(cid:12) (cid:12)
(cid:12) 1 M t (cid:12) 1 t2 different fractions of excitatory neurons, e.g., 10% or 20%,
≈(cid:12) 1+i (cid:12)−
(cid:12)M tn+1 2(tn+1)2 (cid:12) (lower panel of Fig. 8). We exclude stimuli larger than 50%,
(cid:12) Xk=1(cid:16) k k (cid:17)(cid:12) (17) as these cannot be considered a perturbation of the system
(cid:12) (cid:12)
(cid:12) t2 1 t2 (cid:12) 1 2 dynamics. Thoughdifferentstimuliwillgenerally correspond
≈(cid:12) 1− + (cid:12)
s M(cid:16)X(tnk+1)2(cid:17) M2(cid:16)Xtnk+1(cid:17) trootrdainffsemreinttterressypopnrsoedsuocfedthbeystyhseteemxt,eransatlhienpruetleadseepeonfdnseoun-
t2 1 1 2 the intensity of the input itself, the quantity 1−R(t) refers
≈1− − ,
2 (tn+1)2 tn+1 only to the perturbed neurons so that we can compare data
(cid:16)D k E D k E (cid:17) obtainedindifferentstimulations. Weobservethat,indepen-
where the sum runs only on the nodes which have received dentlyoftheinputintensity,themetastablestatesarealways
thesynchronousstimulus and therefore show tn=0. affected by the external perturbation for times that are an
k
Similar plots can be obtained also for different connec- order of magnitude longer than the other dynamical regimes
tivities of the inhibitory hubs, (we made explicit checks for of partial synchronization.
∆=160 and ∆=500): the prebalance regime is again char- Let us point out that an intriguing result of this analysis
acterized by partially synchronousmetastable states and the is that, in order for the system to work in such an efficient
8
6000 1.4 0.465
(a) (b) (c) (d)
0.02
5000 1.2
0.46
i4000 t()0.015 SIik 1 Zt()E
3000 YE Ih 0.455
0.01 0.8
2000
0.45
0.005 0.6
1000
0 0.4 0.445
0
6547 6548 6549 6550 6551 6552 6547 6548 6549 6550 6551 6552 0 100 200 300 400 0 0.005 0.01 0.015 0.02
6000 1.4 0.465
(a) (b) 1.34 finite size, P(k)
0.02
5000 1.2 1.32 finite size, P(4k)
0.46
i34000000 Yt()E00.0.0115 ISIhik0.18 111..22.368 Zt()E0.455 HMF, P(k) (d)
80 100 120
2000 0.45
0.005 0.6 (c)
1000
0 0 0.4 0.445
6525 6526 6527 6528 6529 6530 6525 6526 6527 6528 6529 6530 0 100 200 300 400 0 0.005 0.01 0.015 0.02
t t k YE(t)
FIG. 9. Comparisonofthe dynamics of anetwork withfI =0.1,simulatedthrough both thefinite sizeand theHMFapproach. Inthe
upperfiguresreddataareproducedevolvingthefinitesizeequationswithN =6000neurons,andbluedatatheHMFoneswithM =500
connectivity classes. Thep.d.f.ischaracterizedinbothcasesbyhkEi=100,∆=200,andσ=10. Inthelowerfiguresweshowhowthe
HMF approximation improves as we approach the limitof ki →∞. Again, green (blue) data evolve according to the finite size (HMF)
approach, butforthefinitesizenetworkwehavequadrupledthedegrees (new p.d.f.parameters: hkEi=400, ∆=800,σ=40). Thanks
tothishigherconnectivity,thefinitesizedynamicsbetterconvergestothesameHMFnetwork. Inpanels(a)wecomparetherasterplots.
As in the HMFfigures M is equal to 500, the HMF indices have been opportunely rescaled for a better overlap with the red and green
data. Data below (above) the black line refer to excitatory (inhibitory) nodes. (b) Field YE(t) as a function of time. The comparison
of YI(t) gives similarresults. (c) hISIki as a function of the degree k. In the lower panel, we divided the degrees of green data by 4, to
overlapfinitesizedatawithHMFpoints. Theinsetisazoom ofthehISIkionlyfortheexcitatory neuronsandforallthreesets ofdata.
(d)AttractorsYE(t)vsZE(t).
regime, the network should be closed to the balance state nitegraphswithdifferentaverageconnectivity,whichareap-
where the connectivity ratio between inhibitory and excita- proximated by the same HMF equations. Finally, fixing the
toryneuronsisequaltotheinverseofthecorrespondingfrac- fraction f , we can consider the various dynamical regimes:
I
tions[seeEq. (15)]. Thisresultdependsonthenatureofthe quasi-periodic, synchronouswith balance, and asynchronous.
modelandonthespecificdistributionofnetworkconnections InFig.9weconsiderf =0.1,i.e.,theregimeofpartialsyn-
I
but nonetheless it is an interesting prediction that might be chronization,whenf is(much)smallerthanitsbalancevalue
I
tested experimentally. fB. The blue data are obtained from the simulation of the
I
HMF equations, while the upper panels refer to hk i = 100
E
and the lower panels to hk i=400. Panels (a) represent the
E
VI. BEYOND THE HMF APPROXIMATION raster plots and panels (b) the synaptic fields YE(t), which
in the finitesize approach are defined as
In the last part of this work, we compare the HMF ap- 1
proximation results with the dynamics on a finite network. YE(t)=+YEE(t)−YEI(t) = L yijki− yijki ,
Although a series of papers have already addressed this is- (cid:16)i,jX∈E,E i,jX∈E,I (cid:17)
(18)
sue, wewill show that in thebalance regime theresult isnot
whereListhetotalnumberoflinks, meansthatthe
trivial and further discussions are required. Accordingly, we i,j∈E,I
sum is restricted to the y resources from an inhibitory node
compare the HMF approach, given by Eqs. (13), with the
P
to an excitatory one and similarly for . Panels (c)
simulation of the Eqs. (1)-(7), where the size is N and the i,j∈E,E
degree ki is extracted from PE(k) and PI(k). Clearly we ex- show the time average hISIki of interspiPke interval as a func-
pecttheresulttoimproveforlargeconnectivitiesofthefinite tion of the connectivity and finally in panels (d) we plot the
graph. Thiscanbeverifiedbyconsideringnetworkswithdif- microscopicattractorsYE(t)vsZE(t)[inthefinitesizegraph
ferent average degrees, which are approximated by the same ZE(t) is defined analogously to Eq. (18), using the resources
HMFequations. Inparticular,weintroduceaclassofrandom zij insteadofyij,whileintheHMFsimulationitiscomputed
graphs, where P (k) and P (k) are both Gaussian distribu- similarly to Eqs. (11) and (12), using z(†,∗) instead of y(†,∗)].
E I k k
tions with averages hk i and hk i =3hk i (∆ =2hk i) and Figure 9 shows that HMF and finite size networks display
E I E E
withstandarddeviationσ=hk i/10forbothexcitatoryand a similar behavior and the main effect of finite connectivity
E
inhibitory neurons. By varying the values of hk i we get fi- is to superimpose on the HMF dynamics a noise, which is
E
9
100 (a)
(a)
0.1
0.214
100 (c) (d)
Rt() (b) t()i 00..2238
−10−2 EE
1 (c) ∆ = 200 ∆Y10
h
∆ = 400 (b)
∆ = 800
100 101 t 102 103 0 20 40 60 80 1t00 120 140 160 180 200
6000
(a) (b)
5000 4000
i
4000 2000
i
3000
0
2000 4000 (c) (d)
i
1000
(a) (b) (c) 2000
0
100 102 104 100 102 104 100 102 104
t 0
0 2 4 t 6 60 62 t 64
FIG. 10. Stability of the synchronous state in finite size simula-
FIG. 11. Dynamicsoffinitesizenetworksaftertheapplicationof
tions. In the upper panel we compare the time evolution of the
an external stimulus at time t = 0 to 30% of excitatory neurons,
parameter 1−R(t) (in logarithmic scale) in three different net-
fiwrosrtknsewtwitohrkN(b=lue60d0a0tan)eiusrcohnasraacntderifzIed=by0.2∆5.=T20h0e,phk.dE.fi.=of1t0h0e, σeraanc=hdo1nm0etlawynodcrhkdoswiffeener,pelwnottithiinnhNtihbe=ituo5rpy0p0ef0rrapncaetuinorenolsn(sifn,I∆s(esm=eeil2ot5gh0ae,rilhtekhgEmeniidc=)s.c1aF0loe0r),
σ = 10; the p.d.f. of the second network is doubled (dark green
data, ∆=400,hkEi=200,σ=20); thep.d.f.ofthelastnetwork itnhgeftreommpotr=al0a.veInracgoervraeslupeosnhd∆enYcEeEo(ftt)hieasblaafcuknacrtrioownso,fitnimthee,sltoawretr-
is quadrupled (light green data, ∆ = 800, hkEi = 400, σ = 40). panelsweshowsnapshotsoftherasterplottakenatdifferenttimes
The external stimulus is applied at time t=0 to all neurons. At
thebeginning1−R(t)isalmostnull,becausethesystemistotally for the network with fI = 0.1 [(a) and (b)] and with fI = 0.214
[(c)and(d)]. Thespikeeventsofstimulated(nonstimulated)nodes
synchronous, butthenduetothefinitesizefluctuations itreturns
arerepresentedbylight(dark)bluedots,whiletheorangeandred
to 1, that is the network evolves towards its stable asynchronous
linesare the signals defined inEq. (19). Tobetter overlapsignals
state. Incorrespondenceoftheblackarrows,weshowinthethree
withrasterplots,weshiftedthemaddingapositiveconstant: this
lower figures the microscopic dynamics of the networks through
is why they are not null before the stimulus. After the stimulus,
theirrasterplots. Theindicesbelow(above)theblacklinereferto
signalsturnonsharplyandthendecay. Theindicesbelow(above)
excitatory(inhibitory)neurons.
theblacklinerefertoexcitatory(inhibitory)neurons.
clearly vanishing by increasing the connectivity of the finite
size graph. Analogous conclusions can be found for all the first network, with lower connectivity, the synchronous con-
dynamical regimes where the fraction of inhibitory neurons figuration is for themost part destroyedand thedynamicsis
is much smaller than the balance value. A similar accor- quicklyreturningtotheasynchronousstate, butin theother
dancecanbeobservedalsointheasynchronousregime,when raster plots,where thedegreeis twoor fourtimes larger, the
f > fB. As the network activity is now aperiodic and the networkstillremainsquitesynchronous,showingthatthelife-
I I
fluctuationsaredominantbothforthesynapticfieldsandfor timeoftheperiodicdynamicsdivergeswith theconnectivity.
theattractors, we donot show thecomparison plots. Thesynchronizationofthesystemcanbemeasuredusingthe
A substantial difference emerges in the balance regime, parameter 1−R(t) and plotting it as a function of time (the
where the HMF approach is characterized by a totally syn- totally synchronous stimulus is applied at t = 0). The be-
chronousdynamics. Indeed,inthiscasetheasymptoticstable haviorof 1−R(t)isillustrated in theupperpanelofFig. 10,
state for the finite size system seems to be asynchronous. In underliningthatthelargertheconnectivityinthenetworkis,
particular, the synchronization in themean-field formulation themore stable thesynchronousstate is.
isaconsequenceofaperfectfieldsubtraction,aswehaveseen Eventually, we show that, even if the HMF and finite size
in Sec. IV, but in networks with finite connectivity fluctua- dynamics have different attractors, the conclusions we have
tions dominate and destroy the phases locking inducing the drawn regarding the high performance for input detection in
synchronization. However,thetotallysynchronousstatechar- thisdynamicalregimearestillvalidalsoinfinitesizesamples.
acterizingtheHMFequationsemergesinthefinitesizesystem Forthispurpose,letusconsidertheresponseofthefinitesize
as a metastable state, whose lifetime increases with the net- networktotheapplicationofanexternalstimulus,inanalogy
work connectivity. In Fig. 10 we consider thebalance regime with the results for the HMF approach shown in Sec. V. Let
f =fB =0.25for∆=200: thelowerpanelshowstheraster us synchronize a fraction equal to S = 0.3 of the excitatory
I I
plots of three different networks, taken from time t=100 af- nodes,followingthesameprocedureoftheHMFsimulations.
ter that a stimulus synchronized the whole dynamics. In the Now,thestimulusperturbationcanbemeasuredthroughthe
10
differencebetweentheexcitatoryfieldsproducedbythenon- control role of highly connected neurons, both in input and
synchronizedandsynchronizednodesrespectively,asfollows: inoutputdirection. Weapplyaheterogeneousmean-fieldap-
proachtothefinitesizenetworkdynamics,thatletsusspeed
∆Y (t)=S Ynot−st(t)−(1−S) Yst (t), (19) upnumericalcomputationsandhighlighttheroleofneuronal
EE EE EE
connections distributions.
where Ynot−st(t) and Yst (t) are computed according to the In this modeling approach, hub inhibitory neurons turn
EE EE
first sum in Eq. (18), limiting it to the excitatory neurons outtobestrongly effectivetodrivethenetworksynchroniza-
which were respectively not stimulated and stimulated. We tion, even in the case in which their relative numberis small
also add the normalization factors S and 1−S, as the first with respect to the excitatory component. Indeed, if their
and the second field are produced by a different number of hub character (their connectivity with respect to that of the
neurons: in this way, before the stimulus, ∆Y (t) fluctu- excitatory population) is high, a small fraction of inhibitory
EE
ates around zeros. The temporal evolution of such variable neurons (fI ∼ 0.1) is able to increase network synchroniza-
is plotted with orange and red lines on the raster plots in tion, leading the majority of neural population (actually the
Figs. 11(a)–11(d). As we have periodic signals, in order to excitatory component) to synchronizeat a certain frequency,
better describe their decays, we compute the average value which results to be equal to that of the isolated LIF neu-
for each period and we plot theresults in theupper panel of ron. Aroundthisregime of completesynchronization wefind
Fig. 11. an interesting metastable dynamical phase. In this dynam-
In the initial dynamical regime of partial synchronization ical regime, excitatory neurons spike in clusters of different
(f = 0.1), the perturbation produced by the stimulus de- size, the composition of each cluster changes in time and it
I
cays fast and after few oscillations the system returns to the is dependent on initial conditions. Interestingly, this dynam-
initial configuration [see the raster plots in Figs. 11(a) and icalregime isthemost appropriatefor storinginformation of
11(b)]. In the metastable and in the totally synchronous an external input received by the network. The hub nature
regime, when the inhibitory fraction is around the balance of inhibitory neurons allows us to control the overall popu-
value (f = 0.214 and f = 0.23), the time required to re- lation dynamics even when the network is mainly composed
I I
turn to the original synchronization level is longer and the by excitatory neurons. In fact, in the case of the absence of
signal decay is slower. Then thedecay rate increases again if inhibitory hubs, i.e., when excitatory and inhibitory neurons
weconsidertheasynchronousstatesafter thebalanceregime have the same average connectivity, the totally synchronous
(f = 0.28 for example). These results point out that the regime appears only at high fraction of inhibitory nodes. On
I
long response time of the system to an external stimulus is the other hand, as soon as the connectivity of the inhibitory
an effect which ispresent not only in theHMFapproach but componentgrowswithrespecttotheexcitatoryone(i.e.,the
alsoinfinitenetworks,thoughinthiscasethebalanceregime inhibitory neurons hub character increases), the metastable
is asymptotically asynchronous. dynamical phaseappears for lower and lower inhibitory frac-
tions. Eventually,for veryhigh valuesof f , thedynamicsof
I
the neural population is completely asynchronous. As a re-
sult,thefractionofinhibitoryneuronsandtheirhubcharacter
VII. CONCLUSION AND PERSPECTIVES
are two crucial ingredients for observing a wide range of dy-
namicalphasesanddrivetheoverallnetworksynchronization.
Wehavestudiedtheroleofinhibitoryhubsforthesynchro- In our approach, the coincidence between input and output
nization and input processing of a neural network with short connectivities for each neuron can have strong effects on the
termplasticity. Wehaveusedasimplecorrelationruleforthe overall dynamics. In future steps, we plan to investigate the
networkgeneration,inwhichinputandoutputconnectivities case in which such a correlation is reduced and hub neurons
arethesameforeachneuron: thisallowsustoemphasizethe are identified only by a high input or outputconnectivity.
[1] E. Niedermeyer and F. L. da Silva, Electroencephalog- bridge universitypress, Cambridge, UK,2003), Vol. 12.
raphy: Basic Principles, Clinical Applications, and Re- [9] V. N. Murthy and E. E. Fetz, J. Neurophysiol. 76, 3968
lated Fields (Lippincott Williams & Wilkins, Philadel- (1996).
phia, 2005). [10] B. Ermentrout, NeuralComput. 8, 979 (1996).
[2] S.A.Huettel,A.W.Song, andG.McCarthy,Functional [11] C. B¨orgers and N. Kopell, Neural Comput. 15, 509
Magnetic Resonance Imaging (Sinauer Associates Inc., (2003).
Sunderland,Massachusetts, 2014), Vol. 1. [12] S. Luccioli, S. Olmi, A. Politi, and A. Torcini, Phys.
[3] P. Gerrard and R. Malcolm, Neuropsychiatr.Dis. Treat. Rev. Lett.109, 138103 (2012).
3, 349 (2007). [13] C. van Vreeswijk and H.Sompolinsky,in Computational
[4] W. Klimesch, Int.J. Psychophysiol. 24, 61 (1996). Neuroscience (SpringerUS,Boston,1997), pp.209–213.
[5] A.VonSteinandJ.Sarnthein,Int.J.Psychophysiol.38, [14] V. Capano, H. J. Herrmann, and L. de Arcangelis, Sci.
301 (2000). Rep. 5 (2015).
[6] J. McAuley and C. Marsden, Brain 123, 1545 (2000). [15] P.-A. Salin and J. Bullier, Physiol. Rev.75, 107 (1995).
[7] V.ShustermanandW.C.Troy,Phys.Rev.E77,061911 [16] M. P. van den Heuvel and O. Sporns, Trends Cogn. Sci.
(2008). 17, 683 (2013).
[8] A.Pikovsky,M.Rosenblum, andJ.Kurths,Synchroniza- [17] P. Bonifazi, M. Goldin, M. A. Picardo, I. Jorquera,
tion: A Universal Concept in Nonlinear Sciences (Cam- A. Cattani, G. Bianconi, A. Represa, Y. Ben-Ari, and