Phase-Resetting Curves Determine Synchronization, Phase Locking, and Clustering in Networks of Neural Oscillators

Srisairam Achuthan; Carmen C. Canavier

doi:10.1523/JNEUROSCI.0426-09.2009

Abstract

Networks of model neurons were constructed and their activity was predicted using an iterated map based solely on the phase-resetting curves (PRCs). The predictions were quite accurate provided that the resetting to simultaneous inputs was calculated using the sum of the simultaneously active conductances, obviating the need for weak coupling assumptions. Fully synchronous activity was observed only when the slope of the PRC at a phase of zero, corresponding to spike initiation, was positive. A novel stability criterion was developed and tested for all-to-all networks of identical, identically connected neurons. When the PRC generated using N − 1 simultaneously active inputs becomes too steep, the fully synchronous mode loses stability in a network of N model neurons. Therefore, the stability of synchrony can be lost by increasing the slope of this PRC either by increasing the network size or the strength of the individual synapses. Existence and stability criteria were also developed and tested for the splay mode in which neurons fire sequentially. Finally, N/M synchronous subclusters of M neurons were predicted using the intersection of parameters that supported both between-cluster splay and within-cluster synchrony. Surprisingly, the splay mode between clusters could enforce synchrony on subclusters that were incapable of synchronizing themselves. These results can be used to gain insights into the activity of networks of biological neurons whose PRCs can be measured.

Introduction

We hypothesize that the phase-resetting curve (PRC) contains all information required to predict the activity of a synaptically coupled network of oscillatory neurons. The PRC tabulates transient changes in period produced by an action potential in the presynaptic neuron as a function of the phase at which the perturbation is applied. We determine which aspects of the PRC are critical for generating a fully synchronous mode (Fig. 1A), a splay mode in which neurons fire sequentially (Fig. 1B), and clustering modes that exhibit synchrony within a cluster but phase locking at nonzero phase between clusters (Fig. 1C).

Figure 1.

Firing modes in an all-to-all four-neuron network. Only the most commonly observed, rotationally symmetric modes are shown. A, Fully synchronous mode (one cluster of four neurons). B, Splay mode in which all neurons fire in a sequential manner. C, Two clusters of two neurons in antiphase with each other. The vertical dashed lines indicate the intervals between firing times.

Many examples of synchrony and phase locking can be found in the CNS. Central pattern generators for repetitive motor neuron behavior frequently exhibit approximate splay modes (Mulloney and Hall, 2007; Smarandache and Mulloney, 2008), including antiphase (Stein et al., 1997). Transiently synchronized assemblies of neurons are believed to underlie cognitive functions (Buzsáki, 2006). Pathological synchrony could lead to epileptic seizures (Huguenard and McCormick, 2007; Traub and Jefferys, 1994) and tremor (Hammond et al., 2007), and there is consistent evidence for a reduction of synchronization in schizophrenia (Uhlhaas and Singer, 2006). A better understanding of how the phase resetting of the component oscillators controls network activity may shed light on many biological processes involving neural synchronization. For example, we show how a fully synchronous cluster at one frequency can break into two clusters in antiphase, which doubles the frequency that could be observed using the local field potential or EEG.

We assume only that each component of the network is a pacemaker (Winfree, 1967), that the effect of one perturbation dies out before the next one is received, and that the perturbations received by an oscillator in a closed circuit are similar to those used to generate its PRC. Much theoretical work (Ermentrout and Kopell, 1990, 1991; Van Vreeswijk et al., 1994; Ermentrout, 1996; Li et al., 2003; Netoff et al., 2005a) is based on the assumption of weak coupling, which implies that many cycles of the oscillation are required for synchronization or phase locking. Neural synchrony is frequently achieved rapidly, even within a single cycle (Singer, 1999); in these cases, strong coupling is required.

We propose two ways to extract information from the PRCs: an iterative map and an analytical method. The iterative map uses the PRCs and firing times in the previous cycle to predict firing times in the next (and subsequent) cycle(s) of the network. The analytical method directly calculates the phasic relationships and stability of oscillatory modes that are exhibited by the network from the PRCs. Both methods are tested on networks of model neurons that encompass generic excitability and PRC types. The simple, homogeneous model networks used herein to establish proof of principle do not capture all the complexities of neural dynamics; the point of these examples is that if the appropriate PRCs can be obtained, no further information is required to understand network synchronization properties due to reciprocal pulsatile coupling of intrinsic oscillators [see Nowotny et al. (2008) for other mechanisms of synchronization]. The solutions for homogeneous networks approximate those of mildly heterogeneous networks, and we suggest that modulation of the PRC in such networks may be a viable strategy to promote or discourage synchronization.

Materials and Methods

Simulations.

Two neural models were chosen as exemplars of different excitability types as originally described by Hodgkin (1948). Neurons with Hodgkin's type I excitability exhibit a gradual transition from quiescence to repetitive spiking as the magnitude of an applied depolarizing current step is increased. At the onset of spiking, the frequency is arbitrarily low. Oscillators with type I excitability act as integrators [Izhikevich (2007), Chapter 1] and generally have only advances due to brief excitations and delays due to brief inhibitions (Ermentrout, 1996). Neurons with Hodgkin's type II excitability exhibit an abrupt transition from quiescence to repetitive spiking as the magnitude of an applied depolarizing current step is increased. Spiking cannot be sustained below a threshold frequency. Oscillators with type II excitability act as resonators [Izhikevich (2007), Chapter 1] and can exhibit both advances or delays in response to either excitation or inhibition depending upon the timing of the input (Ermentrout, 1996). The Wang and Buzsáki (WB) model as commonly parameterized is type I. Although the Morris–Lecar (ML) model can be parameterized in either regime, here we have used type II parameters. All networks were comprised of identical, identically connected neurons for simplicity, but the methods can accommodate heterogeneity (Maran and Canavier, 2008).

The differential equations for the Wang and Buzsáki (1996) and Morris and Lecar (1981) networks were simulated using a variable step size implicit fifth-order Runge–Kutta method (Hairer and Wanner, 1991). The current balance equation for each WB model neuron is as follows: where the capacitance C = 1 μF/cm², V is the cell membrane voltage in millivolts, and t is time in milliseconds. The leak current is given by I_L = g_L(V − E_L). The sodium current is given by I_Na = g_Nam_∞³h(V − E_Na). The steady-state activation m_∞ = α_m/(α_m + β_m), where α_m(V) = −0.1(V + 35)/{exp[−0.1(V + 35)] − 1} and β_m(V) = 4exp[−(V + 60)/18]. The rate equation for the inactivation variable h in the expression for sodium current is as follows: where ϕ = 5. (The symbol was changed from φ to ϕ to avoid confusion with the symbol for phase.) The rate constants for the inactivation variable h are given by α_h(V) = 0.07exp[−(V + 58)/20] and β_h(V) = 1/{exp[−0.1(V + 28)] + 1}. The potassium current is given by I_K = g_Kn⁴(V − E_K), where the activation variable n satisfies the following equation: where the rate constants for n are α_n(V) = −0.01(V + 34)/{exp[−0.1(V + 34)] − 1} and β_n(V) = 0.125exp[−(V + 44)/80]. The reversal potentials E_Na, E_K, and E_L were set to 55, −90, and −65 mV, respectively. The maximal sodium (g_Na), potassium (g_K), and leak (g_L) conductances were set to 35, 9, and 0.1 mS/cm², respectively. I_stim is the applied current and was set at 0.5 μA/cm². Unless otherwise stated, the values for the various parameters were equal to those given above.

The synaptic current is given by I_syn = g_syns(V − E_syn), where g_syn is the maximum synaptic conductance and E_syn is equal to −75 mV for inhibitory synaptic connectivity and equal to 0 mV for excitatory synaptic connectivity. The rate of change of the gating variable s is given by the following equation: with T(V_pre) = 1/[1 + exp(−V_pre/2)], where V_pre is the voltage of the presynaptic cell, α = 6.25 ms⁻¹ is the rate constant of the synaptic activation (Bartos et al., 2001), and τ_syn is the synaptic decay time constant and was set to 1.0 ms.

The current balance equation for each ML model neuron is as follows: where the capacitance C = 20 μF/cm², V is the cell membrane voltage in millivolts, and t is time in milliseconds. The calcium current is given by I_Ca = g_Cam_∞(V)(V − E_Ca). The leak current is given by I_L = g_L(V − E_L). The steady-state activation is as follows: where V₁ = −1.2 mV and V₂ = 18 mV. The potassium current is given by I_K = g_Kw(V − E_K). The rate equation for the activation variable w in the expression for the potassium current is as follows: where ϕ was set to 0.04 with the steady-state activation amplitude as follows: and activation rate as follows: with V₃ = 2 mV and V₄ = 30 mV. The reversal potentials E_Ca, E_K, and E_L were set to 120, −84, and −60 mV, respectively. The maximal potassium (g_K) and leak (g_L) conductances were set to 8.0 and 2.0 mS/cm², respectively. For calcium, the maximal conductance (g_Ca) was set to 4.4 mS/cm² and I_stim was set at 100.0 μA/cm². Unless otherwise stated, the values for the various parameters for the type II excitability regime were equal to those given above and were taken from Rinzel and Ermentrout (1998). The synaptic coupling current was the same as in the case of the WB model, except that the synaptic decay time constant τ_syn was set to 10.0 ms.

Generation of the PRC.

PRC computation is illustrated in Figure 2A. The presynaptic model neuron is initialized at its spiking threshold when the coupling is turned on for a single cycle of the presynaptic neuron at different phases in the cycle of the postsynaptic model neuron. The action potential threshold is set to −14 mV in this study because the synaptic conductance trace began to noticeably increase above zero at that potential. Zero phase is defined as an upward crossing of the threshold potential. The PRCs were generated for each neuron by stimulating a single repetitively firing neuron by driving the synaptic conductance in that neuron with a single action potential from an identical neuron. This type of PRC is often called a spike response curve (Acker et al., 2003). The single action potential in the presynaptic model neuron triggers a change in synaptic conductance (Fig. 2A, bottom trace). This serves as a perturbation and is used to generate the PRC. In the open loop condition (i.e., when there is only a single unidirectional coupling) the phase at which a stimulus is received is φ = ts/P_i, where P_i is the intrinsic period and ts is the time between the last action potential in the model neuron for which the PRC is being generated and the action potential initiation in the presynaptic model neuron. The normalized change in the length of the cycle containing the perturbation and the one following it are called the first-order resetting, f₁(φ) = (P₁ − P_i)/P_i, and second-order resetting, f₂(φ) = (P₂ − P_i)/P_i, respectively. P₁ is the length of the cycle containing the perturbation and P₂ is the length of the subsequent cycle (Fig. 2A,B).

Figure 2.

Generation of the phase-resetting curve. A, The voltage waveform represents a regular spiking WB model neuron with intrinsic period P_i. The lower trace corresponds to the postsynaptic conductance resulting from a spike in the presynaptic model neuron. The change in the period P_i as a result of the perturbation received from the stimulus presynaptic action potential at a phase φ = ts/P_i is used to generate the phase-resetting curve. P₁ represents the period of the cycle in which the stimulus is received. The following cycle period is represented by P₂. B, Example of first-order (f₁(φ)) and second-order (f₂(φ)) phase-resetting curves. The parameter values for the WB model neuron considered are as follows: g_syn = 0.35 mS/cm², τ_syn = 1 ms, I_app = 2.0 μA/cm², and otherwise as in Materials and Methods.

Iterative pulse-coupled maps with no predetermined firing order.

