Desynchronizing electrical and sensory coordinated reset neuromodulation

Coordinated reset (CR) stimulation is a desynchronizing stimulation technique based on timely coordinated phase resets of sub-populations of a synchronized neuronal ensemble. It has initially been computationally developed for electrical deep brain stimulation (DBS), to enable an effective desynchronization and unlearning of pathological synchrony and connectivity (anti-kindling). Here we computationally show for ensembles of spiking and bursting model neurons interacting via excitatory and inhibitory adaptive synapses that a phase reset of neuronal populations as well as a desynchronization and an anti-kindling can robustly be achieved by direct electrical stimulation or indirect (synaptically-mediated) excitatory and inhibitory stimulation. Our findings are relevant for DBS as well as for sensory stimulation in neurological disorders characterized by pathological neuronal synchrony. Based on the obtained results, we may expect that the local effects in the vicinity of a depth electrode (realized by direct stimulation of the neurons' somata or stimulation of axon terminals) and the non-local CR effects (realized by stimulation of excitatory or inhibitory efferent fibers) of deep brain CR neuromodulation may be similar or even identical. Furthermore, our results indicate that an effective desynchronization and anti-kindling can even be achieved by non-invasive, sensory CR neuromodulation. We discuss the concept of sensory CR neuromodulation in the context of neurological disorders.


INTRODUCTION
Several brain disorders, such as Parkinson's disease (PD), are characterized by abnormally strong, pathological neuronal synchronization (Lenz et al., 1994;Nini et al., 1995). The standard therapy for medically refractory PD patients is electrical highfrequency (HF) deep brain stimulation (DBS) (Benabid et al., 1991). For this, an electrical HF (>100 Hz) pulse train is delivered to a target brain area through a depth electrode.
The mechanism of action of HF DBS is not yet fully understood. Multiple possible mechanisms may contribute to the therapeutic effects of DBS, depending, e.g., on the neuronal structures being stimulated [for review see (Benabid et al., 2002;Vitek, 2002;McIntyre et al., 2004b)]. Some of these mechanisms take place in the vicinity of the stimulating electrode: for instance, neuronal activity might be blocked because HF DBS changes the activation of voltage-gated currents in the vicinity of the stimulating electrode, in this way inducing a depolarization blockade (Beurrier et al., 2001) or low-amplitude subthreshold oscillations of the membrane potential (Pyragas et al., 2012). Also, neuronal activity near the stimulating electrode might indirectly be inhibited via an excitation of axon terminals connected to neurons by inhibitory synapses (synaptic inhibition) (Dostrovsky et al., 2000). Depolarization blockade, stimulation-induced lowamplitude subthreshold oscillations, and synaptic inhibition are candidate mechanisms underlying the suppression of the somatic firing within the target nucleus. In contrast, independently of the suppression of neuronal activity in the vicinity of the stimulating electrode HF DBS may lead to a HF axonal output via an activation of projection neurons (McIntyre et al., 2004a). Put otherwise, local and non-local effects of HF DBS may differ considerably (McIntyre et al., 2004a).
In a certain number of patients 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, 2001(Tass, , 2003aRosenblum and Pikovsky, 2004;Hauptmann et al., 2005;Popovych et al., 2005Popovych et al., , 2006Pyragas et al., 2007;Popovych and Tass, 2010), which do not significantly influence the natural firing of the individual neurons and only selectively counteract the pathological synchronization of a neuronal target population. Coordinated reset (CR) stimulation (Tass, 2003a,b), considered in this paper, is a robust desynchronizing stimulation technique and aims at a therapeutic reshaping of the synaptic connectivity (Tass and Majtanik, 2006). CR stimulation can be applied in a closed-loop as well as in open-loop setup and requires neither extensive calibration nor technically involved real-time measurements and processing of the pathological neuronal activity as compared to other methods. According to its stimulation protocol, CR stimulation counteracts synchronization in the neuronal target population by splitting the entire population into several sub-populations where the phases of the neuronal oscillations are reset by the stimulation sequentially, i.e., in a timely coordinated manner (Tass, 2003a,b). 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., 2011). Due to the pathologically strong synaptic connectivity, the entire target population runs from the cluster state through a transient characterized by pronounced desynchronization. Accordingly, to exploit this quasi-reactive type of desynchronization, CR stimuli are delivered intermittently (Tass, 2003a,b), for instance, by applying CR in an m:n ON-OFF mode, where m cycles with CR are followed by n cycles without any stimulation (Lysyansky et al., 2011).
In computational studies one may take into account a fundamental property of neuronal tissue, where the weights of the adaptive synapses depend on the underlying neuronal firing via a spike timing-dependent plasticity (STDP) (Gerstner et al., 1996;Markram et al., 1997;Feldman, 2000;Wittenberg and Wang, 2006;Caporale and Dan, 2008). For the models with STDP it has been shown that the CR-induced decrease of the rate of coincident firing leads to a reduction of the mean synaptic weight and, in turn, shifts the network to a healthy model attractor state characterized by desynchronized activity and weak connectivity (Tass and Majtanik, 2006;Hauptmann and Tass, 2007). Accordingly, CR causes an unlearning of both pathological synchronization and pathological connectivity (so-called anti-kindling) (Tass and Majtanik, 2006). CR-induced long-lasting desynchronization has experimentally been verified in rat hippocampal slice rendered epileptic by magnesium withdrawal (Tass et al., 2009). In the present study we consider two models of neuronal networks of spiking Hodgkin-Huxley (HH) and bursting FitzHugh-Rinzel (FHR) neurons with excitatory and inhibitory synaptic coupling which is governed by a STDP rule.
In the computational studies so far, only direct electrical stimulation of the neurons' somata has been taken into account [see, e.g., (Hauptmann and Tass, 2007)]. In contrast, we here also consider indirect, synaptically-mediated CR stimulation. Our motivation is twofold.
1. Electrical deep brain CR stimulation: we study the differences of direct CR effects (mediated by direct somatic stimulation) in the vicinity of the depth electrode as opposed to indirect CR effects (mediated by stimulation of excitatory or inhibitory fibers). The latter stimulation modality may have local and non-local effects depending on that whether afferent axon terminals or efferent fibers are stimulated. 2. Sensory CR stimulation: we study whether it is possible to perform CR stimulation with indirect stimulation only. This is a prerequisite for sensory CR stimulation. The network models with indirect, synaptically-mediated stimulation studied in this paper may be considered as minimal models for sensory stimulation, where in a first step relay stations in the sensory pathway are neglected.
Since the mechanism of action of CR stimulation is based on the phase reset of oscillatory neuronal activity, we may expect that the discussed control technique will be effective for both direct somatic and indirect, synaptically-mediated stimulation. Indeed, the phase reset is a universal phenomenon and can be achieved for a variety of stimulation setups as reported by a great body of modeling and experimental studies, e.g., for hyperpolarizing or depolarizing electrical pulses (Best, 1979;Demir et al., 1997;Tateno and Robinson, 2007;Neiman et al., 2007), excitatory or inhibitory post-synaptic potentials (PSPs) (Perkel et al., 1964;Pinsker, 1977;Lerma and Garcia-Austt, 1985;Jackson et al., 2002;Prinz et al., 2003), sensory stimulation (Givens, 1996;Makeig et al., 2002;Jansen et al., 2003;Ross et al., 2005), and transcranial magnetic stimulation (Van der Werf and Paus, 2006). Below we show that direct electrical (somatic) CR stimulation as well as indirect inhibitory or indirect excitatory CR stimulation can have intriguingly similar effects. The consequences of our results are two-fold. On the one hand the local effects (realized by direct stimulation of the neurons' somata or stimulation of afferent axon terminals) and the non-local effects (realized by stimulation of efferent fibers) of deep brain CR stimulation may be similar or even identical. On the other hand a basic prerequisite for sensory CR stimulation is fulfilled. Hence, the DBS-oriented concept of CR-induced desynchronization (Tass, 2003a,b) and anti-kindling (Tass and Majtanik, 2006) may be extended to a more general CR stimulation concept which can be realized by invasive as well as non-invasive, sensory (e.g., acoustic) stimulation. An effective counteraction of pathological neuronal synchronization by acoustic CR stimulation has recently been reported by Tass et al. (2012).

HODGKIN-HUXLEY SPIKING NEURONS
We consider an ensemble of spiking HH neurons (Hodgkin and Huxley, 1952) [see (Hansel et al., 1993) for model equations and parameter values used here], The constant currents I i are uniformly distributed random numbers: S i (t) denotes the internal synaptic coupling within the neuronal population, and F i (t) defines an external stimulation. The neurons interact via excitatory and inhibitory chemical synapses by means of the PSP s i which is triggered by a spike of neuron i (Gerstner and Kistler, 2002;Izhikevich, 2007) and modeled in the standard way by an additional equation for s i (t) (Golomb and Rinzel, 1993;Terman et al., 2002). The coupling term reads and contains a weighted ensemble average of all PSPs received by neuron i from other neurons in the ensemble. Parameter c ij defines the coupling strength from neuron j to neuron i, and V r,j is the reversal potential taken as V r,j = 20 mV for excitatory coupling and V r,j = −40 mV for inhibitory coupling. The type and strength of neuronal interactions is defined by the spatial profile of coupling considered in the form of a Mexican hat , where σ 1 = 3.5, σ 2 = 2, d ij = d|i − j| is the distance between neurons i and j, and d = 10/(N − 1). To avoid boundary effects, we assume that the neuronal fibers are uniformly distributed within the neuronal population such that the distance |i − j| is replaced by N − |i − j|, if the indices i and j get more than N/2 apart from each other. The positive (negative) values of M ij indicate an excitatory (inhibitory) interaction, which corresponds to closely (distantly) located neurons. Depending on the sign of M ij , the corresponding reversal potential V r,j is taken in coupling (2). In such a way local excitatory connections and distant inhibitory connections are modeled (Dominguez et al., 2006;de la Rocha et al., 2008).
The normalized sum of the PSPs s j (t) in Equation (2) represents the collective synaptic activity of the neuronal population, which yields the local field potential LFP(t) = N −1 N j=1 s j (t) (Buzsaki, 2004;Kelly et al., 2010). The LFP will be used to detect synchronization. Low-amplitude LFP oscillations ( Figure 1B, red curve) are characteristic for a desynchronized regime, for example, for a coupling-and stimulation-free regime [S i = F i = 0 in Equation (1)], where neurons independently fire at different times ( Figure 1B, black dots). Because of different constant currents I i , the individual spiking frequencies {ω i } (average number of spikes per second) are broadly distributed around the mean ω ≈ 70.7 Hz with standard deviation σ ω ≈ 0.6 ( Figure 1B, right plot). For sufficiently strong synaptic coupling the neurons synchronize, their frequencies {ω i } get narrowly distributed around the mean ω ≈ 71.4 Hz with standard deviation σ ω ≈ 0.002 ( Figure 1C, right plot). The neurons fire nearly simultaneously ( Figure 1C, black dots), which results in large-amplitude oscillations of the LFP (Figure 1C, red curve).
The extent of phase synchronization in the neuronal ensemble can also be revealed with the order parameter R(t) defined as (Kuramoto, 1984) where (t) is the mean phase. The quantity ϕ j (t) approximates the phase of neuron j and linearly increases over 2π from one spike to the next (Pikovsky et al., 2001) where t j,k , k = 0, 1, 2, . . ., are the spiking times of neuron j. The synchronization order parameter R(t) is ranging from 1 to 0, which corresponds to perfect in-phase synchronization (all neurons fire at the same time) or its complete absence, respectively. For example, for the desynchronized regime illustrated in Figure 1B, the time-averaged order parameter R(t) ≈ 0.06, whereas R(t) ≈ 0.85 for the synchronized dynamics illustrated in Figure 1C.

FITZHUGH-RINZEL BURSTING NEURONS
We consider an ensemble of N interacting bursting neurons, each neuron being modeled by the FHR system (Rinzel, 1987;Izhikevich, 2001), As before, V i , i = 1, 2, . . . , N is the membrane potential of neuron i, S i (t) denotes the internal synaptic coupling considered in the form (2) with the reversal potential V r,j = 2 (V r,j = −2) for excitatory (inhibitory) coupling, F i (t) defines the external stimulation, and s i is the PSP triggered by a burst of neuron i. The constant currents I i are uniformly distributed random numbers: Consequently, without synaptic interactions and external stimulation, S i = F i = 0 in Equation (5), the individual bursting frequencies ω i (average number of bursts per second) of the FHR neurons (5) are broadly distributed around the mean value ω ≈ 4.87 Hz with standard deviation σ ω ≈ 0.21 ( Figure 2B, right plot). The individual neurons burst independently of each other (Figure 2A and Figure 2B, black dots), which results in low-amplitude fluctuations of the LFP (Figure 2B, red curve). In contrast, strong enough coupling synchronizes the neurons, such that the bursting frequencies get narrowly spread, see Figure 2C (right plot), where {ω i } are distributed around the mean value ω ≈ 3.74 Hz with standard deviation σ ω ≈ 0.08. The neurons fire nearly simultaneously ( Figure 2C, black dots), which results in large-amplitude oscillations of the LFP ( Figure 2C, red curve). The corresponding values of the time-averaged order parameter are R(t) ≈ 0.075 in the desynchronized regime and R(t) ≈ 0.95 in the synchronized regime, respectively.

COORDINATED RESET STIMULATION
The considered neuronal networks (1) and (5) are stimulated with CR stimulation (Tass, 2003a,b). According to the CR stimulation algorithm the stimulation signals are sequentially delivered to M different sub-populations of the neuronal target ensemble via M different stimulation sites ( Figure 3A). We choose M = 4 [see (Tass, 2003a,b)]. The stimulation signals are administered via stimulation sites in a time-coordinated manner such that the next stimulation site is activated with a delay of T s /M after the activation of the preceding stimulation site. The stimulation period T s is optimally chosen close to the mean period of the synchronized neurons. Within one cycle of length T s each stimulation site is activated once.
In this paper we consider two types of CR stimulation, namely, (i) direct electrical stimulation of the neurons' somata and (ii) stimulation of excitatory or inhibitory axons terminals or fibers. The latter indirect stimulation models either synaptically-mediated local or non-local effects of electrical deep brain CR stimulation or the effects of sensory CR stimulation (see Introduction). For direct electrical stimulation brief trains of HF charge-balanced electrical pulses are administered to the neuronal tissue ( Figure 3B). Electrical CR stimulation requires an electrode implantation into the neuronal target population, as, for example, in the case of DBS (Benabid et al., 1991). In fact, CR stimulation has initially been developed for the application to electrical DBS (Tass, 2003a,b). To model direct electrical stimulation, we assume

Frontiers in Human Neuroscience
www.frontiersin.org that the stimulation signals ( Figure 3B) are directly delivered to the somata of the stimulated neurons such that the stimulation term F i (t) in Equations (1) and (5) reads K is the stimulation strength, D(i, j) is the spatial profile of the stimulation strength, decaying with distance between the stimulation site j and neuron i ( Figure 3A). ρ j (t) is an indicator function controlling the activation of the stimulation site j, and P(t) is a HF pulse train ( Figure 3B). We also assume that M = 4 stimulation sites are uniformly placed within the neuronal population ( Figure 3A). For example, for a population of N = 200 neurons the stimulation sites are located at lattice coordinates x j = 25, 75, 125, 175. We consider a quadratic spatial decay profile of the stimulation current, which is characteristic for brain tissue (Richardson et al., 2003). Parameters d = 10/(N − 1) and σ d = 0.8. The HF pulse train P(t) in Equation (6) comprises charge-balanced pulses, where a short positive pulse of unit amplitude and length T p is followed by a longer compensatory negative part of length T n with the inter-pulse interval T p + T n , where t k are times of pulse onsets t k = k(T p + T n ), k = 0, 1, . . .. For example, below we use T p = 1 ms and T n = 9 ms with interpulse interval 10 ms, i.e., the pulse-train frequency is 100 Hz for the electrical stimulation of FHR neurons. For the electrical stimulation of HH neurons we use electrical pulses with parameters T p = 0.4 ms and T n = 1.6 ms, where the pulse-train frequency is again of about one order of magnitude larger than the intrinsic frequency of synchronized neurons. The indicator function ρ j (t) in Equation (6) equals 1 if the jth stimulation site is active at time t and zero otherwise. We apply CR stimulation, where the sequence of stimulation site activation is randomly varied between stimulation cycles of length T s with equal probability. To exploit the transient desynchronization induced by CR stimuli (Tass, 2003a,b), we use a patterned ON-OFF CR stimulation with n = 3 cycles ON and m = 2 cycles OFF protocol (with cycle duration T s ), where the stimulation is delivered only during the ON cycles, and desynchronization is strongest during the OFF cycles, see (Lysyansky et al., 2011).
As a possible alternative to direct electrical CR stimulation, we consider indirect, synaptically-mediated, e.g., sensory CR stimulation. Note, the indirect stimulation also models the stimulation of excitatory or inhibitory axons or fibers induced by electrical deep brain CR stimulation (see above and Introduction). We assume that the stimulation signals arrive at the neuronal target population as PSPs evoked in the dendrites of the stimulated neurons e.g., by external sensory stimulation. We model the stimulation-evoked PSPs by the α-function of the form.
with parameter α = 6M/T s , which are shown in Figure 3C, compare with PSPs evoked by spikes of HH neurons (1) (Figure 1A, black curve) and bursts of FHR neurons (5) (Figure 2A, black curve). Through the paper we use T s = 16 ms (α = 3/2) and T s = 320 ms (α = 3/40) for the sensory stimulation of the HH and FHR neurons, respectively. {t k } are the onset times of the α-stimuli, which correspond to the activation times of the stimulation sites given by the indicator function ρ j (t), and the inter- (1) and (5) reads: As before, K is the stimulation strength, D(i, j) is the spatial profile of the stimulation strength (7). This type of the spatial stimulation profile is characteristic for electrical stimulation of brain tissue (Richardson et al., 2003) and also fits to the shape of the tuning curves in sensory processing (Patuzzi and Robertson, 1988;Gulick et al., 1989;Robertson and Irvine, 1989). P(t) is a periodic α-train of the form (9), see also . The indirect stimulation can be excitatory or inhibitory, and we model these two types of stimulation by considering the reverse potential in Equation (10) V r = 20 mV for HH neurons and V r = 2 mV for FHR neurons to achieve excitatory stimulation effect, or V r = −40 mV for HH neurons and V r = −2 mV for FHR neurons in the case of inhibitory stimulation effect.

CROSS-TRIAL PHASE RESETTING ANALYSIS
The effect of CR stimulation is based on the phase reset of the stimulated neurons (Tass, 2003a,b). In order to reveal whether the considered direct (somatic) (6) or indirect (synapticallymediated) (10) stimulation can effectively reset the neuronal populations (1) and (5), we perform a cross-trial phase resetting analysis, see (Tass, 2003c,d;Krachkovskyi et al., 2006) for further details. For this, a long series of identical brief (single shot) stimuli is delivered to the neuronal population, and we evaluate the impact of the stimuli statistically. For this, all neurons in the ensemble simultaneously receive the stimuli via all four stimulation sites at the same times t k , k = 0, 1, 2, . . . , i.e., the indicator functions ρ j in Equations (6) and (10) are simultaneously set to 1 or 0. As a single direct stimulus we use two pulses of the form (8), with opposite polarity to guarantee charge-balancing. In contrast, a single indirect stimulus is modeled by an α-pulse of the form (9). The stimulus onset times t k are randomized with a relatively long inter-stimulus interval t k+1 − t k = T r + ξ, so that the neurons can relax to their natural dynamics during the poststimulus transient. The parameter T r is fixed T r = 100 ms for the HH ensemble (1) and T r = 2000 ms for FHR ensemble (5), and ξ is a random number being uniformly distributed in the interval Frontiers in Human Neuroscience www.frontiersin.org ξ ∈ [0, 50] ms for the HH ensemble (1) and ξ ∈ [0, 1000] ms for the FHR ensemble (5). Around each stimulus a time window of length T r , is considered, W = [−0.25T r , 0.75T r ], where the stimulus onset time t k is set to zero within this window. The dynamics of the mean phase (3) is analyzed across stimulation trials. To detect stereotyped stimulus responses of the phase dynamics, we calculate the cross-trial (across index k) distribution of the mean phase [ k (t) (mod 2π)], t ∈ W, and the corresponding resetting index.
where L is the number of administered stimuli. The resetting index E ∈ [0, 1] quantifies the extent of the phase reset such that it is large if the phase distribution has a well-localized structure, e.g., a well-pronounced phase peak. If large values of E(t) are observed during the post-stimulus transient time in contrast to small values of E(t) in the pre-stimulus time interval then this means that the considered stimulation, which is administered at a random phase because of the randomized stimulation times t k , causes a phase reset of the neuronal ensemble to a stereotypical phase. The stimulation thus restarts the neuronal phase dynamics from a well-defined initial phase, which is important for an effective desynchronization by CR stimulation (Tass, 2003a,b). The cross-trial phase resetting analysis detects stereotypical features in an ensemble of stimulus responses. Accordingly, for the phase resetting analysis, STDP is turned off and the synaptic weights are considered c ij = 0.5 for the HH neurons (1), (2), and c ij = 0.2 for the FHR neurons (5), (2) to avoid any long-term effects of the stimulation and to make sure that the stimuli hit the synchronized network in its natural, pre-stimulus dynamics.

SPIKE TIMING-DEPENDENT PLASTICITY
We investigate the neuronal dynamics modeled by the ensembles (1) and (5) and controlled by direct electrical or indirect (synaptically-mediated), e.g., sensory CR stimulation for adaptive synapses, where the synaptic weights c ij are considered to be dynamical variables dependent on the timing of the neuronal firing of the pre-and post-synaptic interacting neurons. We equip the considered models with the STDP (Gerstner et al., 1996;Markram et al., 1997;Feldman, 2000;Wittenberg and Wang, 2006;Caporale and Dan, 2008). For STDP the time difference t ij = t i − t j , i.e., the time delay between the nearest onsets of the spikes or bursts of the pre-synaptic neuron j and the post-synaptic neuron i is crucial. Depending on t ij , the synaptic weight c ij is updated in a point process like manner by the increment δ · c ij for excitatory connections and −δ · c ij for inhibitory connections with δ > 0 according to the rule.
Due to STDP [Equation (12), Figure 4, left plot] the synaptic strength of excitatory synapses is potentiated or depressed depending on whether the post-synaptic firing follows or advances the pre-synaptic firing, respectively (Gerstner et al., 1996;Markram et al., 1997;Debanne et al., 1998;Bi and Poo, 1998). For inhibitory or distant synapses the situation can be opposite (Caporale and Dan, 2008), which we consider in this paper. In order to avoid an unbounded growth, the synaptic weights c ij will be confined to the interval c ij ∈ [0, c max ], where c max = 1 for the HH ensemble (1) and c max = 0.5 for the FHR ensemble (5). The neurons do not have self-connections: c ii = 0. For an average "net" effect of STDP on the synaptic weights, one can calculate the integral c(ε) = 1 2ε ε −ε c(ζ)dζ which indicates how fast the synaptic strength will on average be potentiated or depressed if the relative firing time t is uniformly distributed in the interval t ∈ [−ε, ε], see also (Kepecs et al., 2002;Tass and Majtanik, 2006;Hauptmann and Tass, 2007). As follows from the shape of c(ε) (Figure 4, right plot), the synaptic weight of excitatory synapses is up-regulated if pre-and postsynaptic neuron preferentially fire in synchrony ( t is narrowly distributed, small ε), whereas it is down-regulated if they fire in an uncorrelated manner ( t is broadly distributed, large ε), e.g., due to desynchronizing stimulation [see (Tass and Majtanik, 2006;Hauptmann and Tass, 2007)].
STDP contributes to a stabilization of synchronized and desynchronized states, which naturally leads to the emergence of a multi-stability of synchronized and desynchronized states in the neuronal populations with STDP. In the considered neuronal model of bursting FHR neurons (5), for example, a strong initial coupling among the neurons leads to an onset of synchronization, where the neurons within the ensembles fire nearly simultaneously. As discussed above, such a dynamics results in a potentiation of the synaptic weights because of STDP, which, in turn, further supports the synchronized dynamics and leads to its stabilization. On the other hand, weak coupling results in a desynchronized dynamics and, due to STDP, in a further depression of the synaptic weights. In such a way a desynchronized state is also expected to be stable in the neuronal ensemble with STDP. This is illustrated for FHR neurons (5) in red curve), whereas in the latter regime the LFP oscillations are of low-amplitude ( Figure 5B, green curve) which is characteristic for a desynchronized dynamics. In such a way two coexisting stable strongly coupled synchronized and weakly coupled desynchronized states can be observed in the FHR neuronal ensemble (5) with STDP (12). We note that some other stable states may also coexist. Four coexisting stable states, for instance, can be observed for the HH ensemble (1) with STDP (12) for different initial coupling matrices, which is illustrated in Figure 6.

