In this section we introduce a network of coupled phase oscillators as a simple model of a neuronal ensemble subjected to CR stimulation. We consider the well-known Kuramoto system of N all-to-all (or globally) coupled phase oscillators, where each oscillator receives the coupling from the other N − 1 oscillators, and which reflects several main synchronization properties of a number of oscillatory networks (Kuramoto, 1984; Acebrón et al., 2005), (1)θ˙j=ωj+CN∑k=1Nsin(θk−θj)+Sj(t),   j=1, 2, …, N, where θ_j(t) are the phases, ω_j are the natural frequencies, C is the coupling strength in the ensemble, and S_j(t) is the stimulation signal which will be defined below. We consider N = 400 oscillators and coupling C = 0.1. The natural frequencies ω_j are Gaussian distributed with mean ω_mean = π and standard deviation σ_ω = 0.02. For numerical integration we use a Runge–Kutta method of order 5(4) with adaptive step size (Hairer et al., 1993).

Equation (1) governs the time-dependent dynamics of the oscillator phases θ_j(t), t ≥ 0, which, in the absence of stimulation (S_j(t) ≡ 0), have been found to spontaneously synchronize for sufficiently strong coupling (Kuramoto, 1984; Strogatz, 2000). In that case, a large group of oscillators starts to oscillate at the same frequency, and their phases get narrowly distributed by forming a single phase cluster, in this way constituting an in-phase synchronization. The onset of in-phase synchronization in ensemble (Equation 1) is reflected by an increasing amplitude of the mean field, which is given by the first order parameter R₁ defined as Rm=| N−1∑j=1Nexp(imθj) | for m = 1 (Haken, 1983; Kuramoto, 1984; Tass, 1999). Accordingly, large values of R₁, 0 ≤ R₁ ≤ 1, close to 1 correspond to in-phase synchronization of oscillators (Equation 1), whereas a desynchronized state, where the phases θ_j are uniformly distributed on the unit circle, is indicated by small values of all order parameters R_m, m ≥ 1. For example, the time-averaged first order parameter 〈R₁〉 ≈ 0.98 for the considered coupling strength and distribution of the natural frequencies. The above order parameters will thus be used below to characterize the extent of synchronization in ensemble (Equation 1). The order parameters R_m can also reflect the formation of symmetric cluster states in ensemble (Equation 1), where the phases θ_j (mod 2π) are divided into a few synchronized groups (clusters) of equal size being equidistantly spaced on a circle of length 2π. Such a cluster state comprising l clusters is characterized by large values of R_l combined with small values of R_m, 1 ≤ m < l (see, for example, the values of R₁ and R₄ in Figure 2 below for a CR-induced 4-cluster state).

We note here that since model (Equation 1) is dimensionless, all parameters are also considered dimensionless and will, thus, be measured and illustrated in figures below in arbitrary units. We are, however, able to compare the obtained results with realistic setups and provide a mapping between the parameter space of the model and experimentally measured quantities (see section 4).

CR stimulation is delivered to the phase ensemble (Equation 1) in the following setup: The oscillators are assumed to be equidistantly arranged on a 1D lattice of length L = 10 with lattice coordinates x_j = (j − 1)L/(N − 1), j=1, N¯. The stimulation signals are sequentially delivered via N_s stimulation sites, which are equidistantly spaced within the stimulated ensembles with coordinates ck=(k−12)L/Ns, k=1, Ns¯. The control of system (Equation 1) by CR stimulation of strength I is modeled by the following stimulation term (Tass, 2003a,b) in Equation (1): (2)Sj(t)=I∑k=1NsD(xj, k)ρk(t)P(t)cosθj, where P(t) is a HF pulse train of unit amplitude, and ρ_k(t) are the indicator functions such that ρ_k(t) = 1 if the kth stimulation site is active and ρ_k(t) = 0 otherwise. Functions D(x_j, k) define the spatial spread of the stimulation signals in the target population such that the impact of the stimulation signal delivered via the kth stimulating site on the jth oscillator depends on the distance |x_j − c_k| between the oscillator and the stimulation site. Following Richardson et al. (2003) we consider a quadratic spatial profile of the current spread (3)D(xj, k)=11+(xj−ck)2/σ2, where σ defines the spatial decay rate of the stimulation current (Figure 1A).

