Desynchronization boost by non-uniform coordinated reset stimulation in ensembles of pulse-coupled neurons

Several brain diseases are characterized by abnormal neuronal synchronization. Desynchronization of abnormal neural synchrony is theoretically compelling because of the complex dynamical mechanisms involved. We here present a novel type of coordinated reset (CR) stimulation. CR means to deliver phase resetting stimuli at different neuronal sub-populations sequentially, i.e., at times equidistantly distributed in a stimulation cycle. This uniform timing pattern seems to be intuitive and actually applies to the neural network models used for the study of CR so far. CR resets the population to an unstable cluster state from where it passes through a desynchronized transient, eventually resynchronizing if left unperturbed. In contrast, we show that the optimal stimulation times are non-uniform. Using the model of weakly pulse-coupled neurons with phase response curves, we provide an approach that enables to determine optimal stimulation timing patterns that substantially maximize the desynchronized transient time following the application of CR stimulation. This approach includes an optimization search for clusters in a low-dimensional pulse coupled map. As a consequence, model-specific non-uniformly spaced cluster states cause considerably longer desynchronization transients. Intriguingly, such a desynchronization boost with non-uniform CR stimulation can already be achieved by only slight modifications of the uniform CR timing pattern. Our results suggest that the non-uniformness of the stimulation times can be a medically valuable parameter in the calibration procedure for CR stimulation, where the latter has successfully been used in clinical and pre-clinical studies for the treatment of Parkinson's disease and tinnitus.


