Modeling the Network Dynamics of Pulse-Coupled Neurons

2 MB·

by  Sarthak Chandra

Modeling the Network Dynamics of Pulse-Coupled Neurons
Sarthak Chandra,1 David Hathcock,2 Kimberly Crain,3 Thomas M. Antonsen,1 Michelle Girvan,1 and Edward
Ott1
1)University of Maryland, College Park, Maryland 20742, U.S.A.
2)Case Western Reserve University, Cleveland, Ohio 44016, U.S.A.
3)Iowa State University, Ames, Iowa 50011, U.S.A.
We derive a mean-field approximation for the macroscopic dynamics of large networks of pulse-coupled theta
neurons in order to study the effects of different network degree distributions, as well as degree correlations
(assortativity). UsingtheansatzofOttandAntonsen(Chaos,19(2008)037113),weobtainareducedsystem
of ordinary differential equations describing the mean-field dynamics, with significantly lower dimensional-
ity compared with the complete set of dynamical equations for the system. We find that, for sufficiently
large networks and degrees, the dynamical behavior of the reduced system agrees well with that of the full
network. This dimensional reduction allows for an efficient characterization of system phase transitions and
attractors. Fornetworkswithtightlypeakeddegreedistributions,themacroscopicbehaviorcloselyresembles
7
that of fully connected networks previously studied by others. In contrast, networks with scale-free degree
1
distributions exhibit different macroscopic dynamics due to the emergence of degree dependent behavior of
0
2 different oscillators. For nonassortative networks (i.e., networks without degree correlations) we observe the
presence of a synchronously firing phase that can be suppressed by the presence of either assortativity or
n
disassortativity in the network. We show that the results derived here can be used to analyze the effects of
a
J network topology on macroscopic behavior in neuronal networks in a computationally efficient fashion.

1

