## Abstract

We design and analyze the dynamics of a large network of theta neurons, which are idealized type I neurons. The network is heterogeneous in that it includes both inherently spiking and excitable neurons. The coupling is global, via pulselike synapses of adjustable sharpness. Using recently developed analytical methods, we identify all possible asymptotic states that can be exhibited by a mean field variable that captures the network's macroscopic state. These consist of two equilibrium states that reflect partial synchronization in the network and a limit cycle state in which the degree of network synchronization oscillates in time. Our approach also permits a complete bifurcation analysis, which we carry out with respect to parameters that capture the degree of excitability of the neurons, the heterogeneity in the population, and the coupling strength (which can be excitatory or inhibitory). We find that the network typically tends toward the two macroscopic equilibrium states when the neuron's intrinsic dynamics and the network interactions reinforce one another. In contrast, the limit cycle state, bifurcations, and multistability tend to occur when there is competition among these network features. Finally, we show that our results are exhibited by finite network realizations of reasonable size.

## 1.  Introduction

The cortex of the brain is structured into a complex network of many functional neural assemblies (Sherrington, 1906; Hebb, 1949; Harris, 2005). Each assembly typically encompasses a large number of interacting neurons with various dynamical characteristics. Within this network, communication among the neural assemblies generally involves macroscopic signals that arise from the collective behavior of the constituent neurons in each assembly. It has been proposed that functional behavior arises from these interactions, and perceptual representations in the brain are believed to be encoded by the macroscopic spatiotemporal patterns that emerge in these networks of assemblies.

In modeling the brain, a microscopic description of individual neurons and their interactions is important. However, a model for the macroscopic dynamical behavior of large assemblies of neurons is essential for understanding the brain's emergent, collective behavior (Peretto, 1984; Sompolinsky, 1988; Kanamaru & Masatoshi, 2003). In this letter, we construct such a model using the canonical theta neuron of Ermentrout and Kopell (Ermentrout & Kopell, 1986; Ermentrout, 1996). We construct a large heterogeneous network, containing both excitable and spiking neurons, that is globally coupled via smooth pulselike synapses. Then, using the methods of Ott and Antonsen (2008, 2009) and Marvel, Mirollo, and Strogatz (2009) (see also Pikovsky & Rosenblum, 2008, 2011), we derive a low-dimensional “reduced” dynamical system that exhibits asymptotic behavior that is coinciding with that of the macroscopic mean field of the network. We use this reduced system to classify all the asymptotic macroscopic configurations that the network can exhibit.

We show that our network exhibits three fundamental collective states: a partially synchronous rest state (PSR), a partially synchronous spiking state (PSS), and a collective periodic wave state (CPW). In the PSR state, most neurons remain at rest but are excitable, and the macroscopic mean field sits on a stable equilibrium. In the PSS state, the mean field is also on a stable equilibrium, but typically, most individual neurons spike regularly. These states are similar to states that have been called asynchronous (Abbott & van Vreeswijk, 1993; Hansel & Mato, 2001, 2003). We find that they are typically encountered in cooperative networks in which the internal dynamics of the neurons and the interneuronal network interactions reinforce each other. In other parameter regions where the internal dynamics and network interaction are in competition, a collective periodic state, the CPW state, can occur. In the CPW state, the phases of the neurons transiently cluster, and the degree of network coherence waxes and wanes periodically in time in such a way that the macroscopic mean field exhibits a stable limit cycle.

These three macroscopic states can coexist and transition into each other, and here we clarify precisely how this happens using our reduced mean field equation. We obtain a complete bifurcation diagram with respect to parameters that represent the degree of excitability, heterogeneity, and the strength of coupling (both excitatory and inhibitory) within the network. Our model provides a comprehensive description for the asymptotic macroscopic behavior of a large network of heterogeneous theta neurons in the thermodynamic limit.

The remainder of this letter is organized as follows. In section 2, we describe the basic features of our theta neuron network, and in section 3, we derive the mean field reduction using the Ott-Antonsen method (Ott & Antonsen, 2008, 2009; Marvel et al., 2009). Then, in section 4, we use the reduced mean field equation to identify and describe the three possible macroscopic states: PSR, PSS, and CPW. In section 5, we provide a comprehensive bifurcation analysis for the macroscopic dynamics of the network. We summarize and discuss our results in section 6.

## 2.  Microscopic Formulation

The brain contains a huge number of different types of neurons that feature various dynamical characteristics. Here, we construct a mathematically tractable network that encompasses three fundamental characteristics: neuronal excitability, pulselike coupling, and heterogeneity.

### 2.1.  Neuronal Excitability: The Theta Neuron.

A typical neuron at rest begins to spike as a constant injected current exceeds a threshold. Neurons are usually classified into two types based on this behavior (Hodgkin, 1948; Ermentrout, 1996; Izhikevich, 2007). Type I neurons begin to spike at an arbitrarily slow rate, whereas type II neurons spike at a nonzero rate as soon as the threshold is exceeded. Neurophysiologically, excitatory pyramidal neurons are often of type I, and fast-spiking inhibitory interneurons are often of type II (Nowak, Azouz, & Sanchez-Vives, Gray, & McCormick, 2003; Tateno, Harsch, & Robinson, 2004).1

Ermentrout and Kopell (1986; Ermentrout, 1996) showed that near the firing threshold, type I neurons can be represented by a canonical phase model that features a saddle-node bifurcation on an invariant cycle, or a SNIC bifurcation. This model has come to be known as the theta neuron and is given by
2.1
where is a phase variable on the unit circle and is a bifurcation parameter related to an injected current. The SNIC bifurcation occurs for . For , the neuron is attracted to a stable equilibrium which represents the resting state. An unstable equilibrium is also present, representing the threshold. If an external stimulus pushes the neuron's phase across the unstable equilibrium, will move around the circle and approach the resting equilibrium from the other side. When crosses , the neuron is said to have spiked. Thus, for , the neuron is excitable. As the parameter increases, these equilibria merge in a saddle-node bifurcation and disappear, leaving a limit cycle. Consequently, the neuron spikes regularly. This transition is depicted schematically in Figure 1.
Figure 1:

SNIC bifurcation of the theta neuron. For , the neuron is at rest but excitable. For , the neuron spikes regularly. The SNIC bifurcation occurs at . A spike is said to occur when the phase variable crosses .

Figure 1:

SNIC bifurcation of the theta neuron. For , the neuron is at rest but excitable. For , the neuron spikes regularly. The SNIC bifurcation occurs at . A spike is said to occur when the phase variable crosses .

In this letter, we construct our network using only theta neurons. Future work will consider networks containing mixtures of type I and type II neurons.

### 2.2.  Coupling via a Pulselike Synapse.

We consider a network of N theta neurons coupled together via a synaptic current Isyn,
2.2
with . Thus, the synaptic current changes the effective excitability parameter of the jth neuron.
We write Isyn as a collective signal in which each neuron contributes a pulselike synaptic current depending on its phase angle. Thus,
2.3
where , , and an is a normalization constant such that
The parameter n defines the sharpness of the pulselike synaptic current (Ariaratnam & Strogatz, 2001) such that it becomes more and more sharply peaked at as n increases, as shown in Figure 2. The sum in equation 2.3 is over the entire population, and k is the overall synaptic strength for the whole network.
Figure 2:

A plot of the synaptic function for different values of the sharpness parameter n. As n increases from 1 to 9, the shape of the synaptic function becomes more pulselike around the firing state at .

Figure 2:

A plot of the synaptic function for different values of the sharpness parameter n. As n increases from 1 to 9, the shape of the synaptic function becomes more pulselike around the firing state at .

### 2.3.  Heterogeneity.

Neurons in real biological networks exhibit a range of intrinsic excitabilities. To model this, we assume that the parameter for each neuron is different and is drawn at random from a distribution function . Here we assume a Lorentzian distribution,
2.4
where is the center of the distribution and is its half-width at half-maximum. Thus, describes the degree of neuronal heterogeneity in the network. Since always includes both positive and negative 's, our network therefore contains a mixture of both excitable and spontaneously spiking neurons, with the ratio being biased in favor of one or the other depending on the chosen value of .

## 3.  Mean Field Reduction of the Network

In the limit , the system affords a continuum description in which the network of theta neurons can be described by a probability density function (Kuramoto, 1975, 1984). Thus, gives the fraction of oscillators that have phases in and parameters in . The time evolution of F is governed by the continuity equation,
3.1
where , the velocity of a neuron, is the continuum version of equations 2.2 and 2.3:
3.2
In order to explore the collective behavior of this network, we introduce the macroscopic mean field (also known as the order parameter) z(t), which describes the macroscopic coherence of the network and is defined as
3.3

We now derive a reduced dynamical system that exhibits asymptotic behavior that coincides with the asymptotic behavior of the mean field z(t) of the network. Our procedure follows Ott and Antonsen (2008, 2009) and Marvel et al. (2009).

The velocity equation, equation 3.2, can be written in sinusoidally coupled form (Marvel et al., 2009) such that the explicit dependence on the individual oscillator's phase occurs only through the harmonic functions and :
3.4
with
3.5
and
3.6
Here, z is the mean field variable introduced in equation 3.3, n is the sharpness parameter of the synapse as previously described, and H(z, n)=Isyn/k is the rescaled synaptic current given by2
3.7
where
3.8
and
3.9
In these equations, denotes the complex conjugate of z and is the Kronecker delta function on the indices (i, j). Note that is a real-valued function.
Now we adopt the ansatz that the probability density function F can be written as a Fourier expansion in which the coefficients appear as powers of a single yet-to-be-determined complex function :
3.10
This ansatz defines a two-dimensional manifold (parameterized by the real and imaginary parts of ) in the space of all probability density functions. In Ott and Antonsen (2008), the authors showed that this manifold is invariant if and only if satisfies as well as the following differential equation (which is obtained by substituting equation 3.10 into equation 3.1):
3.11
The macroscopic mean field can be written as an integral transform of with as the kernel by substituting equation 3.10 into equation 3.3. This gives
3.12
The integro-differential equations defined by equations 3.11 and 3.12 give the general equation of motion for the asymptotic behavior of the macroscopic mean field z(t). Now, permitting to be complex and analytically continuing into the upper half of the complex plane, and further assuming that is given by the Lorentzian in equation 2.4, the integral in equation 3.12 can be evaluated in closed form using the residue theorem. The result is
where is the center and is the half-width-at-half-maximum of the natural frequency distribution given in equation 2.4.
Finally, by substituting f (see equation 3.5) and h (see equation 3.6) into equations 3.11 and 3.12 and evaluating at the residue, we arrive at the desired reduced dynamical system:
3.13
In accordance with Ott and Antonsen (2009), we find that the attractors of this two-dimensional ordinary differential equation are the attractors of the macroscopic dynamics for the infinite discrete network given by equations 2.2 to 2.4 with .3

## 4.  Macroscopic Dynamics of the Network

By analyzing the reduced mean field equation, equation 3.13, we find that there are three possible asymptotic macroscopic states for the network. Two of these, which we call the partially synchronous rest (PSR) state and the partially synchronous spiking (PSS) state, correspond to equilibria of the macroscopic mean field. The third, which we call the collective periodic wave (CPW) state, corresponds to a limit cycle of the macroscopic mean field.

### 4.1.  Macroscopic Equilibrium States.

In the PSR state, the macroscopic mean field z(t) of the network settles onto a stable node. This state is predominantly (but not exclusively) observed when the distribution of excitability parameters is such that most neurons are in the resting regime (e.g., ), and the neurons are coupled through inhibitory synapses (k<0). Thus, most neurons are inactive, with their phase angles residing near their resting states. Nevertheless, some spiking neurons are present. These come from the tail of the distribution and have a negligible effect on the collective behavior of the network. Figure 3a shows an example of the macroscopic PSR equilibrium, with its local invariant manifolds calculated using the reduced mean field equation, equation 3.13, using , , k=−2, and n=2. (A movie showing both the macroscopic and microscopic behavior of the PSR state in Figure 3 is available in the online supplement.) Figures 3b and c show a trajectory, after the initial transient behavior has been discarded, of the macroscopic mean field z(t) calculated from a network of N=10, 000 theta neurons with the same system parameters. As expected, the slightly noisy trajectory hovers about the predicted equilibrium with fluctuations roughly on the order of .

Figure 3:

Phase portraits of the PSR macroscopic state with , , k=−2, and n=2. (a) The stable node of the PSR state is shown with its local dynamics calculated using the reduced dynamical system (see equation 3.13). (b) A trajectory, with transients removed, of the macroscopic mean field for the PSR state calculated from a network of N=10, 000 theta neurons using equations 2.2 to 2.4. (c) A magnification of the PSR state shown in panel b. The dimensions of the box are x=Re(z(t)): −0.5360 to −0.5300; y=Im(z(t)): −0.8345 to −0.8285. Fluctuations due to finite-size effects are small but visible in this zoomed-in view.

Figure 3:

Phase portraits of the PSR macroscopic state with , , k=−2, and n=2. (a) The stable node of the PSR state is shown with its local dynamics calculated using the reduced dynamical system (see equation 3.13). (b) A trajectory, with transients removed, of the macroscopic mean field for the PSR state calculated from a network of N=10, 000 theta neurons using equations 2.2 to 2.4. (c) A magnification of the PSR state shown in panel b. The dimensions of the box are x=Re(z(t)): −0.5360 to −0.5300; y=Im(z(t)): −0.8345 to −0.8285. Fluctuations due to finite-size effects are small but visible in this zoomed-in view.

In the PSS state, the macroscopic mean field z(t) settles onto a stable focus. The collective behavior of the infinite network is again at rest, but in this case, there is an intrinsic circulation, as the local stability of the equilibrium is given by a pair of complex eigenvalues (with negative real parts). This state occurs predominantly (but not exclusively) when most neurons inherently spike (), with the coupling being either excitatory (k>0) or weakly inhibitory (). Thus, although most neurons are active, the network is partially synchronous and organized such that phase cancellation among the neurons results in a well-defined stationary macroscopic mean.

Figure 4a shows an example of the macroscopic PSS state obtained using the reduced system, equation 3.13, with , , k=2, and n=2. (A movie showing both the macroscopic and microscopic behavior of the PSS state in Figure 4 is available in the online supplement.) As before, panels b and c show the mean field trajectory z(t) of a network of N=10, 000 neurons at the same parameter values. The reduced system once again accurately predicts the ultimate macroscopic state of the network, and finite-size network effects reveal bouts of coherent circulation about the focus.