In this section, we describe a type of map (Canavier et al., 1999) that cannot be easily analyzed because it does not assume a firing pattern, but instead polls the phases of each oscillator to determine which one will fire next. Given the current phase of each oscillator (neuron) and its period, one can easily determine which neuron(s) will fire next. Given the phase-resetting behavior of each neuron in response to the firing of every other neuron, one can then determine how the phases of each neuron are altered by each firing, so one can determine the future sequence of firing for all time, in the absence of noise. The determination of each firing event is considered one map iteration. The map may produce different firing patterns depending upon how it is initialized. We refer to this method as iterative rather than analytical to distinguish it from maps constructed assuming a fixed, predetermined firing order (see next section). The firing intervals can be specified by subtracting the starting phase from the ending phase, adding any resetting presumed to occur during the interval, and multiplying by the intrinsic period. With a single input per cycle, the interval between the firing of a spike (or burst) and the receipt of the next input is called the stimulus interval (ts). In the kth cycle it is defined as follows: where P_i is the intrinsic period of the neuron, φ[k] is the input phase, and f₂(φ[k − 1]) is the second-order reset due to an input in the previous cycle. The time between when the input arrives and the time that the neuron next fires is called the recovery interval. In the kth cycle it is defined as follows: where f₁(φ[k]) is the first-order reset due to an input in the current cycle. The first-order resetting is included in the recovery interval, whereas the second-order resetting is included in the stimulus interval. If multiple inputs arrive during a cycle, the first-order resetting is taken at the time that each input arrives. There are two ways to handle second-order resetting. One is to discard all second-order resetting except the one due to the most recent input. This preserves the assumption that the effect of one input dies out before the next is received and can be justified on the basis that in general only inputs arriving late in the cycle produce significant second-order resetting if the duration of the input is short compared with the cycle period. However, this approach destroys the symmetry required for synchrony to be observed (see Appendix 3) and was only used when summing the second-order resetting caused a problem emulating a splay mode.

The map representing an ideal pulse-coupled system was implemented using code written in C (see nonlinear.c, available at www.jneurosci.org as supplemental material). The map requires the first- and second-order PRCs, the intrinsic period of each neuron, and the initial values of the phase of each neuron as well as the stored value of the saved second-order resetting. The map updates the phases only at the times associated with each episode of neural firing. To determine how to update the phases, first it must be determined which neuron(s) will fire next. This is accomplished by finding the minimum of the expression [P_ij(1 − φ_j)], where P_ij is the intrinsic period of the jth neuron and φ_j is its phase. Once the next firing time was established, the phases of the firing and nonfiring neurons were updated differently. The phase of the firing neuron(s) was first set to zero, then any saved second-order resetting was subtracted from the phase of the firing neuron(s) (see below for the implications of negative phase) and the saved second-order resetting was cleared. Then any first-order resetting of the firing neuron due to any other neurons that were firing simultaneously was also subtracted. The phases of all the nonfiring neurons that the firing neuron(s) projects to were incremented by the appropriate amount based on the time elapsed to the next firing time, and then the first-order resetting due to the firing neurons was subtracted from the incremented phase. The saved second-order resetting was incremented to reflect these firings. This map can be quantified (Canavier et al. 1999) as a modified Winfree (2001) model: where ω is a normalized frequency of one cycle per period, φ_j,sn[k] is the phase at which the jth neuron receives the nth input in the current cycle, φ_j,n[k − 1] is the phase at which the jth neuron received the nth input in the previous cycle, f₁(φ[k]) is the first-order reset due to the input in the current cycle, and f₂(φ[k − 1]) is the second-order reset due to the input in the previous cycle. The first-order resetting is added at the time the input is received as indicated by the δ function. When the phase reaches one it is reset to zero, and as the δ function indicates, the second-order resetting is added at that time.

For sufficiently weak coupling strength, it is appropriate to simply add the resetting due to simultaneous inputs. On the other hand, in a network of strongly pulse-coupled oscillators, the phase resetting due to simultaneous inputs does not add linearly (Fig. 3C1). However, the simultaneously active synaptic conductances do add linearly. To account for this, the PRCs corresponding to each possible number (up to N − 1) of simultaneously active synapses were generated in the open loop configuration (Fig. 3A1) and indexed by the total conductance observed during N simultaneous inputs [f_i(φ, Ng_syn)]. When simultaneous inputs were detected in the iterated map of the closed loop network of four neurons (Fig. 3A2), the amount of resetting corresponding to the correct number of simultaneously active synapses was applied. The PRCs corresponding to the simultaneous inputs up to the maximum observed in the mode shown are illustrated in Figure 3, B1 and C1. The corresponding iterated maps are shown in Figure 3, B2 and C2, with the phase of each oscillator shown in a different line style. In between firing times, the phase is incremented steadily and forms a line segment with a slope of one. The firing time of a neuron is indicated by the vertical lines indicating that the phase has been reset to zero immediately after a phase of one is achieved. The first-order resetting for each neuron is indicated by a step change connecting two line segments with a slope of one indicating an instantaneous reset due to the firing of a presynaptic neuron. In each of these cases the iterative map incorporates both first- and second-order resets, although only the PRCs for first-order resetting are shown. The magnitude of the second-order resets is small and is represented by small negative excursions in Figure 3, B2 and C2. Figure 3C2 shows a mode with two clusters of two, where the red color actually indicates the sum of the second-order resetting due to the previous firing of the two neurons in the other cluster and the first-order resetting due to the other neuron in the same cluster.

Figure 3.

Applying PRCs to iterative prediction of network activity. A1, PRCs for individual oscillators were generated using the open loop topology. A2, Schematic of a four-neuron network of all-to-all identical and identically coupled pulse-coupled oscillators. B1, PRC corresponding to g_syn = 0.01 mS/cm² for a type I neuron with excitatory input. B2, Iterative map output for the splay mode. C1, PRC corresponding to g_syn = 0.08 and 0.16 mS/cm² for a type II neuron with inhibitory input. C2, Iterative map output for the antiphase cluster mode. PRCs are nonlinear with respect to conductance strength. Note that in a network of four oscillators, each oscillator can receive up to three simultaneous inputs. For strong coupling, the resetting does not add linearly as evidenced by the changing shape as the conductance strength is increased to represent a single input [g_syn = 0.01 (in B1) and 0.08 (in C1) mS/cm²], two simultaneous inputs (g_syn = 0.16 mS/cm²), and three simultaneous inputs (data not shown). The phases increment linearly in time except when an input causes an instantaneous reset. In each case both first- and second-order resets (indicated by small negative excursions in B2 and C2) are incorporated.

The examples illustrate that the number of inputs received simultaneously in such a network can vary. In the fully synchronous mode (data not shown) three inputs are received simultaneously. Figure 3B2 shows a splay mode (Fig. 1B) in which a single input is received at a time, corresponding to the PRC (Fig. 3B1) for type I neurons with excitatory synaptic connections with synaptic conductance strength of 0.01 mS/cm². A mode with two clusters in antiphase (Fig. 1C) is shown in Figure 3C2 corresponding to the PRC (Fig. 3C1) for type II neurons with inhibitory synaptic connections. In this mode, a single input is received from the other member of the cluster when the cluster fires, but two simultaneous inputs are received when the other cluster fires. Only two distinct traces are visible because the two traces in each cluster overlap. PRCs with synaptic conductance strength of 0.08 and 0.16 mS/cm² corresponding to a single input and two inputs respectively are shown in Figure 3C1. Note the large difference in shape indicating highly nonlinear summation of the resetting. The worst case scenario for delays is that a delay applied early in the cycle can result in a negative phase at the time a second closely spaced input is received. Since a negative phase does not actually correspond to a point on the limit cycle in this case (Oh and Matveev, 2009), the resetting assigned to the second pulse is arbitrary. If an input was received while the phase of a neuron was negative, we choose to assign the nearest value for resetting, i.e., the resetting at a phase of zero. This can be a source of error near synchrony (see Appendix 3).

In contrast to the iterative pulse-coupled map described in this section, the next several sections will describe analytical methods to determine the existence and stability of the fully synchronous mode, the splay mode, and for phase-locked synchronous clusters. The analytical methods for the determination of stability results require the construction of discrete maps. Discrete maps based on a fixed firing order may be iterated, but it is not necessary to do so because they are analytically tractable; their fixed points and the stability of those fixed points can be calculated using algebraic methods.

Stability criterion for synchrony in an all-to-all network.

In order for a synchronous mode to exist, all oscillators must have the same network period when they receive an input from all their presynaptic neurons at a phase of zero. Since this condition is automatically satisfied in an network of all-to-all identical, identically connected neurons, existence of the fully synchronous mode is guaranteed (see synchronymode.c, available at www.jneurosci.org as supplemental material). However, the only previous stability criterion (Goel and Ermentrout, 2002) for the fully synchronous mode in a pulse-coupled N neuron network assumes that the inputs are always staggered and remain in the same order as synchrony is approached, and assumes that the phase resetting at phases of zero and one is zero. As synchrony is approached in reality, the inputs begin to overlap and arrive simultaneously, and the phase resetting at zero is not generally zero. For strong coupling, the approximation of linear summation of resetting due to multiple inputs becomes poorer as N is increased. Therefore, the PRC must be parameterized by both the phase (φ) and conductance (g_syn), where f_i′(φ, g_syn) refers to the slope of the ith-order PRC at a phase of φ. A stability criterion for an all-to-all cluster of N neurons was developed by extending an existing proof for the case of two neurons (Dror et al., 1999; Oprisan et al., 2004) to the N neuron network. A perturbation is assumed in which only a single neuron is perturbed away from synchrony in a system (Fig. 4) that has been linearized about the fixed point of the map at synchrony. This leaves a subcluster of N − 1 neurons synchronized. The stability result depends critically on the sign and slope of the applicable PRCs and is derived in Appendix 1. In Figure 4, tr₁[k] and tr_N−1[k] represent the recovery intervals in kth cycle for the single neuron and the (N − 1) group of neurons, respectively. Similarly, ts₁[k] and ts_N−1[k] represent the stimulus intervals for the single neuron and the (N − 1) group of neurons, respectively. In general, increasing the conductance does not change the sign of the PRC at zero (corresponding to synchrony); however, it does in general increase the slope. If second-order resetting is ignored, a necessary, but not sufficient, condition for synchrony in the case of two clusters is that all the applicable PRC slopes be the range 0–2. Thus the largest possible slope will likely be the most destabilizing one, hence the choice of the largest possible (N − 1) subcluster for the most destabilizing perturbation, as it will in general produce the slope [f₁′(φ, (N − 1)g_syn)] with the greatest absolute value. However, depending upon how the PRC changes with increasing conductance, other possible perturbations into different size clusters may also need to be examined in a similar manner.

Figure 4.

Presumed firing pattern for stability proof. A perturbation of a single neuron away from the synchronous firing is assumed, and an iterated map based on the perturbation from synchrony on one cycle (ts).

The following criterion includes the effect of second-order resetting, which can be important at synchrony (Oprisan and Canavier, 2001; Oprisan et al., 2004) and produces the following characteristic quadratic equations whose roots (λ) must all have an absolute value less than 1 to guarantee stability: and The roots (λ) are also the eigenvalues of the linear system and matrix M described in Appendix 1. For completeness, there are two characteristic equations because there are two possible ways the firing order can be perturbed, with either the small or the large cluster leading. The leading cluster is assigned a phase of 0 approached from the right and the other a phase of 1 approached from the left. The sum of the first- and second-order resetting is presumed to have a continuous derivative at 0 and 1, but the first- and second-order resetting curves themselves do not have a continuous derivative at zero and at one (Goel and Ermentrout, 2002). In practice, if the coupling is approximately pulsatile in nature and dies out within one cycle of the oscillation, the first-order resetting disappears as a phase of one is approached, from the left, and second-order resetting dies out as a phase of zero is approached from the right. If the magnitude of the resetting is quite small in these regions, then it is reasonable to assume the slope is quite small as well: which reduces the expression above for the eigenvalues to λ = 1 − f₁′(0⁺, g_syn) − f₂′(1⁻, (N − 1)g_syn) and λ = 1 − f₁′(0⁺, (N − 1)g_syn) − f₂′(1⁻, g_syn). However, if the derivative of the sum of the first- and second-order resetting is indeed continuous at 0 and 1, coupled with the assumption that the four quantities given above tend to 0, then a single eigenvalue is obtained: 1 − f₁′(0⁺, g_syn) − f₁′(0⁺, (N − 1)g_syn).

Stability criterion for splay mode.