The CR stimulation signals are short HF pulse trains (i.e., bursts with an intra-burst frequency as for standard HF DBS) sequentially delivered via different stimulation sites (Tass, 2003a,b; Tass et al., 2012b). To define the stimulation signals in our models, we consider a long sequence of rectangular pulses of unit amplitude P(t) with the pulse period T_p = 0.025 (i.e., 40 pulses per time unit) and pulse width T_p/2 = 0.0125. Then the stimulation signal delivered to the neurons via the kth stimulation site during its active period reads ρ_k(t)P(t), where ρ_k(t) controls the switching on and off of the kth stimulation site as defined above. During each stimulation cycle of length T, all N_s stimulation sites are sequentially activated delivering an HF pulse train of length T/N_s, respectively (Figure 1B). For the phase oscillators (Equation 1) the length of the CR stimulation cycle is considered T = 2 which equals the mean period of the stimulation-free synchronized phase ensemble.

In the previous section 2.1 we introduced the phase ensemble (Equation 1) controlled by CR stimulation. For applications, it is important to reveal whether the results for this model presented in sections 3.1 and 3.2 are robust and generic enough to be valid for more realistic neuronal models. Here we present a network of FHN (FitzHugh, 1961; Nagumo et al., 1962) spiking neurons interacting via excitatory chemical synapses, (4) v˙j=vj−13vj3−wj+1+Ij(syn)+Ij(stim),w˙j=εj(vj+0.7−0.8wj), s˙j=2(1−sj)1+exp(−10vj)−sj,   j=1, 2, …, N.

The variable v_j models the membrane potential of a single neuron, w_j is a recovery variable, and the parameter ε_j determines the natural spiking frequency of a neuron (number of spikes per time unit). The synaptic coupling among the neurons is realized via post-synaptic potentials s_j triggered by spikes of neuron j (Terman, 2005; Izhikevich, 2007), which are modeled by additional differential equations for s_j, see Terman et al. (2002) and Terman (2005). Then the synaptic current I^(syn)_j in Equation (4) reads (5)Ij(syn)=C(V−vj)1N∑k=1Nsk, where V is a reversal potential, taken as V = 2 to realize an excitatory coupling, and C defines the coupling strength. We consider N = 400 neurons, coupling strength C = 0.11, and parameters ε_j are Gaussian distributed with mean 0.08 and standard deviation 0.002.

In order to quantify the extent of synchronization in this network, we use the order parameters R_m, m ≥ 1, defined for the phase ensemble, see section 2.1, where phase θ_j of neuron j is approximated based on its spike timings t_{j, k}: θ_j(t) = 2π(t − t_{j, k})/(t_{j, k + 1} − t_{j, k}) + 2πk, t_{j, k} ≤ t < t_{j, k + 1}, k = 0, 1, 2,… (Pikovsky et al., 2001). For example, the time-averaged first order parameter 〈R₁〉 ≈ 0.96 for the considered coupling strength and distribution of parameters ε_j.

The stimulation current I^(stim)_j in Equation (4) of CR stimulation is given by the following expression: (6)Ij(stim)=I∑k=1NsD(xj, k)ρk(t)P(t), where the functions D, ρ_k, and P are the same, as defined in section 2.2 for the ensemble of phase oscillators. We use the stimulation period T = 38 of CR stimulation, which approximates the period of the intrinsic spiking dynamics of the neuronal ensemble (Equations 4, 5). The period of the high-frequency pulse train P(t) is taken T_p = 0.5, and the pulse width T_p/2 = 0.25.

Synchronized firing of bursting neurons has been reported in basal ganglia of MPTP-treated monkeys (Bergman et al., 1994, 1998). In this section we present a model network of coupled adaptive exponential integrate-and-fire (aEIF) bursting neurons (Brette and Gerstner, 2005; Naud et al., 2008; Touboul and Brette, 2008). The model is given by the following equations for the membrane potentials V_j and adaptation currents w_j: (7)CV˙j=−gL(Vj−EL)+gLΔTexp(Vj−VTΔT)−wj+Ij(syn)+Ij(stim)+Ij,  w˙j=(a(Vj−EL)−wj)/τw,   j=1, 2, …, N.