Figure 4:

Phase portraits of the PSS macroscopic state with , , k=2, and n=2. (a) The stable focus of the PSS state is shown with its local dynamics calculated using the reduced model, equation 3.13. (b) A trajectory, with transients removed, of the macroscopic mean field for the PSS state calculated from a network of N=10, 000 theta neurons using equations 2.2 to 2.4. (c) A magnification of the the PSS state shown in panel b. The dimensions of the box are x=Re(z(t)): −0.2815 to −0.2415; y=Im(z(t)): −0.0250 to 0.0150. Fluctuations due to finite-size effects are small but visible in this zoomed-in view.

Figure 4:

Phase portraits of the PSS macroscopic state with , , k=2, and n=2. (a) The stable focus of the PSS state is shown with its local dynamics calculated using the reduced model, equation 3.13. (b) A trajectory, with transients removed, of the macroscopic mean field for the PSS state calculated from a network of N=10, 000 theta neurons using equations 2.2 to 2.4. (c) A magnification of the the PSS state shown in panel b. The dimensions of the box are x=Re(z(t)): −0.2815 to −0.2415; y=Im(z(t)): −0.0250 to 0.0150. Fluctuations due to finite-size effects are small but visible in this zoomed-in view.

Both the PSR and PSS states exhibit stationary behavior in the macroscopic mean field z(t) and reflect partially coherent network configurations. We emphasize here the subtle difference between them: one is a node in the macroscopic mean field and the other is a focus. This observation suggests that transient behavior in the macroscopic mean field z(t) resulting from abrupt perturbations of network parameters should reveal the difference between these two states.

Figure 5 shows time series of the macroscopic mean field z(t) for both the PSR (panels a and b) and the PSS (panels c and d) states. For the PSR state, the system starts with the following parameter set: , , k=−2, and n=2. Then, at t=500, is abruptly switched from −0.2 to −0.5. The new asymptotic state remains a PSR state (with Lyapunov exponents ), but the stable node shifts, and the macroscopic mean field z(t) converges exponentially toward the new asymptotic value. The time series from both the reduced system (see Figure 5a) and a discrete network of 10, 000 neurons (see Figure 5b) clearly demonstrates this exponential convergence.

Figure 5:

Time series of the real part of the macroscopic mean field x=Re(z(t)) showing the very different responses of the PSR and PSS states to a sudden small change in at t=500. (a) The behavior of the reduced equation, equation 3.13. (b) The time series calculated using a network of 10, 000 theta neurons for the PSR state. (c, d) The same for the PSS state. The horizontal dashed lines indicate the asymptotic values of the macroscopic equilibria at the initial and perturbed values. The parameter values are given in the main text.

Figure 5:

Time series of the real part of the macroscopic mean field x=Re(z(t)) showing the very different responses of the PSR and PSS states to a sudden small change in at t=500. (a) The behavior of the reduced equation, equation 3.13. (b) The time series calculated using a network of 10, 000 theta neurons for the PSR state. (c, d) The same for the PSS state. The horizontal dashed lines indicate the asymptotic values of the macroscopic equilibria at the initial and perturbed values. The parameter values are given in the main text.

The results from applying the same procedure to a PSS state (with , k=2, and n=2, and changing from 0.2 to 0.5) are shown in Figures 5c and 5d. In this case, the perturbed PSS state is characterized by a stable focus with a pair of stable complex eigenvalues (). Thus, the transient behavior after the parameter shift exhibits prominent oscillations that do not occur in the PSR case.

The PSR and PSS states are the expected macroscopic states when the network's intrinsic dynamics (as parameterized by the excitability parameter ) and the synaptic interactions (as characterized by the coupling k) are not in competition. The PSR state tends to occur when a majority of the neurons are intrinsically at rest () and the overall coupling is inhibitory (k<0). Conversely, the PSS state tends to occur when a majority of the neurons inherently spike () and the overall coupling is excitatory (k>0). In these cases, the internal dynamics and the network interaction reinforce each other, and the resulting macroscopic dynamics is a simple equilibrium. We show in the following section that a more complicated dynamical state can occur when the intrinsic neuronal dynamics and the network interaction compete with one another.

### 4.2.  Macroscopic Limit Cycle State and Multistability.

In the CPW state, the macroscopic mean field settles onto a stable limit cycle, and z(t) oscillates in time. We adopt this terminology based on previous work in which a collective oscillatory state of this type, when viewed from the microscopic perspective, was described as a wave (Crawford, 1994; Ermentrout, 1998; Bressloff, 1999; Osan, Rubin, & Ermentrout, 2002). Here, we find that this state occurs when most neurons are active () and the synaptic interaction is inhibitory (k<0). The microscopic configuration of the neurons is such that the degree of coherence waxes and wanes in time as the phases of the neurons corral together and spread apart in a periodic manner. Thus, the collective oscillation reflects the interplay between the neurons’ inherent tendency to spike and the suppressive network interaction. Indeed, we show below that the occurrence of CPW states is mediated by Andronov-Hopf and homoclinic bifurcations of the mean field, and thus are emergent properties of the network. In particular, the frequency of the limit cycle for a CPW state is not simply related to the frequencies of the individual neurons. In addition, the macroscopic limit cycle takes on different shapes and sizes for different system parameters, thus indicating different microscopic wave patterns.

A particular example of the CPW state is shown in Figure 6. As before, panel a shows the attractors predicted by equation 3.13 with , , k=−9, and n=2. (A movie showing both the macroscopic and microscopic behavior of the CPW state in Figure 6 is available in the online supplement.) In this case, an attracting CPW limit cycle and a PSR node coexist. Two unstable equilibria are present as well. Panel b shows the asymptotic mean field behavior of a network of 10, 000 neurons, where separate runs with different initial conditions were used to demonstrate the coexistence of the two attractors. Once again, the reduced mean field equation gives an excellent prediction for the asymptotic temporal behavior of the full network.

Figure 6:

(a) Phase portraits of the asymptotic macroscopic states that occur with network parameters , , k=−9, and n=2. (a) The reduced equation, equation 3.13, predicts the coexistence of a stable node (PSR), a saddle equilibrium, an unstable focus, and a stable limit cycle state (CPW) for the macroscopic mean field. (b) The asymptotic macroscopic states exhibited by a finite network with 10, 000 neurons. Two mean field trajectories showing the PSR and CPW states are shown; these were obtained with different initial conditions after transients were discarded. (c) A close-up of a section of the CPW limit cycle. The dimensions of the box are x=Re(z(t)): 0.5050 to 0.6550; y=Im(z(t)): −0.0750 to 0.0750. Fluctuations in the trajectory are due to finite-size effects.

Figure 6:

(a) Phase portraits of the asymptotic macroscopic states that occur with network parameters , , k=−9, and n=2. (a) The reduced equation, equation 3.13, predicts the coexistence of a stable node (PSR), a saddle equilibrium, an unstable focus, and a stable limit cycle state (CPW) for the macroscopic mean field. (b) The asymptotic macroscopic states exhibited by a finite network with 10, 000 neurons. Two mean field trajectories showing the PSR and CPW states are shown; these were obtained with different initial conditions after transients were discarded. (c) A close-up of a section of the CPW limit cycle. The dimensions of the box are x=Re(z(t)): 0.5050 to 0.6550; y=Im(z(t)): −0.0750 to 0.0750. Fluctuations in the trajectory are due to finite-size effects.

