Abstract
Networks of model neurons were constructed and their activity was predicted using an iterated map based solely on the phaseresetting 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 alltoall 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 betweencluster splay and withincluster 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 phaseresetting 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).
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 fifthorder 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^{2}, 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_{Na}m_{∞}^{3}h(V − E_{Na}). The steadystate 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_{K}n^{4}(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^{2}, respectively. I_{stim} is the applied current and was set at 0.5 μA/cm^{2}. Unless otherwise stated, the values for the various parameters were equal to those given above.
The synaptic current is given by I_{syn} = g_{syn}s(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^{−1} 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^{2}, V is the cell membrane voltage in millivolts, and t is time in milliseconds. The calcium current is given by I_{Ca} = g_{Ca}m_{∞}(V)(V − E_{Ca}). The leak current is given by I_{L} = g_{L}(V − E_{L}). The steadystate activation is as follows: where V_{1} = −1.2 mV and V_{2} = 18 mV. The potassium current is given by I_{K} = g_{K}w(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 steadystate activation amplitude as follows: and activation rate as follows: with V_{3} = 2 mV and V_{4} = 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^{2}, respectively. For calcium, the maximal conductance (g_{Ca}) was set to 4.4 mS/cm^{2} and I_{stim} was set at 100.0 μA/cm^{2}. 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 firstorder resetting, f_{1}(φ) = (P_{1} − P_{i})/P_{i}, and secondorder resetting, f_{2}(φ) = (P_{2} − P_{i})/P_{i}, respectively. P_{1} is the length of the cycle containing the perturbation and P_{2} is the length of the subsequent cycle (Fig. 2A,B).
Iterative pulsecoupled 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 phaseresetting 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_{2}(φ[k − 1]) is the secondorder 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_{1}(φ[k]) is the firstorder reset due to an input in the current cycle. The firstorder resetting is included in the recovery interval, whereas the secondorder resetting is included in the stimulus interval. If multiple inputs arrive during a cycle, the firstorder resetting is taken at the time that each input arrives. There are two ways to handle secondorder resetting. One is to discard all secondorder 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 secondorder 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 secondorder resetting caused a problem emulating a splay mode.
The map representing an ideal pulsecoupled 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 secondorder 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 secondorder 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 secondorder resetting was subtracted from the phase of the firing neuron(s) (see below for the implications of negative phase) and the saved secondorder resetting was cleared. Then any firstorder 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 firstorder resetting due to the firing neurons was subtracted from the incremented phase. The saved secondorder 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_{1}(φ[k]) is the firstorder reset due to the input in the current cycle, and f_{2}(φ[k − 1]) is the secondorder reset due to the input in the previous cycle. The firstorder 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 secondorder 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 pulsecoupled 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 firstorder 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 secondorder resets, although only the PRCs for firstorder resetting are shown. The magnitude of the secondorder 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 secondorder resetting due to the previous firing of the two neurons in the other cluster and the firstorder resetting due to the other neuron in the same cluster.
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^{2}. 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^{2} 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 pulsecoupled 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 phaselocked 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 alltoall 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 alltoall 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 pulsecoupled 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 ithorder PRC at a phase of φ. A stability criterion for an alltoall 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_{1}[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_{1}[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 secondorder 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_{1}′(φ, (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.
The following criterion includes the effect of secondorder 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 secondorder resetting is presumed to have a continuous derivative at 0 and 1, but the first and secondorder 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 firstorder resetting disappears as a phase of one is approached, from the left, and secondorder 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_{1}′(0^{+}, g_{syn}) − f_{2}′(1^{−}, (N − 1)g_{syn}) and λ = 1 − f_{1}′(0^{+}, (N − 1)g_{syn}) − f_{2}′(1^{−}, g_{syn}). However, if the derivative of the sum of the first and secondorder 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_{1}′(0^{+}, g_{syn}) − f_{1}′(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 firstorder 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 φ_{1}. 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} >…> φ_{2}> φ_{1} > 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 secondorder 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],…, Δφ_{2}[k], Δφ_{1}[k]]^{T} and S is the following matrix:
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 PRCbased 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 firstorder PRCs at the locking points, indicated by f_{1}′(φ_{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 φ_{1}*, the second at φ_{2}*, 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_{1}(φ_{i−1}*)} and The first interval includes the secondorder resetting from the previous cycle. Subsequent intervals calculate the elapsed time between inputs and account for the firstorder 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 φ_{1}* that makes ts_{1}[∞] = ts_{N}[∞]. Then the values of the φ_{i}* that make ts_{i}[∞] equal to ts_{N}[∞] and ts_{1}[∞] 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 Mdimensional network. Next the existence of a betweencluster splay mode in an N/Mdimensional network is determined by modifying the existence criteria given above for the splay mode such that the first interval includes the firstorder 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_{1}(0, (M − 1)g_{syn}): The stability of the betweencluster 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 withincluster modes.
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 fourneuron 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 pulsecoupled 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 fourneuron 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).
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 secondorder 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 lefthand panels show the PRC for a representative conductance value. The center trace shows the iterative predictions, and the righthand 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).
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 fourneuron 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^{2}, 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 fourneuron 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^{2}, the splay mode was only accurately predicted by the iterative map if only secondorder 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 fourneuron 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 withincluster 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 firstorder 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 fourneuron 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.
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 betweencluster 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 fourneuron 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, withincluster interactions of M neurons are analyzed by studying the synchronous mode for M neurons, and betweencluster 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 withincluster synchrony is predicted to be stable (blue bar), betweencluster 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 withincluster synchrony (blue diamonds) and betweencluster splay (green squares). Note that since we plot the absolute value of the eigenvalues, the blue diamonds in Figure 10B1 corresponding to the withincluster synchronous mode reverse direction when the eigenvalue passes through zero.
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 withincluster 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^{2} to 0.07 mS/cm^{2}). Therefore the cluster mode is predicted to be stable by the analytical approach until 0.07 mS/cm^{2}. Beyond this the withincluster 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 betweencluster mode can somewhat compensate for the destabilizing influence of withincluster 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 withincluster synchrony for the twoneuron network is unstable for all values of g_{syn}. Figure 10A2 has no blue bar corresponding to the parameter region that supports stable withincluster 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 withincluster interaction is >1.0 for g_{syn} values from 0.0 mS/cm^{2} onwards. The prediction that a twoneuron 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 fourneuron network is predicted to be unstable by the analytical approach. However, the observed cluster mode in a fourneuron network is stable throughout the given conductance range. This strongly suggests that the betweencluster interaction contributes to stabilizing the cluster mode. A small betweencluster eigenvalue again seems to compensate for a withincluster 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 oversimplification, the method was successfully applied in networks of four and 12 neurons.
Clustering in a 12neuron network
The above results are not specific to any particular network size. As an illustrative example we use a 12neuron network that is a scaled up version of the fourneuron 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 withincluster synchrony, the green bar predicts stable betweencluster splay mode, and the black bar predicts unstable betweencluster 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^{2}.
Figure 11A1 shows that the two clusters of six mode is analytically predicted to be stable up to g_{syn} = 0.02 mS/cm^{2} (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 withincluster synchrony as conductance is increased. The betweencluster 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 betweencluster splay mode is responsible for limiting the number of clusters that can be observed.
One of the relevant slopes for withincluster synchrony is always f_{1}′(0, g_{syn}). The other is f_{1}′(0, 5g_{syn}) for a cluster of six, f_{1}′(0, 3g_{syn}) for a cluster of four, and f_{1}′(0, 2g_{syn}) for a cluster of three. Withincluster 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 12neuron 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^{2}. 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 phaselocked 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 halfcenter 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 alltoall networks of pulsatile coupled oscillators (Goel and Ermentrout, 2002) did not address secondorder resetting, but only firstorder resetting at phase just after spike initiation (0^{+}) and just before spike initiation (1^{−}). For pulsatile coupling, firstorder resetting at 1^{−} is negligible, but the secondorder resetting at 1^{−} is equal to the firstorder resetting at 0^{+}. This allows the stability criterion at synchrony to collapse to λ = 1 − f_{1}′(0^{+}, g_{syn}) − f_{1}′(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 pulsecoupled network of bursting HindmarshRose 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 betweencluster splay mode is necessary for clusters to form. If withincluster 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 twoneuron 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 pulsecoupled systems enabled by bistability in the twoneuron network, as well as clustering where synchrony in the twoneuron 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 twoneuron 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 phaseclustered states have been mostly found in networks with reciprocal inhibition, except in an alltoall excitatory network in which withincluster 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 alltoall 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 fastspiking 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 fastspiking and stellate cells, are thought to be firing at theta frequency. The betweencluster analysis methods developed here show that in general it is progressively more difficult to find PRCs that support betweencluster 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 PRCbased 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 twoneuron networks to quantitatively assess larger network dynamics (Skinner et al., 2005, Netoff et al., 2005a) or reducing the larger network to a twoneuron 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 steadystate terms from both sides of the equations, we obtain the following: The two equations above constitute a discrete linear system: where Δ[k] = [Δφ_{1}[k], Δφ_{N−1}[k]]^{T} and M is a 2by2 matrix.
The characteristic polynomial corresponding to the matrix M is as follows: For synchrony, the above expression has to be evaluated twice, once at φ_{1}* = 0^{+} and φ_{N−1}* = 1^{−}, and again at φ_{1}* = 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}′(1^{−}, (ig_{syn})) and f_{2}′(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_{1}′(0^{+}, ig_{syn}) = f_{2}′ (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] >…> φ_{2}[k]> φ_{1}[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 firstorder resetting (f_{1}(φ)), which is not indexed by g_{syn} here because only one input is received at a time, and ignoring higherorder 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_{1}(φ_{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} >…> φ̂_{1}. 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 φ_{1}[k + 1] in the stimulus interval ts_{1} 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 φ_{1}[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 steadystate fixed point. After linearizing the PRCs such that f(φ) = f(φ*) + f′(φ*)Δφ and canceling the steadystate 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],…, Δφ_{2}[k], Δφ_{1}[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 pulsecoupled 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 nonnegative. 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 nonnegative. In the case of no simultaneous inputs, this interval is equal to the phase of the next input φ_{i} plus the secondorder resetting f_{2}(φ_{j}) from the previous input, hence −f_{2}(φ_{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 secondorder 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 firstorder 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 secondorder 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 pulsecoupled map based on a strictly memoryless assumption cannot produce the fully synchronous mode or synchronous subclusters because ignoring all secondorder 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 pulsecoupled map to reproduce these modes. If we include all secondorder 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 secondorder resetting can disrupt ability of the pulsecoupled map to reproduce the splay mode in Figure 8C2 when the intervals are short compared with those obtained by including secondorder resetting from each input. For the parameter regime g_{syn} = 0.09–0.17 mS/cm^{2}, the pulsecoupled map was reprogrammed to ignore the secondorder 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 secondorder 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 betweencluster interactions can produce withincluster 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