This proof was first presented by Goel and Ermentrout (2002) as a proof for the stability of synchrony. Here we rederive and reinterpret it as a proof for the stability of the splay mode. First we assume that the oscillators are firing in a sequential order such that no two neurons are synchronized (Fig. 5). As in the iterative map described above, for each neural firing time we add the first-order resetting due to the input (phase after the firing indicated by φ̂), then we advance the phases of the nonfiring neurons to the next firing time. The indices are switched as in Goel and Ermentrout (2002) on each cycle such that the phase of neuron to fire next (the largest phase) is always labeled φ_N−1 and the smallest is labeled φ₁. A variable is not wasted on the firing neuron because it is always 1 before the firing and 0 afterward. If the phase of the oscillators are such that 1 > φ_N−1 > φ_N−2 >…> φ₂> φ₁ > 0, then their phases in the (k + 1)th cycle can be written in terms of their phases in the kth cycle assuming that the phase after an input φ̂_i is a monotonically increasing function of the phase before the input φ_i so that the firing order remains constant. We then assume a fixed point of the mapping, which we indicate by dropping the indices for cycle number and using an asterisk instead: φ_i*. At such a fixed point, each neuron receives inputs at the same phases in each cycle, which corresponds to a periodic mode in the network comprised of the oscillators that generated the PRCs used for the map. We then linearize about this fixed point, ignoring second-order resetting, as explained in Appendix 2, and obtain a set of linear equations that describe how a perturbation vector Δ[k] from the fixed point evolves: where Δ[k] = [Δφ_N−1[k], Δφ_N−2[k],…, Δφ₂[k], Δφ₁[k]]^T and S is the following matrix:

Figure 5.

Mapping that results from assumed splay mode firing pattern: The dashed vertical lines indicate the firing times within an assumed sequential firing order in a network of N neurons. Successive firing intervals denoted here as ts_i (stimulus intervals), for i = 1 to N. For i = 1…N − 1, φ_i[k] and φ̂_i[k] indicate the phase of the nonfiring oscillator in the kth cycle before and after firing.

For a rotationally symmetric system in which each neuron receives exactly the same input, we obtain the single matrix shown above (see Appendix 2). If not, then the matrices corresponding to each firing time until the pattern repeats (N matrices) must be multiplied in the correct order (data not shown). The splay mode of both the PRC-based map and the corresponding periodic mode in the full network is locally stable if all the eigenvalues of S have an absolute value less than one so that the perturbation Δ[k] from the fixed point decays to zero. Note that the matrix S, and therefore the stability, depends only upon the slopes of the first-order PRCs at the locking points, indicated by f₁′(φ_i*).

Existence criterion for splay mode.

For simplicity, we restricted our analytical predictions of existence to symmetric splay modes in which the firing intervals between successive firing of neurons in the circuit were equal. To find the periodic modes after all transients have died out, we use the assumed fixed point of the mapping described in the previous section. The first input for each neuron in Figure 5 is received at φ₁*, the second at φ₂*, and the last at φ_N−1*. P_i is the intrinsic period. The steady firing intervals can be written as follows: For i from 2 to N − 1, ts_i[∞] = P_i {φ_i* − φ_i−1* + f₁(φ_i−1*)} and The first interval includes the second-order resetting from the previous cycle. Subsequent intervals calculate the elapsed time between inputs and account for the first-order resetting due to the input received at the start of the interval. If a value is assumed for φ_N−1*, then all other φ_i* required for equally spaced intervals can be iteratively calculated, by first calculating ts_N and finding the value of φ₁* that makes ts₁[∞] = ts_N[∞]. Then the values of the φ_i* that make ts_i[∞] equal to ts_N[∞] and ts₁[∞] are found. Only sets of φ_i* in which all φ_i* are less than 1 are considered. The equation for ts_N−1[∞] can be used to calculate a value of φ_N−1* that is compared with the originally assumed value. The difference between the assumed and calculated values of φ_N−1* can be plotted, and the zero crossings of this difference identify the unique set φ_i* of which satisfies all periodicity criteria that guarantee the existence of the mode. This procedure for identification of splay modes was implemented in C (see splaymode.c, available at www.jneurosci.org as supplemental material).

Prediction of clustering.

To understand how synchronized clusters and subclusters form, the synchrony within a cluster and the splay phase mode locking between clusters were considered separately (Fig. 6). First the stability of synchrony within a cluster is determined for an M-dimensional network. Next the existence of a between-cluster splay mode in an N/M-dimensional network is determined by modifying the existence criteria given above for the splay mode such that the first interval includes the first-order resetting due to the M − 1 other neurons within the cluster firing simultaneously at a phase of zero (see clustermode.c, available at www.jneurosci.org as supplemental material), producing the resetting with an additional argument that indicates the summed conductance of these inputs, f₁(0, (M − 1)g_syn): The stability of the between-cluster mode is determined by applying the stability criterion at the phases determined by the existence criterion. The prediction is based on the intersection of the predicted between- and within-cluster modes.

Figure 6.

Cluster mode prediction: within- and between-cluster interactions. Left, Schematic of four-neuron network with two clusters of two neurons each. Synchrony is observed within each cluster, and antiphase, which is an example of a splay mode, is observed between clusters. Within-cluster interactions are characterized by determining the stability of synchrony in an isolated two-neuron network (rectangle) with g_syn set to the value for a single synapse. Right, Between-cluster interactions are determined by collapsing the neurons within each rectangle to a single oscillator, and determining the existence and stability of a splay mode between two such oscillators using twice the conductance for a single synapse. This method generalizes to more clusters and to larger clusters.

Results

Iterative maps prove that PRC contains sufficient predictive information

The iterative prediction method throws away all information about the set of coupled nonlinear differential equations except for the PRCs and the intrinsic periods. For illustrative purposes we examined model networks of four identical, identically coupled neurons to test this hypothesis, although the iterative method can accommodate heterogeneity (Maran and Canavier, 2008). In a four-neuron network, some modes that can often be observed (Fig. 1) are the fully synchronous mode (one cluster of four neurons), splay firing mode, and two clusters of two neurons in antiphase. Not all networks exhibit each of these modes, and there are many other possible firing modes. If the iterative map successfully predicts the modes exhibited by the full system of differential equations under the assumption of ideal pulse-coupled oscillators, then we can be certain that the PRCs of the component neurons do indeed contain all the information required to predict the activity of the closed loop network. The output of the iterative map is the number identifying the neuron(s) that just fired and the elapsed time (firing interval) since the last neuron(s) fired. This output was used to quantify the comparison to the solution of the systems of differential equations.

The most complex network activity that we observed was for a four-neuron network of type I WB neurons coupled via inhibition, and the iterative method was quite successful at predicting closed loop network activity using PRCs generated in the open loop. Therefore, the iterative predictions for that system are described first and in detail. Figure 7 plots the predicted and observed firing intervals between successive firings of any neuron(s) in the network as the parameter space of the network was explored by increasing the synaptic conductance from very weak values to rather strong ones. Simultaneous firings are counted as a single event, hence a rotationally symmetric mode such as the fully synchronous mode (open blue diamond), two clusters of antiphase (open red circles), or the ideal splay mode in which all firing intervals are equal (not observed for this network) would contain a single firing interval. At weak coupling strengths, the interval for the two cluster antiphase mode is approximately half of that of the fully synchronous mode, and that of the ideal splay mode would be a quarter of the fully synchronous mode (see firing intervals delimited by vertical dashed lines in Fig. 1).

Figure 7.

Iteratively predicted versus observed firing intervals showing all observed modes. These results are for a four-neuron network with type I excitability and inhibitory synaptic connectivity. A, The firing intervals produced by the iterative map. B, Firing intervals produced by integrating the full system of differential equations. The modes labeled antiphase clusters or near antiphase clusters refer to two clusters in antiphase or near antiphase, whereas the synchronous and nearly synchronous modes refer to a single cluster.

If more than one type of symbol is observed at a given conductance, this indicates bistability. For example, at the weakest conductance values, the fully synchronous mode is bistable with the two clusters of two neurons in antiphase mode. Each mode has a basin of attraction, or set of initial conditions that leads to that mode. In the case of the differential equations, the initial conditions are the state space of the coupled system, whereas in the iterative map, only the phases of each oscillator (and the stored second-order resetting) need be initialized. Different initial conditions can lead to different oscillatory modes. An exhaustive search of initial conditions was not performed, but instead the search focused on initial conditions likely to lead to the three simple modes shown in Figure 1. The analytical predictions (see section on analytical methods) of stable and unstable modes can be used to compute the phase of each neuron at the time that a reference neuron reaches threshold, and these in turn can be used to initialize both the iterative map and the full system of differential equations exactly in the predicted mode. Therefore, every predicted mode was specifically tested for stability. In addition, these conditions were used to test for the presence of the mode at nearby parameter values where the mode was not predicted. Initializing the oscillators exactly at synchrony was not counted as a test of stability—the mode had to be robust to at least a very small perturbation to be considered stable. Initial conditions between the map and the differential equations could be compared by selecting initial conditions along the unperturbed limit cycle at a given phase.

The closed blue diamonds represent nearly synchronous modes in which the firing times of all four neurons are not exact (note that in Fig. 7B a clear bifurcation into slightly different alternating intervals is observed, denoted as near synchrony). The closed red circles represent two clusters of two neurons. The short closed red intervals indicate that the two neurons within each cluster are not exactly in phase with each other, which was denoted as near antiphase cluster. If the variation in the firing interval between any two neurons in the network is <10% of the network period, we label the mode near synchrony. In the near antiphase mode, the same criterion is applied within each cluster, and furthermore the difference between the midpoints of the firing intervals in each cluster must be half the network period. Finally, the closed green squares represent a near splay mode in which the four intervals in Figure 1B are not exactly equal to each other. The periodic modes were captured well by the iterative map, although the near synchronous mode persists at larger conductance values in the full system of differential equations (Fig. 7B) than in the iterative map (Fig. 7A) and there were some small discrepancies in the prediction of the generalized splay mode. The yellow stars represent intervals from complex periodic modes with a period of many cycles, or in some cases aperiodic modes. Even these modes were accurately captured by the iterative map. The correspondence of Figure 7A with Figure 7B argues that, in the network chosen for this example, the open loop PRCs contain all the information required to predict closed loop network activity.

In addition to type I neurons connected via inhibition, we also examined type I neurons connected via excitation and type II neurons connected by either inhibition or excitation. In these networks we only observed the simple modes described in Figure 1, hence we summarize the predictions of the iterative map for only these modes in Figure 8. The left-hand panels show the PRC for a representative conductance value. The center trace shows the iterative predictions, and the right-hand trace shows the observations of the full system of differential equations. Results for networks comprised of type I oscillators connected via reciprocal inhibition are shown in Figure 8A. Figure 8A1 shows a typical PRC for the type I inhibition that was used to generate the iterative map results shown in Figure 7A. Figure 8, A2 and A3, shows expanded views of the low end of the conductance range from Figure 7, A and B, respectively, with only the two simplest modes shown. These modes are those of clusters of two neurons each in antiphase with the other cluster but in exact synchrony within each other (open red circles), and of one cluster of four neurons in exact synchrony with each other. These modes are repeated here to facilitate comparison with the other networks examined. At low conductance values, these two modes are bistable, but first the fully synchronous mode then the two clusters in antiphase mode disappear as conductance is increased in both the iterated map (Fig. 8A2) and the full system of differential equations (Fig. 8A3).

Figure 8.

Iterative map prediction versus the observed network behavior in four-neuron all-to-all networks. Typical PRCs are shown in the leftmost column corresponding to a synaptic conductance strength of 0.10 mS/cm². Iterative map predictions are shown in the center column. The observed network dynamics resulting from integrating the full systems of differential equations are shown on the rightmost column. The open colored marker symbols indicate the firing intervals between spikes specific to a mode. A, WB model neurons (type I excitability) with inhibitory coupling. B, ML model neurons (type II excitability) with excitatory coupling. C, WB model neurons (type I excitability) with excitatory coupling. For g_syn from 0.09–0.17 mS/cm², all second-order resetting except that due to the most recent input was ignored (filled green squares indicate the different implementation of the iterated map). This modification was required (in this regime only) for correct predictions. D, ML model neurons (type II excitability) with inhibitory coupling.