For most regimes in parameter space, the macroscopic behavior of the network is found to exclusively approach just one of the above defined states; there is only a single macroscopic attractor. However, there are significant parameter regions in which the network exhibits multistability, where two or more of these macroscopic states are found to coexist. Indeed, the example shown in Figure 6 is an example of multistability in which both a stable CPW state and a stable PSR state coexist. For parameters within these multistable regions, the network approaches one of the stable macrostates depending on how the neurons in the network are configured initially. We find, based on our bifurcation analysis of the mean field equation (see section 5), that dynamical competition is a necessary ingredient for the emergence of multistability. (A more detailed analysis of the multistable state for a similar but nonautonomous network of theta neurons is reported in So et al. (2013).) (A singular situation occurs with , corresponding to a homogeneous network of neurons. Due to the high degree of symmetry present in this case, the collective behavior consists of many coexisting neutrally stable limit cycles, and the overall macroscopic dynamics can be counterintuitively more complicated than the heterogeneous case studied here. A more detailed analysis of this homogeneous case and other extensions of this work will be reported elsewhere.)

One can also entertain the notion of unstable macroscopic PSR, PSS, or CPW states, as we already mentioned. Although these are not typically observable in the collective behavior of the physical network, we will demonstrate in the next section that they play an important role in mediating the transitions among the three classes of observable (attracting) macroscopic states.

## 5.  Bifurcation Analysis of the Macroscopic States

Having identified the three classes of attractors for the macroscopic mean field z(t), we now turn our attention to the analysis of the bifurcations that they can undergo. Specifically, we identify the bifurcations that occur as the following network parameters are varied: the neurons’ intrinsic excitability parameter , the heterogeneity parameter , and the overall coupling strength k. We consider both excitatory (k>0) and inhibitory (k<0) interaction among the neurons. The bifurcation set will be illustrated in the three-dimensional parameter space defined by , , and k, for fixed values of the synaptic sharpness parameter n. In our examples, we use n=2 and n=9, and our results suggest that the bifurcation scenarios described here are qualitatively robust with respect to n.

We begin by separating the reduced system, equation 3.13, into its real and imaginary parts, where z(t)=x(t)+iy(t):
5.1
Then, by setting the right side of both of these equations equal to zero, we obtain two conditions for the macroscopic equilibria of the network (xe, ye) as a function of the three network parameters:
5.2
Now, instead of solving equations 5.2 for xe and ye given particular values of , , and k, we consider xe, ye, , , and k to be five independent variables and think of equations 5.2 as two constraints that define a three-dimensional submanifold on which the equilibria must reside. Algebraic conditions for the occurrence of a particular kind of bifurcation provide additional constraints, thus defining a lower-dimensional surface (or surfaces) that characterizes the bifurcation of interest.

For a generic codimension 1 bifurcation such as the saddle node (SN) or the Andronov-Hopf (AH) bifurcation, this procedure results in two-dimensional surfaces embedded in the full five-dimensional space. We can visualize these two-dimensional bifurcation sets in the three-dimensional space defined by the network parameters , , and k. In sections 5.1 and 5.2, we examine the saddle node and the Andronov-Hopf bifurcations, respectively. We infer (and numerically verify) that homoclinic bifurcations are present as well. We describe in section 5.3 the transition between the PSR and the PSS states in which a macroscopic equilibrium changes from a node to a focus, or vice versa. We call this a node focus (NF) transition. This transition is not typically classified as a bifurcation in the traditional sense, since the stability of the equilibrium does not change, and additional states are neither created nor destroyed. Nevertheless, it is desirable to know where in parameter space this transition occurs, since the type of equilibrium (i.e., focus or node) can have macroscopic consequences, as illustrated in Figure 5. Collectively, these results lead to an understanding of the various bifurcations and transitions that occur in the attractors of the macroscopic mean field of our network.

The saddle node bifurcation is defined by the condition
5.3
where is the Jacobian of the system given by equations 5.1. Since our reduced equation is two-dimensional, all saddle node bifurcations that occur in our network must necessarily involve PSR states. This is because the creation of a pair of PSS equilibrium states requires at least three dimensions (two corresponding to the complex-conjugate eigenvalues and one along the heteroclinic connection). Note also that the above determinant condition includes the codimension 2 cusp bifurcation when both eigenvalues of J are zero simultaneously. The combination of the three algebraic constraints given in equations 5.2 and 5.3 allows us to solve for , , and k in terms of the remaining 2 degrees of freedom, xe and ye. We then plot the SN bifurcation surface parametrically in by considering all possible values of (xe, ye) within the allowed state space . The SN bifurcation surfaces obtained in this manner are displayed in Figure 7. Panels a and b show the surfaces obtained for synaptic sharpness parameters n=2 and n=9, respectively, and panel c is a magnification of panel a. Note that these figures extend into the unphysical region where . This is done to help the reader visualize the shape of the surfaces, as they are symmetric across .
Figure 7:

The saddle node (SN) bifurcation surfaces in the three-dimensional parameter space for synaptic sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). The rough edges in panel a are due to numerical limitations. Panel c is a magnification of panel a. The black line segments in panels a and c represent paths in parameter space that are keyed to Figures 8c and 11c.

Figure 7:

The saddle node (SN) bifurcation surfaces in the three-dimensional parameter space for synaptic sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). The rough edges in panel a are due to numerical limitations. Panel c is a magnification of panel a. The black line segments in panels a and c represent paths in parameter space that are keyed to Figures 8c and 11c.

The bifurcation set consists of two similar tentlike structures. The edges of the tentlike surfaces correspond to parameter values where a codimension 2 cusp bifurcation occurs. It is notable that these tentlike structures are predominantly (but not exclusively) located in regions where the internal excitability parameter and the coupling strength k are of opposite sign for both excitatory and inhibitory connectivity. This is the dynamically competitive region mentioned above. Furthermore, the similarity between the surfaces in a (for n=2) and b (for n=9) indicate the robustness of our results with respect to the synaptic sharpness parameter n.

Panel a includes two line segments—one parallel to the k axis (with ) and the other parallel to the axis (with ). Panel c is a magnification of panel a in the vicinity of the former, showing how this line segment pierces the SN surfaces. These line segments are keyed to Figures 8c and 11c, discussed below, in order to clarify which macroscopic states exist and which bifurcations occur as parameters are traversed along these lines.

Figure 8:

(a) A two-dimensional slice of the saddle node bifurcation set shown in Figure 7a with fixed at −0.3. Saddle node (SN) bifurcations occur on the solid curves, and these meet at a cusp point. (b) A one-dimensional bifurcation diagram (y=Im(z) versus k) corresponding to the dashed line in panel a at . The solid lines indicate stable PSR (node) or PSS (focus) equilibria as labeled, the broken line is an unstable PSR equilibrium, and the open diamonds indicate node focus (NF) transition points. The circle is the SN bifurcation corresponding to the right SN curve in panel a. The other SN point is so close to an NF point that they cannot be distinguished here and are marked SN/NF. (c) Schematic representation of the sequence of bifurcations (left) and macroscopic states (right) that occur as k traverses the range shown in panel b. This is the same as the vertical line in Figures 7a and c. The shaded interval indicates multistability.