The jth oscillator emits a spike if its membrane potential V_j overcomes some threshold, here we set it to −25 (mV). At this moment the variables (V_j, w_j) are instantaneously reset to the values:  Vj→Vr,wj→wj+b.

The parameters of the model are g_L = 30 nS, E_L = −70.6 mV, V_T = −50.4 mV, Δ_T = 2 mV, τ_w = 40 ms, a = 4 nS, b = 80 pA, V_r = −47.2 mV, and C = 281 pF. The values I_j are randomly chosen according to a Gaussian distribution with mean 780 and σ_I = 1.0 pA. These parameters imply a bursting mode of the oscillators (Touboul and Brette, 2008) with period ≈ 70 ms. The synaptic coupling I^(syn)_j is given by Ij(syn)=K(Vr.p.−Vj)1N∑k=1Nα(t−tk(sp)),α(x)=4x exp(−4x), where K is the coupling strength, t^(sp)_k is the last spike of the kth oscillator, and V_r.p. = −20 mV is a reversal potential. We consider K = 12 nS such that the time-averaged first order parameter 〈R₁〉 ≈ 0.92 calculated as for the FHN ensemble except that the timings t_{j, k} stand for the onsets of bursts. The stimulation current I^(stim)_j of CR stimulation is considered in the form (Equation 6) with parameter T_p = 2 ms. The stimulation period T = 70 ms, and the number of neurons in the ensemble is taken N = 200.

The desynchronizing impact of CR stimulation on the controlled oscillators depends on the stimulation parameters mentioned above, see Lysyansky et al. (2011) where, in particular, the role of the stimulation strength I and the current decay rate σ was investigated. In the present work the main attention is paid to the impact of the number of stimulating sites N_s on desynchronizing effect. CR stimulation can be administered either continuously, where the stimulation cycles described above of length T are applied one after another (without interruption), or in an intermittent way, where a few stimulation cycles are repeatedly followed by a few cycles without stimulation. In this section we consider the former stimulation protocol.

Examples of the time courses of the order parameters R₁ and R₄ before, during, and after 4-site continuous CR stimulation are illustrated in Figure 2A. CR stimulation is administered to an ensemble of strongly synchronized oscillators (Equation 1) (switched on at t = 400) where, because of the strong enough coupling, the phases are narrowly distributed (Figure 2C), and the first order parameter attains a large value R₁ ≈ 0.98 in the pre-stimulus interval (Figure 2A for t < 400). After a short transient the stimulation suppresses the in-phase synchronization marked by small values of R₁ (Figure 2A, black curve for 400 < t < 700). We found that during the stimulation the order parameters slightly fluctuate around their mean values 〈R₁〉 ≈ 0.07, 〈R₂〉 ≈ 0.13, 〈R₃〉 ≈ 0.17, and 〈R₄〉 ≈ 0.55 (Figure 2A, only R₁ and R₄ are depicted by black and red curves, respectively). Therefore, continuous CR stimulation administered via 4 stimulation sites can reliably induce a 4-cluster state such that the phases are grouped into four clusters equidistantly spaced on the unit circle (Figure 2D). The emergence of the cluster state is reflected by large values of R₄ combined with small values of R₁, R₂, and R₃, see section 2.1 and Lysyansky et al. (2011). After the stimulation is switched off (at t = 700) the oscillatory ensemble transiently relaxes from the stimulation-induced cluster state to a desynchronized state where the phases get nearly uniformly distributed (Figure 2E), and which is characterized by small values of both order parameters R₁ and R₄ (Figure 2A). Then the desynchronization is followed by a resynchronization due to the persistent strong coupling. The phases of the stimulation-free oscillators form a single cluster (Figure 2F) and approach the original strongly synchronized state (Figure 2C), and the order parameters increase (Figure 2A for t > 800). The discussed cluster state induced by the continuous CR stimulation is a robust phenomenon (Lysyansky et al., 2011). For example, for the considered 4-site CR stimulation there exists a range of the stimulation strength I where the time-averaged first and fourth order parameters 〈R₁〉 and 〈R₄〉 attain small and large values, respectively, which is characteristic for a four-cluster state, see Figure 2B. The minimal value of 〈R₁〉 obtained for the optimal stimulation strength I_opt is depicted by the blue circle in Figure 2B.