INTRODUCTION
Pathological neuronal synchronization is a hallmark of several neurological disorders like Parkinson's disease (PD), essential tremor, tinnitus, or epilepsy (Lenz et al., 1994;Nini et al., 1995;Mormann et al., 2000;Weisz et al., 2005Weisz et al., , 2007Kane et al., 2009;Schnitzler et al., 2009;Roberts et al., 2010), whereas the neuronal firing is uncorrelated in the normal state (Nini et al., 1995;Wilson et al., 2004) such that the abnormal synchronization is associated with pathology and symptoms (Levy et al., 2000). The standard therapy for medically refractory PD patients is electrical deep brain stimulation (DBS), where a high-frequency (HF >100 Hz) electrical pulse train is administered to target brain areas via chronically implanted depth electrodes (Benabid et al., 1991). HF DBS is found to significantly alter the neuronal activity of the stimulated neurons: The neuronal firing can be suppressed in the vicinity of the stimulation electrode, whereas the neurons are overactivated in the output structures of the stimulated neuronal population such that the pathological firing pattern is changed, see (Deniau et al., 2010) for review.
The mechanism of HF DBS is not fully understood. A modeling study by Wilson et al. (2011) suggests that the HF periodic DBS may induce a chaotic desynchronization, while a desynchronizing impact of a periodic forcing on synchronized populations seems to be a rather general phenomenon (Popovych and Tass, 2011). Also, as shown computationally, the effect of HF DBS strongly depends on the target structures stimulated (Hauptmann and Tass, 2007): Delivering HF DBS (nearly) exclusively to excitatory target structures may cause a desynchronization, whereas a stronger involvement of inhibitory target structures typically causes a pronounced inhibition. In some patients, however, HF DBS may be ineffective or cause side effects (Limousin et al., 1999;Kumar et al., 2003;Volkmann, 2004;Freund, 2005;Rodriguez-Oroz et al., 2005). Accordingly, along the lines of a model-based approach (Tass, 1999) novel stimulation techniques have been developed (Tass, 2001a(Tass, , 2003aRosenblum and Pikovsky, 2004;Hauptmann et al., 2005;Popovych et al., 2005Popovych et al., , 2006Pyragas et al., 2007;Popovych and Tass, 2010), which selectively counteract the pathological synchronization and restore uncorrelated neuronal firing. They are based on either phase resetting or feedback approaches, where from a theoretical standpoint the latter might have a great potential in controlling pathological synchronization, but so far have not been technically realized.
Other methods suggest to stimulate a single oscillator in the network (Nabi and Moehlis, 2011), drive the neurons into a phaseless set in order to achieve desynchronization (Danzl et al., 2009), or focused on the optimization of the standard HF DBS via a closed-loop stimulation setup (Feng et al., 2007a,b). In monkeys rendered parkinsonian with the neurotoxin MPTP Rosin et al. (2011) studied closed-loop DBS under acute conditions. To this end, they delivered a short train (comprising 7 pulses at 130 Hz) through a pair of electrodes located in the Globus pallidus internus (GPi) at a predetermined, fixed latency (80 ms) following each action potential recorded through an electrode placed in the primary motor cortex (M1). This type of stimulation caused a strong decrease of the firing rate of the pallidal neurons together with a pronounced decrease of the oscillatory neuronal activity at tremor frequency (4-7 Hz) and at double tremor frequency (9-15 Hz) along with an amelioration of the MPTP-induced akinesia. After cessation of this type of closed-loop DBS the initial firing pattern reverted back, i.e., pallidal firing rate and pallidal oscillatory activity attained pre-stimulus levels (Rosin et al., 2011). In contrast, standard continuous 130 Hz DBS caused a less pronounced decrease of the pallidal firing rate, the oscillatory neuronal activity and the amelioration of the akinesia (Rosin et al., 2011).
In contrast to the closed-loop DBS tested by Rosin et al. (2011), CR stimulation can be performed in a closed-loop or an openloop mode (Tass, 2003a,b). For instance, an adaptation of the stimulation frequency to the dominant frequency of the pathological neuronal synchronized collective oscillation can increase its efficacy (Tass, 2003a,b), see also (Tass et al., 2009). However, CR stimulation is robust with respect to variations of both stimulation and model parameters as follows from both computational as well as pre-clinical and clinical studies (Tass, 2003b;Tass et al., 2012a,b). More importantly, the goal of the CR approach is fundamentally different. CR stimulation does not aim at a decrease of firing rates and/or an abolishment of oscillatory neuronal activity. Rather, CR stimulation aims at specifically counteracting pathological synchrony by desynchronization (Tass, 2003a,b). This is because neurons have to be active in order to unlearn their pathological synaptic connectivity. In this way sustained long-lasting desynchronization is induced, and, as predicted computationally (Tass and Majtanik, 2006;Hauptmann and Tass, 2007;Tass and Popovych, 2012), therapeutic effects were observed in rat hippocampal slice experiments in the context of epilepsy (Tass et al., 2009) as well as in a clinical proof of concept study in tinnitus patients treated with acoustic CR stimulation (Tass et al., 2012a), where both pathological neuronal synchrony  and pathological effective connectivity (Silchenko et al., 2013) considerably decreased. In addition, in parkinsonian (MPTP) monkeys it was shown that unilateral CR stimulation delivered to the subthalamic nucleus (STN) for only 2 h per day during 5 days leads to significant and sustained therapeutic aftereffects for at least 30 days, while standard 130 Hz DBS has no aftereffects (Tass et al., 2012b). Another motivation why CR stimulation is delivered at frequencies similar to the pathological oscillatory frequency, is that in this case the desynchronizing effect is achieved at favorably small stimulation intensities (Tass, 2003a,b). In fact, as shown computationally, CR stimulation is able to strongly alter neuronal firing rates if delivered at frequencies substantially different to the dominant frequency of the stimulated neuronal population (Lysyansky et al., 2011b). For instance, CR stimulation may effectively activate hypo-or inactive neuronal populations without inducing neuronal synchrony (Lysyansky et al., 2011b).
The CR stimulation (Tass, 2003a,b) is based on the phase reset of oscillatory neuronal activity and has a broad applicability, since the phase reset is a universal phenomenon and can be achieved for a variety of stimulation setups and conditions, see, e.g., Refs. (Winfree, 1977;Paydarfar and Eldridge, 1987;Brandt, 1997;Makeig et al., 2002;Neiman et al., 2007;Thorne et al., 2011). According to its stimulation protocol, CR stimulation counteracts synchronization in the neuronal target population by dividing the entire population into several sub-populations where the phases of the neuronal oscillators within each sub-population are reset by the stimulation sequentially, i.e., in a timely coordinated manner. In this way, the collective neuronal oscillations in the sub-populations get phase shifted with respect to each other, and the total synchronization is replaced by, e.g., a cluster state (Tass, 2003a,b;Lysyansky et al., 2011a). Due to the pathologically strong synaptic connectivity, the entire target population runs from the cluster state through a transient characterized by pronounced desynchronization and finally resynchronizes if left unperturbed. Accordingly, to keep the neuronal ensemble in a desynchronized state, CR stimuli are delivered intermittently (Tass, 2003a,b), for instance, by applying CR in an m:n ON-OFF mode, where a few, say m, cycles with CR are optimally followed by a few, say n, cycles without any stimulation (Lysyansky et al., 2011a), for example, with m = 3 and n = 2 (Tass et al., 2012a,b). Such a stimulation protocol has computationally been found to be effective in inducing transient desynchronization in stimulated neuronal ensembles (Tass, 2003a,b;Lysyansky et al., 2011a).
The post-stimulus transient, where the stimulation-free neurons undergo an unperturbed desynchronized dynamics, plays an important role for the emergence of long-lasting effects of CR stimulation. In computational models taking into account the adaptive synapses governed by spike timing-dependent plasticity (STDP) (Gerstner et al., 1996;Markram et al., 1997;Feldman, 2000;Wittenberg and Wang, 2006;Caporale and Dan, 2008), it has been shown that CR stimulation can lead to a reduction of the mean synaptic weight and, in turn, shift the network to a state characterized by desynchronized activity and weak connectivity (Tass and Majtanik, 2006;Hauptmann and Tass, 2007) which persists after the stimulation is switched off. Modeling shows that CR stimulation is effective for a number of stimulation setups, in particular, for direct somatic stimulation as well as for excitatory or inhibitory synaptically-meditated stimulation which corresponds to stimulation of afferent or efferent fibers . This is particularly important since it has been shown that stimulation of fibers projecting to the STN appears to be responsible for the therapeutic effect of HF DBS delivered through STN electrodes (Gradinaru et al., 2009).
In order to further optimize the therapeutic benefit of CR stimulation, in this paper we investigate the impact of the stimulation parameters and the stimulation protocol on the stimulationinduced desynchronization. In particular, we focus on how the timing of the phase resets of the neuronal sub-populations influences the quality of the stimulation-induced cluster state and the post-stimulation transient. We found that appropriately adapted non-uniform stimulation onsets for the different stimulation sites can divide the phases of the stimulated neurons in such a way that the desynchronized post-stimulation transient gets significantly prolonged, until the population eventually resynchronizes again. To confine the complexity of our analysis, we study phase models without STDP. Phase resetting can be incorporated in such models in a natural way. We here set out to determine an optimal pattern of phase resets for CR stimulation. Put otherwise, we address the question of how to optimally choose the stimulation onsets for the single stimulation sites in CR.
In weakly coupled networks of oscillators, the technique of averaging is often applied to obtain a coupling which involves only the relative phase differences of the interacting oscillators (Ermentrout and Kopell, 1991;Swift et al., 1992;Hoppensteadt and Izhikevich, 1997). In globally coupled systems of identical neurons, these averaged phase models always possess symmetric cluster states, i.e., states, in which m clusters of equal size exist with a phase distance of 2π/m between neighbors (Okuda, 1993). Hence, a natural answer to the question above is to choose the stimulation times such that the phases get uniformly distributed over one period, independently on the type of neurons. In this case target patterns of the stimulation are the symmetric cluster states. However, in non-averaged models, the importance of the coupling between the single neurons becomes apparent. It plays an important role in determining the exact way of how the stimulation should be applied to cause a longer transient. In this work, we use systems of globally pulse-coupled phase oscillators for modeling the dynamics of a neuronal population. In particular, the symmetric cluster states disappear generically, and non-symmetric cluster states become possible candidates as stimulation target states. We propose a method for computing the stimulation times, which resets the system to a suitable cluster state. The timing points of the applied stimuli in these cases are non-uniformly spaced. The desynchronizing post-stimulation transient after such a stimulation turns out to be longer than the corresponding post-stimulation transient after a uniform CR stimulation of the same system. The paper is organized as follows. In section 2.1 we introduce the globally pulse-coupled model that is used to describe the collective dynamics of the neurons. In section 2.2 we study the relevant dynamical properties of the model, i.e., the appearance of synchronization, symmetric clusters and splay states. CR stimulation is introduced to the model in section 2.3, and we derive how one should apply CR stimuli to obtain longer post-stimulation desynchronization transients in section 2.4. The theoretical analysis is illustrated by numerical simulations in section 3. In particular, the robustness of the results to the variation of the stimulation parameters (stimulation intensity and electrode activation time) is studied in section 3.2. The effects of inhomogeneous frequencies is studied in section 3.3.

PULSE-COUPLED PHASE OSCILLATORS
Phase models play a key role in describing the individual dynamics of single oscillators, e.g., oscillatory neurons, see e.g., (Kuramoto, 1984). In particular, a stable periodic dynamics can be modeled by a simple equation for the periodic motion of the phase with frequency ω: ϕ(t) = ϕ (0) + ωt, where ϕ is considered modulo 2π. The direction of the phase is neutrally stable. Therefore, a sufficiently weak temporary perturbation, which does not move the original system far away from the corresponding limit cycle, persists in the phase for all times, while all its other effects die out exponentially fast due to the stability of the limit cycle which corresponds to the periodic motion. In coupled systems, weak interactions can be conceived as perturbations, and the phase reduction can be applied as well (Kuramoto, 1984;Hansel et al., 1993;Hoppensteadt and Izhikevich, 1997;Brown et al., 2004). In fact, phase models are particularly important for studying network dynamics, because many types of synchronization, which are of interest in such models, depend on the relative phases of the units combining the network (Pikovsky et al., 2001).
The effect of perturbations is incorporated into the phase equations by the phase response curve (PRC) (Guckenheimer, 1975;Kuramoto, 1984;Ermentrout, 1996;Winfree, 2001) (see Figure 1). It measures the response of the individual neuron to weak stimuli. We consider the case when it admits a representation in terms of a smooth scalar function Z(ϕ) of the phase. For example, applying a weak current I(t) to a neuron with PRC Z(ϕ) changes its phase dynamics froṁ which describes the phase dynamics in the weak stimulation limit (Kuramoto, 1984). If the perturbation is pulse-like, e.g., a brief electrical stimulation or a synaptic input from another neuron (if synapses are fast), it may be approximated as an instantaneous input which resets the neuron's phase at time t = t 0 of the incoming pulse as following where ϕ t − 0 = lim t↑t 0 ϕ(t) and ϕ t + 0 = lim t↓t 0 ϕ(t). We can formally write the pulse-like perturbation as I(t) = I · δ (t − t 0 ), using the Dirac delta-function (Goel and Ermentrout, 2002). In what follows, we study a system of N identical phase oscillators  (Brown et al., 2004). (B) PRC Z (ϕ) = Z ML (ϕ) of a particular Morris-Lecar type model (Ermentrout, 1996;Sato et al., 2011). which are globally pulse-coupled with weight I = N (Goel and Ermentrout, 2002)

Frontiers in Computational
where t k, are the times where the k-th neuron spikes. For convenience we assume that a single neuron emits a spike if its phase crosses zero (mod 2π). In such a way, if the phase ϕ k of neuron k passes through the spike, where the resetting function μ is defined as In the case that n neurons simultaneously spike at time t = t * , the reset value of neuron ϕ j is taken to be ϕ j t + * = μ n ϕ j t − * , where μ n (ϕ) = μ (μ (. . . μ(ϕ) . . .)) denotes the n-fold superposition of the resetting function. This mapping is preferable to the choice of ϕ j t + * = ϕ j t − * + n N Z ϕ j t − * because it assures a continuous dependence on the initial conditions. More specifically, the effect of the spiking of a cluster of n neurons changes continuously as the neurons split from the cluster, i.e., μ (τ 1 + μ (τ 2 + μ (. . . μ(ϕ) . . .))) → μ n (ϕ) as the interspike intervals τ i → 0. Note that the resetting function for the cluster spikes can also be measured or computed directly in order to achieve more realistic modeling (Achuthan and Canavier, 2009). Applied to our case it would mean using a given measured function μ(ϕ, n) instead of μ n (ϕ) in the case when an n-cluster spikes. However, we restrict our analysis to the function μ n (ϕ) defined as a superposition of individual resetting functions.
Note, that for sufficiently small values of /N the resetting function μ is monotone. This ensures the preservation of the ordering of the phases. We assume that N is sufficiently large and this property holds. For some systems (Goel and Ermentrout, 2002;Brown et al., 2004;Stiefel et al., 2009), it turned out, that the PRC has a small value at the spike moment Z(0) ≈ Z(2π) ≈ 0. For simplicity, we assume that Z(0) = Z(2π) = 0 holds for our model. Let us shortly explain why this is a reasonable approximation for weakly coupled spiking systems. Firstly, the modulus of the PRC must be roughly proportional to the density of isochrons (Guckenheimer, 1975;Winfree, 2001) at the corresponding point of the limit cycle of the full system. This density on the other hand is inversely proportional to the modulus of the vector field, which is large at the spiking point. Therefore, the modulus of the PRC has to be small at the spiking point.
Note, even though this work is focused on systems with pulse coupling (1), the main qualitative message about non-uniform CR stimulation and non-uniform positions of clusters still holds for systems of the forṁ with a smooth, periodic coupling function G, which was proposed in (Winfree, 2001). Systems (1) and (3) were less extensively studied than their averaged versions, which take the formφ where (Ermentrout and Kopell, 1991;van Vreeswijk et al., 1994;Daido, 1997;Kuramoto, 1997). As a result of the averaging, the stability properties of the corresponding solutions of (1) [or (3)] and (4) may differ at the order of O 2 . This is the same magnitude as of the errors made by the phase reduction, and, thus, studying the averaged system suggests itself as a simpler and, presumably, equivalent task. In the next section we show that an important genericity of stationary solutions with distributed phases is overlooked by this choice. In fact, the homogeneous stationary solutions or symmetric clusters in the averaged system (4) correspond to some other, generically non-homogeneous solutions of the original system, whose shape may differ at the order of O ( ), see, for instance, Equations (12) and (13) below. A precise targeting of these solutions by CR can essentially contribute to the efficacy of the desynchronization technique. Since systems (1) and (3) admit non-homogeneous stationary solutions with distributed phases, it is of particular interest to study them in the context of CR.

SYNCHRONIZATION, CLUSTERS, AND STATIONARY SOLUTIONS
In this section we review the dynamical properties of the stimulation-free population (1), which are relevant in the context of desynchronization. In particular, we study the appearance and stability of a synchronized state, cluster states, as well as splay states. We pay special attention to the differences between the mentioned dynamical states of models (1), (3), and (4).

Stability of synchronized solutions
In each of the systems (1), (3) and (4), there exists an in-phase synchronous solution where all neurons are perfectly synchronized The conditions for the stability of the synchronous state are well known for all of the above systems (Goel and Ermentrout, 2002). In particular, for the pulse-coupled system (1), in-phase synchronization is linearly stable iff For (3) the corresponding condition reads For the averaged system (4), the stability condition of the synchronized state is A comparison of conditions (5)-(7) leads to the following relationships: -Vanishing coupling: Conditions (7) and (6) for the averaged and non-averaged systems differ only in the second order of , respectively. Therefore, they coincide in the limit of small , it is positive and concentrated at ϕ = 0, then and the condition for synchronization of the averaged system (7) coincides with the condition for the non-averaged pulsecoupled system (5) provided Z (0) = 0.
Although in this work we will focus on smooth PRCs, it is worth to mention another synchronization effect, which may generically occur for PRCs, which have discontinuous derivatives as a phase of zero is approached from the left and right (Lücken and Yanchuk, 2012): though the synchronous state is locally unstable in this case, the first order parameter R 1 (t) = 1 N k exp (iϕ k (t)) may still approach its maximal value R 1 = 1 due to the appearance of structurally and dynamically stable homoclinic connections to the synchronous state.

Splay states and stationary solutions
Splay states are periodic solutions of system (1), in which all oscillators are spread in a way that the time differences between the subsequent spikes t k + 1, − t k, are always the same Zillmer et al., 2007). Note that the term "splay state" can also be used differently and may, more generally, refer to any state featuring phases which are spread over the periodic interval [0, 2π] (Achuthan and Canavier, 2009). Splay states are important for desynchronization issues, since they possess small order To study splay states in large systems it is useful to consider an equation for the phase distribution density ρ(t, ϕ), since its stationary solution approximates the distribution of the phases of splay states in the limit of large N. For the pulse-coupled system (1) the dynamics of the phase distribution density is governed by the following continuity equation (Ernst et al., 1995;Brown et al., 2004) Its equivalent for the smoothly-coupled model (3) is and for the averaged system (4) one gets Solving (8) for stationary solutions ρ(t, ϕ) = ρ s (ϕ) and taking into account the normalization 2π 0 ρ s (ψ)dψ = 1, we obtain Here ρ s (0) is uniquely determined by the implicit equation obtained from (11) by inserting ϕ = 0 (see Appendix A.1 for details). Thus, (11) describes a unique stationary solution of (8).
For small , this solution can be approximated as whereZ = 1 2π 2π 0 Z (ψ) dψ is the mean value of the PRC. For smoothly-coupled systems (3), one analogously finds that a unique stationary solution of (9) exists, if is not too large (see Appendix A.1 ). Its first-order expansion in reads For the averaged model (4), we find that for any value of the constant distribution densityρ(ϕ) ≡ 1 2π is a stationary solution of (10).
As follows from Equations (12) and (13), phase distributions of the splay states of the pulse-coupled system (1) as well as the non-averaged system (3) deviate from a uniform distribution. For small , the deviations can be estimated in the first order in by (12) and (13), respectively. This is in contrast to splay states of the averaged system (4), which are always uniformly distributed. and Z ML (ϕ) (Figure 1B), respectively. Figure 3 shows that the theoretically obtained stationary phase distribution density ρ s (ϕ) (black curve) convincingly approximates the numerically calculated phase distribution histogram (gray bars) of the splay state in the pulse-coupled system (1) for large number N of oscillators.
The stability of splay states as well as cluster states can be studied by the set of multipliers λ of the corresponding return map, which determine the rate, with which small perturbations from the considered state are growing with time. In particular, if the multiplier with the maximal absolute value is λ max , then a generic perturbation will grow as λ j max with j is the number of spikes of a neuron from the ensemble. Therefore, a stable state corresponds to the case |λ max | < 1 and an unstable to |λ max | > 1. We calculated numerically the maximal multipliers λ max of the splay state of system (1) for two cases: (Figure 4, blue pluses). We observe that the splay states are unstable for all N except N = 3 in case (ii), and their multipliers converge to an asymptotic value λ ∞ > 1 as the number of neurons N increases. Note that the splay state may be stable for synaptic coupling and type-I PRC, as shown in (Achuthan and Canavier, 2009) for a system of a few coupled oscillators. For a more detailed analysis of the stability of splay states in large FIGURE 3 | Comparison of splay state and stationary phase distribution density. Black solid curve depicts the theoretical solution of the stationary phase distribution density ρ s (ϕ) of (8) for the pulse-coupled system (1) with a PRC of Morris-Lecar type Z (ϕ) = Z ML (ϕ) (Figure 1B). Gray bars illustrate the numerically computed and normalized phase distribution histogram of the corresponding splay state for N = 1000 oscillators. Coupling strength = 3.