PHASE RESETTING OF HH AND FHR ENSEMBLES
The phase of the ensembles of spiking HH neurons (1), (2), and bursting FHR neurons (5), (2) can effectively be reset by direct electrical as well as indirect (e.g., sensory) stimulation. To illustrate this, we first consider synchronized bursting FHR neurons (as in Figure 2C) and administer a long series of resetting electrical pulses and sensory stimuli as described in Section. 2.4. For electrical stimulation ( Figure 7A) the distribution density of the mean phases k (t) averaged across stimulation trials is uniform in the pre-stimulus interval (for t < 0). This is due to the randomized stimulation times t k , which also indicates that the stimuli are administered to the neuronal population at random phases. Nevertheless, during the post-stimulus transient (for t > 0) the phase distribution density has a well-pronounced peak which indicates that the stimuli induce a stereotypical response of the collective dynamics of the stimulated ensemble. The direct electrical stimulation can thus induce a phase reset. The same is achieved by indirect excitatory and indirect inhibitory (e.g., sensory) stimulation (Figures 7B,C). The phase resetting index E(t) (Figure 7, left plots, white curves) clearly reflects the stimulationinduced phase reset. It attains small values in the pre-stimulus interval, whereas E(t) increases after the stimulation.
The considered direct electrical and indirect (synapticallymediated) inhibitory or excitatory types of stimulation (modeling either sensory stimulation or electrical stimulation of axons terminals or fibers) can also induce a phase reset in the ensemble of HH spiking neurons (1), (2). We again consider a synchronized ensemble (as in Figure 1C) and administer a long series Frontiers in Human Neuroscience www.frontiersin.org of resetting electrical pulses and sensory stimuli as described in Section 2.4. As shown in Figures 7D-F for the electrical and excitatory and inhibitory sensory stimulations, the cross-trial diagrams show two well-pronounced phase peaks during the poststimulus transient [note the difference to the single-peak crosstrail phase distribution for the FHR neurons (Figures 7A-C)].