As mentioned in the Introduction, in this study we estimate the desynchronizing impact of CR stimulation for different numbers of stimulation sites N_s. For this the stimulation strength I is varied in some interval [0, I_max] whereas the other parameters remain fixed. The minimal (optimal) value 〈R₁〉_opt of the time-averaged first order parameter 〈R₁〉 achieved in the considered interval of I (see Figure 2B) is then used to quantify the optimal effect of CR stimulation for a given number of stimulation sites N_s. As N_s increases, the optimal values 〈R₁〉_opt may demonstrate different behaviors depending on the spatial decay rate σ of the stimulation current, see Figure 3. 〈R₁〉_opt saturates for large N_s, and additional stimulation sites do not further improve the optimal effect of CR stimulation. The saturation levels of 〈R₁〉_opt decisively depend on the values of σ, where small (large) σ implies small (large) values of 〈R₁〉_opt. Moreover, small σ (≲0.5) ensures that the increasing number of stimulation sites N_s up to ≈5–10 sites improves the stimulation outcome as compared to the stimulation delivered via N_s ≈2–3 sites. In contrast, large σ (≳1.25) implies that N_s = 2 sites is the most preferable choice, and larger N_s leads to larger values of 〈R₁〉_opt and to a suboptimal desynchronizing effect of continuous CR stimulation. For intermediate values of σ (≈0.75–1.00) 〈R₁〉_opt does not significantly change with varying number of stimulation sites.

In summary, if the properties of the neuronal tissue imply a spatially selective stimulation profile, corresponding to a small σ, additional stimulation sites can improve the desynchronizing effect of CR stimulation. In contrast, a small number of stimulation sites is the best choice for a broad profile of the spatial current decay (i.e., for large σ).

The suggested protocol of CR stimulation (Tass, 2003a,b) implies an intermittent (ON–OFF) administration of the stimulation, where several cycles of CR stimulation (ON-cycles) are followed by a few rest cycles without stimulation (OFF-cycles). Such an intermittent stimulation allows the oscillators to freely evolve during the OFF-cycles according to their intrinsic unperturbed dynamics. During the post-stimulus transient (OFF cycles) the oscillatory ensemble (Equation 1) relaxes into a desynchronized state followed by resynchronization, both phenomena (transient desynchronization and subsequent resynchronization) being a consequence of the persisting strong coupling in the network (Tass, 2003a,b). This effect is illustrated in Figure 2A, where the post-stimulus desynchronization is reflected by small values of the order parameters. Intermittent stimulation can also decrease the amount of current delivered to the target tissue. This, in turn, can minimize the spread of the stimulation to the neighboring areas and, hence, prevent from possible side effects. In networks with spike timing-dependent plasticity (STDP) (Markram et al., 1997; Bi and Poo, 1998) intermittent CR stimulation causes an anti-kindling effect (Tass and Majtanik, 2006; Hauptmann and Tass, 2007). Anti-kindling is a desynchronization-induced unlearning of pathologically upregulated synaptic connectivity and, in turn, synchrony, i.e., a stabilization of a desynchronized activity in a network, which persists even if the stimulation is completely switched off (Tass and Majtanik, 2006; Hauptmann and Tass, 2009). The time scale of STDP is essentially larger than the m : n timing of the intermittent stimulation and has an order of tens or hundreds periods of the oscillation (Bi and Poo, 1998; Hauptmann and Tass, 2009). In this paper we investigate the impact of CR stimulation on oscillatory networks without STDP.

To quantify the desynchronizing effect of intermittent CR stimulation, we use an approach suggested in Lysyansky et al. (2011). In each rest interval ℐ_k we evaluate the maximal value of the first order parameter such that a set of the maximal values rk=maxt∈ℐkR1(t), k=1, 2, …, k = 1, 2,… is calculated. For example, during the rest intervals the first order parameter R₁(t) can increase and reach its maximal values r_k at the end of the OFF cycles (Figure 4B, red circles). The mean value of r_k, i.e., 〈r〉=Nrp−1∑k=1Nrprk, where N_rp is the number of rest periods, is then used to estimate the effect of the intermittent ON–OFF CR stimulation. In what follows we denote by m and n the number of the ON and OFF cycles of length T of the interchanging stimulation and rest time intervals, respectively, where T is the stimulation period (Figure 4A). If m and n can attain non-integer values, an interesting anti-resonance-like behavior of 〈r〉 vs. m and n can be observed (Lysyansky et al., 2011). In the present paper we consider integer values of m and n only and analyze the impact of the number of stimulation sites N_s on the optimal value of 〈r〉. We also vary the stimulation strength I as well as the stimulation timing, i.e., the values of m and n in order to find optimal parameters for the intermittent CR stimulation.