Figure 8:

(a) A two-dimensional slice of the saddle node bifurcation set shown in Figure 7a with fixed at −0.3. Saddle node (SN) bifurcations occur on the solid curves, and these meet at a cusp point. (b) A one-dimensional bifurcation diagram (y=Im(z) versus k) corresponding to the dashed line in panel a at . The solid lines indicate stable PSR (node) or PSS (focus) equilibria as labeled, the broken line is an unstable PSR equilibrium, and the open diamonds indicate node focus (NF) transition points. The circle is the SN bifurcation corresponding to the right SN curve in panel a. The other SN point is so close to an NF point that they cannot be distinguished here and are marked SN/NF. (c) Schematic representation of the sequence of bifurcations (left) and macroscopic states (right) that occur as k traverses the range shown in panel b. This is the same as the vertical line in Figures 7a and c. The shaded interval indicates multistability.

Figure 8a shows a two-dimensional slice through the n=2 tent at . A typical fold structure with two saddle node curves meeting at a codimension 2 cusp point is seen. Panel b shows the one-dimensional bifurcation diagram, plotting y=Im(z) versus k, that results from following k along the line (dashed line in panel a; this is the same as the vertical line segment in Figures 7a and 7c). This diagram shows how the equilibrium solutions evolve as k increases from zero. Initially there is an attracting PSS state. This changes into a PSR state at the NF transition point indicated by the open diamond at k=0.1028. As k increases further, this PSR state gradually migrates toward higher values of y. Then a SN bifurcation and a NF transition point occur in rapid succession at k=0.9067 and k=0.9075, respectively (these points are not resolvable at the resolution shown in the figure and are therefore marked “SN/NF”). This SN bifurcation creates new stable and unstable PSR states in a separate region of state space (near y=−0.08), and at the NF point, the stable PSR changes into a stable PSS state. As k increases further, this stable PSS state persists, while the unstable PSR state migrates toward smaller values of y and collides with the coexisting stable PSR state. These annihilate each other via the SN bifurcation at k=1.1237. This sequence of events is shown schematically collapsed onto the k axis in panel c, so that it can be compared to the vertical line segment in Figures 7a and 7c. The shaded region indicates an interval in which more than one stable attractor exists. Note also that an NF transition occurs at a negative value of k (−0.5697) that is not visible in panel b.

### 5.2.  Andronov-Hopf Bifurcation.

The Andronov-Hopf bifurcation is defined, for our two-dimensional system, by two conditions,
5.4
Equations 5.4 combined with equations 5.2 give three equations for five unknowns, with the additional constraint that det[J] must be greater than zero. Proceeding as before, we obtain two-dimensional parametric plots of the AH bifurcation surface, shown in Figure 9. In this case, there is qualitative similarity between the shapes for n=2 (panel a) and n=9 (panel b) cases, but there are quantitative differences in the location of the surfaces. Note also the line segment included in panel a (with ); this is the same as the horizontal line segment that appears in Figure 7a and is keyed to Figure 11c.
Figure 9:

The Andronov-Hopf (AH) bifurcation surface in the three-dimensional parameter space for sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). The black line segment in panel a is the same as the horizontal line segment that appears in Figure 7a and represents a path in parameter space that is keyed to Figure 11c.

Figure 9:

The Andronov-Hopf (AH) bifurcation surface in the three-dimensional parameter space for sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). The black line segment in panel a is the same as the horizontal line segment that appears in Figure 7a and represents a path in parameter space that is keyed to Figure 11c.

The result is a tube or funnel-shaped surface that opens and flattens out on one side. The funnel emanates from the regime of large inhibitory coupling and less heterogeneity with (i.e., most neurons are very close to their SNIC bifurcations), and then opens up and flattens out for increasing values of (i.e., greater dominance of spiking neurons). As in the case of the SN bifurcation, the surface occurs most prominently where there is dynamic competition within the network. However, in this case, the surface exists only where the competition is specifically between predominantly active neurons and inhibitory network interaction ( and k<0).

Figure 10a shows the two-dimensional bifurcation diagram that results from slicing through the n=2 AH and SN surfaces at k=−9. The two SN curves again meet at a cusp, and the AH curve intersects the left SN curve at a codimension 2 Bogdanov-Takens (BT) point. The dotted rectangular region shown in panel a is magnified in panel b, making it easier to see the AH curve, as well as the homoclinic (HC) bifurcation curve that also emerges from the BT point (Kuznetsov, 2004). We identified the latter curve numerically.4

Figure 10:

(a) Superimposed two-dimensional slices of both the SN (see Figure 7a) and AH (see Figure 9a) bifurcation sets at k=−9. The saddle node (SN) curves meet at a cusp point, and a Bogdanov-Takens (BT) point (triangle) occurs on the left SN curve. The dotted rectangular region in panel a is magnified in panel b, showing the Andronov-Hopf (AH) and homoclinic (HC) bifurcation curves. The HC curve was interpolated from the points indicated by the circles, which were found numerically.

Figure 10:

(a) Superimposed two-dimensional slices of both the SN (see Figure 7a) and AH (see Figure 9a) bifurcation sets at k=−9. The saddle node (SN) curves meet at a cusp point, and a Bogdanov-Takens (BT) point (triangle) occurs on the left SN curve. The dotted rectangular region in panel a is magnified in panel b, showing the Andronov-Hopf (AH) and homoclinic (HC) bifurcation curves. The HC curve was interpolated from the points indicated by the circles, which were found numerically.

To further clarify the identity of the macroscopic network states, Figure 11a shows the one-dimensional bifurcation diagram (in this case, x=Re(z) versus ) obtained by varying along the line (dashed lines in Figure 10; these are the same as the horizontal lines in Figures 7a and 9a). Here, the heavy solid lines represent stable equilibria. The lower equilibrium branch corresponds to the PSR state, and it persists until it collides with an unstable PSR state in a saddle node bifurcation. Moving along the upper stable equilibrium with decreasing , the network exhibits the PSS state before encountering the AH bifurcation, which is supercritical. At this point, the equilibrium loses stability, and an attracting limit cycle emerges: the CPW state. The amplitude of this limit cycle subsequently increases until it collides with the unstable (uPSR) equilibrium in a homoclinic bifurcation.

Figure 11:

(a) A one-dimensional bifurcation diagram (x=Re(z) versus ) corresponding to the dashed lines in Figure 10 at . The solid (dashed) curves are stable (unstable) PSR and PSS equilibria as indicated. SN denotes a saddle node bifurcation. The solid black circles indicate the maximum and minimum values of x for the CPW limit cycle that emerges from the supercritical Andronov-Hopf (AH) point. The existence of this limit cycle also transitions at the indicated homoclinic (HC) bifurcation. (b) A magnification of the region near the node focus (NF) point, showing the other SN bifurcation. (c) Schematic representation of the sequence of bifurcations (top) and macroscopic states (bottom) that occur as traverses the range shown in panel a; this is the same as the horizontal lines in Figures 7a and 9a. The shaded interval indicates multistability.

Figure 11:

(a) A one-dimensional bifurcation diagram (x=Re(z) versus ) corresponding to the dashed lines in Figure 10 at . The solid (dashed) curves are stable (unstable) PSR and PSS equilibria as indicated. SN denotes a saddle node bifurcation. The solid black circles indicate the maximum and minimum values of x for the CPW limit cycle that emerges from the supercritical Andronov-Hopf (AH) point. The existence of this limit cycle also transitions at the indicated homoclinic (HC) bifurcation. (b) A magnification of the region near the node focus (NF) point, showing the other SN bifurcation. (c) Schematic representation of the sequence of bifurcations (top) and macroscopic states (bottom) that occur as traverses the range shown in panel a; this is the same as the horizontal lines in Figures 7a and 9a. The shaded interval indicates multistability.

Figure 11b shows a magnification of the vicinity of the SN/NF point in panel a, showing both the SN and NF points distinctly. This SN point corresponds to the left SN curve in Figure 10a and in this case, leads to the creation of two unstable PSR states.

Finally, the sequence of events in Figure 11a is shown schematically collapsed onto the axis in panel c, so that it can be compared to the horizontal line segments shown in Figures 7a and 9a.

### 5.3.  PSR to PSS Transition.

In Figures 8 and 11, we indicated that the stable equilibria depicted there make transitions between the PSR and PSS states. These transitions occur when the equilibrium changes from a node (PSR state; real eigenvalues) to a focus (PSS state; complex eigenvalues), or vice versa. This change is not normally considered a bifurcation since there is no change in stability or creation of other states. Nevertheless, it is possible to identify surfaces in the parameter space that correspond to this transition. We call this transition the node focus (NF) transition.

The NF transition occurs when the discriminant of the characteristic equation of the Jacobian equals zero, thus signifying the presence of equilibria with real eigenvalues of multiplicity 2:
5.5
To identify the transition surface, we proceed as before by directly plotting the two-dimensional parametric surface in the three-dimensional parameter space using the three algebraic constraints given in equations 5.2 and 5.5. The result is shown in Figure 12. As before, panels a and b show the surfaces for n=2 and n=9, respectively, and panel c is a magnification of a. The horizontal and vertical lines of Figures 7a and 9a, which are linked to Figures 8c and 11c, are included.
Figure 12:

The node focus (NF) surface in the three-dimensional parameter space for sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). (c) A magnification of panel a. The black line segments in panels a and c represent paths in parameter space that are keyed to Figures 8c and 11c.

Figure 12:

The node focus (NF) surface in the three-dimensional parameter space for sharpness parameter (a) n=2 and (b) n=9. To aid in visualization, the figures extend into the unphysical region where (the surfaces are symmetric across ). (c) A magnification of panel a. The black line segments in panels a and c represent paths in parameter space that are keyed to Figures 8c and 11c.

Figure 12 reveals two surfaces: a lower surface with an internal pleat somewhat like a fortune cookie and an upper folded surface like the nose cone of an airplane. The PSS state occurs in the region in the far upper right corner of panels a and b. Here, the network dynamics is cooperative in that predominantly spiking neurons () interact via excitatory synapses (k>0), leading to an active network. In contrast, in the far lower left corner of these figures, predominantly resting neurons () interact cooperatively via inhibitory synapses (k<0), and the network exhibits the quiescent PSR state. The lower fortune cookie–like surface marks a NF boundary between these two regions, but the SN and AH bifurcations discussed above occur nearby as well. Interestingly, however, the upper nose cone surface encloses another region of PSR states. Within the nose cone, networks consist of predominantly resting but excitable neurons interacting via weak excitatory synapses. In this case, the resting states of most neurons are relatively far from their thresholds (), so that weak synaptic excitation is not sufficient to cause most neurons to fire. Thus, the network exhibits the PSR state.

It is important to note that the NF transition applies to both stable and unstable equilibria. In the one-dimensional bifurcation diagram of Figure 8b (with , ), the network transverses through both the lower and upper parts of the nose cone as k increases from 0 to 2.5. The intersection with the lower, fortune cookie–like surface occurs for a negative value of k, as indicated in Figure 8c. All these NF transitions involve stable equilibria.

In contrast, the PSR–PSS transition in Figure 11, visible in panel b, involves an unstable equilibrium. This transition corresponds to traversing the fortune cookie–like surface along the horizontal line shown in Figure 12c (at and k=−9; same as the horizontal lines in Figures 8 and 11). The NF transition occurs as the repeller (unstable PSR) that was created at the nearby saddle node bifurcation (at ) changes into an unstable focus (unstable PSS) at . Although these unstable PSR and PSS states are unobservable in the full network, the stabilization of the PSS state through the AH bifurcation at requires the preexistence of this unstable PSS state.

## 6.  Summary and Discussion

Using the well-known theta neuron model, we constructed a heterogeneous network containing a mixture of at-rest but excitable neurons as well as spontaneously spiking neurons. These were globally coupled together through pulselike synaptic interactions whose shape depends on a sharpness parameter n. We applied a recently developed reduction technique to derive a low-dimensional dynamical equation that completely describes the asymptotic behavior of the network's mean field in the thermodynamic limit of large network size. By analyzing this reduced system, we identified not only all possible asymptotic states of the mean field, but also the bifurcations that occur as three network parameters are varied. We also showed that the predicted behavior derived in this way is exhibited by finite networks of 10, 000 neurons.

Of course, real neurons are not theta neurons. We caution that real type I neurons can be expected to be well approximated by theta neurons only near the onset of spiking. Future work will examine the extent to which our results carry over to networks of more realistic type I neurons. Nevertheless, the theta neuron model does capture important aspects of the generic behavior of type I neurons, including pyramidal neurons. These make up the majority of neurons in the mammalian brain (e.g., approximately 80% in the hippocampus and 70% in the cortex, across many species; Feldman, 1984; Chen & Dzakpasu, 2010). Also, real neuronal networks are clearly not globally coupled. But since our analysis focuses on the dynamics of the macroscopic mean field variable defined in equation 3.3, our results are best interpreted as providing an idealized mesoscopic description of neuronal dynamics (i.e., describing the behavior of a region of brain tissue that is small but nevertheless contains many neurons and many connections). Notable strengths of our network structure are that it includes heterogeneity in the internal dynamics of the individual neurons (i.e., both resting and spiking neurons) and that the coupling is via pulselike synapses of tunable sharpness that approximate real IPSPs without being delta functions. It is also important to note that our approach provides exact results for the asymptotic mean field behavior, although only in the limit of large system size. Nevertheless, we showed that the features and dynamics we derived are in fact exhibited by reasonably sized finite network instantiations.

We found that the asymptotic mean field of our network exhibits only three possible states: two corresponding to equilibrium solutions, the PSR and PSS states, and one limit cycle solution, the CPW state. The most obvious solution is perhaps the PSR state, in which the network consists of predominantly resting neurons that inhibit each other. In the most extreme case, the neurons form a single stationary cluster, and the mean field assumes a constant value. However, PSR and PSS states can also be realized such that all or a significant portion of the neuronal population spike regularly. Despite this spiking activity, the mean field for such partially synchronized states remains constant. These states are similar to asynchonous states that have been described by others (e.g., Abbott & van Vreeswijk, 1993; Hansel & Mato, 2001, 2003). For example, Abbott and van Vreeswijk (1993) defined the asynchronous state by the requirement that the total input into one neuron from all the others be constant, and they argued that the existence of these states justifies the use of firing rate models.

