Modeling the Network Dynamics of Pulse-Coupled Neurons
Sarthak Chandra, David Hathcock, Kimberly Crain, Thomas M. Antonsen, Michelle Girvan, Edward Ott
MModeling the Network Dynamics of Pulse-Coupled Neurons
Sarthak Chandra, David Hathcock, Kimberly Crain, Thomas M. Antonsen, Michelle Girvan, and EdwardOtt University of Maryland, College Park, Maryland 20742, U.S.A. Case Western Reserve University, Cleveland, Ohio 44016, U.S.A. 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 thetaneurons in order to study the effects of different network degree distributions, as well as degree correlations(assortativity). Using the ansatz of Ott and Antonsen (
Chaos , (2008) 037113), we obtain a reduced systemof 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 sufficientlylarge networks and degrees, the dynamical behavior of the reduced system agrees well with that of the fullnetwork. This dimensional reduction allows for an efficient characterization of system phase transitions andattractors. For networks with tightly peaked degree distributions, the macroscopic behavior closely resemblesthat of fully connected networks previously studied by others. In contrast, networks with scale-free degreedistributions exhibit different macroscopic dynamics due to the emergence of degree dependent behavior ofdifferent oscillators. For nonassortative networks (i.e., networks without degree correlations) we observe thepresence of a synchronously firing phase that can be suppressed by the presence of either assortativity ordisassortativity in the network. We show that the results derived here can be used to analyze the effects ofnetwork topology on macroscopic behavior in neuronal networks in a computationally efficient fashion. In April 2013, the U.S. President announced ‘TheBrain Initiative,’ an extensive, long range planof scientific research on human brain function.Computer modeling of brain neural dynamics isan important component of this long-term overalleffort. A barrier to such modeling is the practicallimit on computer resources given the enormousnumber of neurons in the human brain ( ∼ ).Our work addresses this problem by developinga method for obtaining low dimensional macro-scopic descriptions for functional groups consist-ing of many neurons. Specifically, we formulatea mean-field approximation to investigate macro-scopic network effects on the dynamics of largesystems of pulse-coupled neurons and use theansatz of Ott and Antonsen to derive a reducedsystem of ordinary differential equations describ-ing the dynamics. We find that solutions of thereduced system agree with those of the full net-work. This dimensional reduction allows for moreefficient characterization of system phase transi-tions and attractors. Our results show the utilityof these dimensional reduction techniques for an-alyzing the effects of network topology on macro-scopic behavior in neuronal networks. I. INTRODUCTION
Networks of coupled oscillators have been shown tohave a wide variety of biological, physical and engineeringapplications . In modelling the dynamics of such net-works, simulating the microscopic behavior at each nodecan be a computationally intensive task, especially whenthe network is extremely large. In this regard, we note that the dimension reduction analyses in Refs. hasrecently proved to be very effective and has been used toderive the macroscopic behavior of large systems of cou-pled dynamical units in a variety of settings .In particular, Refs. consider network with globallycoupled neurons and use these dimension reduction tech-niques to analyze the macroscopic behavior of the sys-tems.In 1986, Ermentrout and Kopell introduced the thetaneuron model. Their work, along with later studies byErmentrout and by Izhikevich , established the appli-cability of the theta neuron model for studying networksof Class I excitable neurons (as defined by Hodgkin, i.e., those neurons whose activity lies near the transitionbetween a resting state and a state of periodic spiking,and can exhibit spiking with arbitrarily low frequencies).Previous studies modeling networks of thetaneurons have generally been restricted toparticular classes of network topologies. In this paperwe study the macroscopic dynamics of networks of pulsecoupled theta neurons on networks with fairly generaltopologies including arbitrary degree distributions andcorrelations between the degrees of nodes at oppositeends of a link, resulting in so-called ‘assortativity’ or‘disassortativity’ . Assortativity (disassortativity)occurs when network nodes connect preferentially tothose with similar (different) degrees. We note that,studies have shown the biological relevance ofassortativity. Motivated by the results of Restrepo andOtt on networks of Kuramoto oscillators, we use amean field approach in conjunction with the analyticaltechniques developed by Ott and Antonsen to studythe behavior of pulse coupled theta neurons on networkswith arbitrary degree distributions and assortativity.We obtain a reduced system of equations describing themean-field dynamics of the system, with lower dimen- a r X i v : . [ n li n . C D ] J a n sionality compared with the complete set of dynamicalequations for the system. This allows us to examinethe behavior of the network under various conditionsin a computationally efficient fashion. We primarilyuse the example of a scale free degree distribution asan application of the obtained dynamical equationsfor the order parameter and observe the existence of apartially resting phase, an asynchronously firing phase,and a synchronously firing phase that is sensitive tothe presence of assortativity or disassortativity in thenetwork. We also demonstrate that, in contrast tonetworks with sharply peaked degree distributions,networks with scale-free degree distributions exhibitdifferent macroscopic dynamics due to the emergence ofdegree dependent behavior of different oscillators.The remainder of this paper is organized as follows.In Sec. II we describe the model of pulse coupled thetaneurons used on an arbitrary network. In Sec. III setupa mean field description of the behavior on the network,and then (Sec. IV) show how the methods developedby Ott and Antonsen can be used to write a low di-mensional set of equations describing the dynamics of themean field order parameter. In Sec. V we then use thislow dimensional system to describe the behavior of thesystem under different parameters and network topolo-gies. Section VI concludes the paper with further discus-sion and summary of the main result. II. THE MODEL
The theta neuron model encodes the dynamics of asingle neuron in isolation as follows,˙ θ = (1 − cos θ ) + (1 + cos θ ) η, (1)where θ represents the neuron’s state and the parameter η specifies its excitability. The dynamics can be visual-ized as a point traveling around the unit circle (Fig. 1).A neuronal spike is said to occur each time the phase an-gle of the neuron, θ , crosses the leftmost point at θ = π .When η <
0, there are two zeros of the right hand sideof Eq. (1), representing a stable rest state (solid circlein Fig. 1(a)) and an unstable equilibrium (open circle inFig. 1(a)). Thus, starting from a typical initial condition,the state of the neuron goes towards the stable equilib-rium at the rest state represented by the filled circle. Aresting neuron will spike if an external force pushes itsstate (i.e. the angle θ ) from the rest state past the unsta-ble equilibrium (termed as the ‘spiking threshold’). As η is increased above 0, the neuron exhibits a Saddle Nodebifurcation on an Invariant Cycle (SNIC). In this casethere are no fixed points (i.e. no zeros of the right handside of Eq. (1)), and the neuron now fires periodically,as shown in Fig. 1(c). Note that the neuron does notmove at the same rate along the entire circle, and maygo faster or slower around θ = π dependent on whether η is less than or greater than 1, respectively (eq. see theplot of (1 − cos θ ) versus time in Fig. 1(c)). Spiking ThresholdRest state
Spike (a) (b) (c)
FIG. 1. The dynamics of the theta neuron undergo an SNIC(Saddle node on an Invariant Cycle) bifurcation at η = 0. Fornegative η the neuron lies in a rest state, with a threshold forexcitation, and for positive η the oscillator undergoes periodicspiking. The theta neuron model can be extended from a singleneuron in isolation to networks of neurons. We consider asystem of N theta neurons coupled together in a generalnetwork via pulse-like synaptic signals, I i , to each neuron i : ˙ θ i = (1 − cos θ i ) + (1 + cos θ i )[ η i + I i ] , (2) I i = K (cid:104) k (cid:105) N (cid:88) j =1 A ij P n ( θ j ) , (3)where A ij is the adjacency matrix of a network; A ij = 1if there is a directed edge from node j to node i , and A ij = 0 otherwise. The average degree is then given by (cid:104) k (cid:105) = (cid:80) i,j A ij /N . P n ( θ ) = d n (1 − cos θ ) n represents thepulse-like synapse, whose sharpness is controlled by theinteger parameter n . The normalization constant d n isdetermined so that (cid:82) π P n dθ = 2 π . Note that in the caseof a fully connected network, where A ij = 1 for all i and j , this model reduces to that of Luke et al. III. MEAN FIELD FORMULATION
We consider the limit of many neurons, N (cid:29)
1, and as-sume the network is randomly generated from a given de-gree distribution P ( k ) (normalized such that (cid:80) k P ( k ) = N ), where k , the node degree, represents a two-vector ofthe in-degree and the out-degree, ( k in , k out ). Addition-ally, we consider an assortativity function a ( k (cid:48) → k ),which specifies the probability of a link from a node ofdegree k (cid:48) to one of degree k . In this N → ∞ limit, we as-sume that the state of the neurons can be represented bya continuous probability distribution, f ( θ, η | k , t ), suchthat f ( θ, η | k , t ) dθdη is the probability that a node ofdegree k has an excitability parameter in the range[ η, η + dη ] and a phase angle in the range [ θ, θ + dθ ] at time t . Since we are assuming that the excitability parametersdo not vary with time, we define g ( η | k ) = (cid:82) f dθ , which isthe time independent distribution of the excitability pa-rameters η i in the network for a randomly chosen nodeof degree k .In order to describe the synchronization behavior ofthis system, we define the order parameter to be , R ( t ) = 1 N N (cid:88) j =1 e iθ j . (4)As in previous work by Restrepo and Ott , we hypothe-size that in networks with large nodal degrees, the orderparameter can be well approximated via a mean field or-der parameter, defined by a continuum version of Eq.(4),¯ R ( t ) = 1 N (cid:88) k (cid:48) P ( k (cid:48) ) (cid:90) (cid:90) f ( θ (cid:48) , η (cid:48) | k (cid:48) , t ) e iθ (cid:48) dθ (cid:48) dη (cid:48) . (5)Additionally, the distribution f is constrained by the con-tinuity equation, ∂f∂t + ∂∂θ ( v θ f ) = 0 , (6)where v θ is the continuous version of the right hand sideof Eq. (2), v θ = (1 − cos θ ) + (1 + cos θ ) (cid:34) η + d n K (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) ) a ( k (cid:48) → k ) × (cid:90) (cid:90) f ( θ (cid:48) , η (cid:48) | k (cid:48) , t )(1 − cos θ (cid:48) ) n dθ (cid:48) dη (cid:48) (cid:35) . (7) IV. DIMENSION REDUCTION
Employing the dimensional reduction method of Ottand Antonsen , and following its previous applicationto the theta neuron , we assume that f is given by theFourier expansion, (8) f ( θ, η | k , t ) = g ( η | k )2 π (cid:40) ∞ (cid:88) p =1 (cid:2) b ( η, k , t ) p e − ipθ + b ∗ ( η, k , t ) p e ipθ (cid:3)(cid:41) . We then use the binomial theorem to expand the pulsefunction P n ( θ ) using(1 − cos θ ) n = A + n (cid:88) p =1 A p [ e ipθ + e − ipθ ] , (9)where A p = n (cid:88) j,m =0 δ j − m,p Q jm , (10) and Q jm = ( − j n !2 j m ! ( n − j )! ( j − m )! . (11)If we now assume a Lorentz distribution of the excitabil-ity parameters, g ( η | k ) = 1 π ∆( k )[ η − η ( k )] + ∆ ( k ) , (12)we obtain (cid:90) (cid:90) f ( θ (cid:48) , η (cid:48) | k , t ) e ipθ dθ (cid:48) dη (cid:48) = ˆ b ( k , t ) p , p > , p = 0ˆ b ∗ ( k , t ) | p | , p < , (13)with ˆ b ( k , t ) ≡ b ( η ( k ) + i ∆( k ) , k , t ). This now allows usto rewrite v θ in terms of ˆ b ( k , t ) as v θ = ge iθ + h + g ∗ e − iθ , (14)where g = −
12 (1 − η − K (cid:104) k (cid:105) H n ( k , t )) , h = 1 + η + K (cid:104) k (cid:105) H n ( k , t ) , (15)and H n ( k , t ) = d n (cid:88) k (cid:48) (cid:40) P ( k (cid:48) ) a ( k (cid:48) → k ) × (cid:34) A + n (cid:88) p =1 A p (ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p ) (cid:35)(cid:41) . (16)Substituting the phase velocity Eq. (14) and the Ott-Antonsen ansatz Eq. (8) into the continuity equation(6), we find that b ( η, k , t ) satisfies: ∂b∂t = i ( gb + hb + g ∗ ) . (17)Inserting the forms for g and h from Eq. (15) and (16)into this expression, and evaluating each quantity at thepole, η = η ( k ) + i ∆( k ), we obtain a reduced system ofequations for ˆ b ( k , t ) describing the mean field dynamicsof the neuronal network, ∂ ˆ b ( k , t ) ∂t = − i (ˆ b ( k , t ) − b ( k , t ) + 1) (cid:40) − ∆( k )+ iη ( k ) + id n K (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) ) a ( k (cid:48) → k ) × (cid:34) A + n (cid:88) p =1 A p (ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p ) (cid:35)(cid:41) . (18)The mean field order parameter, ¯ R ( t ) can now be writ-ten in terms of ˆ b ( k , t ). Using the assumed form for f ( θ, η | k , t ), we can evaluate the integrals in Eq. (5) usingCauchy’s residue theorem to obtain¯ R ( t ) = 1 N (cid:88) k P ( k )ˆ b ( k , t ) . (19)For the discussion in this paper, we will restrict theassortativity function to be of the form used previouslyby Restrepo and Ott a ( k (cid:48) → k ) = h ( a k (cid:48) → k ) , (20)where h ( x ) = min(max( x, ,
1) is defined to ensure that a ( k (cid:48) → k ) is a valid probability (i.e. 0 ≤ a ( k (cid:48) → k ) ≤ a k (cid:48) → k = k (cid:48) out k in N (cid:104) k (cid:105) (cid:20) c (cid:18) k (cid:48) in − (cid:104) k (cid:105) k (cid:48) out (cid:19) (cid:18) k out − (cid:104) k (cid:105) k in (cid:19)(cid:21) , (21)where c is a parameter used to vary the network assorta-tivity (with c > c < c = 0), the probability offorming a link between two nodes is simply proportionalto the out-degree of the source node and the in-degree ofthe target node.The in-out Pearson assortativity coefficient, r , is astatistic used to characterize the overall assortativity ofa network, and is defined as r = (cid:80) e [( k (cid:48) in − (cid:104) k (cid:105) )( k out − (cid:104) k (cid:105) )] (cid:112)(cid:80) e ( k (cid:48) in − (cid:104) k (cid:105) ) (cid:112)(cid:80) e ( k out − (cid:104) k (cid:105) ) , (22)where (cid:80) e is the sum over all edges connecting a nodeof degree k (cid:48) to a node of degree k . Assuming that a ( k (cid:48) → k ) = a k (cid:48) → k , and that the in and out degreedistributions are independent, we can relate the assorta-tivity coefficient to the parameter c as r = c (cid:104) k (cid:105) (cid:113) ( (cid:104) k in (cid:105) − (cid:104) k (cid:105) )( (cid:104) k out (cid:105) − (cid:104) k (cid:105) ) , (23)which can be seen by noting that the sum of a quan-tity Q ( k , k (cid:48) ), defined on each edge connecting a nodeof degree k (cid:48) to a node of degree k , over edges in ourmean field formulation would be given by (cid:80) e 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) ).The expression for the assortativity coefficient as afunction of c , Eq. (23), is unbounded, while the Pearsonassortativity is by definition bounded between − c , theassortativity function given in Eq. (21) is not a proba-bility. However, for the network parameters used in ournumerical example below, we find that Eq. (23) is veryaccurate for | c |≤ .
5, corresponding to an assortativityrange, | r | (cid:46) . g ( η | k ) ≡ g ( η )) and the ˆ b ’s are given k independent identical initial con-ditions, ˆ b ( k , ≡ ˆ b (0), then there are a few notable casesin which particular degree distributions and our chosenassortativity function Eq. (21) allow for further dimen-sional reduction. For networks with a delta-function de-gree distribution, P ( k ) = δ k in ,k δ k out ,k , the Eq. (18) re-duces to a single equation describing the mean field dy-namics, (24) ∂ ˆ b ( t ) ∂t = − i (ˆ b ( t ) − b ( t ) + 1) (cid:40) − ∆ + iη + id n K (cid:34) A + n (cid:88) p =1 A p (ˆ b ( t ) p + ˆ b ∗ ( t ) p ) (cid:35)(cid:41) . We note that this equation is identical to earlier resultsfor a fully connected network . Thus, networks withonly a single allowed degree have identical asymptoticdynamics to a fully connected network. This result isconsistent with analogous results by Barlev et al for anetwork of Kuramoto oscillators. More generally, if thenetwork has fixed in-degree, P ( k ) = P ( k out ) δ k in ,k , thesystem is similarly reduced to the single dynamical equa-tion, Eq. (24). On the other hand, if the out-degree isfixed, P ( k ) = P ( k in ) δ k out ,k , then dynamics of ˆ b ( k , t ) isindependent of k out , further reducing the dimensionalityof the problem. Reduction efficiency
Equation (18) represents a reduction of the originalsystem of N theta neurons to a system with as manyequations as there are k values in the support of the de-gree distribution P ( k ). We denote this quantity by M k ,which, in the case of independent in and out-degree dis-tributions, is equal to M in × M out , where M in and M out are the number of possible in-degrees and out-degrees re-spectively. In general, simulating the full network, Eq.(2), requires O ( N ) floating point operations per timestep. Using the form of the assortativity function givenin Eq. (21) the sum over k (cid:48) in the reduced system ofequations can be split into two sums, each independentof k , k in N (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) ) k (cid:48) out A + c k out − (cid:104) k (cid:105) N (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) )( k (cid:48) in −(cid:104) k (cid:105) ) A . (25)where A = A + (cid:80) np =1 A p (cid:16) ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p (cid:17) . Sincethe two sums in Eq. (25) are independent of k , eachmust be calculated only once per simulation iteration.Thus, simulating the reduced system Eq. (18) only re-quires O ( M k ) floating point operations per time step — M k operations performed once for each of these two sumsand M k operations for each of the ˆ b ( k , t ) equations. Inmany cases, M k (cid:28) N , so that simulating Eq. (18) is FIG. 2. The effect of varying levels of interpolation on thecalculated results for the trajectories of ¯ R ( t ) in the complexplane starting from an initial condition of ¯ R ( t ) = 0 and endingat a fixed point attractor for K = 3 in a network with neutralassortativity, with η = − .
1. Calculation ofthe order parameter dynamics is robust to a large range inthe level of interpolation. Using as few as 10% of the totalavailable degrees and interpolating the remaining 90% giveresults close to the calculation without interpolation. In therest of this paper we employ a 10% interpolation level in allour mean field calculations. The black arc is a segment of theunit circle | ¯ R ( t ) | = 1. significantly more efficient that simulating the full net-work. Furthermore, if c is set to 0, which is the case ofnetworks with neutral assortativity, then ˆ b ( k , t ) will haveno dependence on k out , and hence the overall problemis reduced to M in independent equations, allowing evengreater computational efficiency.Since ˆ b ( k (cid:48) , t ), P ( k (cid:48) ), and a ( k (cid:48) → k ) are each smoothlyvarying functions, we can achieve further dimensional re-duction by interpolating the summand in Eq. (18) us-ing a coarse-grained grid of k values. In particular, Eq.(18) is not solved for ˆ b ( k , t ) for all of the M k values of k , but only for the small subset of k values that lie onthe coarse-grained grid in k -space The summands on theright hand side of Eq. (18) at k values not on the gridare approximated by a bilinear interpolation of the valuesat the surrounding chosen k values. To perform the bi-linear interpolation, we first interpolate linearly betweenneighboring grid values in one direction. The value ofthe summand at a given k value is then approximatedby linearly interpolating in the other direction betweenvalues estimated with the previous linear interpolation.We find that using as few as 10% of the network degreesyield very accurate results, while an even coarser inter-polation still produces the same qualitative behavior ascan be seen in Fig. 2. 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 P ( k ) = k < k min Ak − γ if k min ≤ k < k max k max ≤ k. (26)The exponent of the power law distribution, γ , was set to3, and k min and k max were set to 750 and 2000, respec-tively. As mentioned earlier, the normalization constant A is chosen to make (cid:80) k P ( k ) = N . We will also setthe parameter n controlling the sharpness of the synap-tic pulse to 2 for all examples considered, and will use aninterpolation level of 10% for all calculations using thereduced system of equations for the mean field theory(cf. Fig. 2).From numerical simulations of the reduced equations,(18), we find that the long term dynamics of the order pa-rameter can be broadly classified into one of three phases– (1) the partially resting (PR) phase; (2) the asyn-chronously firing (AF) phase; and (3) the synchronouslyfiring (SF) phase. The PR phase and the AF phase ap-pear as fixed points in the dynamics of the order param-eter, whereas the SF phase appears as a limit cycle of theorder parameter. A. Fixed points
As a particular example to illustrate the differenttypes of fixed points, we look at a network with neu-tral assortativity ( c = 0) having excitability parametersdistributed according to a Lorentzian distribution withmean η = − . | ¯ R | = 1. In this phase, most of the indi-vidual neurons in the network are independently in theirresting states, in a fashion similar to Fig. 1(a). Thiscorresponds to the case of K = 1 in Fig. 3(a), in whichthe fixed point is located near the edge of the unit cir-cle marked in black. Further, the time series of a fewrandomly chosen neurons (Fig. 3(b)) demonstrates thatalmost all of the neurons are in a resting state. Whilethere may be a small number of neurons that are in thespiking phase due to the spread in the distribution ofvalues of excitability parameters, η , these do not haveany significant effect on the full order parameter of thesystem.As we increase the coupling constant K , the systemtransitions to the asynchronously firing (AF) phase, inwhich the order parameter goes to a fixed point locatednear the center of the unit circle. In this phase, most ofthe individual neurons in the network are asynchronously FIG. 3. (a): Fixed points of R ( t ) observed in networks with neutral assortativity, η = − .
1, for three values ofthe 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 randomlychosen neurons demonstrates that in the PR phase almost all neurons are in a resting state, and as the system approaches theAF state, more nodes transition to an oscillating, excited state.The thick dashed line corresponds to the position of the fixedpoint of the order parameter for the corresponding value of K . firing, in a fashion similar to Fig. 1(c). This can beseen in the case of K = 6 in Fig. 3(c) which showsthat almost all of the neurons are in a recurrent spikingstate. Note that even though the neurons are spiking asynchronously , i.e., their firing times are independentof one another , the fixed point of the order parameteris not at ¯ R = 0. This is because the angular velocityof an individual neuron is not constant along the circle,thus in the average over the ensemble of neurons a biasis present towards the direction for which the angularvelocity of neurons is minimized. As discussed in Sec. II,this may occur at either θ = 0 or at θ = π , dependent onhow large the excitability parameter is for the neuron.We now examine the transition from the PR phase tothe AF phase. Microscopically, in the PR phase, almostall of the neurons are individually in a resting phase,whereas in the AF phase almost all neurons are in thespiking state. To examine the behavior at an interme-diate point, we look at the fixed point for the case of K = 3, as shown in Fig. 3(c). At this intermediate valueof the coupling constant, a fraction of the neurons are inthe spiking state. In particular, the nodes that begin tospike first are those which have larger in-degrees. Thisis demonstrated in Fig. 4, in which we examine ˆ b ( k ) atthe fixed point for K = 3. Since we are looking at a net-work with neutral assortativity ( c = 0), Eq. (25) impliesthat the sum only depends on the out-degree through acommon multiplicative factor. Thus ˆ b is only plotted asa function of k in . Analogously, for the fixed point of thedynamics on the full network, the range of degrees from k min to k max is divided uniformly into several intervals,and for each interval we find a partial order parameter,calculated such that the average in Eq. (4) is only per-formed over those nodes whose in-degree lie within thatinterval, i.e., R ( k in , t ) = 1 ||N || (cid:88) j ∈N e iθ j , (27)where N is the set of nodes having an in-degree within FIG. 4. Comparison of | ˆ b ( k in ) | from the reduced systemof equations and the time average of | R ( k in ) | from the fullsystem, Eq. (27), for a network with neutral assortativity( c = 0), η = −
2, and ∆ = 0 . K = 3. The dynamics un-der these parameters were simulated in a network with 5000nodes, and the network was allowed to relax to a fixed point.Nodes were divided into classes according to their in-degree tocalculate the time averaged effective order parameter for eachclass, which is shown in blue, with the error bars denoting theroot mean squared time fluctuation of the order parameter forthat class. The time fluctuations are due to the finite numberof nodes in each class. (See text for details.) one of the intervals of the range of degrees, ||N || is thenumber of nodes in the set, and k in is the average in-degree of nodes within that set.In addition, we find that the transition from the PRphase to the AF phase occurs via a hysteretic processmediated by saddle node bifurcations. To illustrate this,we evolved the dynamics of the full network in a stepwise fashion by increasing the coupling constant K insmall increments of 0.2, and allowing the system to relaxto an equilibrium before the next increment (Fig. 5(a)).We also compare this with the analogous hysteresis curveobserved for the evolution of the system dynamics on anErd˝os-R´enyi network having the same size and averagedegree as the scale free network being considered (Fig.5(b)). While the hysteretic region begins at around thesame value of the coupling constant, K , for both networktopologies, we find that for the case of the Erd˝os-R´enyinetwork, which has a sharply peaked degree distribution,the range in K that allows hysteresis (3 (cid:46) K (cid:46) .
25) issignificantly larger than the corresponding range for thenetwork with the scale free degree distribution (3 . (cid:46) K (cid:46) ∂ ˆ b ( k , t ) /∂t = 0 for the fixed points, we find that theequilibrium ˆ b ( k ) satisfy,ˆ b ± ( k ) = 1 ± z ( k )1 ∓ z ( k ) , (28)where (29) iz ( k ) = − ∆ + iη + id n K (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) ) a ( k (cid:48) → k ) × (cid:34) A + n (cid:88) p =1 A p (ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p ) (cid:35) , and the sign is chosen to ensure | ˆ b ( k ) |≤
1. Using our formof the assortativity function Eq. (21), we may again splitthe above sum into two parts as in Eq. (25). Thus wemay rewrite Eq. (29) as iz ( k ) = − ∆ + iη + ik in X + i ( k out − (cid:104) k (cid:105) ) Y, (30)where X and Y are given by, (31a) X = d n KN (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) ) k (cid:48) out (cid:34) A + n (cid:88) p =1 A p (ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p ) (cid:35) (31b) Y = d n KN (cid:104) k (cid:105) (cid:88) k (cid:48) P ( k (cid:48) )( k (cid:48) in − (cid:104) k (cid:105) ) (cid:34) A + n (cid:88) p =1 A p (ˆ b ( k (cid:48) , t ) p + ˆ b ∗ ( k (cid:48) , t ) p ) (cid:35) . These simplifications allow for efficient calculation of thesystem fixed points. Choosing initial values, X and Y , we calculate the associated z ( k ) and ˆ b ( k ) using Eq. (30)and Eq. (28), and then recalculate new values, X and Y using Eq. (31). For fixed points of the reduced equa-tions δX = X − X and δY = Y − Y are both zero. Wecalculate δX and δY for several different initial valuesat regularly spaced intervals for X and Y , and iden-tify the fixed points as the points where δX = δY = 0.The interpolation procedure described earlier can also beapplied to this calculation to further increase efficiency.For the nonassortative case ( c = 0), Y = 0 always, soidentifying the fixed points in this case only requires cal-culating the variation in the single parameter X . We usethis method to evaluate the fixed points of the reducedequations for the range of K over which hysteresis wasobserved, and find close agreement between the results ofthis fixed point analysis and the direct evolution of thefull network (Fig. 5). B. Limit Cycles
As a representative example of limit cycles of ¯ R ( t ), weconsider a network with neutral assortativity with ex-citability parameters η distributed as a Lorentzian withmean η = 10 .
75 and width ∆ = 0 .
5, and with a couplingconstant K = −
9. In the SF phase, the order parametergoes to a limit cycle in the complex plane. In this phase,a majority of the neurons are synchronously in a spikingstate. Plots for such limit cycles are shown in Fig. 6,in which we plot the trajectory of the order parameterin the complex plane (after removing transients) for anetwork with the scale free degree distribution given inEq. (26) (blue solid curve), a corresponding Erd˝os-R´enyinetwork having a Poissonian degree distribution (greendashed curve), and a regular network having a delta func-tion degree distribution (i.e. P ( k ) = δ k in ,k δ k out ,k ) (reddotted curve), each having the same average degree. Asseen earlier in Eq. (24), a network with a delta functiondegree distribution has mean field dynamics identical tothose of a fully connected network, and the correspond-ing limit cycle in Fig. 6 is identical to the limit cycleobtained at these parameters for the fully connected net-work by Luke et al. In comparison with the limit cy-cles that are observed for the case of the regular networkor the Erd˝os-R´enyi network, the limit cycles in networkswith scale free degree distributions are diminished in size,due to the large variation in nodal behavior as a functionof degree. Nodes with smaller in-degrees were observedto predominantly be in the spiking phase, with high syn-chronization and a larger limit cycle for the partial orderparameter, whereas nodes with larger in-degrees were inthe resting phase. Due to this differentiation of behaviorwith degree, the averaged full order parameter exhibits alimit cycle that is somewhat reduced in size when com-pared with the results for a fully connected network byLuke et al . However, we see that the limit cycles forthe Erd˝os-R´enyi network are similar in shape and struc-ture to the limit cycles obtained for the regular network, FIG. 5. A sweeping value of K was used to observe the change in phase from the PR state to the AF state. Hysteresis wasobserved on the network with a scale free degree distribution (a) as well as a corresponding Erd˝os-R´enyi network having thesame size and the same average degree (b). For the full network, at each value of K the mean of the order parameter afterignoring the transients have been marked as triangles. A close match is observed with the fixed points as computed from meanfield equations directly (see text for details). Hysteresis is observed for 3 . (cid:46) K (cid:46) (cid:46) K (cid:46) .
25 in the corresponding Erd˝os-R´enyi network (Note the difference in scales for the x-axis in both plots).An apparent crossing of the fixed point curve is seen in (b), which is an artifact of the non-self-intersecting ¯ R curve lying inthe two dimensional complex space, which has been projected onto the real axis in this plot. as would be expected in accordance with the discussionin Sec. IV, since the Poissonian degree distribution forthe Erd˝os-R´enyi network is sharply peaked about the av-erage degree and hence cannot admit a large variationof behavior with nodal degree. As the average degree, (cid:104) k (cid:105) increases, the red and green curves converge becausethe Poisson degree distribution appropriate for an Erd˝os-R´enyi network approaches a delta function. C. Effect of Assortativity
We now consider the effect of assortativity on the limitcycle dynamics of the order parameter in the network.While limit cycle behavior exists in networks with neu-tral assortativity ( c = 0), introduction of assortativity ordissasortativity in the network can cause the limit cycleattractor to transform to a fixed point attractor (AF likestate) via a Hopf bifurcation. This is demonstrated inFig. 7, in which we show that varying c away from zero to ± . r ≈ ± . N induced noise as seenfrom the size of the clouds surrounding the fixed pointposition calculated from the reduced system. VI. CONCLUSION
Using a mean field approximation, in conjunction withthe Ott-Antonsen ansatz, we obtained a reduced systemof equations that successfully model the macroscopic or-der parameter dynamics of a large network of theta neu-rons. This reduced system of equations allows us to ex-amine the effects of varying the network parameters andthe network topology (in terms of degree distributions,as well as degree correlations) in a computationally effi-cient fashion. The order parameter of the network is usedfor describing the macroscopic behavior of the network oftheta neurons, whose attractors can be of various types.In particular, we find resting states, asynchronously fir-ing states and synchronously firing states, the first twoof which appear as a fixed point for the order param-eter (Fig. 3), while the third appears as a limit cyclefor the order parameter (Fig. 6). We also used the re-duced system of equations to observe the effect of varyingthe assortativity in the system and demonstrated that asynchronously firing phase was only present for networkswith neutral or small assortativity, and the addition ofmoderate amounts of assortativity or disassortativity tothe network causes the system to go to an asynchronouslyfiring state instead (Fig. 7). Further, for networks withscale free degree distributions, we find that nodes withdifferent values of their degrees admit a large variation ofbehavior (Fig. 4), a phenomenon not possible in networkswith all-to-all connectivity. In all cases close agreementwas observed between the order parameter dynamics aspredicted by the reduced system of equations (Eq. 18),
FIG. 6. Comparison of the limit cycle attractor for ¯ R ( t ) inthe complex plane across varying degree distributions in anetwork with neutral assortativity ( c = 0) with η = 10 . . K = −
9. The scale free network (blue solidcurve) has a degree distribution according to Eq. (26),the Erd˝os-R´enyi network (green dashed curve) has a Poisso-nian degree distribution, and the regular network (red dottedcurve) has a delta function degree distribution. The blackcircle is the unit circle | ¯ R | =1 and as calculated by evolution of the full system of equa-tions Eq. (2). ACKNOWLEDGEMENTS
This work was supported by the Army Research Of-fice under Grant No. W911NF-12-1-0101, and by theNational Science Foundation under Grant No. PHY-1461089. D. C. Michaels, E. P. Matyas, and J. Jalife, “Mechanisms ofsinoatrial pacemaker synchronization: a new hypothesis.” Circu-lation Research , 704–714 (1987). K. Wiesenfeld, P. Colet, and S. H. Strogatz, “Frequency lock-ing in josephson arrays: connection with the kuramoto model,”Physical Review E , 1563 (1998). I. Z. Kiss, Y. Zhai, and J. L. Hudson, “Emerging coherencein a population of chemical oscillators,” Science , 1676–1678(2002). A. E. Motter, S. A. Myers, M. Anghel, and T. Nishikawa, “Spon-taneous synchrony in power-grid networks,” Nature Physics ,191–197 (2013). B. A. Carreras, V. E. Lynch, I. Dobson, and D. E. Newman,“Complex dynamics of blackouts in power transmission systems,”Chaos , 643–652 (2004). L. Glass and S. A. Kauffman, “The logical analysis of continuous,non-linear biochemical control networks,” Journal of TheoreticalBiology , 103–129 (1973). M. Aldana and P. Cluzel, “A natural class of robust networks,”Proceedings of the National Academy of Sciences , 8710–8714(2003). T. B. Luke, E. Barreto, and P. So, “Complete classification ofthe macroscopic behavior of a heterogeneous network of thetaneurons,” Neural Computation , 3207–3234 (2013). M. M. Abdulrehem and E. Ott, “Low dimensional description ofpedestrian-induced oscillation of the millennium bridge,” Chaos , 013129 (2009). E. Montbri´o, D. Paz´o, and A. Roxin, “Macroscopic descriptionfor networks of spiking neurons,” Physical Review X , 021028(2015). D. Paz´o and E. Montbri´o, “Low-dimensional dynamics of popula-tions of pulse-coupled oscillators,” Physical Review X , 011009(2014). C. R. Laing, “Derivation of a neural field model from a networkof theta neurons,” Physical Review E , 010901 (2014). Z. Lu, K. Klein-Carde˜na, S. Lee, T. M. Antonsen, M. Girvan,and E. Ott, “Resynchronization of circadian oscillators and theeast-west asymmetry of jet-lag,” Chaos , 094811 (2016). E. Ott and T. M. Antonsen, “Low dimensional behavior of largesystems of globally coupled oscillators,” Chaos , 037113 (2008). E. Ott and T. M. Antonsen, “Long time evolution of phase os-cillator systems,” Chaos , 023117 (2009). E. Ott, B. R. Hunt, and T. M. Antonsen Jr, “Comment onLong time evolution of phase oscillator systems [Chaos 19, 023117(2009)],” Chaos , 025112 (2011). J. G. Restrepo and E. Ott, “Mean-field theory of assortativenetworks of phase oscillators,” Europhysics Letters , 60006(2014). E. A. Martens, E. Barreto, S. Strogatz, E. Ott, P. So, and T. An-tonsen, “Exact results for the kuramoto model with a bimodalfrequency distribution,” Physical Review E , 026204 (2009). G. Barlev, T. M. Antonsen, and E. Ott, “The dynamics of net-work coupled phase oscillators: An ensemble approach,” Chaos , 025103 (2011). P. S. Skardal, J. G. Restrepo, and E. Ott, “Frequency assorta-tivity can induce chaos in oscillator networks,” Physical ReviewE , 060902 (2015). D. Paz´o and E. Montbri´o, “From quasiperiodic partial synchro-nization to collective chaos in populations of inhibitory neuronswith delay,” Physical Review Letters , 238101 (2016). J. Roulet and G. B. Mindlin, “Average activity of excitatory andinhibitory neural populations,” Chaos , 093104 (2016). G. B. Ermentrout and N. Kopell, “Parabolic bursting in an ex-citable system coupled with a slow oscillation,” SIAM Journal onApplied Mathematics , 233–253 (1986). B. Ermentrout, “Type i membranes, phase resetting curves, andsynchrony,” Neural Computation , 979–1001 (1996). E. M. Izhikevich, “Class 1 neural excitability, conventionalsynapses, weakly connected networks, and mathematical foun-dations of pulse-coupled models,” IEEE Transactions on NeuralNetworks , 499–507 (1999). A. L. Hodgkin, “The local electric changes associated with repeti-tive action in a non-medullated axon,” The Journal of Physiology , 165 (1948). C. B¨orgers and N. Kopell, “Synchronization in networks of exci-tatory and inhibitory neurons with sparse, random connectivity,”Neural computation , 509–538 (2003). M. E. Newman, “Assortative mixing in networks,” Physical re-view letters , 208701 (2002). P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey,V. J. Wedeen, and O. Sporns, “Mapping the structural core ofhuman cerebral cortex,” PLoS Biology , e159 (2008). S. Bialonski and K. Lehnertz, “Assortative mixing in functionalbrain networks during epileptic seizures,” Chaos , 033139(2013). E. Barzegaran, A. Joudaki, M. Jalili, A. O. Rossetti, R. S. Frack-owiak, and M. G. Knyazeva, “Properties of functional brainnetworks correlate with frequency of psychogenic non-epilepticseizures,” Frontiers in Human Neuroscience (2012). W. de Haan, Y. A. Pijnenburg, R. L. Strijers, Y. van der Made,W. M. van der Flier, P. Scheltens, and C. J. Stam, “Func- FIG. 7. For the parameters η = 4, ∆ = 0 . K = − .
8, in a network with neutral assortativity ( c = 0), the system lies inan SF state (as in (b)). Varying the assortativity in either direction ( c = ± .
5, corresponding to r ≈ ± . tional neural network analysis in frontotemporal dementia andalzheimer’s disease using eeg and graph theory,” BMC Neuro-science , 1 (2009). S. Teller, C. Granell, M. De Domenico, J. Soriano, S. G´omez, andA. Arenas, “Emergence of assortative mixing between clusters ofcultured neurons,” PLoS Computational Biology , e1003796(2014). Some authors, such as Restrepo and Ott define the order pa-rameter differently so as to be weighted with the out-degree ateach node, i.e., R ( t ) = (cid:80) Ni =1 (cid:80) Nj =1 A ij e iθ j / (cid:16)(cid:80) Ni =1 (cid:80) Nj =1 A ij (cid:17) . J. G. Foster, D. V. Foster, P. Grassberger, and M. Paczuski,“Edge direction and the structure of networks,” Proceedings ofthe National Academy of Sciences , 10815–10820 (2010). For another, often useful, definition of a coefficient quantitativelycharacterizing the assortativity or disassortativity of a networksee Ref. . This definition of asynchronous spiking is consistent with remarksby other authors , wherein asynchronous states have beendefined as states in which at each neuron the term coupling itto the other neurons in the network is independent of time, as isobserved in the cases of fixed points. J. G. Restrepo, E. Ott, and B. R. Hunt, “Approximating thelargest eigenvalue of network adjacency matrices,” Physical Re-view E , 056119 (2007). L. Abbott and C. van Vreeswijk, “Asynchronous states in net-works of pulse-coupled oscillators,” Physical Review E , 1483(1993). D. Hansel and G. Mato, “Existence and stability of persistentstates in large neuronal networks,” Physical Review Letters86