In Figure 5A we illustrate the dynamics of 〈r〉 vs. the stimulation intensity I for fixed number m = 3 of the ON cycles and different numbers n of the OFF cycles. As in the case of the continuous stimulation protocol (Figure 2B), there exists an optimal stimulation strength I_opt = I_opt(m, n), where 〈r〉 attains a minimal value 〈r〉_opt (Figure 5A, red circles). Increasing n results in larger values of 〈r〉 and 〈r〉_opt since longer rest intervals imply more time for resynchronization and, hence, larger maximal values of R₁ occurring during the rest intervals. We follow the optimal values 〈r〉_opt as the number of stimulation sites N_s increases and find that 〈r〉_opt saturates if N_s gets large enough (Figure 5B). This phenomenon is tightly related to the saturation observed for the continuous CR stimulation (Figure 3).

Since we are interested in the longest possible post-stimulus desynchronized transient as discussed above, we are looking for the stimulation parameters leading to the maximal admissible length n_max of the rest intervals. For this we first investigate the dynamics of 〈r〉_opt induced by the intermittent CR stimulation vs. the length of the rest intervals n for different numbers of stimulation sites N_s, see Figure 6. As expected (see also Figure 5A), 〈r〉_opt = 〈r〉_opt(n) increases as the rest intervals get longer. At some n = n_max + 1, 〈r〉_opt exceeds the predefined threshold 〈r〉_opt = 0.5 (Figure 6, dashed line). We thus consider the value n = n_max (Figure 6, empty circles) as the maximal admissible length of the rest intervals for the stimulation strength I ∈ [0, I_max]. n_max can thus serve as a criterion for an optimal parameter choice for the intermittent CR stimulation since large n_max guarantees long stimulation-free OFF intervals. We found that an increasing number of stimulation sites N_s can prolong the admissible rest periods, i.e., n_max increases, and thus improve the desynchronizing impact of the intermittent CR stimulation if the latter is spatially rather selective, i.e., the stimulation current is narrowly distributed around the stimulation site within the stimulated population (Figure 6A for σ = 0.5). For a broad spatial spread of the stimulation current the situation, however, is exactly opposite, i.e., n_max decreases for larger N_s (Figure 6B for σ = 2.0) such that the longest rest periods can be achieved for the smallest number of stimulation sites N_s = 2.

The minimal values 〈r〉_opt of the first order parameter illustrated in Figure 6 are attained at the optimal stimulation strengths I = I_opt ∈ [0, I_max] (Figure 5A) shown in Figure 7 vs. the number n of the OFF cycles in the rest intervals. Accordingly, one has to increase the stimulation strength in order to reach the possible minimal value of the order parameter 〈r〉_opt during the longer post-stimulus rest periods if the number of stimulation sites N_s is fixed. Nevertheless, the optimal intermittent CR stimulation with the maximal admissible length n_max of the rest periods (Figure 6) can be realized for smaller optimal stimulation strength I_opt if the number of stimulation sites is increased (Figure 7, empty circles). Note, this effect is well pronounced for small σ (Figure 7A) in contrast to the case of large σ (Figure 7B).

We have shown above that the desynchronizing impact of the continuous as well as the intermittent CR stimulation essentially depends on the spatial decay rate of the stimulation current σ as the number of stimulation sites varies (Figures 3, 6). We summarize our findings in Figure 8. If the number of stimulation sites N_s increases, the maximal length n_max of the admissible rest periods can either increase or decrease depending on the values of σ, see Figure 8A. The spatially selective stimulation with small σ (Figure 8A, black curve for σ = 0.5) implies relatively short rest intervals for the stimulation delivered via small number of sites. Additional stimulation sites strongly elongate the admissible rest intervals, which, however, saturate for large N_s. In contrast, the intermittent CR stimulation with a broad spatial spread of the stimulation signal (large σ) admits the longest rest periods for a small number of stimulation sites (Figure 8A, red curve for σ = 2). Adding further stimulation sites can only worsen the stimulation outcome in this case.