Symmetric clusters
In this paper we consider strongly synchronized neuronal ensembles, where the splay state is expected to be unstable. When an external desynchronization technique would be able to move the system in a vicinity of such a state, the achieved order parameter would be very low for a relatively long time. In the case, when the number of stimulation sites is naturally limited to a low number, a natural substitute for the target state of CR stimulation is a cluster state. For the sake of simplicity we will consider symmetric cluster states consisting of m clusters, each of them containing N/m neurons. Within each cluster, the neurons are synchronized and have the same phase, whereas the phases of different clusters are shifted with respect to each other. For systems (1) and (3) there exists a unique stationary, symmetric m-cluster state, at least for moderate coupling strength. The phases ψ j of individual clusters are not equidistantly distributed, i.e., in general |ψ j+1 (t) − ψ j (t)| = 2π/m. In contrast, in the averaged systems of the form (4), equidistantly distributed m-clusters with |ψ j+1 (t) − ψ j (t)| = 2π/m always exist. In Appendix B, we explain how one can determine multipliers of symmetric cluster states for system (1). We have computed them for the cases Figures 4A,B). It is convenient to distinguish between tangential and transverse multipliers, which correspond to perturbations within the invariant cluster space and perturbations which disperse the single clusters by destroying the perfect synchrony within the clusters, respectively. For case (i), and 1 ≤ m ≤ 25, all m-cluster states are unstable. In case (ii), there is a single value m = 3 for which the cluster state is stable. The tangential multipliers seem to asymptote towards a limit with increasing m (Figure 4, red circles). If the transverse multipliers are smaller than the tangential, as it is the case for both PRCs (Figure 4, red markers), one expects that the perturbed dynamics stays near the invariant cluster space as long as the linear prediction is valid. Detailed analysis of cluster states for different PRCs has been performed in (Ashwin and Swift, 1992;Hansel et al., 1993;Okuda, 1993;Chandrasekaran et al., 2011;Lücken and Yanchuk, 2012).