In contrast to the bistable system in Figure 8A which also produced complex modes for some initial conditions (Fig. 7), the differential equations describing the remaining three networks that were examined produced only a single mode no matter what initial conditions were given, and this mode corresponded to the mode predicted by the iterated map in every case. For a four-neuron network comprised of oscillators with type II excitability and mutual excitation, only the fully synchronous mode was predicted (Fig. 8B2) and observed (Fig. 8B3). This mode retained stability until the conductance became so strong (3.75 mS/cm², data not shown) that the assumptions were violated because the oscillators became stuck at a depolarized fixed point. Both PRCs that supported the fully synchronous mode (Fig. 8A1,B1) had an initial positive slope (a phase of zero corresponds to synchrony). In the four-neuron network comprising oscillators with type I excitability coupled by mutual excitation, only a splay mode was predicted (Fig. 8C2) and observed (Fig. 8C3). Within a certain conductance regime i.e., for g_syn from 0.09 to 0.17 mS/cm², the splay mode was only accurately predicted by the iterative map if only second-order resetting due to the last input in each cycle was considered. In this regime, the network period ranged between 6.34 ms to 2.2 ms, considerably shorter than the intrinsic period of ∼30 ms, and the assumption that each neuron return to the limit cycle between inputs is therefore severely challenged, especially since the times between inputs are reduced to within 1.0–2.0 ms while the synaptic time constant is in the same range, 1 ms. Here, the PRC saturates at large conductance values (Prinz et al., 2003; Oprisan et al., 2004) and the splay mode is never observed to lose stability.

Finally, in the four-neuron network comprised of oscillators with type II excitability with mutual inhibition, the only mode predicted (Fig. 8D2) and observed (Fig. 8D3) was a mode composed of two clusters of two neurons that were synchronized with the other neuron in its cluster but in antiphase with the neurons in the other cluster. Neither of the PRCs with an initial negative slope (Fig. 8C1,D1) could support the fully synchronous mode, and yet in the final case shown, the two neurons within each cluster were synchronized. The result is surprising because a single isolated cluster of two cannot synchronize itself (data not shown), hence the interaction between clusters must somehow stabilize within-cluster synchronization. The good correspondence between the predicted and actual modes in all cases suggests that for these model networks, all information required to predict network activity is indeed contained within the PRCs. The determinant of whether the coupling can be approximated as pulsatile is the duration of the synaptic input compared with the network period. The network period for type I neurons ranged from 2 to 31 ms compared with a synaptic time constant of 1 ms. The network period for type II neurons ranged from 82 to 86 ms compared with a synaptic time constant of 10 ms.

Analytical criteria also prove that PRC contains sufficient predictive information

Having shown using iterative maps that the PRCs in different cases indeed contain sufficient information to predict network activity, we next used analytical existence and stability criteria to gain insight into why certain modes were exhibited by individual networks at the given parameter settings whereas others were not. A stability criterion for the fully synchronous mode, as well as existence and stability criteria for a symmetric splay mode, is given in Materials and Methods. Each of the four networks was analyzed throughout the parameter regimes given in Figure 8 to determine the parameter settings at which the synchronous and splay modes would be predicted to exist and be stable. Then, as described in Materials and Methods, the existence of two clusters of two was predicted by the intersection of parameter values that support synchrony in a cluster of two neurons as well as the antiphase mode between two such clusters connected via a double synapse, with the slight modification as described in Materials and Methods that the stimulus intervals contain the first-order resetting due to the other neuron in the cluster.

Figure 9 shows the results of the analytical predictions side by side with the observations. For the first case of a four-neuron network with type I excitability and mutual inhibition, the analytical predictions (Fig. 9A2) perform somewhat better than the iterative one (Fig. 8A2) in that the point at which the fully synchronous mode loses stability is better predicted by the analytical method (compare with Fig. 9A3) for reasons described in Appendix 3. The two methods give the same performance for the two clusters of antiphase mode, which is explained in more detail below. This is an informative example because the two modes shown continue to exist but lose stability, and the analytical method correctly predicts the point at which stability is lost. The next two examples also show correct performance but are not as illuminating because, as described above, a loss of stability is not observed. As the synaptic conductance strength is increased the synchronous mode and the splay mode, corresponding to type II neurons and type I neurons respectively; remain stable. This is because the PRCs saturate in both cases. However, the existence criteria correctly predicted the firing intervals in every case. The fully synchronous mode is correctly predicted both analytically (Fig. 9B2) and iteratively (Fig. 8B2) in every case (compare with Fig. 9B3) for the network with type II excitability and mutual excitation. Similarly, the splay mode is correctly predicted both analytically (Fig. 9C2) and iteratively (Fig. 8C2) in every case (compare with Fig. 9C3) for the network with type I excitability and mutual excitation.

Figure 9.

Analytical prediction versus the observed network behavior in four-neuron networks. Typical PRCs are shown in the leftmost column corresponding to a synaptic conductance strength of 0.10 mS/cm². Analytical predictions are shown in the center column. The observed modes are shown once again on the rightmost column. The filled circles on the PRCs indicate the phase at which inputs are received by each neuron of the network. The firing intervals between spikes specific to each mode are indicated just as in Figure 8. The black open circles in D2 are used to indicate that the cluster mode is not in agreement with the observed mode (D3) with regard to stability. This is analyzed further in Figure 10.

The phase at which the inputs are received within a cycle are denoted by filled circles on the PRCs. The analytical approach (Fig. 9C2) predicts the observed mode (Fig. 9C3) quite well. The slopes of the PRC at the indicated points on the PRC are −3.39, 0.741, and 0.736, respectively. A key observation that is relevant to predicting clusters, which can only exist if the splay mode between clusters exists and is stable, is that the inputs phase in stable splay modes (dark filled circles) often, but not always, fall in regions of positive slope on the PRCs (Fig. 9C1) that support the splay mode. It is worth noting that in no case was a mode predicted but not observed, for example the networks in Figure 9, A, B, and D, in all cases predicted the splay mode to be unstable, and it was not observed.

One failure of the analytical (Fig. 9D2), but not the iterative (Fig. 8D2), method is observed in the case of two clusters in antiphase (Fig. 9D3) in the network with type II excitability and mutual inhibition. The black open circles in Figure 9D2 are used to indicate that this cluster mode is in agreement with the observed mode (Fig. 9D3) with respect to existence but not stability. Thus, although the analytical approach that treats within- and between-cluster interactions separately does not provide a complete understanding of cluster formation, it nonetheless provides significant boundaries on when cluster formation might be expected (see next section).

Further analysis of clustering in a four-neuron network

A positive PRC slope (Fig. 9A1,D1) at a phase of P_N/2P_i (where P_N is the network period and P_i is the intrinsic period of the component oscillators) determines the stability of the antiphase cluster mode. Clustering is examined by treating a synchronous subcluster of M neurons as a single oscillator with the same period and PRC as a single component neuron, then finding the splay modes for a network of N/M oscillators. In a network of N neurons, within-cluster interactions of M neurons are analyzed by studying the synchronous mode for M neurons, and between-cluster interactions of N/M clusters are analyzed by studying the splay mode for a network of N/M neurons. Figure 10 illustrates the cluster prediction method. The top traces show bars that correspond to the parameter regions in which within-cluster synchrony is predicted to be stable (blue bar), between-cluster splay is predicted to exist and be stable (green bar), and in which the cluster mode is actually observed (red bar). The lower traces show the largest eigenvalue for the stability prediction for within-cluster synchrony (blue diamonds) and between-cluster splay (green squares). Note that since we plot the absolute value of the eigenvalues, the blue diamonds in Figure 10B1 corresponding to the within-cluster synchronous mode reverse direction when the eigenvalue passes through zero.

Figure 10.

A1, A2, Qualitative cluster mode prediction. When the cluster mode is observed, predictions are made analytically using within- and between-cluster interactions. The vertical dotted line indicates the intersection region where the within- and between-cluster interactions agree with the observed cluster mode (red line). For the four-neuron network, the blue line indicates stable synchronous mode in a two-neuron network corresponding to the within-cluster interactions. The green line indicates the stability of splay mode (two clusters of two neurons) corresponding to the between-cluster interactions. B1, B2, The absolute value of the largest eigenvalue (|λ_max|) for the synchronous mode in a two-neuron network (blue diamonds) and for the two clusters of two-neuron mode (green squares) are plotted against the synaptic conductance strength). A value less than 1.0 indicates stable mode and a value greater than 1.0 indicates unstable mode.

First we analyze the cluster mode exhibited in networks of type I neurons coupled by mutual inhibition shown in Figures 8, A2 and A3, and 9, A2 and A3. The intersection of the parameter values encompassed by both the green and blue bars is indicated by the vertical dotted line and predicts the range of the observed values (red bar) reasonably well. Fig. 10B1 shows that the absolute value of the largest eigenvalue (|λ_max|) corresponding to both the within-cluster and the between cluster is <1.0 and is to the left of the vertical dotted line (for g_syn values from 0.0 mS/cm² to 0.07 mS/cm²). Therefore the cluster mode is predicted to be stable by the analytical approach until 0.07 mS/cm². Beyond this the within-cluster synchrony loses stability, and as a result the cluster mode is predicted to lose stability. Note that the red bar extends slightly to the right of the dashed line, indicating that the cluster mode persists for slightly stronger values than predicted, indicating that the small maximal eigenvalues for the between-cluster mode can somewhat compensate for the destabilizing influence of within-cluster eigenvalues with absolute value greater than one. This effect is much more striking in the clusters observed in networks of type II neurons coupled by mutual inhibition corresponding to Figure 8, D2 and D3. Figure 10A2 indicates that the analytical approach does not predict the stability of the observed cluster mode since the within-cluster synchrony for the two-neuron network is unstable for all values of g_syn. Figure 10A2 has no blue bar corresponding to the parameter region that supports stable within-cluster synchrony, hence no vertical line indicating the intersection of the blue and green bars. This is because Figure 10B2 shows that the absolute value of the largest eigenvalue (|λ_max|) corresponding to the within-cluster interaction is >1.0 for g_syn values from 0.0 mS/cm² onwards. The prediction that a two-neuron cluster is unstable is correct as tested in the full system of differential equations (data not shown). As a result, the cluster mode in a four-neuron network is predicted to be unstable by the analytical approach. However, the observed cluster mode in a four-neuron network is stable throughout the given conductance range. This strongly suggests that the between-cluster interaction contributes to stabilizing the cluster mode. A small between-cluster eigenvalue again seems to compensate for a within-cluster eigenvalue greater than one. In contrast to the idea that neurons cannot synchronize if they do not possess the correct PRC slope at zero, the splay mode between clusters was sufficient to synchronize subclusters that were not capable of synchronizing themselves. This is a striking and somewhat counterintuitive finding, with implications for synchronization under both physiological and pathological conditions. Despite the obvious over-simplification, the method was successfully applied in networks of four and 12 neurons.

Clustering in a 12-neuron network

The above results are not specific to any particular network size. As an illustrative example we use a 12-neuron network that is a scaled up version of the four-neuron network with inhibitory inputs to the type I neurons illustrated in Figures 7, 8A, and 9A. In this network, we predicted and observed multistability between the fully synchronous mode (data not shown), two clusters of six (Fig. 11A), and three clusters of four (Fig. 11B). Modes with more than three clusters (Fig. 11C) were found but were always unstable. In Figure 11, the blue bar shows the parameter range that the analytical method predicts will support the within-cluster synchrony, the green bar predicts stable between-cluster splay mode, and the black bar predicts unstable between-cluster splay mode. The intersection of the blue and green bars (demarcated by vertical dashed line) is the analytical prediction of a clustered mode. The red bar shows the parameter regime in which a cluster solution was observed in the full system of differential equations and the orange bar shows the parameter regime in which the iterated map produced a cluster mode. The phases at which each cluster receives input from the others are shown (open squares) on the PRC on the right side of each panel for the case corresponding to g_syn = 0.01 mS/cm².

Figure 11.