As follows from Figures 5A, 7 the optimal value 〈r〉_opt is attained for larger optimal stimulation strength I_opt if the number n of the OFF cycles increases. It is thus important to clarify how much stimulation an individual oscillator receives for the optimal intermittent CR stimulation as in Figure 8A. We therefore calculate the effective amount of the administered stimulation Ieff=0.5Imm+n1Ns×N∑k=1Ns∑j=1ND(xj, k) which is the time and ensemble average of the stimulation signal (Equation 2). In Figure 8B the values of I_eff are plotted vs. the number of stimulation sites N_s for the maximal admissible length of the rest interval n_max from Figure 8A. For small σ (spatially selective stimulation) the desynchronizing effect of the intermittent CR stimulation is significantly improved by increasing N_s, where the length of the admissible rest periods is increased (Figure 8A, black curve), and the amount of required stimulation is decreased (Figure 8B, black curve). For the large σ the situation is opposite, see Figures 8A,B (red curves). However, due to the saturation effect, the improvement of the stimulation impact for the considered small σ = 0.5 becomes less effective if the number of stimulation sites gets larger than N_s ≈ 5 − 7 sites (Figures 8A,B, black curves). It is thus unreasonable to strongly increase the number of stimulation sites, which can either induce no significant improvements of the desynchronizing impact of CR stimulation, or even worsen it if the spatial stimulation profile is broad enough.

Since the best intermittent CR stimulation is leading to the longest OFF periods achieved for a small stimulation current, we introduce a measure Q of the quality of the intermittent CR stimulation, which is proportional to the maximal admissible length of the rest intervals given by n_max and inversely proportional to the effective amount of the stimulation I_eff and the number of stimulation sites N_s, Q = n_max/(I_eff × N_s). This function estimates the efficacy of CR stimulation “per each stimulation site” such that for similar values of n_max/I_eff more preferable is the stimulation setup with a smaller number of sites. We found that for fixed spatial decay rate of the stimulation signal from a range σ ∈ [0.25, 2.5], Q achieves a clear maximum Q_max = Q_max(σ) for an optimal numbers of stimulation sites N_{s, opt} depicted by a black circle in Figure 9, respectively. Ranges of N_s are indicated by colored stripes where the quality Q of the stimulation does not deviate for more than 30% from its best value Q_max. In contrast, the white regions correspond to parameter values of clearly sub-optimal CR stimulation and should be avoided. This diagram supports our conclusion that an optimal CR stimulation has to be administered via a small number of sites if the stimulation current is broadly spread in the tissue, i.e., if σ is large, whereas for small σ (spatially selective stimulation) a larger number of stimulation sites leads to a better desynchronizing impact of CR stimulation.

For the FHN model (Equations 4, 5) controlled by CR stimulation (Equation 6) we use the same evaluation techniques applied to the phase oscillators (Equation 1) and described in section 3.1. We thus calculate the optimal (minimal) values 〈R₁〉_opt of the time-averaged first order parameter 〈R₁〉 obtained under variation of the stimulation strength I ∈ [0, 10]. The effect of continuous CR stimulation on the FHN network is illustrated in Figure 10 vs. the number of stimulation sites N_s. Different curves in the plot represent the stimulation-induced minimal values of the first order parameter 〈R₁〉_opt for different values of the spatial decay rate σ. This diagram is very similar to that in Figure 3 obtained for the phase oscillators. Indeed, 〈R₁〉_opt saturates for large number of stimulation sites. The saturation levels of 〈R₁〉_opt depend on the values of σ as for the phase oscillators, such that for small (large) σ CR stimulation leads to a small (large) order parameter which saturates for large N_s.