MODELING DEEP BRAIN STIMULATION
Strong enough electrical stimuli or synaptic input can reset the phase of a neuron in such a way that its oscillation restarts after the stimulation from a definite phase (Winfree, 1977;Best, 1979;Tass, 1999;Popovych and Tass, 2012). The general mechanism of phase resetting can be understood as follows. In a stimulated neuron, a stable steady (e.g., hyperpolarized) state appears, which is approached during the stimulation. When the stimulation terminates, the steady state disappears, and the system relaxes back to the limit cycle with an asymptotic phase ϕ s determined by the isochron on which the stimulation specific steady state was located. However, for the phase resetting procedure which we propose, the exact value of ϕ s is not essential. Therefore we include the stimulus in the simplest way which provides the model with a qualitatively adequate stimulus response. The phase model with the incorporated stimulus readṡ where I j (t) corresponds to the stimulus applied to the j-th neuron. We assume, that there are m stimulation sites, each stimulating a distinct group of N/m neurons (Figure 5). These stimulation sites can either be active or not, that means where I j,i are stimulation time intervals and I is the stimulation intensity. For two neurons ϕ j and ϕ k from the same group we have S j = S k and I j (t) ≡ I k (t).
The dynamics of uncoupled neurons under stimulation is described asφ If IZ(ϕ) < 0 for some ϕ, and if the intensity I is of sufficiently large magnitude, there appears a pair of fixed points, stable ϕ s and unstable ϕ u , satisfying ω + IZ(ϕ) = 0. If only one such pair exists, the neuron will approach the stable fixed point ϕ s after some time of stimulation and stay there until the stimulation terminates (see Figure 6). In such a situation, we call ϕ s the resetting point.
The stimulation described by (14) aims at establishing a distribution of phases of the neuronal ensemble that prolongs the post-stimulus transient as much as possible before the ensemble synchronizes again. In principle, the strategy is to establish a state as close as possible to some stationary desynchronized state. With a given number m of stimulation sites which influence equally large groups of N/m neurons, the target states for the control are restricted to m-cluster configurations with clusters of equal size. We call a target pattern the state, which is intended to be realized at the end of the stimulation. A series of successive activations and deactivations of the stimulation sites is called stimulation sequence, and a time interval during which the resetting stimuli are delivered at all m sites is called stimulation cycle. For the averaged system (4), the target pattern consists of equidistant clusters ψ j = 2π m j + ϕ s , 1 ≤ j ≤ m, such that the last stimulation-induced cluster is located at the resetting point ψ m = ϕ s at the end of the stimulation cycle. For systems of type (1) or (3), the target pattern is a cluster state ψ j = 2π m j + ϕ s + O( ), which is in general not equidistantly distributed. Among all possible stimulation sequences, we restrict our considerations to those where each stimulation site is activated once per stimulation cycle. The activity of the j-th stimulation site is confined to the time interval S j = t j , t j + τ , where t 1 , . . . , t m are the onset times within the stimulation cycle, and τ is the stimulation duration, which is the same for all stimulation sites. In practice, stimulation sequences have to be administered repeatedly after the system recovers to some undesired level of synchronization (Tass, 2003a).