We have also described the subtle difference between the PSR and PSS states—that the former is a node and the latter is a focus. The consequences of this distinction can be observed in Figure 5, where it is shown that the mean field responses to perturbations are markedly different. When a PSR state is perturbed, the mean field relaxes to the equilibrium directly. In contrast, when a PSS state is perturbed, the mean field displays decaying oscillations. In this latter case, as the equilibrium microscopic configuration is approached, the neurons alternate between bouts of scattering and clumping, or desynchronization and resynchronization, in a manner such that each bout is less severe than the preceding one.

Similar microscopic dynamics underlie the CPW state, except that for the CPW state, the alternation between the de- and resynchronizing bouts persists indefinitely. Consequently, the asymptotic mean field approaches a limit cycle. This state is related to, but is more general than, the synchronous state described by Wang and Buzsáki (1996). The latter synchronous state occurs for homogeneous (or very weakly heterogeneous) networks when the phases of most neurons lock, as in the partially synchronous state shown in their Figure 5C. In this state, almost all neurons fire together. Thus, the order parameter of such a network exhibits a CPW state with a constant order parameter magnitude very close to one and a frequency of oscillation identical to that of an individual neuron. In contrast, the CPW states that we describe in this letter are more general in that typically, the neurons do not phase-lock. Instead, they form clusters of periodically varying coherence, as is quantified by the oscillating order parameter magnitude. The animation provided in the online supplement makes this point clear. It is important to note that the frequency of oscillation exhibited by the mean field while on a CPW limit cycle is an emergent property of the network as a whole. This frequency is therefore not related to the frequencies of the individual neurons in a simple manner. Indeed, in parameter space regions near the Bogdanov-Takens bifurcation, the CPW frequency can be very small. Such low-frequency rhythms have been reported in EEG (Avella Gonzalez et al., 2012), and slow population activities have been observed in hippocampal slice studies (Ho, Strüber, Bartos, Zhang, & Skinner, 2012).

These findings have implications for the interpretation of experimental measurements of network dynamics. For example, raster diagrams are often used to identify the presence or absence of synchrony. Our results suggest that such diagrams may not be sufficient to differentiate between the three macroscopic states that we have identified here.

The reduced equation, equation 3.13, also gives us the ability to fully analyze the relevant transitions among these macroscopic states using the degree of excitability, heterogeneity, and coupling strength (both excitatory and inhibitory) as bifurcation parameters. The results for two generic codimension 1 bifurcations, the saddle node and Andronov-Hopf bifurcations, were summarized visually as two-dimensional surfaces in Figures 7 and 9, respectively. We also examined the node focus transition similarly (see Figure 12). These results are exact only in the limit of large system size. In finite network realizations, fluctuations can build up when the network is near a bifurcation. Hence, the occurrence of these bifurcations can be expected to be blurred somewhat in the parameter space. Nevertheless, we found that the bifurcations we describe here can be reliably identified in networks of 10, 000 theta neurons. In this letter, we did not investigate smaller networks.

By examining these bifurcation surfaces, we found that the dynamical interplay between the internal neuronal dynamics, as parameterized by , and the interneuronal coupling, parameterized by k, appears to play an important role in determining the complexity of the possible macroscopic mean field dynamics of the network. In particular, when the network is preferentially cooperative, such that and k are of the same sign, the network tends to settle into one of the macroscopic equilibrium states. In contrast, when and k are of opposite sign, richer dynamics are seen, including the CPW state, bifurcations of the macroscopic mean field, and multistability.

Finally, we note that the results we reported were for two specific values of the sharpness parameter: n=2 and n=9. Further analysis of the effect of increasing this parameter shows no qualitative change in our results. Specifically, we found no changes to the three classes of macroscopic states and no significant qualitative changes to the bifurcation surfaces for values of n up to n=15. While we suspect that there may be important differences that arise as , we note that low values of the sharpness parameter are actually better approximations of real postsynaptic potentials.

Our results provide general predictions for the macroscopic dynamics of an idealized large network of type I neurons. In a related study (So et al., 2013), a similar but nonautonomous theta neuron network in which was made to oscillate in time was investigated. Varying this parameter such that the SNIC bifurcation is repeatedly traversed is a common method of modeling bursting neurons, specifically parabolic bursters. That letter builds on the results reported here and shows that more complicated behavior, including macroscopic quasiperiodicity, chaos, multistability, and final-state uncertainty, can occur in the nonautonomous case.

Future work will extend our results to include separate and different interacting populations of neurons, including type II neurons, in order to model the excitatory and inhibitory and layered structure of the brain's neuronal networks.

## References