Cluster predictions for three networks of 12 type I neurons with inhibitory inputs. The model neurons were type I WB neuron coupled by inhibition. PRCs are shown for g_syn = 0.01 mS/cm². The brown squares on each PRC indicate the phases at which each cluster receives inputs from the other clusters (one, two, and three, respectively) in the network. The observed results are indicated by red horizontal lines and the results predicted by the iterative method are indicated by orange horizontal lines. A, Two clusters of six neurons. A1, Two clusters of size six were observed to be stable up to g_syn = 0.02 mS/cm². Within-cluster synchrony is stable until g_syn = 0.02 mS/cm². The splay mode captures the between-cluster interaction and remains stable beyond g_syn = 0.10 mS/cm². The dashed vertical line indicates that the two clusters of six neurons are predicted to be stable until g_syn = 0.02 mS/cm² by the analytical approach. A2, PRC. B, Three clusters of four neurons. B1, Three clusters of size four were observed to be stable up to g_syn = 0.04 mS/cm². Within-cluster synchronous mode is stable until g_syn = 0.03 mS/cm². The splay mode captures the between-cluster interaction and is stable for a range of conductance values. The vertical line indicated that the three clusters of four neurons are predicted to be stable until g_syn = 0.03 mS/cm² by the analytical approach. B2, PRC. C, Four clusters of size three. C1, Four clusters of size three were not observed. Within-cluster synchronous mode is stable until g_syn = 0.04 mS/cm². The between-cluster splay mode exists but is not stable (indicated by black bar) for the given range of conductance values. Four clusters of three neurons are predicted to be unstable by the analytical approach. C2, PRC.

Figure 11A1 shows that the two clusters of six mode is analytically predicted to be stable up to g_syn = 0.02 mS/cm² (vertical dashed line), which agrees with observation (red bar). The iterated map (orange bar) correctly predicted the existence of the mode during part of the parameter region but was incorrect on the borderline of the observed mode. Similar results were obtained for three clusters of four (Fig. 11B1). The loss of stability is due to the loss of within-cluster synchrony as conductance is increased. The between-cluster splay mode remains stable at very large conductances, presumably due to saturation of inhibition such that the IPSP essentially clamps the membrane potential at the synaptic reversal potential. For four clusters of three, both the iterated map and the analytical method correctly predicted that this mode would not be observed in the full system of differential equations. The loss of stability of the between-cluster splay mode is responsible for limiting the number of clusters that can be observed.

One of the relevant slopes for within-cluster synchrony is always f₁′(0, g_syn). The other is f₁′(0, 5g_syn) for a cluster of six, f₁′(0, 3g_syn) for a cluster of four, and f₁′(0, 2g_syn) for a cluster of three. Within-cluster synchrony loses stability at g_syn = g_max/(M − 1), where g_max is the total value for the conductance of the M − 1 cluster when stability is lost. Since the slope of the PRC decreases as the conductance decreases, the stability of the six neuron cluster implies the stability of smaller clusters, just as the stability of the 12-neuron cluster (data not shown) implies the stability of all smaller clusters. On the other hand, the stability criterion for the splay mode is complex, but the slope corresponding to the last input received within the cycle appears more frequently and may be more critical. Empirically, we note that stability can be lost when this last input falls in a more steeply negative region of the PRC as in Figure 11C2 compared with Figure 11, A2 and B2. Increases in the largest eigenvalue indicate a destabilizing trend. The largest eigenvalues corresponding to Figure 11, A2, B2, and C2, follow such a destabilizing trend, and are 0.834, 0.973, and 1.009, respectively, for the case corresponding to g_syn = 0.01 mS/cm². Another empirical observation (data not shown) is that it is destabilizing for multiple input phases to fall in the negative slope region.

Discussion

Applications to real neural networks

We have derived general results for generic phase-locked modes based on the slopes of the PRCs at the respective locking points. Although these results were derived for idealized networks of identical and identically connected neurons, our intent is to provide clues about network synchronization behavior in more general networks as described in the next three sections. The network components are not required to be single neurons, but can be generalized to half-center oscillators that rely on inhibitory rebound or to more complex circuits such as gamma modules consisting of multiple neurons (Tort et al., 2007). The analysis can also be generalized to bursting neurons (Luo et al., 2004; Oprisan et al., 2004; Sieling et al., 2009) except that the burst is used as the perturbation instead of a spike as in this study. These methods can be generalized to circuits containing noise (Mancilla et al., 2007); moreover, if the number and stability of the fixed points identified by the stability criteria are robust to slight shifts of the PRC, then the identified modes should be robust to PRC measurement error and other sources of noise, as confirmed both in iterated maps with added noise and in hybrid circuits constructed with the dynamic clamp (Sieling et al., 2009).

Further generalization can be realized by considering heterogeneous networks as perturbed versions of homogeneous, identical, identically connected networks. The modes predicted for all of the homogeneous circuits used herein are robust to the introduction of some heterogeneity (data not shown). Thus PRCs that would lead to exact phase locking in a homogeneous network are predicted to lead to locking points in mildly heterogeneous networks, with phase lags clustered about the exact locking in the idealized circuit. The iterated map requires one PRC per input per oscillator, and there is no requirement that they be identical, so the circuits it can emulate are limited only by the availability of PRCs (Maran and Canavier, 2008).

Implications for synchrony

The analytical approach (Fig. 9A2,B2) builds on the observation that positive slopes of the PRC at a phase of zero corresponding to spike initiation promote synchrony. This implies that at sufficiently weak conductance strengths, type I neurons tend to synchronize with inhibition but not excitation, whereas type II neurons synchronize with excitation but not inhibition. For weak coupling, a positive slope guarantees synchrony (Ermentrout and Kopell, 1990, 1991; Ermentrout, 2002). Previous work on stability in all-to-all networks of pulsatile coupled oscillators (Goel and Ermentrout, 2002) did not address second-order resetting, but only first-order resetting at phase just after spike initiation (0⁺) and just before spike initiation (1⁻). For pulsatile coupling, first-order resetting at 1⁻ is negligible, but the second-order resetting at 1⁻ is equal to the first-order resetting at 0⁺. This allows the stability criterion at synchrony to collapse to λ = 1 − f₁′(0⁺, g_syn) − f₁′(0⁺, (N − 1)g_syn).

All possible perturbations (clusters of N − j and j) should technically be considered, but we found that the criterion above is generally the most severe. The simple form of the eigenvalue λ provides valuable intuition since it is stated in terms of the relevant PRC slopes. For example, it is clear that if both slopes are negative, the eigenvalue will be greater than one and therefore full synchrony in an N neuron cluster will be unstable. Also, if both slopes are positive but their sum is greater than two, the eigenvalue will be less than negative one, and synchrony will again be unstable. Changing the magnitude of the conductance at zero usually does not change the sign of the initial slope. Therefore, if the initial slopes are positive, any manipulation that increases their magnitude sufficiently will destabilize synchrony. These manipulations include increasing the size of the cluster or increasing the individual synaptic conductances, consistent with previous modeling studies of a pulse-coupled network of bursting Hindmarsh-Rose neurons (Belykh et al., 2005) in which stability of the fully synchronous state depended only on the total input each neuron received. For heterogeneous circuits, the near synchronous state will tend to be destabilized when the absolute value of the eigenvalue, calculated using the appropriate PRCs, exceeds one. The appropriate PRCs to use for each neuron represent (1) the effect of that neuron on the rest of the cluster, which may be approximated by its effect on a representative neuron within the cluster, and (2) the simultaneous effects of the all other neurons in the cluster on a specific neuron. The dynamic clamp protocol (Sharp et al., 1993a,b) may be helpful in approximating the latter.

Implications for clustering

The existence and stability of the between-cluster splay mode is necessary for clusters to form. If within-cluster synchrony is also stable, then clusters will be observed, but if it is not, it is still possible that clusters could form (Fig. 10). Previously, Jeong and Gutkin (2007) hypothesized that bistability in a two-neuron system between synchrony and antiphase implies the ability to form two synchronous clusters in a larger network. Another modeling study (Ernst et al., 1995) showed the formation of stable clusters in pulse-coupled systems enabled by bistability in the two-neuron network, as well as clustering where synchrony in the two-neuron network was unstable. Only the latter clusters were susceptible to temporary disruption by noise. We extend these previous results by setting specific criteria in terms of the PRCs for bistability. For example, the PRC for type I inhibition (Fig. 8A1) has an initial region of positive slope that extends beyond a phase of 0.5, so at least for weak values of conductance, the fully synchronous mode is bistable with the antiphase mode in a two-neuron network (data not shown), which does lead to clustering in the map (Fig. 8A2), the analytical prediction (Fig. 9A2) and the full system of differential equations (Fig. 8A3). We also show that making the connections within clusters stronger by either increasing the conductance of the individual synapses or by increasing the size of the cluster, can destabilize clustering just as it does the fully synchronous state.

Stable phase-clustered states have been mostly found in networks with reciprocal inhibition, except in an all-to-all excitatory network in which within-cluster connections were weakened to reduce the network to a bidirectional ring with all adjacent neurons in antiphase, producing synchronous clusters comprised of nonadjacent neurons (Li et al., 2003). We were able to identify one all-to-all homogeneous excitatory network that supported bistability between antiphase and synchrony in a network of four homogeneous neurons, and also supported clustering (supplemental Fig. S1, available at www.jneurosci.org as supplemental material). Therefore the PRC, and not the sign of the reciprocal coupling, determines whether synchronous clusters will form.

In large network simulations (Pervouchine et al., 2006), with many entorhinal stellates and fast-spiking cells, the stellates break up into two clusters, each firing at a theta rhythm, with a population rhythm twice as fast in the beta range. The two clusters are invoked to explain a beta peak seen by Cunningham et al. (2004) in slices of entorhinal cortex in which the main cells involved, the fast-spiking and stellate cells, are thought to be firing at theta frequency. The between-cluster analysis methods developed here show that in general it is progressively more difficult to find PRCs that support between-cluster modes with N/M larger than 2 (Fig. 11), which likely explains why at most two clusters are observed.

Summary of major results

Two systematic PRC-based approaches have been proposed: in one the interactions are presumed to be weak and in the other they are presumed to be pulsatile (Ermentrout and Chow, 2002). We have chosen to pursue the pulsatile approach because neural synchrony is frequently achieved rapidly, even within a single cycle (Singer, 1999), implying strong coupling in which the relative phases of the oscillators change quickly with respect to their absolute phase. For strong coupling, phase resetting does not scale linearly as synaptic conductance is increased (Fig. 3A3). Our key extension of the pulsatile coupling method was to use a lookup table for simultaneously received inputs that accounts for nonlinear scaling of PRCs in the strong coupling regime. This allows larger networks that receive multiple inputs per cycle to be analyzed directly rather than using two-neuron networks to quantitatively assess larger network dynamics (Skinner et al., 2005, Netoff et al., 2005a) or reducing the larger network to a two-neuron network (Pervouchine et al., 2006). The methods that we propose could be broadly applicable because they can be used to gain insight into networks composed of biological neurons, for which the PRCs can be determined directly. Only strong coupling methods can predict the loss of stability as the coupling becomes strong, or as the net coupling received by a given neuron due to increasing network size becomes larger. The major limitation is the assumption of pulsatile coupling: if the inputs are so close together the trajectory has not returned close to the original cycle before each input is received then the assumption is violated.

In sum, we predict that only certain PRC shapes can support a single synchronous cluster, that only certain PRC shapes can support a splay mode, that a subset of those shapes which supports both modes can robustly support multiple clusters, and that as the coupling strength is increased, clusters tend to break up. These are the types of very general predictions that may be testable even in realistic brain networks for which we will have access only to partial information.

Appendix 1: derivation of the stability criterion for the fully synchronous mode

From Figure 4, the periodicity criteria are as follows: Substituting for the recovery intervals (tr) and stimulus intervals (ts) in the above equations, we have the following: Rewriting the above equations, we have the following: In the above two equations we assume a perturbation Δφ_n[k] about the fixed point φ_n* so that φ[k] = φ_n* + Δφ_n[k], n = 1 or N − 1. After linearizing the PRCs such that f(φ) = f(φ*) + f′(φ*) Δφ and canceling the steady-state terms from both sides of the equations, we obtain the following: The two equations above constitute a discrete linear system: where Δ[k] = [Δφ₁[k], Δφ_N−1[k]]^T and M is a 2-by-2 matrix.