STIMULATION-INDUCED STATIONARY M-CLUSTER STATES
Since the stimulation target pattern has to be established as precisely as possible, one has to take into account the influences of the coupling among neurons on the stimulation-induced pattern. We now describe how appropriate stimulation sequences can be found when we restrict ourselves to the stimulation timing S j = t j , t j + τ , j = 1, . . . , m, mentioned above. The long-term desynchronizing effects of such a stimulation sequence will be discussed in the next section 3. For the brevity of notations, let us introduce the resetting map of the stimulated system (14) by where ϕ ∈ [0, 2π] N is the state of the system at the onset of the stimulation at t = t on def = min{t 1 , . . . , t m }, and t = (t 1 , . . . , t m ) is the vector of the onset times of the stimulation sites. (t; ϕ) describes the N-dimensional state of the system at the end of the stimulation. If the duration τ and magnitude of the stimulation intensity I are large enough (see section 3.2), each neuron ϕ k of the j-th group is reset to the collective cluster phase ϕ k t j + τ ≈ ϕ s at the offset time t j + τ of the j-th stimulation site and continues to evolve in a cluster of a common phase which we denote by ψ j (t; ϕ). The dependence of ψ j on initial conditions ϕ follows from the fact that this cluster may still be influenced by the other clusters and neurons. In practice, the state ϕ before stimulation is unknown. To ensure full control over the resulting cluster state (t; ϕ), it should be independent of the initial state ϕ. To this aim we have to assume that the resetting mechanism results in an accurate reset of the stimulated neurons to a determined phase ϕ = ϕ s . Furthermore, we have to ensure that the reset neurons are not affected by neurons, which have not yet been reset in the current stimulation cycle. Both conditions can be fulfilled by choosing a sufficiently large duration τ and intensity I for the reset. Thus, we can approximately identify the state (t; ϕ) with the lower dimensional cluster state (t) = (ψ 1 (t off ) , . . . , ψ m (t off )), where t off = τ + max{t 1 , . . . , t m } denotes the offset time of the entire stimulation sequence, and moreover, (t) is independent on ϕ. Now, let * = ψ * 1 , . . . , ψ * m denote the phases of the stationary cluster state which serves as the stimulation target pattern and can be obtained as described in Appendix B. Then, the problem is to find a solution t of (t) = * .
In this study we do not aim to provide a general algorithm to solve this equation. Moreover, due to discontinuities that are caused by the pulse-like interactions of the ensemble, the function (t) is only piecewise smooth (see Figure 8). We approach the solution of (17) numerically by starting from the uniformly distributed stimulation sequence t 0 as an initial guess. Then we apply the minimization Nelder-Mead simplex search algorithm (Lagarias et al., 1998), which is implemented in the MATLAB function   The time elapsed between t* j and the preceding stimulation onset is denoted by t * j . Note the non-uniformity of the stimulation. See also the corresponding illustration in Figure 7. fminsearch, to minimize (t) − * . Table 1 and Figure 7 illustrate an example of the computed stimulation sequence for the case of four stimulation sites m = 4 and a Morris-Lecar type of PRC Z(ϕ) = Z ML (ϕ) (see Figure 1). For the parameter values given in caption to Table 1, the target pattern has been computed as well as the stimulation sequence t. The optimal stimulation sequence deviates by up to ∼ 4% from the uniformly distributed one where t * j is the same for all j = 2, . . . , m. Figure 7 illustrates the corresponding switching times of the stimulation contacts. It also shows the spiking times of the obtained clusters after the stimulation. In order to find the stimulation sequence, the discontinuous system (17) has been solved. Figure 8 illustrates the emerging types of discontinuities by plotting t * 1 , t * 2 , t * 3 , t 4 − * versus t 4 for fixed t * 1 , t * 2 , and t * 3 .

ADVANTAGES TO UNIFORM STIMULATION
Equidistant clusters are stationary if the coupling depends only on the phase differences as in system (4). In the non-averaged systems (1) Table 1.

FIGURE 8 | Illustration of the discontinuity of t
, t 4 − * is plotted versus t 4 for fixed t * 1 , t * 2 , and t * 3 . Parameters as in Table 1. The discontinuity occurs at such a value of t 4 = t d that leads to ψ 4 t − off = 2π, i.e., the onset of the post-stimulation spike of the cluster j = 4 just coincides with t off , see Figure 7. If t 4 > t d , the impact of the spike of cluster ψ 4 on the cluster ψ 2 has to be taken into account when calculating the resulting cluster positions.

FIGURE 9 | Advantages of the stimulation-induced target patterns consisting of stationary cluster states [plots (B), (D)] compared to those induced by conventional CR stimulation with equidistant stimulation times [plots (A), (C)].
The time courses of the first order parameter R 1 (t) (red curves) and the fourth order parameter R 4 (t) (blue curves) of the neuronal ensemble (14)  with PRCs Z = Z ML and Z = Z H (see Figure 1) are presented in Figure 9 where plots (A,C) illustrate the effect of the uniform CR stimulation and (B,D) are related to the non-uniform CR stimulation. Time courses of the first order parameter R 1 (t) (red curves) and the fourth order parameter R 4 (t) (blue curves), defined as R n = 1 N k exp (inϕ k (t)) ∈ [0, 1], are shown. Large values of the first order parameter are indicative of an in-phase synchronization. On the other hand, for approximately equidistantly distributed m-cluster states high values of the m−th order parameter are combined with low values of all order parameters with lower indices. We use these properties to detect synchronization and the discussed slightly non-uniform cluster states.
All simulations in Figure 9 are started at t = 0 with the neurons' phases randomly distributed in [0, 2π]. For both PRCs, Z = Z ML and Z = Z H , we observe a steady increase of R 1 (t) (red curves) in the pre-stimulation epoch t t on , which indicates the onset of in-phase synchronization of the entire ensemble. This process is significantly faster in the case Z = Z H (Figures 9C,D) than for Z = Z ML (Figures 9A,B) due to Z H (0) Z ML (0) 0 which strengthens the linear stability of the synchronous solution in the case Z = Z H (see section 2.2). When CR stimulation is turned on at t on = 250, the stimulated phases are successively caught at ϕ s and released when the corresponding stimulation site is deactivated. At the end of the stimulation cycle at t off ≈ 266.75, four clusters are established, which are well distributed in [0, 2π]. This leads to a low value of R 1 (t off ) ≈ 0 (A: R 1 (t off ) ≈ 0.027, B: R 1 (t off ) ≈ 0.020, C: R 1 (t off ) ≈ 0.000, D: R 1 (t off ) ≈ 0.044) and a high value of R 4 (t off ) ≈ 1 (A: R 4 (t off ) ≈ 0.987, B: R 4 (t off ) ≈ 0.994, C: R 4 (t off ) ≈ 1.000, D: R 4 (t off ) ≈ 0.969). All simulations show a post-stimulation transient before the system resynchronizes again and the first order parameter R 1 (t) approaches unity. For both considered PRCs, the advantages of the method proposed in section 2.4 to establish nonequidistant clusters is substantial. For the uniform CR stimulation (Figures 9A,C) the post-stimulation desynchronization transient is of approximately the same duration as the initial transient in the pre-stimulation epoch, when starting from a random distribution of the phases. The post-stimulation transient is significantly prolonged by the non-uniform CR stimulation in both cases, for Z = Z ML (Figure 9B) by doubling the transient duration and for Z = Z H (Figure 9D) by tripling the duration of the desynchronization transient. Note that small-scale oscillations of the order parameters originate from the discontinuities of the system's trajectory, which occur whenever a cluster crosses the firing threshold.