Abbott
,
L. F.
, &
van Vreeswijk
,
C.
(
1993
).
Asynchronous states in networks of pulse-coupled oscillators
.
Phys. Rev. E
,
48
,
1483
1490
.
Abdulrehem
,
M. M.
, &
Ott
,
E.
(
2009
).
Low dimensional description of pedestrian-induced oscillation of the Millennium Bridge
.
Chaos
,
19
,
013129
.
Abrams
,
D. M.
,
Mirollo
,
R.
,
Strogatz
,
S. H.
, &
Wiley
,
D. A.
(
2008
).
Solvable model for chimera states of coupled oscillators
.
Phys. Rev. Letts
.,
101
,
084103
.
Alonso
,
L. M.
, &
Mindlin
,
G. B.
(
2011
).
Average dynamics of a driven set of globally coupled excitable units
.
Chaos
,
21
,
023102
.
Ariaratnam
,
J. T.
, &
Strogatz
,
S. H.
(
2001
).
Asynchronous states in networks of pulse-coupled oscillators
.
Phys. Rev. E
,
48
,
1483
1490
.
Avella Gonzalez
,
O. J.
,
van Aerde
K. I.
,
van Elburg
,
R.A.J.
,
Poil
,
S-S
,
Mansvelder
,
H. D.
,
,
K.
, et al
(
2012
).
External drive to inhibitory cells induces alternating episodes of high- and low-amplitude oscillations
.
PLoS Computational Biology
,
8
(
8
),
e1002666
.
Bressloff
,
P. C.
(
1999
).
Synaptically generated wave propagation in excitable neural media
.
Phys. Rev. Letts.
,
82
,
2979
2982
.
Chen
,
X.
, &
Dzakpasu
,
R.
(
2010
).
Observed network dynamics from altering the balance between excitatory and inhibitory neurons in cultured networks
.
Phys. Rev. E
,
82
,
31907
.
Crawford
,
J. D.
(
1994
).
Amplitude expansions for instabilities in populations of globally-coupled oscillators
.
J. Stat. Phys.
,
74
,
1047
1084
.
Daido
,
H.
(
1996
).
Onset of cooperative entrainment in limit-cycle oscillations with uniform all-to-all interactions: Bifurcation of the order function
.
Physica D
,
91
,
24
66
.
Ermentrout
,
G. B.
(
1996
).
Type I membranes, phase resetting curves, and synchrony
.
Neural Comput.
,
8
,
979
1001
.
Ermentrout
,
G. B.
(
1998
).
The analysis of synaptically generated traveling waves
.
J. Comp. Neurosci
,
5
,
191
208
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1986
).
Parabolic bursting in an excitable system coupled with a slow oscillation
.
SIAM J. Appl. Math.
,
46
,
233
253
.
Feldman
,
M. L.
(
1984
).
Cellular components of the cerebral cortex
.
New York
:
Plenum Press
.
Hansel
,
D.
, &
Mato
,
G.
(
2001
).
Existence and stability of persistent states in large neuronal networks
.
Phys. Rev. Letts.
,
86
,
4175
4178
.
Hansel
,
D.
, &
Mato
,
G.
(
2003
).
Asynchronous states and the emergence of synchrony in large networks of interacting excitatory and inhibitory neurons
.
Neural Comput.
,
15
,
1
56
.
Harris
,
K. D.
(
2005
).
Neural signatures of cell assembly organization
.
Nature Reviews Neuroscience
,
6
,
399
407
.
Hebb
,
D. O.
(
1949
).
The organization of behavior: A neuropsychological theory
.
New York
:
Wiley
.
Ho
,
E.C.Y
,
Strüber
,
M.
,
Bartos
,
M.
,
Zhang
,
L.
, &
Skinner
,
F. K.
(
2012
).
Inhibitory networks of fast-spiking interneurons generate slow population activities due to excitatory fluctuations and network multistability
.
Journal of Neuroscience
,
32
(
29
),
9931
9946
.
Hodgkin
,
A. L.
(
1948
).
The local electric charges associated with repetitive action in non-modulated axons
.
J. Physiol.
,
107
,
165
181
.
Izhikevich
,
E.
(
2007
).
Dynamical systems in neuroscience: The geometry of excitability and bursting
.
Cambridge, MA
:
MIT Press
.
Kanamaru
,
T.
, &
Masatoshi
,
S.
(
2003
).
Analysis of globally connected active rotators with excitatory and inhibitory connections using Fokker-Planck equation
.
Physical Review E
,
67
,
031916
.
Kuramoto
,
Y.
(
1975
).
Lecture notes in physics
. In
H. Araki
(Ed.),
International Symposium on Mathematical Problems in Theoretical Physics
,
39
.
Berlin
:
Springer-Verlag
.
Kuramoto
,
Y.
(
1984
).
Chemical oscillations, waves and turbulence
.
New York
:
Springer
.
Kuznetsov
,
Y. A.
(
2004
).
Elements of applied bifurcation theory
(3rd ed.).
New York
:
Springer
.
Martens
,
E. A.
,
Barreto
,
E.
,
Strogatz
,
S. H.
,
Ott
,
E.
,
So
,
P.
, &
Antonsen
,
T. M.
(
2009
)
Exact results for the Kuramoto model with a bimodal frequency distribution
.
Physical Review E
,
79
,
026204
.
Marvel
,
S. A.
,
Mirollo
,
R. E.
, &
Strogatz
,
S. H.
(
2009
).
Identical phase oscillators with global sinusoidal coupling evolve by Mobius group action
.
Chaos
,
19
,
043104
.
Marvel
S. A.
, &
Strogatz
,
S. H.
(
2009
).
Invariant submanifold for series arrays of Josephson junctions
.
Chaos
,
19
,
013132
.
Montbrió
,
E.
, &
Pazó
,
D.
(
2011
).
Shear diversity prevents collective synchronization
.
Phys. Rev. Letts.
,
106
,
254101
Nowak
,
L. G.
,
Azouz
,
R.
,
Sanchez-Vives
,
M. V.
,
Gray
,
C. M.
, &
McCormick
,
D. A.
(
2003
).
Electrophysiological classes of cat primary visual cortical neurons in vivo as revealed by quantitative analyses
.
J. Neurophysiol.
,
89
,
1541
1566
.
Omel'chenko
,
O. E.
, &
Wulfrum
,
M.
(
2012
).
Nonuniversal transitions to synchrony in the Sakaguchi-Kuramoto model
.
Phys. Rev. Letts.
,
109
,
164101
.
Osan
,
R.
,
Rubin
,
J.
, &
Ermentrout
,
G. B.
(
2002
).
Regular traveling waves in a one-dimensional network of theta neurons
.
SIAM J. Appl. Math.
,
62
,
1197
1221
.
Ott
,
E.
, &
Antonsen
,
T. M.
(
2008
).
Low dimensional behavior of large systems of globally coupled oscillators
.
Chaos
,
18
,
037113
.
Ott
,
E.
, &
Antonsen
,
T. M.
(
2009
).
Long time evolution of phase oscillator systems
.
Chaos
,
19
,
023117
.
Peretto
,
P.
(
1984
).
Collective properties of neural networks: A statistical physics approach
.
Biol. Cybern.
,
50
,
51
62
.
Pikovsky
,
A.
, &
Rosenblum
,
M.
(
2008
).
Partially integrable dynamics of hierarchical populations of coupled oscillators
.
Phys. Rev. Lett.
,
101
,
264103
.
Pikovsky
,
A.
, &
Rosenblum
,
M.
(
2011
).
Dynamics of heterogeneous oscillator ensembles in terms of collective variables
.
Physica. D
,
240
,
872
881
.
Sherrington
,
C. S.
(
1906
).
The integrative action of the nervous system
.
New York
:
Scribners
.
Skinner
,
F. K.
(
2013
).
Moving beyond type I and type II neuron types [v1; ref status: indexed, http://f1000r.es/w7]
.
F1000Research
,
2
:
19. doi:10.3410/f1000research.2-19.v1
.
So
,
P.
, &
Barreto
,
E.
(
2011
).
Generating macroscopic chaos in a network of globally coupled phase oscillators
.
Chaos
,
21
,
033127
.
So
,
P.
,
Luke
,
T.
, &
Barreto
,
E.
(
2013
).
Networks of theta neurons with time-varying excitability: Macroscopic chaos, multistability, and final-state uncertainty
.
Physica D
.
doi:10.1016/j.physd.2013.04.2009
Sompolinsky
,
H.
(
1988
).
Statistical mechanics of neural networks
.
Physics Today
,
41
,
70
82
.
Tateno
,
T.
,
Harsch
,
A.
, &
Robinson
,
H.P.C.
(
2004
).
Threshold firing frequency-current relationships of neurons in rat somatosensory cortex: Type 1 and Type 2 dynamics
.
J. Neurophysiol.
,
92
,
2283
2294
.
Wang
,
X.-J.
, &
Buzsáki
,
G.
(
1996
).
Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model
.
J. Neurosci.
,
16
(
20
),
6402
6413
.

## Notes

1

Most theoretical studies consider only these two stereotypical behaviors, but see Skinner (2013).

2

In equation 3.7, the qth powers of the order parameter z and their complex conjugates appear. More generally, these terms should be replaced by the Daido moments (Daido, 1996). However, with the the choice of the Lorentzian distribution for , we have zq=zq for , and for q<0.

3

This method of analysis has been applied to other coupled networks of similar form (see, e.g., Pikovsky & Rosenblum, 2008, 2011; Abrams, Mirollo, Strogatz, & Wiley, 2008; Marvel & Strogatz, 2009; Martens, Barreto, Strogatz, Ott, So, & Antonsen, 2009; Abdulrehem & Ott, 2009; So & Barreto, 2011; Montbrió & Pazó, 2011; Alonso & Mindlin, 2011; Omel'chenko & Wulfrum, 2012; So, Luke, & Barreto, 2013).

4

Because the HC bifurcation is a global bifurcation, it cannot be specified by a set of simple algebraic conditions.