The characteristic polynomial corresponding to the matrix M is as follows: For synchrony, the above expression has to be evaluated twice, once at φ₁* = 0⁺ and φ_N−1* = 1⁻, and again at φ₁* = 1⁻ and φ_N−1* = 0⁺. This accounts for two possible perturbations from synchrony, one in which neuron group N − 1 leads and one in which neuron 1 leads (Fig. 4). Often, f₁′(1⁻, (ig_syn)) and f₂′(0⁺, (ig_syn)) are near 0, for i = 1 and N − 1, because the PRC is flat in that region, thus the expression reduces to a single eigenvalue for each evaluation: f₁′(0⁺, ig_syn) = f₂′ (1⁻, ig_syn) for i = 1 and N − 1 because they measure the same quantity, namely, the resetting in the next cycle when the perturbation occurs at the time of action potential initiation. Thus, we obtain the following eigenvalue twice:

Appendix 2: derivation of the stability criterion for the splay mode

Suppose that the oscillators are firing in a sequential order. Let the phases of the oscillators immediately before oscillator N fires be the following: 1 > φ_N−1[k] > φ_N−2[k] >…> φ₂[k]> φ₁[k] > 0. The phase of the firing oscillator φ_N is taken to be 1 in the kth cycle and the oscillators are reindexed on every cycle so that only N − 1 variables are required. Considering only first-order resetting (f₁(φ)), which is not indexed by g_syn here because only one input is received at a time, and ignoring higher-order resetting, we obtain expressions for the phase immediately after a perturbation (φ̂) in terms of the phase immediately before the perturbation (φ): φ̂_N = 0 and φ̂_i = φ_i − f₁(φ_i) for i = 1,…,N − 1. Assuming that the phase transition function is monotonically increasing (a key assumption), we have that φ̂_N−1 > φ̂_N−2 >…> φ̂₁. The next step is to advance all of the phases until just before the next neuron fires, but keeping the numbering such that the neuron that just fired as oscillator N in cycle k is now oscillator 1 in cycle k + 1, oscillator 1 is now 2, and so on (Fig. 5). The phase φ₁[k + 1] in the stimulus interval ts₁ is equal to the normalized recovery interval 1 − φ̂_N−1[k], so we can calculate it as follows: Each phase can then be updated to the next firing time by adding φ₁[k + 1] to the phase after the previous input and increasing the index by one as in Figure 5. In the above equations we can assume a perturbation Δφ_i[k] about a fixed point φ_i*φ_i[k] = φ_i* + Δφ_i[k] for i = 1,…N − 1; where φ_i* is the steady-state fixed point. After linearizing the PRCs such that f(φ) = f(φ*) + f′(φ*)Δφ and canceling the steady-state terms from both sides of the equations, we obtain the following linear system: .............. These equations can be written as a single matrix equation as follows: where Δ[k] = [Δφ_N−1[k], Δφ_N−2[k],…, Δφ₂[k], Δφ₁[k]]^T and S is the following matrix:

The eigenvalues of the above matrix determine the stability of the splay mode. Since this matrix describes the linearization of a discrete map at the fixed point of the mapping that corresponds to a splay firing pattern, the eigenvalues of the matrix determine not only the stability of the linear system, but also whether a perturbation of the fixed point would grow or dissipate, with implications for the stability of the corresponding firing pattern in both the map and the full system of differential equations. The mode is predicted to be stable only if the absolute value of all eigenvalues is <1.

Note that if the intervals between firings are not identical and/or the neurons are not identical, then it is not sufficient to consider what happens when a single neuron fires, but instead N sequential firings must be considered in which N neurons receive N − 1 inputs each, so N square matrices of dimension N − 1 containing N − 1 PRC slopes in the format shown above, for a total of N(N − 1) distinct slopes, must be multiplied by each other to obtain the correct linear system to determine stability in this case. Also, the intrinsic periods cannot be ignored as they were in the above derivation. The resultant matrix will still have N − 1 eigenvalues.

Appendix 3: technical limitations of the pulse-coupled maps

Causality

An input cannot advance the next spike to before the input is received. In practice, this implies that the recovery interval between when an input is received and when the next spike is generated must be non-negative. If inputs were allowed to sum linearly in a nonlinear regime, causality violations of this type can occur, but the lookup table method prevents this type of causality violation. The stimulus interval between when a neuron fires and when it receives its next input must also be non-negative. In the case of no simultaneous inputs, this interval is equal to the phase of the next input φ_i plus the second-order resetting f₂(φ_j) from the previous input, hence −f₂(φ_j) must be less than φ_i.

Strengths and weaknesses of the iterated map with no presumed firing order

The strength of the iterative map method described in Materials and Methods is that a single map can reproduce all firing orders, whereas by definition a map based on a presumed firing order cannot (Pervouchine et al., 2006). Maps similar to the iterative map have been used in the past to study network dynamics (Canavier et al., 1997, 1999; Netoff et al., 2005b; Maran and Canavier, 2008; Oh and Matveev, 2009; Sieling et al., 2009). The iterated map method could predict two clusters of two in antiphase (Fig. 8D2) in a case in which analytical methods failed (Fig. 9D2), and it easily handles multistability and aperiodic modes (Fig. 7). However, the fixed points of this map cannot easily be determined a priori due to the nonlinear step of polling which neuron(s) will fire next. Therefore, stability criteria based on the slopes of the PRC cannot be derived from this map, which limits the intuitive understanding of synchronization and phase locking that can be derived from this method. There are also some subtle problems with simultaneous inputs that caused a poorer prediction of full synchronization (Fig. 8A2) compared with the analytical method (Fig. 9A2) in one case.

The subtle problems at synchronization include the possibility of a negative phase immediately after spiking when second-order resetting is a delay (positive resetting in our convention). In Figure 3C2, for example, a small dip into negative phases can be observed right after the phase of each neuron is reset from one to zero to indicate that the neuron has fired. Since a negative phase does not actually correspond to a point on the limit cycle in this case (Oh and Matveev, 2009), the resetting assigned to the second pulse is arbitrary. As described in Materials and Methods, if an input was received while the phase of a neuron was negative, we choose to assign the nearest value for resetting, i.e., the resetting at a phase of zero. Thus, this can introduce error near synchrony.

A similar problem occurs in the convergence to synchrony of two inputs to a given neuron, as in the two clusters of two case. When two inputs are received in succession, all of the resetting due to the first input is assumed to be complete by the time the second input is received, but if the two inputs arrive close together this may not be a good approximation. If the pulses are simultaneous, the resetting due to two pulses is f(φ,2 g_syn), whereas if the two pulses are separated by a normalized delay τ the total resetting is f(φ) + f(φ − f(φ) + τ) (Oprisan and Canavier, 2003). In the limit as τ approaches zero, the resetting due to two pulses approaches f(φ) + f(φ − f(φ)), causing a discontinuity at φ for simultaneous pulses unless f(φ, 2g_syn) = f(φ, g_syn) + f(φ − f(φ, g_syn), g_syn), which is not necessarily true in general, another reason that the map may not always converge to synchrony as desired.

The assumption that the resetting be complete by the time the next input is received means that the system is memoryless. However, if the interval between two inputs contains a threshold event such as an action potential, then the interval is comprised of two parts: the recovery interval between when the input is received and the action potential, and the stimulus interval between the action potential and the next input. In this case, the first-order resetting due to the first input can be taken in the first interval and the second order in the second interval without violating the assumption of a return to the limit cycle between inputs. The second-order resetting due to inputs that do not immediately precede an action potential are dropped under a strict memoryless assumption. However, in a system with more than two neurons, the pulse-coupled map based on a strictly memoryless assumption cannot produce the fully synchronous mode or synchronous subclusters because ignoring all second-order resetting except that due to the last input received in a cycle very slightly disrupts the symmetry required to achieve exact synchrony. This is only significant because the theoretical proofs presented in the paper only deal with exact synchrony, hence we would like for the pulse-coupled map to reproduce these modes. If we include all second-order resetting, the memoryless assumption is somewhat relaxed, as it must be in the full system of differential equations and easily achieve exact synchrony if the stability criteria are fulfilled. On the other hand, including all second-order resetting can disrupt ability of the pulse-coupled map to reproduce the splay mode in Figure 8C2 when the intervals are short compared with those obtained by including second-order resetting from each input. For the parameter regime g_syn = 0.09–0.17 mS/cm², the pulse-coupled map was reprogrammed to ignore the second-order contribution of all but the most recent input to successfully reproduce the splay mode in Figure 8C3. Since inputs occur with small separation in the splay mode, this is a severe test of the condition that all second-order resetting be complete by the time the next input is received. A final weakness of the iterated map is the tradeoff between the strict memoryless assumption and the symmetry requirement of synchrony. Mismatches between the performance of the full system of differential equations can be observed in either case when the assumptions are pushed too hard, but very good agreement is observed in a wide variety of cases over large parameter regimes.

Strengths and weaknesses of the analysis of maps with a presumed firing order

The strength of the analytical prediction method is that it allows us to write an explicit criterion for stability of the presumed firing pattern in terms of the slopes of the PRC that can be applied to the fixed points of the map, which are easily identified. The major intuition gained for the synchronous mode that is derived is only possible using the analytical method. As noted above, the analytical method (Fig. 9A2) outperformed the iterative map method in predicting when synchrony would break down (Fig. 8A2) in one case. Although the analytical method also allows us to develop intuition about how synchronous clusters can form, the iterative map method (Fig. 8A2,D2) was better able to identify which parameter sets support this type of clustering than the analytical method (Fig. 9A2,D2). The iterative map cannot provide insight on why this type of clustering occurs, though, some more analytical work is required to determine under what conditions between-cluster interactions can produce within-cluster synchrony in a cluster that cannot synchronize on its own.

Footnotes

This work was supported by National Institutes of Health Grant NS54281 under the Collaborative Research in Computational Neuroscience program.
Correspondence should be addressed to Srisairam Achuthan, Neuroscience Center of Excellence, Louisiana State University Health Sciences Center, 2020 Gravier Street, Suite D, New Orleans, LA 70112. sachut{at}lsuhsc.edu

References