FIGURE 7 | Cross-trial diagrams of the distribution densities of the mean phase k (t) averaged across stimulation trials for (A)-(C) synchronized FHR neurons (5), (2), and (D)-(F) synchronized HH neurons (1), (2) without STDP.
Note, for the phase resetting analysis, STDP is turned off to avoid any long-term effects of the stimulation and to guarantee that after the post-stimulus transients the network has re-established its natural, pre-stimulus dynamics. In plots ( Figure 1C for HH neurons and in Figure 2C for FHR neurons. This indicates that two stereotypical responses can be induced by the electrical and sensory stimulation in the HH neuronal ensemble. Nevertheless, the time dynamics of the resetting index E(t) supports the stimulation-induced phase reset (Figures 7D-F, white curves). E(t) is small in the pre-stimulus interval (the cross-trials phase distribution is uniform), whereas E(t) attains large values in the post-stimulus interval.

REWIRING OF NEURONAL ENSEMBLES WITH STDP BY CR STIMULATION
The direct electrical as well as the indirect (synaptically-mediated, e.g., sensory) CR stimulation can have a well-pronounced desynchronizing effect on the stimulated neuronal population. Due to STDP CR stimulation does not only influence the collective dynamics, but also the pattern of the synaptic couplings among the neurons, i.e., it may induce a rewiring of the stimulated neuronal ensemble: the control of the collective dynamics enables a control of the connectivity. To illustrate the above effects of CR stimulation, we first consider the FHR neurons (5), (2) with STDP (12) in a strongly coupled and synchronized regime, see Figure 5.
Weak direct electrical as well as weak indirect (e.g., sensory) stimulation has no noticeable impact on the mean coupling of the stimulated neuronal ensemble, as expected (Figure 8, blue curves). The stimulated neurons remain strongly coupled and, Frontiers in Human Neuroscience www.frontiersin.org thus, synchronized during as well as after the epoch of CR stimulation. However, for moderate stimulation intensities the time course of the mean coupling C(t) in the stimulated neuronal population is different (Figure 8, green curves): it gradually decays during the stimulation and, finally, approaches a small magnitude during the post-stimulation transient. Accordingly, both direct electrical and indirect (e.g., sensory) CR stimulation shifts the neuronal ensemble from a strongly coupled and synchronized regime to a stable weakly coupled and desynchronized state. These qualitatively very different states of the FHR neuronal ensemble have been identified as stable in our analysis above (see Figure 5). The desynchronizing and rewiring impact of CR stimulation is a robust effect and preserved within a large range of the stimulation strength (Figure 8, red curves). However, for strong stimulation, both direct electrical and indirect (e.g., sensory), the synaptic weights among the stimulated neurons hardly deviate during the stimulation from those of the initially strongly coupled regime. Nevertheless, CR stimulation can still induce a long-lasting desynchronizing and rewiring effect where the coupling weights approach small values during the post-stimulation transient (Figure 8, red curves), and the stimulation-free neuronal ensemble eventually stabilizes in a weakly coupled and desynchronized regime observed in Figure 5.
Analogous arguments are applicable to the network of HH neurons (1), (2) with STDP (12). Indeed, weak direct electrical as well as indirect (e.g., sensory) CR stimulation do not induce any long-lasting change in the coupling (Figure 9, blue curves), whereas a stronger CR stimulation significantly suppresses the mean coupling C(t) during the stimulation and leads to a longlasting rewiring and desynchronization (Figure 9, green curves), see also Figure 6. If the stimulation gets too strong, it can still significantly suppress the synaptic weights during the stimulation, but it does not cause any long-lasting effect (Figure 9, red curves). The neuronal population returns to the strongly coupled and synchronized regime (as in the pre-stimulation epoch) after the stimulation is switched off. The latter effect differs from that of the FHR neurons (Figure 8). However, also for the FHR ensemble, too strong CR stimulation fails to induce a long-lasting rewiring and desynchronization.
For both considered neuronal models with STDP there exists a broad range of the stimulation strength where the corresponding direct electrical or indirect (e.g., sensory) CR stimulation is effective, i.e., where it induces long-lasting rewiring and desynchronization. For illustration, we register the mean coupling C(t) in the stimulated ensemble during the stimulation period C on , i.e., in the stimulation ON condition, and then at the end of the poststimulus transient C off , i.e., in the stimulation OFF condition. These two quantities are depicted for FHR neurons in Figure 10 versus the stimulation intensity K by green dashed and green solid curves, respectively. As follows, CR stimulation is effective for K ∈ (K 1 , K 2 ) (Figure 10, vertical dashed lines) where the neuronal ensemble converges during the post-stimulation transient to a weakly coupled and desynchronized regime (Figure 10, solid green and red curves). As CR stimulation gets stronger, it induces a weaker suppression of the coupling during stimulation, i.e., the mean synaptic weight C on increases with K (Figure 10, green dashed curve). However, the long-lasting effect of CR stimulation is still preserved until C on reaches a critical value and CR stimulation becomes ineffective for K > K 2 (Figure 10, solid green and red curves). For the HH neurons with STDP there also exists an interval of the stimulation intensity K ∈ (K 1 , K 2 ), where the direct electrical as well as the indirect (e.g., sensory) CR stimulation is effective, i.e., it shifts the neuronal population to a weakly-coupled and desynchronized regime (Figure 11, solid green and red curves). The difference can be observed for large stimulation strength K, where CR stimulation still induces a significant rewiring of the neuronal ensemble only during the stimulation, i.e., C on is small (Figure 11, green dashed curves), while it does not lead to longlasting changes in coupling and collective neuronal dynamics for K > K 2 .

MIXED STIMULATION TYPE AND RANDOM NETWORK TOPOLOGY
This section serves to illustrate that the effects of CR stimulation can also be observed for a qualitatively different network topology as well as for mixed excitatory-inhibitory stimulation. To illustrate that the observed long-lasting effects of CR stimulation do not depend on the particular type of the network topology, we deliver CR stimulation to a network comprising FHR model neurons as studied above [Equations (2), (5)], but with random connections, i.e. without any particular underlying topological structure (Figure 12). For this, we randomly and uniformly redistribute the parameters M ij of the coupling profile in the interval [−p, 1 − p], where 0 ≤ p ≤ 1 defines the fraction of the inhibitory coupling (M ij < 0) in the neuronal population. For the sake of illustration, we consider a FHR neuronal network (2), (5) with p = 0.2 which corresponds to 20% of inhibitory connections. Similar results can also be obtained for smaller p, e.g., for purely excitatory coupling where p = 0. Such a network with STDP also exhibits a multistability of synchronized and desynchronized states where, depending on the initial conditions, a stable synchronized state characterized by large-amplitude oscillations of LFP ( Figure 12B, red curve) can be observed. Being captured in that dynamical attractor state, the considered STDP rule (12) results in a potentiation of the excitatory connections and a depression of the inhibitory connections ( Figure 12C, left plot). CR stimulation can effectively desynchronize such an ensemble with random topology. We illustrate this for an indirect (synaptically-mediated) stimulation, where a fraction of the population of randomly selected neurons receives an inhibitory stimulation, and the rest of the ensemble is stimulated excitatorily. Put otherwise, the CR stimuli affect a randomly selected portion of the neurons in an inhibitory manner and the rest in an excitatory manner. Irrespectively of the ratio of neurons being stimulated excitatorily vs. inhibitorily, CR stimulation shifts the neuronal ensemble from a strongly coupled and synchronized regime to a weakly coupled and desynchronized regime ( Figure 12A). The latter is characterized by low-amplitude oscillations of the LFP ( Figure 12B, green curve) and depressed excitatory and potentiated inhibitory connections ( Figure 12C, right plot). This state is stable and preserved if the stimulation is switched off, which constitutes a long-lasting desynchronizing and rewiring effect of CR stimulation.

DISCUSSION
The mechanism of action of CR stimulation essentially relies on the ability to perform a phase reset of the oscillatory dynamics of neuronal populations (Tass, 2003a,b). To achieve the desynchronizing effect of CR stimulation, phase resetting stimuli are sequentially delivered at different sites, i.e., to different neuronal subpopulations. Since the subpopulations are reset at different times, their phases get shifted with respect to each other (Tass, 2003a,b). In this computational paper we show that a phase reset of neuronal populations can robustly and in an intriguingly similar way be achieved with (i) direct electrical stimulation, (ii) indirect (synaptically-mediated) excitatory stimulation, or (iii) indirect (synaptically-mediated) inhibitory stimulation (Figure 7). Furthermore, we show that the CR-induced Frontiers in Human Neuroscience www.frontiersin.org desynchronization and anti-kindling (unlearning of pathological synchrony and connectivity) can be achieved no matter which stimulation modality-(i), (ii), or (iii)-we use. For all three stimulation modalities the CR mechanism robustly works within a broad range of the stimulation strength. We illustrate these effects on two very different network models with spiking and bursting neurons, respectively.
We note that the observed long-lasting effects of CR stimulation do not depend on the particular type of the coupling topology. In fact, CR stimulation robustly works also in a network with random topology (Figure 12).
An exclusively excitatory fiber tract was selected as target for standard HF DBS in a case with spinocerebellar ataxia type 2 and severe tremor (Freund et al., 2007). According to a computational study, the beneficial effect of this particular target might be a consequence of the pure excitatory stimulation (Hauptmann and Tass, 2007). However, typically both excitatory and inhibitory structures are affected by stimulation. In general, pure excitatory as well as pure inhibitory stimulation are limiting cases, where in clinical applications mixed excitatory-inhibitory stimulation is more realistic. On the single neuron level the phase reset is a fundamental and robust dynamical mechanism, which can be achieved by hyperpolarizing or depolarizing electrical pulses (Best, 1979;Demir et al., 1997;Neiman et al., 2007;Tateno and Robinson, 2007) as well as by excitatory or inhibitory PSPs (Perkel et al., 1964;Pinsker, 1977;Lerma and Garcia-Austt, 1985;Jackson et al., 2002;Prinz et al., 2003). Accordingly, on a population level a timely coordinated phase reset in different sub-populations can be robustly achieved, no matter which portion of randomly selected neurons undergoes a phase reset via excitatory or inhibitory stimulation (Figure 12). We have chosen a mixed stimulation set-up, where a portion (ranging from 0 − 100%) of randomly selected neurons receives inhibitory stimulation. Remarkably, CR robustly causes an anti-kindling over the entire inhibitory range, i.e. from 0 to 100% (Figure 12).
Our results have four major consequences: 1. Local effects of HF DBS (in the vicinity of the stimulation electrode) and non-local effects of HF DBS caused by an activation of fibers may differ considerably (McIntyre et al., 2004b). The local effects may be inhibitory due to direct electrical effects of HF DBS on voltage-gated currents leading to a depolarization blockade (Beurrier et al., 2001) or to low-amplitude subthreshold oscillations of the membrane potential (Pyragas et al., 2012). In addition, inhibitory local effects may arise from an excitation of axon terminals connected to neurons by inhibitory synapses (synaptic inhibition) (Dostrovsky et al., 2000). In contrast, the non-local effects may be excitatory due to a HF axonal output via an activation of excitatory projection neurons (McIntyre et al., 2004a). In contrast, our results indicate that the local effects (mediated by direct stimulation of somata or stimulation of inhibitory axon terminals) and the non-local effects (mediated by stimulation of excitatory fibers) of deep brain CR stimulation may be the same, so that the beneficial effect, the desynchronization, is achieved both locally (in the vicinity of the electrode) and non-locally (via projecting fibers). 2. From a theoretical standpoint our results indicate that the initially DBS oriented concept of direct electrical brain CR stimulation may be extended to a more general concept of CR neuromodulation comprising both invasive and non-invasive (e.g., sensory) applications. In fact, acoustic CR neuromodulation has already been applied to counteract chronic subjective tinnitus and the underlying pathological neuronal synchronization . Subjective tinnitus is a percept of sound without an objectively identifiable sound source, which often emerges as consequence of a hearing loss (Lockwood et al., 2002;Eggermont, 2003;Moller, 2003;Weisz et al., 2005). Pathologically increased delta band activity is observed in cortical regions deprived of afferent input (Steriade, 2006). Correspondingly, pathologically enhanced neuronal synchronization was found, e.g., in auditory brain areas Frontiers in Human Neuroscience www.frontiersin.org of tinnitus patients (Llinas et al., 1999;Weisz et al., 2005Weisz et al., , 2007. To perform an acoustic CR neuromodulation for the treatment of chronic subjective tinnitus, the tonotopic organization of the central auditory system was used, and electrical stimulation bursts applied to different brain sites were replaced by acoustically delivered tones of different pitch . This type of treatment turned out to significantly counteract both tinnitus symptoms and the underlying pathological neuronal synchronization . In a prospective, randomized, single blind, placebo-controlled trial in 63 patients with chronic tonal tinnitus it was shown that CR treatment was safe and well-tolerated and resulted in a highly significant decrease of tinnitus loudness and symptoms as measured by VAS and TQ scores . Furthermore, as shown by means of EEG recordings, acoustic CR stimulation significantly decreased the tinnitus-related abnormal neuronal synchronization processes . 3. For strong stimulation in our models we observed a nontrivial relationship between the effects during CR stimulation and after turning off CR stimulation: in the FHR network direct electrical CR stimulation and indirect (e.g., sensory) CR stimulation can lead to a long-lasting rewiring and desynchronization without any noticeable acute effect (Figures 8  and 10), whereas for the HH network a pronounced desynchronization during CR stimulation may not guarantee any long-lasting desynchronization (Figures 9 and 11). In fact, in our models we found an optimal range of the stimulation intensity, within which CR stimulation effectively induces long-lasting desynchronization and anti-kindling. Previously, we reported on such an optimal parameter range already for simpler models of coupled phase oscillators and FitzHugh-Nagumo spiking neurons both without STDP (Lysyansky et al., 2011). In this paper we generalize our findings to more complicated and more realistic neuronal models and coupling topologies comprising STDP as well as for the three stimulation modalities (being direct electrical, indirect excitatory, and indirect inhibitory). Our findings indicate that in the case of strong stimulation desynchronizing effects achieved during CR stimulation need not be predictive of long-lasting desynchronizing effects, persisting after stimulus offset. 4. For our models we may compare the robustness of the direct electrical CR stimulation to that of the indirect (e.g., sensory) excitatory and indirect (e.g., sensory) inhibitory CR stimulation. For this we compare the relative lengths of the parameter K intervals K R = (K 2 − K 1 )/K 1 related to effective CR stimulation, see Figures 10 and 11. We found that K R ≈ 32, K R ≈ 29, and K R ≈ 27 for the direct electrical, indirect (e.g., sensory) excitatory, and indirect (e.g., sensory) inhibitory CR stimulation for the HH neurons, respectively (Figure 11). For the FHR neurons the corresponding values read K R ≈ 38, K R ≈ 635, and K R ≈ 624 (Figure 10). We conclude that in our models the indirect (synaptically-mediated, e.g., sensory) CR stimulation is at least as robust as the direct electrical CR stimulation. In some cases, however, it may be even much more robust than the direct electrical CR stimulation concerning an induction of long-lasting desynchronizing effects. These results indicate that it might be a fruitful approach to counteract pathological neuronal synchronization by desynchronizing non-invasive, sensory CR neuromodulation [see ].