EFFECTS OF DIFFERENT STIMULATION INTENSITIES AND DURATION
An important question is how the stimulation intensity I and the stimulation duration τ influence the desynchronization transient. Figure 10 shows results of numerical simulations for the PRC Z = Z ML . An increase of the stimulation duration, τ ∈ [0, 10], leads to an increase in the desynchronized transient equally for both, uniform ( Figure 10A) and optimized, non-uniform stimulation timing ( Figure 10B). This is indicated by longer intervals of decreased order parameter R 1 after the stimulation between t on and t off . For the various stimulation lengths, we have t on = 200 and t off ∈ [205,217]. Beyond τ ≈ 4.5, the effect of the uniform timing does not enhance, while by the optimized protocol it increases further until about τ = τ c ≈ 6.6. This value corresponds to the duration which assures independence of the stimulation outcome on the system's state prior to stimulation (see section 2.4).
Similarly, small magnitudes of the stimulation intensity 1 |I| 2 yield a prolonged transient for both protocols (Figures 10C,D).
Beyond |I| = 2, no further increase can be observed for the uniform timing ( Figure 10C), but the transient after the nonuniform protocol continues to grow at least until |I| = 6 ( Figure 10D).  In (A,B), t off ranges from t off ≈ 205 to t off ≈ 217, and in (C,D), the stimulation interval is constant. In all cases the initial phases were drawn randomly from a uniform distribution on [0, 2π]. Other parameters: = 0.5, N = 240.

ROBUSTNESS TO VARIATIONS OF NATURAL FREQUENCIES
In the above approach we assumed that the neurons are identical. In more realistic situations, the parameters of individual neurons can vary. In order to test the robustness of the proposed resetting technique with respect to parameter changes, which break the symmetry of the system, we consider an ensemble (1) with non-identical natural frequencies ω j , j = 1, . . . , N, e.g., randomly chosen from a uniform distribution in [1 − ω, 1 + ω]. The results of simulations are shown in Figure 11. We present them for the PRC Z = Z ML (Figure 1), but qualitatively similar results have been obtained for Z = Z H as well. It turns out, that a significant prolongation of the post-stimulation desynchronization transient can be observed for a range of ω. Indeed, one observes a clear difference between the post-stimulation behavior of the order parameters R 1 and R 4 for the suggested CR stimulation with non-uniform stimulation timing (Figures 11C,D) and that for the conventional CR stimulation with equidistant stimulation times (Figures 11A,B). This effect of the optimized CR stimulation, however, decreases for a broader distribution of the natural frequencies. Nevertheless, our calculations suggest that the optimization procedure of CR stimulation can robustly improve its desynchronizing impact on neuronal populations exhibiting undesired synchronization. It is well known that the broadening of the frequency range induces a desynchronizing transition (Kuramoto, 1984) such that the system with a very broad frequency range does not synchronize even in the absence of stimulation. This transition occurs at larger width of the frequency distribution if the coupling strength is increased. In our illustration with coupling strength = 0.5 (Figure 11) the desynchronization transition and the loss of the advantageous effect of the proposed stimulation technique take place at approximately the same frequency mismatch. The suggested non-uniform CR stimulation is more effective than the uniform one for a range of ω which supports synchronized dynamics.
All results in the sections 3.2 and 3.3 hold within a range of the coupling strength, 0 < 3.5 for Z = Z ML and 0 < 5.0 for Z = Z H (data not shown). In these ranges, with larger , the re-synchronization after the stimulation happens faster, but the optimized protocol is still superior to the uniform. This means that the length of the transient following the non-uniform stimulation exceeds the one after the uniform by a similar factor as in the presented case for = 0.5. Moreover, the results are independent of the population size N if it is sufficiently large. In practice this is already the case for N ≥ 100, i.e., if each cluster contains 25 neurons.