↵
1. Acker CD,
2. Kopell N,
3. White JA
(2003) Synchronization of strongly coupled excitatory neurons: relating network behavior to biophysics. J Comput Neurosci 15:71–90.
OpenUrl CrossRef PubMed
↵
1. Bartos M,
2. Vida I,
3. Frotscher M,
4. Geiger JRP,
5. Jonas P
(2001) Rapid signaling at inhibitory synapses in a dentate gyrus interneuron network. J Neurosci 21:2687–2698.
OpenUrl Abstract/FREE Full Text
↵
1. Belykh I,
2. de Lange E,
3. Hasler M
(2005) Synchronization of bursting neurons: what matters in the network topology. Phys Rev Lett 94:188101.
OpenUrl CrossRef PubMed
↵
1. Buzsáki G
(2006) Rhythms of the brain (Oxford UP, New York).
↵
1. Canavier CC,
2. Butera RJ,
3. Dror RO,
4. Baxter DA,
5. Clark JW,
6. Byrne JH
(1997) Phase response characteristics of model neurons determine which patterns are expressed in a ring circuit model of gait generation. Biol Cybern 77:367–380.
OpenUrl CrossRef PubMed
↵
1. Canavier CC,
2. Baxter DA,
3. Clark JW,
4. Byrne JH
(1999) Control of multistability in ring circuits of oscillators. Biol Cybern 80:87–102.
OpenUrl PubMed
↵
1. Cunningham MO,
2. Pervouchine DD,
3. Kopell NJ,
4. Whittington MA
(2004) Cellular and network mechanisms of slow activity (<1 Hz) in the entorhinal cortex. Soc Neurosci Abstr 30:638–9.
OpenUrl
↵
1. Dror RO,
2. Canavier CC,
3. Butera RJ,
4. Clark JW,
5. Byrne JH
(1999) A mathematical criterion based on phase response curves for stability in a ring of coupled oscillators. Biol Cybern 80:11–23.
OpenUrl CrossRef PubMed
↵
1. Ermentrout B
(1996) Type I membranes, phase resetting curves and synchrony. Neural Comput 8:979–1001.
OpenUrl CrossRef PubMed
↵
1. Ermentrout B
(2002) Simulating, analyzing, and animating dynamical systems (SIAM, Philadelphia).
↵
1. Ermentrout GB,
2. Chow CC
(2002) Modeling neural oscillations. Physiol Behav 77:629–633.
OpenUrl CrossRef PubMed
↵
1. Ermentrout GB,
2. Kopell NJ
(1990) Oscillator death in systems of coupled neural oscillators. SIAM J Appl Math 50:125–146.
OpenUrl CrossRef
↵
1. Ermentrout GB,
2. Kopell NJ
(1991) Multiple pulse interactions and averaging in coupled neural oscillators. J Math Biol 29:195–217.
OpenUrl CrossRef
↵
1. Ernst U,
2. Pawelzik K,
3. Geisel T
(1995) Synchronization induced by temporal delays in pulse-coupled oscillators. Phys Rev Lett 74:1570–1573.
OpenUrl CrossRef PubMed
↵
1. Goel P,
2. Ermentrout B
(2002) Synchrony, stability, and firing patterns in pulse-coupled oscillators. Physica D 163:191–216.
OpenUrl CrossRef
↵
1. Hairer E,
2. Wanner G
(1991) Stiff and differential-algebraic problems. Springer series in computational mathematics, Solving ordinary differential equations II (Springer, Berlin), Vol 14.
↵
1. Hammond C,
2. Bergman H,
3. Brown P
(2007) Pathological synchronization in Parkinson's disease: networks, models and treatments. Trends Neurosci 30:357–364.
OpenUrl CrossRef PubMed
↵
1. Hodgkin AL
(1948) The local electric changes associated with repetitive action in a non-medullated axon. J Physiol 107:165–181.
OpenUrl FREE Full Text
↵
1. Huguenard JR,
2. McCormick DA
(2007) Thalamic synchrony and dynamic regulation of global forebrain oscillations. Trends Neurosci 30:350–356.
OpenUrl CrossRef PubMed
↵
1. Izhikevich EM
(2007) Dynamical systems in neuroscience: the geometry of excitability and bursting (MIT, Cambridge, MA).
↵
1. Jeong HY,
2. Gutkin B
(2007) Synchrony of neuronal oscillations controlled by GABAergic reversal potentials. Neural Comput 19:706–729.
OpenUrl CrossRef PubMed
↵
1. Li YX,
2. Wang YQ,
3. Miura R
(2003) Clustering in small networks of excitatory neurons with heterogenous coupling strengths. J Comput Neurosci 14:139–159.
OpenUrl CrossRef PubMed
↵
1. Luo C,
2. Clark JW Jr.,
3. Canavier CC,
4. Baxter DA,
5. Byrne JH
(2004) Multimodal behavior in a four neuron ring circuit: mode switching. IEEE Trans Biomed Eng 51:205–218.
OpenUrl CrossRef PubMed
↵
1. Mancilla JG,
2. Lewis TJ,
3. Pinto DJ,
4. Rinzel J,
5. Connors BW
(2007) Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. J Neurosci 27:2058–2073.
OpenUrl Abstract/FREE Full Text
↵
1. Maran SK,
2. Canavier CC
(2008) Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks in which firing order is not always preserved. J Comput Neurosci 24:37–55.
OpenUrl CrossRef PubMed
↵
1. Morris C,
2. Lecar H
(1981) Voltage oscillations in the barnacle giant muscle fiber. Biophys J 35:193–213.
OpenUrl CrossRef PubMed
↵
1. Mulloney B,
2. Hall WM
(2007) Local and intersegmental interactions of coordinating neurons and local circuits in the swimmeret system. J Neurophysiol 98:405–413.
OpenUrl Abstract/FREE Full Text
↵
1. Netoff TI,
2. Acker CD,
3. Bettencourt JC,
4. White JA
(2005a) Beyond two-cell networks: experimental measurement of neuronal responses to multiple synaptic inputs. J Comput Neurosci 18:287–295.
OpenUrl CrossRef PubMed
↵
1. Netoff TI,
2. Banks MI,
3. Dorval AD,
4. Acker CD,
5. Haas JS,
6. Kopell N,
7. White JA
(2005b) Synchronization in hybrid neuronal networks of the hippocampal formation. J Neurophysiol 93:1197–1208.
OpenUrl Abstract/FREE Full Text
↵
1. Nowotny T,
2. Huerta R,
3. Rabinovich MI
(2008) Neuronal synchrony: peculiarity and generality. Chaos 18:037119.
↵
1. Oh M,
2. Matveev V
(2009) Loss of phase-locking in non-weakly coupled inhibitory networks of type-I neurons. J Comput Neurosci 26:303–320.
OpenUrl CrossRef PubMed
↵
1. Oprisan SA,
2. Canavier CC
(2001) Stability analysis of rings of pulse-coupled oscillators: the effect of phase resetting in the second cycle after the pulse is important at synchrony and for long pulses. Differ Equ Dyn Syst 9:243–258.
OpenUrl
↵
1. Oprisan SA,
2. Canavier CC
(2003) Stability analysis of entrainment by two periodic inputs with a fixed delay. Neurocomputing 52–54:59–63.
OpenUrl
↵
1. Oprisan SA,
2. Prinz AA,
3. Canavier CC
(2004) Phase resetting and phase locking in hybrid circuits of one model and one biological neuron. Biophys J 87:2283–2298.
OpenUrl CrossRef PubMed
↵
1. Pervouchine DD,
2. Netoff TI,
3. Rotstein HG,
4. White JA,
5. Cunningham MO,
6. Whittington MA,
7. Kopell NJ
(2006) Low dimensional maps encoding dynamics in entorhinal cortex and hippocampus. Neural Comput 18:2617–2650.
OpenUrl CrossRef PubMed
↵
1. Prinz AA,
2. Thirumalai V,
3. Marder E
(2003) The functional consequences of changes in the strength and duration of synaptic inputs to oscillatory neurons. J Neurosci 23:943–954.
OpenUrl Abstract/FREE Full Text
↵
1. Koch C,
2. Segev I
1. Rinzel J,
2. Ermentrout GB
(1998) in Methods in neuronal modeling from ions to networks, Analysis of neural excitability and oscillations, eds Koch C, Segev I (MIT, Cambridge, MA).
↵
1. Sharp AA,
2. O'Neil MB,
3. Abbott LF,
4. Marder E
(1993a) The dynamic clamp—artificial conductances in biological neurons. Trends Neurosci 16:389–394.
OpenUrl CrossRef PubMed
↵
1. Sharp AA,
2. O'Neil MB,
3. Abbott LF,
4. Marder E
(1993b) Dynamic clamp—computer-generated conductances in real neurons. J Neurophysiol 69:992–995.
OpenUrl Abstract/FREE Full Text
↵
1. Sieling FH,
2. Canavier CC,
3. Prinz AA
(2009) Predictions of phase-locking in excitatory hybrid networks: excitation does not promote phase-locking in pattern generating networks as reliably as inhibition. J Neurophysiol, in press.
↵
1. Singer W
(1999) Neural synchrony: a versatile code for the definition of relations? Neuron 24:49–65.
OpenUrl CrossRef PubMed
↵
1. Skinner FK,
2. Bazzazi H,
3. Campbell SA
(2005) Two-cell to N-cell heterogeneous, inhibitory networks: precise linking of multi-stable and coherent properties. J Comput Neurosci 18:343–352.
OpenUrl CrossRef PubMed
↵
1. Smarandache CR,
2. Mulloney B
(2008) Gradients in strength of synapses in a circuit that coordinates limb movements. Soc Neurosci Abstr 34:371–3.
OpenUrl
↵
1. Stein SG,
2. Grillner S,
3. Selverston AI,
4. Stuart DG
, eds (1997) Neurons, networks, and motor behavior (MIT, Cambridge, MA).
↵
1. Tort ABL,
2. Rotstein HG,
3. Dugladze T,
4. Gloveli T,
5. Kopell NJ
(2007) On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus. Proc Natl Acad Sci U S A 104:13490–13495.
OpenUrl Abstract/FREE Full Text
↵
1. Traub RD,
2. Jefferys JG
(1994) Are there unifying principles underlying the generation of epileptic afterdischarges in vitro? Prog Brain Res 102:383–394.
OpenUrl PubMed
↵
1. Uhlhaas PJ,
2. Singer W
(2006) Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52:155–168.
OpenUrl CrossRef PubMed
↵
1. Van Vreeswijk C,
2. Abbott LF,
3. Ermentrout GB
(1994) Inhibition, not excitation, synchronizes coupled neurons. J Comput Neurosci 1:313–321.
OpenUrl CrossRef PubMed
↵
1. Wang XJ,
2. Buzsáki G
(1996) Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J Neurosci 16:6402–6413.
OpenUrl Abstract/FREE Full Text
↵
1. Winfree AT
(1967) Biological rhythms and the behavior of populations of coupled oscillators. J Theor Biol 16:15–42.
OpenUrl CrossRef PubMed
↵
1. Winfree AT
(2001) The geometry of biological time (Springer, New York).

In this issue

View Full Page PDF

Citation Tools

Respond to this article

Request Permissions

Cited By...

Articles

Show more Articles

Behavioral/Systems/Cognitive

Show more Behavioral/Systems/Cognitive

[1] ↵
Acker CD,
Kopell N,
White JA
(2003) Synchronization of strongly coupled excitatory neurons: relating network behavior to biophysics. J Comput Neurosci 15:71–90.
OpenUrl CrossRef PubMed

[2] Acker CD,

[3] Kopell N,

[4] White JA

[5] ↵
Bartos M,
Vida I,
Frotscher M,
Geiger JRP,
Jonas P
(2001) Rapid signaling at inhibitory synapses in a dentate gyrus interneuron network. J Neurosci 21:2687–2698.
OpenUrl Abstract/FREE Full Text

[6] Bartos M,

[7] Vida I,

[8] Frotscher M,

[9] Geiger JRP,

[10] Jonas P

[11] ↵
Belykh I,
de Lange E,
Hasler M
(2005) Synchronization of bursting neurons: what matters in the network topology. Phys Rev Lett 94:188101.
OpenUrl CrossRef PubMed

[12] Belykh I,

[13] de Lange E,

[14] Hasler M

[15] ↵
Buzsáki G
(2006) Rhythms of the brain (Oxford UP, New York).

[16] Buzsáki G

[17] ↵
Canavier CC,
Butera RJ,
Dror RO,
Baxter DA,
Clark JW,
Byrne JH
(1997) Phase response characteristics of model neurons determine which patterns are expressed in a ring circuit model of gait generation. Biol Cybern 77:367–380.
OpenUrl CrossRef PubMed

[18] Canavier CC,

[19] Butera RJ,

[20] Dror RO,

[21] Baxter DA,

[22] Clark JW,

[23] Byrne JH

[24] ↵
Canavier CC,
Baxter DA,
Clark JW,
Byrne JH
(1999) Control of multistability in ring circuits of oscillators. Biol Cybern 80:87–102.
OpenUrl PubMed

[25] Canavier CC,

[26] Baxter DA,

[27] Clark JW,

[28] Byrne JH

[29] ↵
Cunningham MO,
Pervouchine DD,
Kopell NJ,
Whittington MA
(2004) Cellular and network mechanisms of slow activity (<1 Hz) in the entorhinal cortex. Soc Neurosci Abstr 30:638–9.
OpenUrl

[30] Cunningham MO,

[31] Pervouchine DD,

[32] Kopell NJ,

[33] Whittington MA

[34] ↵
Dror RO,
Canavier CC,
Butera RJ,
Clark JW,
Byrne JH
(1999) A mathematical criterion based on phase response curves for stability in a ring of coupled oscillators. Biol Cybern 80:11–23.
OpenUrl CrossRef PubMed

[35] Dror RO,

[36] Canavier CC,

[37] Butera RJ,

[38] Clark JW,

[39] Byrne JH