Also for the intermittent ON–OFF CR stimulation of the FHN ensemble the impact of the stimulation can be quantified in the same way as described in section 3.2 for the network of phase oscillators. In particular, Figure 11A illustrates the influence of the number of stimulation sites N_s on the maximal admissible length n_max of the rest intervals. As for the phase oscillators (Figure 8A), a spatially selective (i.e., with small σ) intermittent CR stimulation is less effective for small number of stimulation sites N_s (n_max is small), whereas the length of the rest intervals can significantly be prolonged if N_s increases (Figure 11A, black curve for σ = 0.5). Simultaneously, the amount of the stimulation current I_eff received on average by a single neuron in the network decays with increasing N_s (Figure 11B, black curve). Therefore, a selective intermittent CR stimulation optimally delivered to a neuronal network via a large number of stimulation sites enables long stimulation-free desynchronization time intervals with a minimal amount of the administered stimulation current. In contrast, for a broad spatial spread of the stimulation current in the neuronal tissue (large σ) the intermittent CR stimulation has to be optimally administered via a small number on stimulation sites, e.g., N_s = 2. Only for such a stimulation protocol CR stimulation can lead to the longest possible rest periods achieved for the smallest possible delivered stimulation current (Figure 11, red curves for σ = 2). In this case, a larger number of stimulation sites can worsen the desired desynchronizing impact of CR stimulation.

Note, that for the FHN neuronal ensemble (Equations 4, 5) controlled by the intermittent CR stimulation (Equation 6) we have observed the same behavior of the order parameter 〈r〉_opt and the optimal stimulation strength I_opt vs. the number of the OFF cycles n in the rest intervals (figures are not shown) as for the phase ensemble, see Figures 6, 7. Together with the results illustrated in the above Figures 10, 11 this indicates the robustness of the reported phenomena, which can equally be found for phase oscillators and for spiking neurons interacting via excitatory synapses.

When continuous CR stimulation is administered to synchronized aEIF bursting neurons (Equation 7), the in-phase synchronization is replaced by a cluster state (Figures 12A,B). The bursting neuronal discharges get organized in a few sub-groups (clusters) corresponding to the number of stimulation sites, and the neurons within the same cluster fire simultaneously with equidistant time shift between neighboring clusters. The same dynamics is observed for the phase oscillators (Figure 2) and FHN neurons. The synchronization order parameter of the aEIF neuronal ensemble is suppressed from 〈R₁〉 ≈ 0.92 in the in-phase synchronized regime to 〈R₁〉 ≈ 0.014 in the cluster state for the parameter values used in Figures 12A,B. In the latter regime the other order parameters are 〈R₂〉 ≈ 0.063, 〈R₃〉 ≈ 0.088, and 〈R₄〉 ≈ 0.766, which is a clear indication of the 4-cluster state that can be seen in Figures 12A,B.

The spatial decay rate σ of the stimulation current and the number N_s of stimulation sites can significantly influence the desynchronizing impact of permanent CR stimulation as measured by the values of the first order parameter. In particular, for large σ the time-averaged first order parameter 〈R₁〉 grows as N_s increases, whereas the opposite is true for small σ, see Figure 12C. This effect revealed for the considered bursting neurons is again in perfect agreement with the results obtained for the phase oscillators (Figure 3) and FHN spiking neurons (Figure 10).

Also the intermittent ON–OFF CR stimulation of the aEIF neuronal ensemble can become more effective for a large number of stimulation sites. However, this requires spatially selective stimulation, i.e., small σ (Figure 13A, black diamonds). In this case the maximally admissible number of OFF intervals n_max increases with the number of stimulation sites, which permits longer time intervals with desynchronized dynamics of the stimulation-free neurons. We illustrate this effect on the synchronized dynamics for σ = 0.5 and different numbers of stimulation sites in Figures 13B,C. The suboptimal number of stimulation sites N_s = 2 (Figure 13C) implies shorter OFF-intervals and a faster re-increase of both order parameter R₁ and amplitude of the ensemble mean field 〈V〉 during the OFF-intervals in comparison to the stimulation delivered via a larger number N_s = 5 of stimulation sites (Figure 13B). The situation is opposite for a broad spatial spread of the stimulation current in the neuronal tissue (large σ), where the intermittent CR stimulation can optimally be administered via a small number of stimulation sites only (Figure 13A, red circles). In this case, a larger number of stimulation sites can significantly shorten the admissible length of the rest intervals. The above effects of the intermittent CR stimulation for the aEIF bursting neurons (Equation 7) are again in good agreement with those observed for the phase oscillators (Figure 8) and FHN spiking neurons (Figure 11), which further confirms the generality of the presented results.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.