DISCUSSION
A number of pulsatile stimulation techniques have been developed which enable to directly shift a synchronized neuronal population into a desynchronized state, irrespective of the initial state at which the stimulus is delivered (Tass, 2001a(Tass, ,b, 2002a. However, less favorably, these techniques require careful calibration of the stimulation parameters and their continuous adaptation to varying model parameters. To overcome this limitation and provide a desynchronizing stimulation technique which is robust and does not require time-consuming or technically involved calibration procedures, CR's indirect approach to desynchronization was developed (Tass, 2003a,b): Inducing a cluster state by means of time-shifted phase resetting stimuli delivered to different neuronal sub-populations can robustly be achieved and does not require relevant calibration (Tass, 2003a,b). The cluster states, in turn, are relevant since they lead to long post-stimulus desynchronized transients (Lysyansky et al., 2011a), and in the presence of STDP (Gerstner et al., 1996;Markram et al., 1997) the related decrease of the rate of coincidences induces an antikindling (Tass and Majtanik, 2006;Hauptmann and Tass, 2007;Tass and Popovych, 2012). Neither in preclinical nor in clinical studies adverse effects of CR stimulation have been observed (Tass et al., 2009(Tass et al., , 2012a. We considered model networks of weakly pulse-coupled neurons with phase resetting curves and compared them to averaged models, where the phase dynamics depends only on the phase differences between the oscillators. Whereas the latter models are better analytically tractable and attained a great attention in the literature (Strogatz, 2000;Winfree, 2001;Acebrón et al., 2005), they neglect some important information about stationary states of the original systems. In particular, the stationary splay and cluster states are not uniform for the pulse-coupled networks, contrary to those of the averaged models. These non-uniformly distributed cluster states can serve as target states for CR stimulation in ensembles of pulse-coupled neurons. We have found that the optimal stimulation sequence should be slightly non-uniform in order to approach the non-uniform cluster state at the end of the stimulation.
We have shown that the phase response curves of the stimulated neurons determine the phase distribution densities of splay and cluster states, which, in turn, influence the timing of the stimulation sequences. The proposed non-uniform stimulation sequences result in significant improvements of the stimulation outcome and lead to several times longer post-stimulation transients in comparison to the equally spaced stimulation sequences. Intriguingly, modifications of the stimulation timing points of only a few percent (e.g., 4%, see Figure 7) actually double or even triple the duration of the post-stimulation desynchronization transient (Figure 9). The proposed approach takes into account and compensates for the interactions among neurons during the stimulation.
We also showed that the discussed stimulation protocol is robust with respect to variation of the natural frequencies, stimulation parameters, and coupling strength. It can lead to a prolonged transient for a range of non-identical frequencies of the single oscillators. One can expect that the non-uniform stimulation technique can be superior to a series of equally timed stimulations in more diverse and realistic setups, where transmission delays and coupling functions are heterogeneous and this should be confirmed in further studies. Moreover, since the mechanism of the discussed desynchronizing method is based on the phase reset of the neuronal oscillations, which is a universal phenomenon of the neuronal dynamics (Winfree, 1977;Best, 1979;Tass, 1999), and the timing of the optimal stimulation sequence is determined by the phase response curves, the presented approach is generic and can be applied to other neuronal models and other stimulation-induced target states. In particular, to models, which employ PRCs of second, or higher order (Oprisan et al., 2004;Achuthan and Canavier, 2009). For our purpose, we restricted the investigation to the system (1), since it is one of the simplest models, which already possesses enough features illustrating our main finding, that the optimal stimulation timing is non-uniform. In the framework of the considered model one can, however, incorporate the PRCs measured either experimentally or obtained by detailed modeling of the globus pallidus and STN (Schultheiss et al., 2010;Farries and Wilson, 2012b,a) which are possible target regions for CR deep brain stimulation Tass et al., 2012b), and where a change in PRC structure might contribute to disease-related changes in synchronous activity. Another target for non-invasive acoustic CR neuromodulation adapted for the treatment of tinnitus Tass et al., 2012a) is the auditory cortex where a phase reset can be achieved by different types of auditory stimuli (Brandt, 1997;Thorne et al., 2011).
As yet, only slight modifications of the τ/m timing of CR stimulation (with onset times t 1 , t 1 + τ/m, t 1 + 2τ/m, . . . , t 1 + (m − 1)τ/m) have been investigated in the Kuramoto model as well as in an ensemble of synaptically coupled FitzHugh-Nagumo oscillators modeling spiking neurons in the context of M:N ON-OFF CR stimulation, where M cycles with CR are followed by N cycles without stimulation (Lysyansky et al., 2011a). For that stimulation protocol a τ/m timing of CR stimulation was used. However, in one variant of this protocol the M-th stimulation cycle was prematurely terminated at the break time t off , where t off < τ. An optimal choice of t off (e.g., t off = 0.44τ for a particular set of model parameters tested) caused a pronounced increase of the desynchronization transient, i.e., an increase of the time elapsing to resynchronization by a factor of approx. 2. In contrast, inappropriate values of t off induced a decrease (e.g., a halving) of the resynchronization time.
To study post-stimulus desynchronization transients of two phase oscillators with time-delayed coupling subject to coincident, but phase shifted stimulation the transient time T tr (defined as the time it takes a trajectory after stimulus offset to permanently enter into an ε vicinity of the stable phase-locked state) was computed (Krachkovskyi et al., 2006). For vanishing time delays the phase space of that two-oscillator model is simple and the optimal phase shift in the coupling term puts the system's phase difference onto an unstable fixed point (Figure 9a in Krachkovskyi et al., 2006). In contrast, for non-vanishing time delay the phase space gets considerably more complex, and the optimal phase shift puts the system onto a particular point in phase space where the system gets trapped by a stable manifold, leading to a particularly high transient time T tr . Incorporating time delays into the coupling of the model studied in this paper will certainly increase the complexity of its phase space. Given the results by Ref. (Krachkovskyi et al., 2006), we expect that such timed delays may have an impact on the resynchronization transient and, hence, on the optimal timing of CR stimulation.
Another important direction for future analysis comes from the fact that biological networks typically comprise neurons of different kind. Consequently, our approach has to be extended to mixed populations, containing neurons of different type by, for example, including inhibitory neurons found in human STN (Levesque and Parent, 2005). In this work we considered a simple all-to-all coupling topology which provides an easy way to obtain a synchronized neuronal dynamics serving as a model for the dynamical regimes encountered in Parkinsonian patients and monkeys (Nini et al., 1995;Levy et al., 2000). This was also motivated by the reported strong functional connectivity of tremor-related neurons in STN (Amtage et al., 2009) and anatomical intranuclear connectivity as follows from experimental and modeling studies (Iwahori, 1978;Kita et al., 1983;Willshaw, 1998, 2004;Shen and Johnson, 2006), but see Ref. (Wilson et al., 2004). The considered weak coupling is supported by the observed gradual decay and recovery of pathological oscillations at the onset and offset of DBS (Kang and Lowery, 2011). CR stimulation has been shown to work for other coupling topologies and stimulation setups, e.g., for sensory stimulation Tass and Popovych, 2012;Tass et al., 2012a). For further details of the effects of CR stimulation, more realistic coupling topologies and sophisticated neuronal models as well as connections to other neural populations within the basal ganglia and cortical brain areas have to be considered. Also, the spatial spread of the stimulation current, 3D effects as well as optimized electrode geometries (see e.g., Buhlmann et al., 2011) have to be taken into account in future studies. Pursuing such studies, our ultimate goal is to come up with CR sequences which enable to further minimize the stimulation current for DBS. This might contribute to a decrease of the rate of side effects caused by stimulation spread to neighboring brain areas. By the same token, this might enable considerably smaller geometries of the implantable pulse generator due to a significant reduction of battery size.

A.1 PULSE COUPLED SYSTEMS
In this section we show how a stationary solution ρ (t, ϕ) = ρ s (ϕ) for (8) can be found and show that it is unique. It follows from (8) that this solution satisfies with some constant C. Hence .
The constant C can be determined from the normalization condition Combining (A1) and (A2), we obtain the expression (11) for ρ s (ϕ). Evaluating (11) at ϕ = 0 and taking into account Z(0) = 0, the equation for ρ s (0) reads To investigate the resolvability of this equation we consider the function which is well-defined and smooth on x ∈ (0, ζ) with ζ = ∞ in case Z(ϕ) ≥ 0, and ζ = −ω min ψ Z(ψ) otherwise. It is easy to see that F (0) = 1 2π > 0. Furthermore, we will show that The conditions (i) and (ii) imply that F (x) has a unique root, which corresponds to a unique solution ρ s (ϕ) of (8). In order to show (i), we treat two cases separately. First, if ζ < ∞, then the smoothness of Z (ψ) leads to and, therefore, lim x→ζ F (x) = −ζ. In the second case, when ζ = ∞, we use the following estimates: Since Z(ψ) is smooth and Z(0) = 0, one can always find a large enough constant L > 0 such that Hence lim x → ∞ F (x) = −∞. Now let us show that (ii) holds. The second derivative of F (x) can directly be calculated and reads The sign of F coincides with the sign of F 1 . We apply the Cauchy-Schwarz's inequality to obtain 2π 0 g 2 fdψ Thus, F 1 ≤ 0 holds and claim (ii) is proven. We have shown that the equation for the phase density for the pulse-coupled system (8) has a unique stationary solution. Note that smallness of was not assumed here. The first order approximation in of this density takes the form (12).