[40] ↵
Ermentrout B
(1996) Type I membranes, phase resetting curves and synchrony. Neural Comput 8:979–1001.
OpenUrl CrossRef PubMed

[41] Ermentrout B

[42] ↵
Ermentrout B
(2002) Simulating, analyzing, and animating dynamical systems (SIAM, Philadelphia).

[43] Ermentrout B

[44] ↵
Ermentrout GB,
Chow CC
(2002) Modeling neural oscillations. Physiol Behav 77:629–633.
OpenUrl CrossRef PubMed

[45] Ermentrout GB,

[46] Chow CC

[47] ↵
Ermentrout GB,
Kopell NJ
(1990) Oscillator death in systems of coupled neural oscillators. SIAM J Appl Math 50:125–146.
OpenUrl CrossRef

[48] Ermentrout GB,

[49] Kopell NJ

[50] ↵
Ermentrout GB,
Kopell NJ
(1991) Multiple pulse interactions and averaging in coupled neural oscillators. J Math Biol 29:195–217.
OpenUrl CrossRef

[51] Ermentrout GB,

[52] Kopell NJ

[53] ↵
Ernst U,
Pawelzik K,
Geisel T
(1995) Synchronization induced by temporal delays in pulse-coupled oscillators. Phys Rev Lett 74:1570–1573.
OpenUrl CrossRef PubMed

[54] Ernst U,

[55] Pawelzik K,

[56] Geisel T

[57] ↵
Goel P,
Ermentrout B
(2002) Synchrony, stability, and firing patterns in pulse-coupled oscillators. Physica D 163:191–216.
OpenUrl CrossRef

[58] Goel P,

[59] Ermentrout B

[60] ↵
Hairer E,
Wanner G
(1991) Stiff and differential-algebraic problems. Springer series in computational mathematics, Solving ordinary differential equations II (Springer, Berlin), Vol 14.

[61] Hairer E,

[62] Wanner G

[63] ↵
Hammond C,
Bergman H,
Brown P
(2007) Pathological synchronization in Parkinson's disease: networks, models and treatments. Trends Neurosci 30:357–364.
OpenUrl CrossRef PubMed

[64] Hammond C,

[65] Bergman H,

[66] Brown P

[67] ↵
Hodgkin AL
(1948) The local electric changes associated with repetitive action in a non-medullated axon. J Physiol 107:165–181.
OpenUrl FREE Full Text

[68] Hodgkin AL

[69] ↵
Huguenard JR,
McCormick DA
(2007) Thalamic synchrony and dynamic regulation of global forebrain oscillations. Trends Neurosci 30:350–356.
OpenUrl CrossRef PubMed

[70] Huguenard JR,

[71] McCormick DA

[72] ↵
Izhikevich EM
(2007) Dynamical systems in neuroscience: the geometry of excitability and bursting (MIT, Cambridge, MA).

[73] Izhikevich EM

[74] ↵
Jeong HY,
Gutkin B
(2007) Synchrony of neuronal oscillations controlled by GABAergic reversal potentials. Neural Comput 19:706–729.
OpenUrl CrossRef PubMed

[75] Jeong HY,

[76] Gutkin B

[77] ↵
Li YX,
Wang YQ,
Miura R
(2003) Clustering in small networks of excitatory neurons with heterogenous coupling strengths. J Comput Neurosci 14:139–159.
OpenUrl CrossRef PubMed

[78] Li YX,

[79] Wang YQ,

[80] Miura R

[81] ↵
Luo C,
Clark JW Jr.,
Canavier CC,
Baxter DA,
Byrne JH
(2004) Multimodal behavior in a four neuron ring circuit: mode switching. IEEE Trans Biomed Eng 51:205–218.
OpenUrl CrossRef PubMed

[82] Luo C,

[83] Clark JW Jr.,

[84] Canavier CC,

[85] Baxter DA,

[86] Byrne JH

[87] ↵
Mancilla JG,
Lewis TJ,
Pinto DJ,
Rinzel J,
Connors BW
(2007) Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. J Neurosci 27:2058–2073.
OpenUrl Abstract/FREE Full Text

[88] Mancilla JG,

[89] Lewis TJ,

[90] Pinto DJ,

[91] Rinzel J,

[92] Connors BW

[93] ↵
Maran SK,
Canavier CC
(2008) Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks in which firing order is not always preserved. J Comput Neurosci 24:37–55.
OpenUrl CrossRef PubMed

[94] Maran SK,

[95] Canavier CC

[96] ↵
Morris C,
Lecar H
(1981) Voltage oscillations in the barnacle giant muscle fiber. Biophys J 35:193–213.
OpenUrl CrossRef PubMed

[97] Morris C,

[98] Lecar H

[99] ↵
Mulloney B,
Hall WM
(2007) Local and intersegmental interactions of coordinating neurons and local circuits in the swimmeret system. J Neurophysiol 98:405–413.
OpenUrl Abstract/FREE Full Text

[100] Mulloney B,

[101] Hall WM

[102] ↵
Netoff TI,
Acker CD,
Bettencourt JC,
White JA
(2005a) Beyond two-cell networks: experimental measurement of neuronal responses to multiple synaptic inputs. J Comput Neurosci 18:287–295.
OpenUrl CrossRef PubMed

[103] Netoff TI,

[104] Acker CD,

[105] Bettencourt JC,

[106] White JA

[107] ↵
Netoff TI,
Banks MI,
Dorval AD,
Acker CD,
Haas JS,
Kopell N,
White JA
(2005b) Synchronization in hybrid neuronal networks of the hippocampal formation. J Neurophysiol 93:1197–1208.
OpenUrl Abstract/FREE Full Text

[108] Netoff TI,

[109] Banks MI,

[110] Dorval AD,

[111] Acker CD,

[112] Haas JS,

[113] Kopell N,

[114] White JA

[115] ↵
Nowotny T,
Huerta R,
Rabinovich MI
(2008) Neuronal synchrony: peculiarity and generality. Chaos 18:037119.

[116] Nowotny T,

[117] Huerta R,

[118] Rabinovich MI

[119] ↵
Oh M,
Matveev V
(2009) Loss of phase-locking in non-weakly coupled inhibitory networks of type-I neurons. J Comput Neurosci 26:303–320.
OpenUrl CrossRef PubMed

[120] Oh M,

[121] Matveev V

[122] ↵
Oprisan SA,
Canavier CC
(2001) Stability analysis of rings of pulse-coupled oscillators: the effect of phase resetting in the second cycle after the pulse is important at synchrony and for long pulses. Differ Equ Dyn Syst 9:243–258.
OpenUrl

[123] Oprisan SA,

[124] Canavier CC

[125] ↵
Oprisan SA,
Canavier CC
(2003) Stability analysis of entrainment by two periodic inputs with a fixed delay. Neurocomputing 52–54:59–63.
OpenUrl

[126] Oprisan SA,

[127] Canavier CC

[128] ↵
Oprisan SA,
Prinz AA,
Canavier CC
(2004) Phase resetting and phase locking in hybrid circuits of one model and one biological neuron. Biophys J 87:2283–2298.
OpenUrl CrossRef PubMed

[129] Oprisan SA,

[130] Prinz AA,

[131] Canavier CC

[132] ↵
Pervouchine DD,
Netoff TI,
Rotstein HG,
White JA,
Cunningham MO,
Whittington MA,
Kopell NJ
(2006) Low dimensional maps encoding dynamics in entorhinal cortex and hippocampus. Neural Comput 18:2617–2650.
OpenUrl CrossRef PubMed

[133] Pervouchine DD,

[134] Netoff TI,

[135] Rotstein HG,

[136] White JA,

[137] Cunningham MO,

[138] Whittington MA,

[139] Kopell NJ

[140] ↵
Prinz AA,
Thirumalai V,
Marder E
(2003) The functional consequences of changes in the strength and duration of synaptic inputs to oscillatory neurons. J Neurosci 23:943–954.
OpenUrl Abstract/FREE Full Text

[141] Prinz AA,

[142] Thirumalai V,

[143] Marder E

[144] ↵
Koch C,
Segev I
Rinzel J,
Ermentrout GB
(1998) in Methods in neuronal modeling from ions to networks, Analysis of neural excitability and oscillations, eds Koch C, Segev I (MIT, Cambridge, MA).

[145] Koch C,

[146] Segev I

[147] Rinzel J,

[148] Ermentrout GB

[149] ↵
Sharp AA,
O'Neil MB,
Abbott LF,
Marder E
(1993a) The dynamic clamp—artificial conductances in biological neurons. Trends Neurosci 16:389–394.
OpenUrl CrossRef PubMed

[150] Sharp AA,

[151] O'Neil MB,

[152] Abbott LF,

[153] Marder E

[154] ↵
Sharp AA,
O'Neil MB,
Abbott LF,
Marder E
(1993b) Dynamic clamp—computer-generated conductances in real neurons. J Neurophysiol 69:992–995.
OpenUrl Abstract/FREE Full Text

[155] Sharp AA,

[156] O'Neil MB,

[157] Abbott LF,

[158] Marder E

[159] ↵
Sieling FH,
Canavier CC,
Prinz AA
(2009) Predictions of phase-locking in excitatory hybrid networks: excitation does not promote phase-locking in pattern generating networks as reliably as inhibition. J Neurophysiol, in press.

[160] Sieling FH,

[161] Canavier CC,

[162] Prinz AA

[163] ↵
Singer W
(1999) Neural synchrony: a versatile code for the definition of relations? Neuron 24:49–65.
OpenUrl CrossRef PubMed

[164] Singer W

[165] ↵
Skinner FK,
Bazzazi H,
Campbell SA
(2005) Two-cell to N-cell heterogeneous, inhibitory networks: precise linking of multi-stable and coherent properties. J Comput Neurosci 18:343–352.
OpenUrl CrossRef PubMed

[166] Skinner FK,

[167] Bazzazi H,

[168] Campbell SA

[169] ↵
Smarandache CR,
Mulloney B
(2008) Gradients in strength of synapses in a circuit that coordinates limb movements. Soc Neurosci Abstr 34:371–3.
OpenUrl

[170] Smarandache CR,

[171] Mulloney B

[172] ↵
Stein SG,
Grillner S,
Selverston AI,
Stuart DG
, eds (1997) Neurons, networks, and motor behavior (MIT, Cambridge, MA).

[173] Stein SG,

[174] Grillner S,

[175] Selverston AI,

[176] Stuart DG

[177] ↵
Tort ABL,
Rotstein HG,
Dugladze T,
Gloveli T,
Kopell NJ
(2007) On the formation of gamma-coherent cell assemblies by oriens lacunosum-moleculare interneurons in the hippocampus. Proc Natl Acad Sci U S A 104:13490–13495.
OpenUrl Abstract/FREE Full Text

[178] Tort ABL,

[179] Rotstein HG,

[180] Dugladze T,

[181] Gloveli T,

[182] Kopell NJ

[183] ↵
Traub RD,
Jefferys JG
(1994) Are there unifying principles underlying the generation of epileptic afterdischarges in vitro? Prog Brain Res 102:383–394.
OpenUrl PubMed

[184] Traub RD,

[185] Jefferys JG

[186] ↵
Uhlhaas PJ,
Singer W
(2006) Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron 52:155–168.
OpenUrl CrossRef PubMed

[187] Uhlhaas PJ,

[188] Singer W

[189] ↵
Van Vreeswijk C,
Abbott LF,
Ermentrout GB
(1994) Inhibition, not excitation, synchronizes coupled neurons. J Comput Neurosci 1:313–321.
OpenUrl CrossRef PubMed

[190] Van Vreeswijk C,

[191] Abbott LF,

[192] Ermentrout GB

[193] ↵
Wang XJ,
Buzsáki G
(1996) Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J Neurosci 16:6402–6413.
OpenUrl Abstract/FREE Full Text

[194] Wang XJ,

[195] Buzsáki G

[196] ↵
Winfree AT
(1967) Biological rhythms and the behavior of populations of coupled oscillators. J Theor Biol 16:15–42.
OpenUrl CrossRef PubMed

[197] Winfree AT

[198] ↵
Winfree AT
(2001) The geometry of biological time (Springer, New York).

[199] Winfree AT

Main menu

User menu

Search