## 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}

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.

*I*as a collective signal in which each neuron contributes a pulselike synaptic current depending on its phase angle. Thus, where , , and

_{syn}*a*is a normalization constant such that The parameter

_{n}*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.

### 2.3. Heterogeneity.

## 3. Mean Field Reduction of the Network

*F*is governed by the continuity equation, where , the velocity of a neuron, is the continuum version of equations 2.2 and 2.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).

*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*)=

*I*/

_{syn}*k*is the rescaled synaptic current given by

^{2}where and 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.

*F*can be written as a Fourier expansion in which the coefficients appear as powers of a single yet-to-be-determined complex function : 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): 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

*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.

*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: 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 .

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.

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.

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.

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*.

*z*(

*t*)=

*x*(

*t*)+

*iy*(

*t*): 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 (

*x*,

_{e}*y*) as a function of the three network parameters: Now, instead of solving equations 5.2 for

_{e}*x*and

_{e}*y*given particular values of , , and

_{e}*k*, we consider

*x*,

_{e}*y*, , , and

_{e}*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.

### 5.1. Saddle Node Bifurcation.

*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,

*x*and

_{e}*y*. We then plot the SN bifurcation surface parametrically in by considering all possible values of (

_{e}*x*,

_{e}*y*) 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

_{e}*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 .

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 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.

**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.

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}

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.

### 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.

*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 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

## Notes

^{1}

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

^{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.