A.2 SMOOTHLY COUPLED SYSTEMS (3) AND THEIR AVERAGED COUNTERPART (4)
To unify notations, we write the dynamics for the single neurons in the limit of large population aṡ The phase density equation reads and its stationary solution ρ s (ϕ) fulfills with some constant C depending on ρ s (·). Solving for ρ s (ϕ) and integration over [0, 2π] with respect to ϕ gives an expression for C We see that the stationary solution is generally aperturbation of the constant state. Without going into details, we remark that for small the right hand side of (A4) is a contraction with respect to ρ s (·) on all sets

B. STABILITY OF SYMMETRIC CLUSTER STATES
In this section, we present a method for calculating characteristic multipliers of the symmetric cluster states of system (1) (see section 2.2.3). To this end, we introduce a discrete time map (B1) which describes the system's evolution in between two subsequent spikes. The main idea of our analysis is a separation of the effects of tangential and transversal perturbations to the cluster state. This allows us to derive the expressions (B3) and (B5) for the characteristic multipliers which can be evaluated easily. Finally, we explain how these multipliers behave asymptotically as N → ∞.
Let us consider a time moment t 0 , when some neuron of the ensemble (1) just emitted a spike. We choose its index to be N and all other indices in such a way, that the phases are ordered as 2π ≥ ϕ 1 (t + 0 ) ≥ · · · ≥ ϕ N (t + 0 ) = 0. Until the next spiking event, all neurons advance with the same phase velocity ω. Without loss of generality, we can assume ω = 1. Then the next neuron ϕ 1 reaches the spiking threshold ϕ 1 (t − 1 ) = 2π at time t 1 = t 0 + 2π − ϕ 1 (t + 0 ). The other neurons ϕ k are now located at . After the spike, the first neuron is reset to ϕ 1 (t + 1 ) = 0 and all others to Because of the monotonicity of the resetting map μ, the order of the phases is preserved as 2π ≥ ϕ 2 (t + 1 ) ≥ · · · ≥ ϕ N (t + 1 ) ≥ ϕ 1 (t + 1 ) = 0. We define the firing map as F : T N → T N , T = [0, 2π] / {0 ∼ 2π} (the points 0 and 2π are identified as the same point) componentwise by and F N (ϕ) := 0. It captures the threshold crossing and spiking of the oscillator ϕ 1 and shifts the indices k → k − 1. This map takes the ordered phases after some spiking event 2π ≥ ϕ 1 ≥ · · · ≥ ϕ N = 0 and maps them to ordered phases after the successive spiking event 2π ≥ F 1 (ϕ) ≥ · · · ≥ F N (ϕ) = 0. Since we assume the resetting function μ(ϕ) to be smooth on T, F is also smooth on T N . The N-fold iteration of the firing map R := F • · · · • F = F N is called the return map. Given an initial state ϕ with ϕ 1 ≥ · · · ≥ ϕ N , the map R(ϕ) returns the state after each neuron has fired once. In order to simplify notations, we adopt double indices, ϕ ,j := ϕ ( −1)n+j , to address the j-th member of the -th cluster. A stationary, symmetric m-cluster state is a fixed point ϕ * = F n (ϕ * ) for the n = N/m-th iterate of the firing map, which satisfies ϕ * ,j = ψ with 1 ≤ ≤ m, 1 ≤ j ≤ n, and ordered cluster positions 2π > ψ 1 > · · · > ψ m = 0. For such a symmetric cluster state, the time τ elapsed between two successive threshold crossings of clusters is the same. This time τ can be determined as the smallest solution of the equation G m−1 τ (0) + τ = 2π. Here the function G τ (x) := μ n (x + τ) describes the phase resetting of a neuron located at position x + τ induced by a cluster spike. The value of G m−1 τ (0) = G τ • G τ • · · · • G τ (0) is the position of a cluster initially located at ψ = 0 after m − 1 other clusters have fired. The positions of other clusters are then given as ψ = G m− τ (0). The linear stability of the cluster state ϕ * can be determined by its characteristic multipliers with respect to R, i.e., the eigenvalues of DR (ϕ * ). In the following we calculate these multipliers. Using double indices for F in the same way as for ϕ, we have F n ,j ϕ * = μ n (ψ +1 + 2π − ψ 1 ) = ψ , where the index + 1 is considered modulo m within the range from 1 to m. For a perturbed cluster state, the phase of each neuron can be written as ϕ ,j = ψ j + η ,j . To study the dynamics of the perturbations, we introduce the following subspaces W = m k = 1 η k,j = η k, i , for 1 ≤ i, j ≤ n , V = k = η k, j = 0, for 1 ≤ j ≤ n ∩ η ,n = 0 , 1 ≤ ≤ m. Elements of W correspond to perturbations, which rule the relative motions of the clusters. Correspondingly, the subspace V contain perturbations, which split the -th cluster. Since the perturbations from V cannot split any other clusters, we have This implies DF n ϕ * V ⊂ V −1 ⊕ W, and DF n ϕ * W ⊂ W.
By m applications of the map DF n = DF n (ϕ * ), we obtain the following sequence Hence, for DR (ϕ * ) = (DF n (ϕ * )) m , we have the skew structure DR (ϕ * ) V ⊂ W ⊕ V and DR (ϕ * ) W ⊂ W. This implies that the matrix DR (ϕ * ) has a block triangular form in an appropriate basis and the spectrum of eigenvalues of DR (ϕ * ) is the union of the spectra, which are restricted to the subspaces W and V . The eigenvalues belonging to the subspace V , can be found by considering all possible perturbations in V . These perturbations can be described as a linear combination of the following basis vectors η i k,j = 1, for k = , 1 ≤ j ≤ i, 0, else, for 1 ≤ i ≤ n − 1. Each element η i corresponds to a different way to split the -th cluster into two subclusters. It is possible to show that the multiplier corresponding to η i is given as the change of the distance of these two subclusters during one return.
Calculating the perturbed dynamics step by step (from one cluster firing to the next) and linearizing in η i leads to DR ϕ * η i = λη i + w with some w ∈ W and the eigenvalue λ independent on and i.
The perturbations in W can be described by the following basis vectors η k,j = 1, for k = , 1 ≤ j ≤ n 0, else.
Each element η corresponds to a perturbation of the -th cluster, which does not split the cluster. Calculations (for brevity, we omit the details here) of the dynamics of η show that the eigenvalues of DR(ϕ * ) restricted to W are given as the m-th powers of the eigenvalues of the matrix This means that tangential multipliers of the symmetric m-cluster state are given as λ m where λ are the solutions of det(A − λI) = 0.