InApril2013, theU.S.Presidentannounced‘The that the dimension reduction analyses in Refs.14–16 has
]
D Brain Initiative,’ an extensive, long range plan recentlyprovedtobeveryeffectiveandhasbeenusedto
C of scientific research on human brain function. derive the macroscopic behavior of large systems of cou-
Computer modeling of brain neural dynamics is pled dynamical units in a variety of settings10,11,17–22.
.
n an important component of this long-term overall In particular, Refs.8,10–12 consider network with globally
i
l effort. A barrier to such modeling is the practical coupledneuronsandusethesedimensionreductiontech-
n
limit on computer resources given the enormous niques to analyze the macroscopic behavior of the sys-
[
number of neurons in the human brain (∼ 1011). tems.

1 Our work addresses this problem by developing In1986,ErmentroutandKopell23introducedthetheta
v a method for obtaining low dimensional macro- neuron model. Their work, along with later studies by
2
scopic descriptions for functional groups consist- Ermentrout24 and by Izhikevich25, established the appli-
1
2 ing of many neurons. Specifically, we formulate cability of the theta neuron model for studying networks
0 a mean-field approximation to investigate macro- of Class I excitable neurons (as defined by Hodgkin,26
0 scopic network effects on the dynamics of large i.e., those neurons whose activity lies near the transition
1. systems of pulse-coupled neurons and use the between a resting state and a state of periodic spiking,
0 ansatz of Ott and Antonsen to derive a reduced andcanexhibitspikingwitharbitrarilylowfrequencies).
7 system of ordinary differential equations describ-
Previous studies modeling networks of theta
1 ing the dynamics. We find that solutions of the
neurons8,18,22,27 have generally been restricted to
: reduced system agree with those of the full net-
v particular classes of network topologies. In this paper
work. Thisdimensionalreductionallowsformore
i we study the macroscopic dynamics of networks of pulse
X efficient characterization of system phase transi-
coupled theta neurons on networks with fairly general
r tions and attractors. Our results show the utility
a topologies including arbitrary degree distributions and
of these dimensional reduction techniques for an-
correlations between the degrees of nodes at opposite
alyzing the effects of network topology on macro-
ends of a link, resulting in so-called ‘assortativity’ or
scopic behavior in neuronal networks.
‘disassortativity’28. Assortativity (disassortativity)
occurs when network nodes connect preferentially to
those with similar (different) degrees. We note that,
studies29–33 have shown the biological relevance of
I. INTRODUCTION assortativity. Motivated by the results of Restrepo and
Ott17 on networks of Kuramoto oscillators, we use a
Networks of coupled oscillators have been shown to mean field approach in conjunction with the analytical
haveawidevarietyofbiological,physicalandengineering techniquesdevelopedbyOttandAntonsen14–16 tostudy
applications1–13. In modelling the dynamics of such net- the behavior of pulse coupled theta neurons on networks
works, simulating the microscopic behavior at each node with arbitrary degree distributions and assortativity.
can be a computationally intensive task, especially when We obtain a reduced system of equations describing the
the network is extremely large. In this regard, we note mean-field dynamics of the system, with lower dimen-
2
sionality compared with the complete set of dynamical (a) (b) (c)
equations for the system. This allows us to examine Spiking
the behavior of the network under various conditions Threshold
in a computationally efficient fashion. We primarily Spike
use the example of a scale free degree distribution as Rest state
an application of the obtained dynamical equations
for the order parameter and observe the existence of a
partially resting phase, an asynchronously firing phase,
and a synchronously firing phase that is sensitive to
the presence of assortativity or disassortativity in the FIG. 1. The dynamics of the theta neuron undergo an SNIC
network. We also demonstrate that, in contrast to (SaddlenodeonanInvariantCycle)bifurcationatη=0. For
negativeη theneuronliesinareststate,withathresholdfor
networks with sharply peaked degree distributions,
excitation,andforpositiveηtheoscillatorundergoesperiodic
networks with scale-free degree distributions exhibit
spiking.
different macroscopic dynamics due to the emergence of
degree dependent behavior of different oscillators.
The remainder of this paper is organized as follows.
Thethetaneuronmodelcanbeextendedfromasingle
In Sec. II we describe the model of pulse coupled theta
neuroninisolationtonetworksofneurons. Weconsidera
neurons used on an arbitrary network. In Sec. III setup
system of N theta neurons coupled together in a general
a mean field description of the behavior on the network,
networkviapulse-likesynapticsignals,I ,toeachneuron
i
and then (Sec. IV) show how the methods developed
i:
by Ott and Antonsen14–16 can be used to write a low di-
mensionalsetofequationsdescribingthedynamicsofthe
mean field order parameter. In Sec. V we then use this θ˙ =(1−cosθ )+(1+cosθ )[η +I ], (2)
i i i i i
low dimensional system to describe the behavior of the
N
system under different parameters and network topolo- K (cid:88)
I = A P (θ ), (3)
gies. SectionVIconcludesthepaperwithfurtherdiscus- i (cid:104)k(cid:105) ij n j
j=1
sion and summary of the main result.
where A is the adjacency matrix of a network; A =1
ij ij
if there is a directed edge from node j to node i, and
II. THE MODEL
A = 0 otherwise. The average degree is then given by
ij
(cid:104)k(cid:105)=(cid:80) A /N. P (θ)=d (1−cosθ)n represents the
i,j ij n n
The theta neuron model encodes the dynamics of a pulse-like synapse, whose sharpness is controlled by the
single neuron in isolation as follows, integer parameter n. The normalization constant d is
n
(cid:82)2π
θ˙ =(1−cosθ)+(1+cosθ)η, (1) determinedsothat 0 Pndθ =2π. Notethatinthecase
of a fully connected network, where A =1 for all i and
ij
where θ represents the neuron’s state and the parameter j, this model reduces to that of Luke et al.8
η specifies its excitability. The dynamics can be visual-
ized as a point traveling around the unit circle (Fig. 1).
Aneuronalspikeissaidtooccureachtimethephasean- III. MEAN FIELD FORMULATION
gle of the neuron, θ, crosses the leftmost point at θ =π.
When η < 0, there are two zeros of the right hand side Weconsiderthelimitofmanyneurons,N (cid:29)1,andas-
of Eq. (1), representing a stable rest state (solid circle sumethenetworkisrandomlygeneratedfromagivende-
(cid:80)
in Fig. 1(a)) and an unstable equilibrium (open circle in greedistributionP(k)(normalizedsuchthat P(k)=
k
Fig. 1(a)). Thus,startingfromatypicalinitialcondition, N), where k, the node degree, represents a two-vector of
the state of the neuron goes towards the stable equilib- the in-degree and the out-degree, (k ,k ). Addition-
in out
rium at the rest state represented by the filled circle. A ally, we consider an assortativity function a(k(cid:48) → k),
resting neuron will spike if an external force pushes its which specifies the probability of a link from a node of
state(i.e. theangleθ)fromthereststatepasttheunsta- degreek(cid:48) tooneofdegreek. InthisN →∞limit,weas-
ble equilibrium (termed as the ‘spiking threshold’). As η sumethatthestateoftheneuronscanberepresentedby
is increased above 0, the neuron exhibits a Saddle Node a continuous probability distribution, f(θ,η|k,t), such
bifurcation on an Invariant Cycle (SNIC). In this case that f(θ,η|k,t)dθdη is the probability that a node of
there are no fixed points (i.e. no zeros of the right hand degree k has an excitability parameter in the range
side of Eq. (1)), and the neuron now fires periodically, [η,η+dη]andaphaseangleintherange[θ,θ+dθ]attime
as shown in Fig. 1(c). Note that the neuron does not t. Sinceweareassumingthattheexcitabilityparameters
(cid:82)
move at the same rate along the entire circle, and may donotvarywithtime,wedefineg(η|k)= fdθ,whichis
go faster or slower around θ = π dependent on whether the time independent distribution of the excitability pa-
η is less than or greater than 1, respectively (eq. see the rameters η in the network for a randomly chosen node
i
plot of (1−cosθ) versus time in Fig. 1(c)). of degree k.
3
In order to describe the synchronization behavior of and
this system, we define the order parameter to be34,
(−1)jn!
Q = . (11)
jm 2jm!(n−j)!(j−m)!
N
1 (cid:88)
R(t)= eiθj. (4) If we now assume a Lorentz distribution of the excitabil-
N
j=1 ity parameters,
AsinpreviousworkbyRestrepoandOtt17,wehypothe-
1 ∆(k)
size that in networks with large nodal degrees, the order g(η|k)= , (12)
π[η−η (k)]2+∆2(k)
parameter can be well approximated via a mean field or- 0
der parameter, defined by a continuum version of Eq.
we obtain
(4),
R¯(t)= 1 (cid:88)P(k(cid:48))(cid:90) (cid:90) f(θ(cid:48),η(cid:48)|k(cid:48),t)eiθ(cid:48)dθ(cid:48)dη(cid:48). (5) (cid:90) (cid:90) f(θ(cid:48),η(cid:48)|k,t)eipθdθ(cid:48)dη(cid:48) =ˆb1(,k,t)p, pp=>00
N
k(cid:48) ˆb∗(k,t)|p|, p<0,
(13)
Additionally,thedistributionf isconstrainedbythecon-
with ˆb(k,t) ≡ b(η (k)+i∆(k),k,t). This now allows us
tinuity equation, 0
to rewrite v in terms ofˆb(k,t) as
θ
∂f ∂
+ (v f)=0, (6)
∂t ∂θ θ v =geiθ+h+g∗e−iθ, (14)
θ
where v is the continuous version of the right hand side
θ where
of Eq. (2),
(cid:34) 1 K K
g =− (1−η− H (k,t)), h=1+η+ H (k,t),
v =(1−cosθ)+(1+cosθ) η 2 (cid:104)k(cid:105) n (cid:104)k(cid:105) n
θ
(15)
K (cid:88) and
+d P(k(cid:48))a(k(cid:48) →k)
n(cid:104)k(cid:105) (cid:40)
k(cid:48) (cid:88)
(cid:90) (cid:90) (cid:35) Hn(k,t)=dn P(k(cid:48))a(k(cid:48) →k)
× f(θ(cid:48),η(cid:48)|k(cid:48),t)(1−cosθ(cid:48))ndθ(cid:48)dη(cid:48) . k(cid:48)
(cid:34) n (cid:35)(cid:41)
× A +(cid:88)A (ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p) .
(7) 0 p
p=1
(16)
IV. DIMENSION REDUCTION
Substituting the phase velocity Eq. (14) and the Ott-
Antonsen ansatz Eq. (8) into the continuity equation
Employing the dimensional reduction method of Ott
andAntonsen14–16,andfollowingitspreviousapplication (6), we find that b(η,k,t) satisfies:
to the theta neuron8, we assume that f is given by the
Fourier expansion,
∂b
=i(gb2+hb+g∗). (17)
(cid:40) ∂t
g(η|k)
f(θ,η|k,t)= 1
Insertingtheformsforg andhfromEq. (15)and(16)
2π
(8) into this expression, and evaluating each quantity at the
∞ (cid:41)
+(cid:88)(cid:2)b(η,k,t)pe−ipθ+b∗(η,k,t)peipθ(cid:3) . pole, η = η0(k)+i∆(k), we obtain a reduced system of
equations for ˆb(k,t) describing the mean field dynamics
p=1
of the neuronal network,
We then use the binomial theorem to expand the pulse
(cid:40)
function P (θ) using ∂ˆb(k,t) (ˆb(k,t)−1)2 (ˆb(k,t)+1)2
n =−i + −∆(k)
∂t 2 2
n
(cid:88)
(1−cosθ)n =A0+ Ap[eipθ+e−ipθ], (9) +iη (k)+id K (cid:88)P(k(cid:48))a(k(cid:48) →k)
p=1 0 n(cid:104)k(cid:105)
k(cid:48)
where (cid:34) n (cid:35)(cid:41)
× A +(cid:88)A (ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p) .
n 0 p
(cid:88)
A = δ Q , (10) p=1
p j−2m,p jm
(18)
j,m=0
4
The mean field order parameter, R¯(t) can now be writ- and theˆb’s are given k independent identical initial con-
ten in terms of ˆb(k,t). Using the assumed form for ditions,ˆb(k,0)≡ˆb(0), thenthereareafewnotablecases
f(θ,η|k,t),wecanevaluatetheintegralsinEq. (5)using in which particular degree distributions and our chosen
Cauchy’s residue theorem to obtain assortativity function Eq. (21) allow for further dimen-
sional reduction. For networks with a delta-function de-
R¯(t)= 1 (cid:88)P(k)ˆb(k,t). (19) gree distribution, P(k) = δkin,kδkout,k, the Eq. (18) re-
N duces to a single equation describing the mean field dy-
k
namics,
For the discussion in this paper, we will restrict the
assortativity function to be of the form used previously ∂ˆb(t) (ˆb(t)−1)2 (ˆb(t)+1)2 (cid:40)
by Restrepo and Ott17 =−i + −∆+iη
∂t 2 2 0
(24)
a(k(cid:48) →k)=h(ak(cid:48)→k), (20) (cid:34) n (cid:35)(cid:41)
+id K A +(cid:88)A (ˆb(t)p+ˆb∗(t)p) .
n 0 p
where h(x)=min(max(x,0),1) is defined to ensure that
p=1
a(k(cid:48) →k)isavalidprobability(i.e. 0≤a(k(cid:48) →k)≤1),
and We note that this equation is identical to earlier results
k(cid:48) k (cid:20) (cid:18)k(cid:48) −(cid:104)k(cid:105)(cid:19)(cid:18)k −(cid:104)k(cid:105)(cid:19)(cid:21) for a fully connected network8. Thus, networks with
a = out in 1+c in out , only a single allowed degree have identical asymptotic
k(cid:48)→k N(cid:104)k(cid:105) k(cid:48) k
out in dynamics to a fully connected network. This result is
(21) consistent with analogous results by Barlev et al19 for a
wherecisaparameterusedtovarythenetworkassorta-
network of Kuramoto oscillators. More generally, if the
tivity(withc>0andc<0correspondingtoassortative
network has fixed in-degree, P(k) = P(k )δ , the
and disassortative networks, respectively). In networks out kin,k
systemissimilarlyreducedtothesingledynamicalequa-
with neutral assortativity (c = 0), the probability of
tion, Eq. (24). On the other hand, if the out-degree is
forming a link between two nodes is simply proportional fixed, P(k) = P(k )δ , then dynamics of ˆb(k,t) is
totheout-degreeofthesourcenodeandthein-degreeof in kout,k
independent of k , further reducing the dimensionality
the target node. out
of the problem.
The in-out Pearson assortativity coefficient, r, is a
statistic used to characterize the overall assortativity of
a network, and is defined35 as
Reduction efficiency
(cid:80) [(k(cid:48) −(cid:104)k(cid:105))(k −(cid:104)k(cid:105))]
r = e in out , (22)
(cid:112)(cid:80) (k(cid:48) −(cid:104)k(cid:105))2(cid:112)(cid:80) (k −(cid:104)k(cid:105))2 Equation (18) represents a reduction of the original
e in e out system of N theta neurons to a system with as many
where (cid:80) is the sum over all edges connecting a node equations as there are k values in the support of the de-
e
of degree k(cid:48) to a node of degree k36. Assuming that gree distribution P(k). We denote this quantity by Mk,
a(k(cid:48) → k) = a , and that the in and out degree which, in the case of independent in and out-degree dis-
k(cid:48)→k
distributions are independent, we can relate the assorta- tributions, is equal to Min×Mout, where Min and Mout
tivity coefficient to the parameter c as arethenumberofpossiblein-degreesandout-degreesre-
spectively. In general, simulating the full network, Eq.
c (cid:113) (2), requires O(N2) floating point operations per time
r = ((cid:104)k2 (cid:105)−(cid:104)k(cid:105)2)((cid:104)k2 (cid:105)−(cid:104)k(cid:105)2), (23)
(cid:104)k(cid:105)2 in out step. Using the form of the assortativity function given
in Eq. (21) the sum over k(cid:48) in the reduced system of
which can be seen by noting that the sum of a quan-
equations can be split into two sums, each independent
tity Q(k,k(cid:48)), defined on each edge connecting a node
of k,
of degree k(cid:48) to a node of degree k, over edges in our
mean field formulation would be given by (cid:80) Q(k,k(cid:48))=
(cid:80)k(cid:80)k(cid:48)P(k(cid:48))a(k(cid:48) →k)P(k)Q(k,k(cid:48)). e kin (cid:88)P(k(cid:48))k(cid:48) A+ckout−(cid:104)k(cid:105)(cid:88)P(k(cid:48))(k(cid:48) −(cid:104)k(cid:105))A.
The expression for the assortativity coefficient as a N(cid:104)k(cid:105) out N(cid:104)k(cid:105) in
k(cid:48) k(cid:48)
function of c, Eq. (23), is unbounded, while the Pearson
(25)
assortativityisbydefinitionboundedbetween−1and1. (cid:16) (cid:17)
where A = A +(cid:80)n A ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p . Since
Thisdifferencearisesbecause, forsufficientlylargec, the 0 p=1 p
assortativity function given in Eq. (21) is not a proba- the two sums in Eq. (25) are independent of k, each
bility. However, for the network parameters used in our must be calculated only once per simulation iteration.
numerical example below, we find that Eq. (23) is very Thus, simulating the reduced system Eq. (18) only re-
accurate for |c|≤ 2.5, corresponding to an assortativity quires O(Mk) floating point operations per time step —
range, |r|(cid:46)0.198. Mk operationsperformedonceforeachofthesetwosums
If we assume the excitability parameters are drawn and M operations for each of the ˆb(k,t) equations. In
k
from a degree independent distribution (g(η|k) ≡ g(η)) many cases, M (cid:28) N2, so that simulating Eq. (18) is
k
5
V. NUMERICAL SIMULATIONS AND RESULTS
In the following examples, we consider a directed net-
work of N = 5000 nodes, with in and out degrees cho-
sen from independent, identical heavy-tailed distribu-
tions given by

0 if k <k
 min
P(k)= Ak−γ if k ≤k <k (26)
min max
0 if k ≤k.
max
Theexponentofthepowerlawdistribution,γ,wassetto
3, and k and k were set to 750 and 2000, respec-
min max
tively. As mentioned earlier, the normalization constant
(cid:80)
A is chosen to make P(k) = N. We will also set
k
the parameter n controlling the sharpness of the synap-
ticpulseto2forallexamplesconsidered,andwillusean
FIG. 2. The effect of varying levels of interpolation on the
interpolation level of 10% for all calculations using the
calculated results for the trajectories of R¯(t) in the complex
reduced system of equations for the mean field theory
planestartingfromaninitialconditionofR¯(t)=0andending
(cf. Fig. 2).
atafixedpointattractorforK =3inanetworkwithneutral
From numerical simulations of the reduced equations,
assortativity, with η = −2 and ∆ = 0.1. Calculation of
0
the order parameter dynamics is robust to a large range in (18),wefindthatthelongtermdynamicsoftheorderpa-
the level of interpolation. Using as few as 10% of the total rametercanbebroadlyclassifiedintooneofthreephases
available degrees and interpolating the remaining 90% give – (1) the partially resting (PR) phase; (2) the asyn-
results close to the calculation without interpolation. In the chronously firing (AF) phase; and (3) the synchronously
rest of this paper we employ a 10% interpolation level in all firing (SF) phase. The PR phase and the AF phase ap-
ourmeanfieldcalculations. Theblackarcisasegmentofthe pear as fixed points in the dynamics of the order param-
unit circle |R¯(t)|=1.
eter,whereastheSFphaseappearsasalimitcycleofthe
order parameter.
A. Fixed points
significantly more efficient that simulating the full net-
work. Furthermore, if c is set to 0, which is the case of
As a particular example to illustrate the different
networkswithneutralassortativity,thenˆb(k,t)willhave types of fixed points, we look at a network with neu-
no dependence on kout, and hence the overall problem tral assortativity (c = 0) having excitability parameters
is reduced to Min independent equations, allowing even distributed according to a Lorentzian distribution with
greater computational efficiency. mean η =−2 and width ∆=0.1 (Fig. 3).
0
When the network is in the PR phase, the order pa-
Sinceˆb(k(cid:48),t), P(k(cid:48)), and a(k(cid:48) →k) are each smoothly rameter goes to a fixed point that lies near the edge of
varyingfunctions,wecanachievefurtherdimensionalre- the unit circle |R¯|= 1. In this phase, most of the indi-
duction by interpolating the summand in Eq. (18) us- vidual neurons in the network are independently in their
ing a coarse-grained grid of k values. In particular, Eq. resting states, in a fashion similar to Fig. 1(a). This
(18) is not solved for ˆb(k,t) for all of the M values of corresponds to the case of K = 1 in Fig. 3(a), in which
k
k, but only for the small subset of k values that lie on the fixed point is located near the edge of the unit cir-
thecoarse-grainedgridink-spaceThesummandsonthe cle marked in black. Further, the time series of a few
right hand side of Eq. (18) at k values not on the grid randomly chosen neurons (Fig. 3(b)) demonstrates that
areapproximatedbyabilinearinterpolationofthevalues almost all of the neurons are in a resting state. While
at the surrounding chosen k values. To perform the bi- there may be a small number of neurons that are in the
linear interpolation, we first interpolate linearly between spiking phase due to the spread in the distribution of
neighboring grid values in one direction. The value of values of excitability parameters, η, these do not have
the summand at a given k value is then approximated any significant effect on the full order parameter of the
by linearly interpolating in the other direction between system.
values estimated with the previous linear interpolation. As we increase the coupling constant K, the system
We find that using as few as 10% of the network degrees transitions to the asynchronously firing (AF) phase, in
yield very accurate results, while an even coarser inter- which the order parameter goes to a fixed point located
polation still produces the same qualitative behavior as near the center of the unit circle. In this phase, most of
can be seen in Fig. 2. theindividualneuronsinthenetworkareasynchronously
6
FIG. 3. (a): Fixed points of R(t) observed in networks with neutral assortativity, η = −2 and ∆ = 0.1, for three values of
0
the coupling strength K. Fixed points in the PR state (K =1) and the AF state (K =6) are marked in the complex plane.
The fixed point at an intermediate value of K is also marked.(b),(c),(d): Time series of the cosine of the phase of 5 randomly
chosenneuronsdemonstratesthatinthePRphasealmostallneuronsareinarestingstate,andasthesystemapproachesthe
AF state, more nodes transition to an oscillating, excited state.The thick dashed line corresponds to the position of the fixed
point of the order parameter for the corresponding value of K.
firing, in a fashion similar to Fig. 1(c). This can be
seen in the case of K = 6 in Fig. 3(c) which shows
that almost all of the neurons are in a recurrent spiking
state. Note that even though the neurons are spiking
asynchronously, i.e., their firing times are independent
of one another37, the fixed point of the order parameter
is not at R¯ = 0. This is because the angular velocity
of an individual neuron is not constant along the circle,
thus in the average over the ensemble of neurons a bias
is present towards the direction for which the angular
velocityofneuronsisminimized. AsdiscussedinSec. II,
thismayoccurateitherθ =0oratθ =π, dependenton
how large the excitability parameter is for the neuron.
We now examine the transition from the PR phase to
the AF phase. Microscopically, in the PR phase, almost
all of the neurons are individually in a resting phase,
whereas in the AF phase almost all neurons are in the
spiking state. To examine the behavior at an interme- FIG. 4. Comparison of |ˆb(k )| from the reduced system
in
diate point, we look at the fixed point for the case of
of equations and the time average of |R(k )| from the full
in
K =3, asshowninFig. 3(c). Atthisintermediatevalue
system, Eq. (27), for a network with neutral assortativity
of the coupling constant, a fraction of the neurons are in (c=0), η =−2, and ∆=0.1 at K =3. The dynamics un-
0
the spiking state. In particular, the nodes that begin to der these parameters were simulated in a network with 5000
spike first are those which have larger in-degrees. This nodes,andthenetworkwasallowedtorelaxtoafixedpoint.
is demonstrated in Fig. 4, in which we examine ˆb(k) at Nodesweredividedintoclassesaccordingtotheirin-degreeto
the fixed point for K =3. Since we are looking at a net- calculatethetimeaveragedeffectiveorderparameterforeach
class,whichisshowninblue,withtheerrorbarsdenotingthe
work with neutral assortativity (c=0), Eq. (25) implies
rootmeansquaredtimefluctuationoftheorderparameterfor
that the sum only depends on the out-degree through a
thatclass. Thetimefluctuationsareduetothefinitenumber
common multiplicative factor. Thus ˆb is only plotted as
of nodes in each class. (See text for details.)
a function of k . Analogously, for the fixed point of the
in
dynamics on the full network, the range of degrees from
k to k is divided uniformly into several intervals,
min max one of the intervals of the range of degrees, ||N|| is the
and for each interval we find a partial order parameter,
number of nodes in the set, and k is the average in-
calculated such that the average in Eq. (4) is only per- in
degree of nodes within that set.
formed over those nodes whose in-degree lie within that
interval, i.e., In addition, we find that the transition from the PR
phase to the AF phase occurs via a hysteretic process
R(k ,t)= 1 (cid:88)eiθj, (27) mediated by saddle node bifurcations. To illustrate this,
in ||N|| we evolved the dynamics of the full network in a step
j∈N
wise fashion by increasing the coupling constant K in
where N is the set of nodes having an in-degree within smallincrementsof0.2, andallowingthesystemtorelax
7
to an equilibrium before the next increment (Fig. 5(a)). we calculate the associated z(k) andˆb(k) using Eq. (30)
Wealsocomparethiswiththeanalogoushysteresiscurve and Eq. (28), and then recalculate new values, X and
1
observed for the evolution of the system dynamics on an Y using Eq. (31). For fixed points of the reduced equa-
1
Erd˝os-R´enyi network having the same size and average tionsδX =X −X andδY =Y −Y arebothzero. We
1 0 1 0
degree as the scale free network being considered (Fig. calculate δX and δY for several different initial values
5(b)). While the hysteretic region begins at around the at regularly spaced intervals for X and Y , and iden-
0 0
samevalueofthecouplingconstant,K,forbothnetwork tify the fixed points as the points where δX = δY = 0.
topologies, we find that for the case of the Erd˝os-R´enyi Theinterpolationproceduredescribedearliercanalsobe
network,whichhasasharplypeakeddegreedistribution, applied to this calculation to further increase efficiency.
the range in K that allows hysteresis (3 (cid:46) K (cid:46) 7.25) is For the nonassortative case (c = 0), Y = 0 always, so
significantly larger than the corresponding range for the identifying the fixed points in this case only requires cal-
network with the scale free degree distribution (3.25 (cid:46) culatingthevariationinthesingleparameterX. Weuse
K (cid:46)4). this method to evaluate the fixed points of the reduced
To compare with the simulation of the dynamics on equations for the range of K over which hysteresis was
the full network, we also calculate the fixed points of the observed,andfindcloseagreementbetweentheresultsof
mean field equations Eq. (18). While the fixed points this fixed point analysis and the direct evolution of the
cannot be readily determined analytically, we can effi- full network (Fig. 5).
ciently compute them via a numerical calculation. Set-
ting ∂ˆb(k,t)/∂t=0 for the fixed points, we find that the
equilibriumˆb(k) satisfy, B. Limit Cycles
As a representative example of limit cycles of R¯(t), we
1±z(k)
ˆb±(k)= 1∓z(k), (28) consider a network with neutral assortativity with ex-
citability parameters η distributed as a Lorentzian with
where meanη0 =10.75andwidth∆=0.5,andwithacoupling
constant K =−9. In the SF phase, the order parameter
goes to a limit cycle in the complex plane. In this phase,
K (cid:88)
iz2(k)=−∆+iη0+idn(cid:104)k(cid:105) P(k(cid:48))a(k(cid:48) →k) a majority of the neurons are synchronously in a spiking
k(cid:48) (29) state. Plots for such limit cycles are shown in Fig. 6,
(cid:34) n (cid:35) in which we plot the trajectory of the order parameter
× A0+(cid:88)Ap(ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p) , in the complex plane (after removing transients) for a
p=1 network with the scale free degree distribution given in
Eq. (26)(bluesolidcurve),acorrespondingErd˝os-R´enyi
andthesignischosentoensure|ˆb(k)|≤1. Usingourform network having a Poissonian degree distribution (green
oftheassortativityfunctionEq. (21),wemayagainsplit dashedcurve),andaregularnetworkhavingadeltafunc-
the above sum into two parts as in Eq. (25). Thus we tion degree distribution (i.e. P(k) = δ δ ) (red
may rewrite Eq. (29) as dotted curve), each having the same avekrina,gkekdoeugt,rkee. As
seen earlier in Eq. (24), a network with a delta function
degree distribution has mean field dynamics identical to
iz2(k)=−∆+iη +ik X+i(k −(cid:104)k(cid:105))Y, (30)
0 in out those of a fully connected network, and the correspond-
ing limit cycle in Fig. 6 is identical to the limit cycle
where X and Y are given by,
obtainedattheseparametersforthefullyconnectednet-
(cid:34) work by Luke et al.8 In comparison with the limit cy-
X =d K (cid:88)P(k(cid:48))k(cid:48) A clesthatareobservedforthecaseoftheregularnetwork
nN(cid:104)k(cid:105)2 out 0 or the Erd˝os-R´enyi network, the limit cycles in networks
k(cid:48) (31a)
n (cid:35) withscalefreedegreedistributionsarediminishedinsize,
+(cid:88)A (ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p) duetothelargevariationinnodalbehaviorasafunction
p
of degree. Nodes with smaller in-degrees were observed
p=1
topredominantlybeinthespikingphase, withhighsyn-
(cid:34) chronizationandalargerlimitcycleforthepartialorder
K (cid:88)
Y =d P(k(cid:48))(k(cid:48) −(cid:104)k(cid:105)) A parameter, whereas nodes with larger in-degrees were in
nN(cid:104)k(cid:105)2 in 0 the resting phase. Due to this differentiation of behavior
k(cid:48) (31b)
n (cid:35) withdegree,theaveragedfullorderparameterexhibitsa
+(cid:88)A (ˆb(k(cid:48),t)p+ˆb∗(k(cid:48),t)p) . limit cycle that is somewhat reduced in size when com-
p
pared with the results for a fully connected network by
p=1
Luke et al8. However, we see that the limit cycles for
Thesesimplificationsallowforefficientcalculationofthe the Erd˝os-R´enyi network are similar in shape and struc-
system fixed points. Choosing initial values, X and Y , ture to the limit cycles obtained for the regular network,
0 0
8
FIG. 5. A sweeping value of K was used to observe the change in phase from the PR state to the AF state. Hysteresis was
observed on the network with a scale free degree distribution (a) as well as a corresponding Erdo˝s-R´enyi network having the
same size and the same average degree (b). For the full network, at each value of K the mean of the order parameter after
ignoringthetransientshavebeenmarkedastriangles. Aclosematchisobservedwiththefixedpointsascomputedfrommean
field equations directly (see text for details). Hysteresis is observed for 3.25 (cid:46) K (cid:46) 4 in the scale free network (a), and is
observedfor3(cid:46)K (cid:46)7.25inthecorrespondingErdo˝s-R´enyinetwork(Notethedifferenceinscalesforthex-axisinbothplots).
An apparent crossing of the fixed point curve is seen in (b), which is an artifact of the non-self-intersecting R¯ curve lying in
the two dimensional complex space, which has been projected onto the real axis in this plot.
as would be expected in accordance with the discussion VI. CONCLUSION
in Sec. IV, since the Poissonian degree distribution for
theErd˝os-R´enyinetworkissharplypeakedabouttheav-
Usingameanfieldapproximation,inconjunctionwith
erage degree and hence cannot admit a large variation
the Ott-Antonsen ansatz, we obtained a reduced system
of behavior with nodal degree. As the average degree,
of equations that successfully model the macroscopic or-
(cid:104)k(cid:105) increases, the red and green curves converge because
der parameter dynamics of a large network of theta neu-
thePoissondegreedistributionappropriateforanErd˝os-
rons. This reduced system of equations allows us to ex-
R´enyi network approaches a delta function.
amine the effects of varying the network parameters and
the network topology (in terms of degree distributions,
as well as degree correlations) in a computationally effi-
cientfashion. Theorderparameterofthenetworkisused
fordescribingthemacroscopicbehaviorofthenetworkof
theta neurons, whose attractors can be of various types.
C. Effect of Assortativity In particular, we find resting states, asynchronously fir-
ing states and synchronously firing states, the first two
of which appear as a fixed point for the order param-
Wenowconsidertheeffectofassortativityonthelimit eter (Fig. 3), while the third appears as a limit cycle
cycle dynamics of the order parameter in the network. for the order parameter (Fig. 6). We also used the re-
While limit cycle behavior exists in networks with neu- ducedsystemofequationstoobservetheeffectofvarying
tralassortativity(c=0),introductionofassortativityor the assortativity in the system and demonstrated that a
dissasortativity in the network can cause the limit cycle synchronouslyfiringphasewasonlypresentfornetworks
attractortotransformtoafixedpointattractor(AFlike with neutral or small assortativity, and the addition of
state) via a Hopf bifurcation. This is demonstrated in moderate amounts of assortativity or disassortativity to
Fig. 7,inwhichweshowthatvaryingcawayfromzeroto thenetworkcausesthesystemtogotoanasynchronously
±2.5 (corresponding to Pearson assortativity coefficients firing state instead (Fig. 7). Further, for networks with
of r ≈±0.198) is sufficient to cause the Hopf bifurcation scale free degree distributions, we find that nodes with
andsendthesystemtoafixedpointattractor. Thefixed differentvaluesoftheirdegreesadmitalargevariationof
pointsfortheorderparametersinthesenetworksexhibit behavior(Fig. 4),aphenomenonnotpossibleinnetworks
relativelylargeamountsoffiniteN inducednoiseasseen with all-to-all connectivity. In all cases close agreement
from the size of the clouds surrounding the fixed point was observed between the order parameter dynamics as
position calculated from the reduced system. predicted by the reduced system of equations (Eq. 18),
9
8T. B. Luke, E. Barreto, and P. So, “Complete classification of
the macroscopic behavior of a heterogeneous network of theta
neurons,”NeuralComputation25,3207–3234(2013).
9M.M.AbdulrehemandE.Ott,“Lowdimensionaldescriptionof
pedestrian-inducedoscillationofthemillenniumbridge,”Chaos
19,013129(2009).
10E.Montbri´o, D.Pazo´, andA.Roxin,“Macroscopicdescription
for networks of spiking neurons,” Physical Review X 5, 021028
(2015).
11D.Paz´oandE.Montbri´o,“Low-dimensionaldynamicsofpopula-
tionsofpulse-coupledoscillators,”PhysicalReviewX4,011009
(2014).
12C.R.Laing,“Derivationofaneuralfieldmodelfromanetwork
ofthetaneurons,”PhysicalReviewE90,010901(2014).
13Z. Lu, K. Klein-Carden˜a, S. Lee, T. M. Antonsen, M. Girvan,
and E. Ott, “Resynchronization of circadian oscillators and the
east-westasymmetryofjet-lag,”Chaos26,094811(2016).
14E.OttandT.M.Antonsen,“Lowdimensionalbehavioroflarge
systemsofgloballycoupledoscillators,”Chaos18,037113(2008).
15E. Ott and T. M. Antonsen, “Long time evolution of phase os-
cillatorsystems,”Chaos19,023117(2009).
16E. Ott, B. R. Hunt, and T. M. Antonsen Jr, “Comment on
FIG. 6. Comparison of the limit cycle attractor for R¯(t) in Longtimeevolutionofphaseoscillatorsystems[Chaos19,023117
(2009)],”Chaos21,025112(2011).
the complex plane across varying degree distributions in a
17J. G. Restrepo and E. Ott, “Mean-field theory of assortative
network with neutral assortativity (c = 0) with η = 10.75,
0 networks of phase oscillators,” Europhysics Letters 107, 60006
∆ = 0.5 and K = −9. The scale free network (blue solid
(2014).
curve) has a degree distribution according to Eq. (26), 18E.A.Martens,E.Barreto,S.Strogatz,E.Ott,P.So, andT.An-
the Erdo˝s-R´enyi network (green dashed curve) has a Poisso- tonsen, “Exact results for the kuramoto model with a bimodal
niandegreedistribution,andtheregularnetwork(reddotted frequencydistribution,”PhysicalReviewE79,026204(2009).
curve) has a delta function degree distribution. The black 19G.Barlev,T.M.Antonsen, andE.Ott,“Thedynamicsofnet-
circle is the unit circle |R¯|=1 work coupled phase oscillators: An ensemble approach,” Chaos
21,025103(2011).
20P. S. Skardal, J. G. Restrepo, and E. Ott, “Frequency assorta-
tivitycaninducechaosinoscillatornetworks,”PhysicalReview
andascalculatedbyevolutionofthefullsystemofequa- E91,060902(2015).
tions Eq. (2). 21D. Paz´o and E. Montbri´o, “From quasiperiodic partial synchro-
nizationto collective chaos inpopulationsofinhibitory neurons
withdelay,”PhysicalReviewLetters116,238101(2016).
22J.RouletandG.B.Mindlin,“Averageactivityofexcitatoryand
inhibitoryneuralpopulations,”Chaos26,093104(2016).
ACKNOWLEDGEMENTS
23G. B. Ermentrout and N. Kopell, “Parabolic bursting in an ex-
citablesystemcoupledwithaslowoscillation,”SIAMJournalon
This work was supported by the Army Research Of- AppliedMathematics46,233–253(1986).
fice under Grant No. W911NF-12-1-0101, and by the 24B.Ermentrout,“Typeimembranes,phaseresettingcurves,and
synchrony,”NeuralComputation8,979–1001(1996).
National Science Foundation under Grant No. PHY-
25E. M. Izhikevich, “Class 1 neural excitability, conventional
1461089. synapses, weakly connected networks, and mathematical foun-
dationsofpulse-coupledmodels,”IEEETransactionsonNeural
Networks10,499–507(1999).
1D. C. Michaels, E. P. Matyas, and J. Jalife, “Mechanisms of 26A.L.Hodgkin,“Thelocalelectricchangesassociatedwithrepeti-
sinoatrialpacemakersynchronization: anewhypothesis.”Circu- tiveactioninanon-medullatedaxon,”TheJournalofPhysiology
lationResearch61,704–714(1987). 107,165(1948).
2K. Wiesenfeld, P. Colet, and S. H. Strogatz, “Frequency lock- 27C.B¨orgersandN.Kopell,“Synchronizationinnetworksofexci-
ing in josephson arrays: connection with the kuramoto model,” tatoryandinhibitoryneuronswithsparse,randomconnectivity,”
PhysicalReviewE57,1563(1998). Neuralcomputation15,509–538(2003).
3I. Z. Kiss, Y. Zhai, and J. L. Hudson, “Emerging coherence 28M. E. Newman, “Assortative mixing in networks,” Physical re-
inapopulationofchemicaloscillators,”Science296,1676–1678 viewletters89,208701(2002).
(2002). 29P.Hagmann,L.Cammoun,X.Gigandet,R.Meuli,C.J.Honey,
4A.E.Motter,S.A.Myers,M.Anghel, andT.Nishikawa,“Spon- V. J. Wedeen, and O. Sporns, “Mapping the structural core of
taneous synchrony in power-grid networks,” Nature Physics 9, humancerebralcortex,”PLoSBiology6,e159(2008).
191–197(2013). 30S.BialonskiandK.Lehnertz,“Assortativemixinginfunctional
5B. A. Carreras, V. E. Lynch, I. Dobson, and D. E. Newman, brain networks during epileptic seizures,” Chaos 23, 033139
“Complexdynamicsofblackoutsinpowertransmissionsystems,” (2013).
Chaos14,643–652(2004). 31E.Barzegaran,A.Joudaki,M.Jalili,A.O.Rossetti,R.S.Frack-
6L.GlassandS.A.Kauffman,“Thelogicalanalysisofcontinuous, owiak, and M. G. Knyazeva, “Properties of functional brain
non-linearbiochemicalcontrolnetworks,”JournalofTheoretical networks correlate with frequency of psychogenic non-epileptic
Biology39,103–129(1973). seizures,”FrontiersinHumanNeuroscience6(2012).
7M.AldanaandP.Cluzel,“Anaturalclassofrobustnetworks,” 32W.deHaan,Y.A.Pijnenburg,R.L.Strijers,Y.vanderMade,
ProceedingsoftheNationalAcademyofSciences100,8710–8714 W. M. van der Flier, P. Scheltens, and C. J. Stam, “Func-
(2003).
10
FIG. 7. For the parameters η =4, ∆=0.5 and K =−4.8, in a network with neutral assortativity (c=0), the system lies in
0
an SF state (as in (b)). Varying the assortativity in either direction (c=±2.5, corresponding to r≈±0.198) causes the limit
cycles to be replaced by a fixed point instead (as in (a),(c)). The behavior predicted by the mean field theory is in agreement
with the simulations of the full network.
tional neural network analysis in frontotemporal dementia and 36Foranother,oftenuseful,definitionofacoefficientquantitatively
alzheimer’s disease using eeg and graph theory,” BMC Neuro- characterizing the assortativity or disassortativity of a network
science10,1(2009). seeRef.38.
33S.Teller,C.Granell,M.DeDomenico,J.Soriano,S.G´omez, and 37Thisdefinitionofasynchronousspikingisconsistentwithremarks
A.Arenas,“Emergenceofassortativemixingbetweenclustersof by other authors39,40, wherein asynchronous states have been
cultured neurons,” PLoS Computational Biology 10, e1003796 defined as states in which at each neuron the term coupling it
(2014). totheotherneuronsinthenetworkisindependentoftime,asis
34Some authors, such as Restrepo and Ott17 define the order pa- observedinthecasesoffixedpoints.
rameter differently so as to be weighted with the out-degree at 38J. G. Restrepo, E. Ott, and B. R. Hunt, “Approximating the
eachnode,i.e.,R(t)=(cid:80)Ni=1(cid:80)Nj=1Aijeiθj/(cid:16)(cid:80)Ni=1(cid:80)Nj=1Aij(cid:17). largest eigenvalue of network adjacency matrices,” Physical Re-
viewE76,056119(2007).
35J. G. Foster, D. V. Foster, P. Grassberger, and M. Paczuski, 39L. Abbott and C. van Vreeswijk, “Asynchronous states in net-
“Edge direction and the structure of networks,” Proceedings of
worksofpulse-coupledoscillators,”PhysicalReviewE48,1483
theNationalAcademyofSciences107,10815–10820(2010).
(1993).
40D. Hansel and G. Mato, “Existence and stability of persistent
states in large neuronal networks,” Physical Review Letters 86,
4175(2001).