How stimulation frequency and intensity impact on the long-lasting effects of coordinated reset stimulation

Several brain diseases are characterized by abnormally strong neuronal synchrony. Coordinated Reset (CR) stimulation was computationally designed to specifically counteract abnormal neuronal synchronization processes by desynchronization. In the presence of spike-timing-dependent plasticity (STDP) this may lead to a decrease of synaptic excitatory weights and ultimately to an anti-kindling, i.e. unlearning of abnormal synaptic connectivity and abnormal neuronal synchrony. The long-lasting desynchronizing impact of CR stimulation has been verified in pre-clinical and clinical proof of concept studies. However, as yet it is unclear how to optimally choose the CR stimulation frequency, i.e. the repetition rate at which the CR stimuli are delivered. This work presents the first computational study on the dependence of the acute and long-term outcome on the CR stimulation frequency in neuronal networks with STDP. For this purpose, CR stimulation was applied with Rapidly Varying Sequences (RVS) as well as with Slowly Varying Sequences (SVS) in a wide range of stimulation frequencies and intensities. Our findings demonstrate that acute desynchronization, achieved during stimulation, does not necessarily lead to long-term desynchronization after cessation of stimulation. By comparing the long-term effects of the two different CR protocols, the RVS CR stimulation turned out to be more robust against variations of the stimulation frequency. However, SVS CR stimulation can obtain stronger anti-kindling effects. We revealed specific parameter ranges that are favorable for long-term desynchronization. For instance, RVS CR stimulation at weak intensities and with stimulation frequencies in the range of the neuronal firing rates turned out to be effective and robust, in particular, if no closed loop adaptation of stimulation parameters is (technically) available. From a clinical standpoint, this may be relevant in the context of both invasive as well as non-invasive CR stimulation.

A number of computational studies were dedicated on desynchronizing synchronized ensembles or networks of oscillators or neurons [57][58][59][60][61][62]. The clinical need for stimulation techniques that cause desynchronization irrespective of the network's initial state [63], thereby being reasonably robust against variations of system parameters and, hence, not requiring time-consuming calibration, motivated the design of Coordinated Reset (CR) stimulation [64,65]. CR stimuli aim at disrupting in-phase synchronized neuronal populations by delivering phase resetting stimuli typically equidistantly in time, separated by time differences T s /N s , where T s is the duration of a stimulation cycle, and N s is the number of active stimulation sites [64,65]. The spatiotemporal sequence by which all stimulation sites are activated exactly once in a CR stimulation cycle is called the stimulation site sequence, or briefly sequence. Taking into account STDP [27][28][29][30] in oscillatory neural networks qualitatively changed the scope of the desynchronization approach: Computationally, it was shown that CR stimulation reduces the rate of coincident firing and, mediated by STDP, also decreases the average synaptic weight, ultimately preventing the network from generating abnormally increased synchrony [33]. This anti-kindling, i.e., unlearning of abnormally strong synaptic connectivity and of excessive neuronal synchrony, causes long-lasting sustained effects that persist cessation of stimulation [33, 43-45, 66, 67]. As shown computationally, anti-kindling can robustly be achieved in networks with plastic excitatory and inhibitory synapses, no matter whether CR stimulation is administered directly to the soma or through synapses [45,68]. In line with these computational findings, long-lasting CR-induced desynchronization and/or therapeutic effects were accomplished with invasive as well as non-invasive stimulation modalities. Longlasting desynchronization was induced by electrical CR stimulation in rat hippocampal slices rendered epileptic by magnesium withdrawal [69]. Electrical CR deep brain stimulation (DBS) caused long-lasting therapeutic after-effects in parkinsonian non-human primates [70,71]. Bilateral therapeutic after-effects for at least 30 days were caused by unilateral CR stimulation delivered to the subthalamic nucleus (STN) of parkinsonian MPTP monkeys for only 2 h per day during 5 consecutive days [70]. In contrast, standard permanent high-frequency deep brain stimulation did not induce any sustained after-effects [70], see also [72]. In patients with Parkinson's disease electrical CR-DBS delivered to the STN caused a significant and cumulative reduction of abnormal beta band oscillations along with a significant improvement of motor function [73]. Non-invasive, acoustic CR stimulation was developed for the treatment of patients suffering from chronic subjective tinnitus [68,74]. In a proof of concept-study acoustic CR stimulation caused a statistically and clinically significant sustained reduction of tinnitus symptoms [74][75][76] together with a concomitant decrease of abnormal neuronal synchrony [74,77], abnormal effective connectivity [78] as well as abnormal cross-frequency coupling [79] within a tinnitus-related network of brain areas.
So far, the pre-clinical [74,80] and clinical [70,73] proof of concept studies for invasive and non-invasive CR stimulation were driven by computationally derived hypotheses and predictions. Theoretically predicted phenomena and mechanisms, such as long-lasting stimulation effects [33,43,45,66], cumulative stimulation effects [67], and improvement by weak stimulus intensity [81] were verified based on dedicated theory-driven study protocols for pre-clinical and clinical proof of concepts [70,73,74,80].
We here set out to investigate the impact of the CR stimulation frequency and intensity on the effects during stimulus delivery (so-called acute effects), on transient effects emerging directly after cessation of stimulation (so-called acute after-effects), and on effects outlasting cessation of stimulation (so-called sustained after-effects). The ultimate goal of this study is to improve the calibration of CR stimulation, in particular, by providing computationally generated predictions that can be tested in subsequent pre-clinical and clinical studies. The computational study presented here is organized around three hypotheses: Hypothesis #1: Due to the inherently periodic structure of CR stimulation the relation between CR stimulation frequency and the spontaneous neuronal firing rates (prior to stimulation) matters. Periodic delivery of CR stimuli with fixed sequence basically constitutes a time-shifted entrainment of the different neuronal subpopulations [64,65]. A particular closed loop embodiment of CR stimulation, periodic stimulation with demand-controlled length of highfrequency pulse train, is basically a time-shifted entrainment of the different neuronal subpopulations with stimulus intensities adapted to the amount of undesired synchrony [64,65]. Accordingly, the duration of a stimulation cycle was selected to be reasonably close to the mean period of the synchronized neuronal oscillation [64,65]. In STDP-free networks of Kuramoto [82] and FitzHugh-Nagumo [83] model neurons the impact of CR stimulation intensity and frequency on the desynchronizing outcome of CR was studied in detail [81].
Hypothesis #2: Different embodiments of CR stimulation may differ with respect to effect size and robustness. In a series of computational studies [45,46,64,65,81,84,85] and in all pre-clinical [69,70] and clinical studies [73][74][75][76][77][78] performed so far, CR was applied either with fixed sequences or rapidly varying sequences (RVS), where the sequence was randomly varied from cycle to cycle. In a recent computational study, it was shown that at intermediate stimulation intensities the CRinduced anti-kindling effect may significantly be improved by CR with slowly varying sequences (SVS), i.e. by appropriate repetition of the sequence with occasional random switching to the next sequence [84]. However, this study was not performed for a larger range of CR stimulation frequencies. By definition, SVS CR stimulation features significantly more periodicity of the stimulus pattern. Accordingly, the dependence of resonance and/or anti-resonance effects on the CR stimulation frequency might be more pronounced for SVS CR as opposed to RVS CR.
Hypothesis #3: Pronounced acute effects might provide a necessary, but not sufficient condition for pronounced sustained after-effects. In a pre-clinical study in Parkinsonian monkeys with CR-DBS delivered at an optimal and a less favorable intensity, it was shown that long and pronounced acute therapeutic after-effects coincide with long-lasting, sustained after-effects [74]. However, according to computational studies the relationship between acute after-effects and sustained long-lasting effects might be more involved, at least for particular parameter combinations [84].
Related to these hypotheses, to assess the robustness of CR stimulation against initial network conditions we performed our numerical simulations for different network initializations, respectively. In this study we did not systematically vary the stimulation duration. Rather, based on a pre-series of numerical simulations, we here used a fixed stimulation duration that is reasonably short, but nevertheless enabled to robustly achieve an anti-kindling for properly selected values of stimulation frequency and intensity. In fact, our goal was to find stimulation parameters enabling short, but notwithstanding effective CR stimulation. Keeping the stimulation duration at moderate levels may be beneficial for applying the CR approach to different invasive as well as non-invasive stimulation modalities. For instance, standard DBS, i.e. permanent electrical high-frequency pulse train stimulation delivered to dedicated target areas through implanted depth electrodes, used for the treatment of, e.g. Parkinson's disease [86][87][88] may cause side effects. If side effects are caused by stimulation of non-target tissue, they may be reduced by adapting the spatial extent of the current spread to the target's anatomical borders by appropriate electrode designs as introduced, e.g. by [89][90][91], in particular, to spatially tailor stimuli by means of directional DBS [92][93][94][95][96][97]. However, some side effects may at least partly be caused by stimulating the target region itself [98,99]. Accordingly, no matter how precisely stimuli are delivered to DBS targets, the amount of stimulation should be decreased as much as possible. Another example refers to non-invasive applications of CR. In general, non-invasive CR stimulation requires the patients' compliance to actually pursue treatment prescriptions. Obviously, patients might prefer shorter treatment sessions.
To come up with favorable combinations of stimulation parameters, in our numerical analysis we used different data analysis methods, e.g. macroscopic measures assessing the average amount of the population's synchrony and synaptic connectivity. These measures are appropriate to demonstrate relevant stimulation effects, such as stimulation-induced transitions from pronounced neuronal synchrony to desynchronized states.
In summary, in this paper we first explain the computational model and analysis methods. We then apply RVS CR stimulation in a wide parameter range of stimulation frequencies and intensities. We repeat the same analysis for SVS CR stimulation and investigate the differential characteristics of RVS CR and SVS CR with respect to efficacy and robustness. Finally, we analyze the relationship between stimulation-induced acute effects and after-effects. Our results provide the foundation for the development of novel control techniques that will be the topic of a forthcoming study.

The Hodgkin-Huxley Spiking Neuron Model
In this study we use the conductance-based Hodgkin-Huxley neuron model [100] for the description of an ensemble of spiking neurons. The set of equations and parameters are (see also [101,102]): The variable V i ,with i = 1,. . .,N, describes the membrane potential of neuron i, and: The total number of neurons is N = 200, while g Na = 120 mS/cm 2 , g K = 36 mS/cm 2 , g l = 0.3 mS/ cm 2 are the maximum conductance per unit area for the sodium, potassium and leak currents respectively. The constants V Na = 50 mV, V K = −77 mV and V l = −54.4 mV refer to the sodium, potassium and leak reversal potentials respectively. C is the constant membrane capacitance (C = 1 μF/cm 2 ), I i the constant depolarizing current injected into neuron i, determining the intrinsic firing rate of the uncoupled neurons. For the realization of different initial networks, we used random initial conditions drawn from uniform distributions, i.e. initial values of the neural synaptic weights c ij are picked from a normal distribution N(μ c = 0.5 μA/cm 2 , σ c = 0.01 μA/cm 2 ) (see also [45,84] for more details). Hence, in this setup the neurons are not identical. The S i term refers to the internal synaptic input of the neurons within the network to neuron i, while F i represents the current induced in neuron i by the external CR stimulation.

Network and neuron coupling description
The  [27,103] and modelled using an additional equation (see also [104,105]): Initially we draw s i 2 [0,1] (randomly from a uniform distribution) and then, the coupling term S i from Eq (1) (see [45]) contains a weighted ensemble average of all postsynaptic currents received by neuron i from the other neurons in the network and is described by the following term: where c ij is the synaptic coupling strength from neuron j to neuron i and V r,j is the reversal potential of the synaptic coupling (20 mV for excitatory and -40 mV for inhibitory coupling).
In accordance with previous studies [45,84,85] the inhibitory reversal potential was set to −40 mV. The latter makes neurons' more susceptible to input, e.g. stimuli. We performed the same sets of simulations for a subset of stimulation parameters with a different value of the inhibitory reversal potential, -80 mV as in [106] instead of -40 mV. In this way, we obtained very similar results. There are no neuronal self-connections within the network (c ii = 0 mS/cm 2 ). The term M ij , which describes the spatial profile of coupling between neurons i and j, is given by: It has the form of a Mexican hat [107][108][109] and defines the strength and type of neuronal interaction: strong short-range excitatory (M ij > 0) and weak long-range inhibitory interactions (M ij < 0). Here d ij = d|i − j| is the distance between neurons i and j, while determines the distance on the lattice between two neighboring neurons within the ensemble, d 0 is the length of the neuronal chain (d 0 = 10), σ 1 = 3.5, and σ 2 = 2.0. In order to limit boundary effects, we consider that the neurons are distributed in such a way that the distance d ij is taken as: d Á min(|i − j|,N − |i -j|) when the i,j > N/2.

Spike-timing-dependent plasticity
We follow the concepts described in [28,29], regarding the synaptic coupling strengths change dependence on the precise timing of pre-and post-synaptic spikes. Hence, we consider all the synaptic weights c ij to be dynamic variables that depend on the time difference (Δt ij ) between the onset of the spikes of the post-synaptic neuron i and pre-synaptic neuron j (denoted by t i and t j respectively). Then the STDP rule for the change of synaptic weights is given by [29,45]: where β 1 = 1, β 2 = 16, γ 1 = 0.12, γ 2 = 0.15, τ = 14 ms and δ = 0.002. According to the value of Δt ij , the synaptic weight c ij is updated in an event-like manner, i.e. we add or subtract an increment δ Á Δc ij for excitatory or inhibitory connections respectively, with learning rate δ > 0 every time a neuron spikes. Furthermore, we restrict the values of c ij on the interval [0,1] mS/ cm 2 for both excitatory and inhibitory synapses, ensuring in this way that their strengthening or weakening remains bounded. The maximal inhibitory synaptic weight c max was set to be 1 in all our stimulations. However, a more detailed investigation about the effect and variation of this value was performed in [84] where when increasing c max of the inhibitory neurons no significant impact was observed regarding (de)synchronization effects accompanied with a lower average network connectivity. The time window of the plasticity is adjusted with respect to the intrinsic firing rate of the neuron population in order to exhibit multistability, as also discussed in [45]. There, different time-windows (via different choices of parameters) were selected for the STDP for two different neuron models, i.e. one with bursting neurons (FitzHugh-Rinzel) and one for spiking neurons (Hodgkin-Huxley). In our simulations, the STDP tends to simply stabilize the ongoing ensemble evolution and does not, by itself, (de-)synchronize the network. The parameters were, in general, chosen such that the ratio Dt ij g 1;2 t is normalized, and the plasticity takes place within a time interval associated with the spiking period of the individual neurons. We analyzed two additional cases for small variation of the plasticity time-window (τ = 12 and τ = 16) and obtained very similar general effects. The selected fixed value τ = 14, used throughout the entire study, also allows us to compare our results with previously published studies.

Coordinated reset stimulation
Coordinated Reset (CR) stimulation was applied to the neuronal ensemble of N spiking Hodgkin-Huxley neurons. This was done sequentially via N s (= 4 in this study) equidistantly spaced stimulation sites [64]: one stimulation site was active during T s /N s , while the other stimulation sites were inactive during that period. After that another stimulation site was active during the next T s /N s period. All N s stimulation sites were stimulated exactly once within one stimulation ON-cycle. Therefore, the duration of each ON-cycle is T s (in ms). The spatiotemporal activation of stimulation sites is represented by the indicator functions ρ k (t) (k {1,. . .,N}): ( The stimulation signals induced single brief excitatory post-synaptic currents. The evoked time-dependent normalized conductances of the postsynaptic membranes are represented by α-functions given in [102]: Here τ stim = T s /(6N s ) denotes the time-to-peak of G stim , and t k is the onset of the k th activation of the stimulation site. Note that the period (or frequency) through the τ stim parameter of the CR stimulation determines not only the random onset timing of each single signal but also its temporal duration. The spatial spread of the induced excitatory postsynaptic currents in the network is defined by a quadratic spatial decay profile (see [102] for more details) given as a function of the difference in index of neuron i and the index x k of the neuron at stimulation site k: with d the lattice distance between two neighboring neurons as defined in Eq (5) and σ d = 0.8 the spatial decay rate of the stimulation current. Thus, the total stimulation current from Eq (1) is expressed by the following equations: where V r = 20 mV denotes the excitatory reversal potential, V i the membrane potential of neuron i, K the stimulation intensity, and ρ, G, D are given by Eqs (7), (8) and (9) respectively. For the RVS CR stimulation, sequences are randomly chosen from a set of N s ! (= 24) different possible sequences during each ON-cycle (an example is shown in Fig 1A for CR stimulation period T s = 10 ms for the first 90 ms of an activated CR period). Each newly drawn sequence is indicated by a different color and lasts exactly one ON cycle. On the other hand, for the SVS-l CR stimulation, one first randomly picks a sequence and repeats it l times before switching to another one, as shown by the example in Fig  Stimulation frequency, intensity and anti-kindling 84]). Depending on the T s value, more (or less) ON-cycles may be administered within a fixed time interval. In this panel, the total time spans up to two completed ON-and OFF cycles (up to~125 ms in this case) and the color changes at each new sequence.

Macroscopic measurements and statistical tools
The synaptic weights, being affected by the STDP and the different intrinsic periods of the neurons, change dynamically in time. One efficient way to measure the strength of the coupling within the neuronal population at time t is given by the following synaptic weight (averaged over the neuron population): where M ij is defined in Eq (4) and sgn is the sign-function. Furthermore, one may additionally measure the degree of the neuronal synchronization within the ensemble, using the order parameter [82,110]: is a linear approximation of the phase of neuron j between its m th and (m + 1) th spikes at spiking times t j,m and t j,m+1 . R(t) is influenced by the synaptic weights, as the latter are time dependent due to the STDP. The order parameter R measures the extent of phase synchronization in the neuronal ensemble and takes values between 0 (absence of in-phase synchronization) and 1 (perfect in-phase synchronization). In our numerical calculations, we estimate C av [see Eq (11)] and R av . The latter quantity is averaged over the last 100 Á T s . Whenever we plot the order parameter versus time, we determine the moving average <R> over a time window of 400 Á T s , because of the presence of strong fluctuations. For the statistical description and analysis of the non-Gaussian distributed C av and R av data (n = 11 samples), we use the median as well as the Inter-Quartile Range (IQR) [111]. The IQR measures the statistical dispersion, namely the width of the middle 50% of the distribution and is represented by the box in a boxplot. It is also used to determine outliers in the data: an outlier falls more than 1.5 times IQR below the 25% quartile or more than 1.5 times IQR above the 75% quartile. Selecting the appropriate sample size is a complex issue (see e.g. http://www.itl.nist.gov/div898/handbook/index.htm), especially when the standard deviation is unknown. Following the steps described in (http://www.itl.nist.gov/div898/handbook/ prc/section2/prc222.htm), we use the formula to get a first rough estimation of the number of measurements (n) to be included in our sample, where α s refers to the risk of rejecting a true hypothesis, and β s is the risk of accepting a false null hypothesis when a particular value of the alternative hypothesis is true, s d the unknown standard deviation, δ s the confidence interval, and x the values from student's t-distribution. Using 11 samples, as minimum sample size, one is able to reach quite small p-values, much smaller than the significance level a s = 0.05.

Simulation description
In this study, for each initial network of N = 200 non-identical-neurons and parameter conditions (or simply "network"), we apply RVS and SVS CR signals (different per network). For each network, the initial conditions for each neuron were randomly drawn from random distributions as given in the Hodgkin-Huxley Spiking Neuron Model subsection. We start the simulation with an equilibration phase, which lasts 2 s. Later on, we evolve the network under the influence of STDP (which will be present until the end of the simulation). We then integrate the network for 60 s with STDP without any external stimulation yet, where a rewiring of the connections takes place, resulting in a strongly synchronized state. Next, we apply CR stimulation for 128 s (resetting the starting time to t = 0 s). During this CR-on period three stimulation ON-cycles (the stimulation is on) alternated with two OFF-cycles (the stimulation is off) as in the example stimulation signals shown in Fig 2. Each ON-and OFF-cycle lasts T s . After 128 s the CR stimulation ceases permanently and we continue the evolution of the CR-off period for extra 128 s. In order to probe and chart the CR stimulation intensity and frequency parameter space, we restrict the CR stimulation intensity to values in the interval (K 2 [0.20,. . .,0.50]). This particular choice is based on our previous experience and numerical studies (see e.g. [45,84]) where it was found that weaker intensities were not able to sufficiently desynchronize the neuron ensemble while larger intensities did not significantly improve (or sometimes even worsen) the outcome of RVS and SVS CR stimulation signals. We then set an initial-central value for the CR stimulation period (that defines the initial/starting frequency) which in principle is selected close to the intrinsic firing rate of the strongly synchronized network. In this case, and before applying the CR stimulation, the intrinsic firing rate of the network is~71 Hz which corresponds to T s % 14 ms. Hence, we begin with the CR stimulation period T 0 = 16 ms The box frames depict the middle 50%, the upper and lower whiskers the first and last 25% respectively while the outliers (black dots) are set as 1.5 times the length of the box (above/below). There is no statistically significant difference between the data sets at which gives an initial stimulation frequency f 0 = 1/T 0 (in a similar manner just like in [45,84] and adjusted to a value close to the intrinsic one). Then we define such a period interval [T smin ; T smax ] in ms (T s : integer) that allows us to create an "approximately" equidistant grid on the frequency space: f stim 2 [25%f 0 ,. . .,175%f 0 ]. This initial T 0 − value is also well studied for different types of signal patterns aiming to optimize the CR effect with the use of different type of CR stimulation sequences (see e.g. [84]). Then, we define the ratio (%) of CR sequence frequency per ON-cycle (f stim ) over the frequency of the reference stimulation frequency (f 0 = 62.5 Hz, T 0 = 16 ms) as r 0 = (f stim /f 0 ) Á 100 and we end up in studying the intensity and frequency-ratio (K,r 0 ) -parameter space. In Table 1, we show the conversion between the stimulation frequency-ratio and period. For comparison reasons, we also give the corresponding ratios r int (%) of CR stimulation frequency per ON-cycle (f stim ) over the frequency of the intrinsic firing rate of the network frequency (f int = 71.4 Hz, T int = 14 ms) without any external stimulation.

Impact of CR stimulation duration and signals on different initial networks
Before presenting the core of our findings, let us first start by discussing how the RVS CR stimulation duration affects the long-lasting anti-kindling of different initial randomly chosen networks. In Fig 2, we show the evolution of the mean synaptic weight C av as a function of time for different total CR-on time durations: t = 64 s (Fig 2A), t = 128 s (Fig 2B), and t = 256 s ( Fig  2C). 128 s is the standard CR-on period used throughout the paper. The CR stimulation intensity is K = 0.20, and the period T s = 10 ms. A general trend appears in this sequence of panels, i.e. the longer the CR stimulation lasts, less spread of the C av regarding the long-lasting antikindling effect is observed after stimulation offset. This is shown in Fig 2D with boxplots. The last box (corresponding to t = 256 s of total CR-on period) has no outliers and shows a more "uniform" long-lasting effect (as shown in Fig 2C) for all 11 network initializations, not only during the CR-on period but also afterwards during the CR-off period. However, there is no statistically significant decrease of the median of the C av from t = 64 s to t = 128 s (right-sided Wilcoxon rank sum test [112], p = 0.0209, 5% significance level). Moreover, the median value of the C av does not change significantly between t = 128 s ( Fig 2B) and t = 256 s (Fig 2C, bothsided Wilcoxon rank sum test, p = 0.8955). Hence, the intermediate stimulation duration t = 128 s provides fairly good results. Furthermore, for considerably larger stimulation durations the anti-kindling is typically, but not always more pronounced. From a clinical standpoint, it is desirable to achieve reasonably pronounced stimulation effects without excessive stimulation durations. Accordingly, in this computational study we choose t = 128 s as total CR-on time, and t = 256 s as total CR-on/off time. For the different simulations, we use different random initial networks and CR signals. For the sake of generality, we do not pick any optimal combination of random initial network and RVS CR stimulation signal that would induce a favorable or biased behavior. This is to assess whether CR effects are robust with respect to different initial conditions. Fig 2B shows a typical example where 11 different random stimulation signals where applied to 11 different initial networks during the CR-on period, with CR stimulation intensity K = 0.20 and stimulation period T s = 10 ms. The CR-on/off period lasts 128 ms respectively. During the CR-on period the mean synaptic weights C av evolve in a similar manner for all networks, with little deviations between the different curves. They reach approximately the same small value at the end of the CR-on period. The latter corresponds to weak excitatory synaptic connectivity and, in most cases in this paper, to globally well-desynchronized states. However, the post-stimulation dynamics of C av may be quite diverse. Some networks retain their weak average connectivity while others, like network 2 and 9 ( Fig 2B) relapse back to states with strong synaptic connectivity. Next, we study what happens if we fix the CR stimulation signal for the 11 different initial networks (Fig 2E). The results are similar to Fig 2B: The outcome at the end of the CR-on period is fairly uniform, while the post-stimulation dynamics of C av is diverse. Replacing one random external stimulation signal by another one may improve the long-term outcome in some cases (e.g. network 8 -green dotted line), but worsen the outcome in others (e.g. network 3 -blue solid line). These plots indicate that both the random initialization of the network and the different stimulation signals during the CR-on period impact on the final outcome at the end of the CR-off period in a complex manner.

Impact of RVS CR stimulation intensity and frequency on acute effects
Next, we investigate how stimulation intensity and stimulation frequency impact on the mean synaptic weight and synchronization at the end of the RVS CR-on period. Fig 3A shows the median of the mean synaptic weight C av , and Fig 3B of the order parameter R av (averaged over the last 100 Á T s ) as a function of stimulation intensity K and stimulation frequency f stim . The color bars show the median values which were calculated from 11 different random initial network configurations. Overall, at the end of the RVS CR-on period we observe a weak excitatory coupling. In other words, CR stimulation shifts the networks' couplings towards more inhibition, the inhibitory couplings get stronger, and desynchronized states emerge for most of the (K,r 0 ) pairs, except for the two columns at f stim = 25%f 0 (T s = 64 ms) and f stim = 145%f 0 (T s = 11 ms). For the former frequency, CR stimulation fails to weaken both the inter-neural connectivity and synchrony, whereas for the latter frequency CR down-regulates synaptic connectivity, but elevated levels of synchrony persist. Fig 3C and 3D show their Inter-Quartile-Range (IQR) respectively, which gives a measure of the data dispersion around these median values.
All IQR values being close to zero indicate that the middle 50% of the distribution are very close to the median value.  Fig 4C and 4D display the corresponding IQRs, showing that the dispersion around the median of the C av results is very small in large parts of the parameter plane.

Impact of RVS CR stimulation intensity and frequency on sustained aftereffects
In contrast, small IQRs are found only for small R av , in regions with strong desynchronization. Fig 4A and 4B display two main bands in the (K,r 0 ) − parameter space associated with small dispersion: The first band is aligned along the horizontal axis, for weak stimulation intensities (i.e. K = 0.20 and K = 0.25) and stimulation frequencies greater than 40% of the standard f 0 corresponding to a stimulation period of T 0 = 16 ms. The second band runs along the vertical stimulation intensity K axis, and for relatively high frequencies, i.e. for f stim = 160%f 0 (T s = 10 ms) and f stim = 175%f 0 (T s = 9 ms) which correspond to~155% and~140% of the firing rate of the synchronized neurons, respectively. For these (bottom-horizontal and right-hand-sidevertical bands) the dispersion around the median values is quite small for both C av and R av (Fig 4C and 4D). In addition, the vertical stripe at the reference frequency value f 0 ("100%", T s = 16 ms), studied in [84], but with a t = 64 s CR-on period, is also associated with robust long-lasting anti-kindling and desynchronization for all CR stimulation intensity values K. Another region with similar characteristics lies at the center of Fig 4A and 4B

for intermediate stimulation intensity and frequency values.
At a first glance, among those two bands in Fig 4A and 4B, where dark color dominates suggesting long-lasting anti-kindling after cessation of CR stimulation, the horizontal band seems especially intriguing. Along the lines of our model analysis, the horizontal band corresponds to pronounced desynchronizing outcome at favorably weak CR stimulation intensities within Stimulation frequency, intensity and anti-kindling a range of stimulation frequencies. However, we have to keep in mind that the discrete grid is not very dense. Hence, in order to investigate whether this conclusion is justified, we calculated C av and R av for all the integer period T s values for K = 0.20, ranging from f stim = 175%f 0 (T s = 10 ms) to f stim = 40%f 0 (T s = 40 ms). Fig 5 shows this fine-grained analysis. The boxplot for C av is shown in Fig 5A, and for the R av in Fig 5B. Note, in this figure the horizontal axis shows the CR stimulation period instead of the frequency. And it is sorted from larger to smaller values for an easier comparison between the two representations. The red and green dots indicate the reference stimulation period T 0 = 16 ms and intrinsic firing rate period T int = 14 ms ms respectively. For T s 2 [9,. . .,24 ms] we observe robust anti-kindling effects. In contrast, for T s 2 [25,. . .,28 ms] many networks tend to be in a synchronized state, while for T s 2 [29,. . .,38 ms] the anti-kindling is found to be robust again, before finally reaching the largest T s value where the CR stimulation signals are not effective at all. In summary, at weak stimulation intensities favorable stimulation outcomes are achieved within wide ranges of the stimulation frequency. For further analyses of stimulation induced effects observed in particular ranges of the stimulation intensity/frequency parameter plane, we refer to the Supporting Information. For particular stimulation parameters, similar acute effects, as assessed with macroscopic quantities R av and C av , may lead to qualitatively different results. Neither prominent features of the connectivity matrix nor the dynamical states of the individually stimulated subpopulations at the end of the CR-on period enabled us to predict the long-term outcome (see Supporting Information). Furthermore, this analysis revealed that CR may be effective without causing side-effects that are time-locked to the individual stimuli (see Supporting Information). Stimulation frequency, intensity and anti-kindling

Impact of SVS CR stimulation intensity and frequency on stimulation effects
Next, we address the robustness of the long-lasting anti-kindling achieved by SVS CR stimulation in the (K,r 0 ) − parameter plane. We use SVS-100 CR stimulation, where the random switching occurs after 100 repetitions of the CR sequence (for motivation see [84]). In Fig 6 we show the total outcome of C av and R av , obtained by delivering SVS-100 CR to the same 11 initial networks as in Figs 3 and 4 and varying the CR stimulation frequency and intensity. Let us compare these results with the results for RVS CR (Fig 3A and 3B). Regarding the medians of the C av , both RVS and SVS CRs (Fig 3A vs Fig 6A) overall the parameter dependence outcomes are similar, where the outcome plots of SVS CR (Fig 6) contain more vertical stripes, associated with greater outcome variability. Let us consider some of the differences between RVS CR and SVS CR: For low intensity (K = 0.20) and high frequencies f stim = 175%f 0 , 160%f 0 (corresponding to T s = 9 ms, 10 ms respectively) SVS-100 does neither cause pronounced acute desynchronizing effects nor sustained long-lasting effects. For low CR frequency 25%f 0 (corresponding to T s = 64 ms, leftmost column) it requires even stronger intensities to induce an anti-kindling compared to RVS (Fig 3A). Regarding the median of R av (Fig 3B vs Fig 6B) for almost all (K,r 0 ) Stimulation frequency, intensity and anti-kindling − parameters the networks are shifted to a desynchronized state at the end of the CR-on period, with only a few exceptions, in particular (K,r 0 ) = (0.20, 175%f 0 ), (0.20, 160%f 0 ) and for 25%f 0 . Moreover, for the frequency f stim = 145%f 0 the SVS CR stimulation achieves more pronounced anti-kindling effects (at the end of the CR-on period) for all intensities K compared to the RVS CR stimulation.
In Fig 6C and 6D we present the outcome for the medians of C av and R av at the end of the CR-off period for SVS-100 CR. The main differences compared to RVS CR (Fig 4A and 4B) are the 'stripes' at f stim = 115%f 0 and, in particular, at f stim = 145%f 0 where SVS-100 neither reduces C av nor R av for all K-values. Moreover, also for the lowest intensity value K = 0.20 and frequencies f stim = 175%f 0 , 160%f 0 , 25%f 0 no anti-kindling is achieved. However, there is a substantial overlap of the (K, f stim ) − parameter range where both RVS and SVS-100 CR lead to long-lasting anti-kindling, mainly for high frequencies f stim = 175%f 0 , 160%f 0 for K ! 0.25 as well as for 40%f 0 ≲ f stim ≲ 100%f 0 and a wide range of K-values. Interestingly, whenever SVS CR stimulation causes an anti-kindling, the long-term effects on the connectivity are particularly robust, irrespective of different network initializations and parameters.
Let us now investigate a denser T s period sample for the weakest intensity K = 0.20, with the same format as in Fig 5 for RVS CR stimulation. Boxplots of C av (Fig 7A) and R av (Fig 7B) at the end of the CR-off period show that SVS CR stimulation at this weak intensity is overall less efficient in inducing long-lasting anti-kindling effects compared to RVS CR (Fig 5). In particular, there is no distinct range of T s periods where SVS CR causes a pronounced anti-kindling. However, for a few values of T s for the long-term outcome for SVS is stronger than for RVS, e.g. for T s = 15 ms, 16 ms. Fig 8 enables us to display the stimulation's global performance in a more concise manner. Namely, it shows the dependence of stimulation outcome on CR stimulation intensity and frequency for the time-averaged mean synaptic weights C av (Fig 8A and 8E) and time-averaged order parameter R av (Fig 8B and 8F), both at the end of the CR-off period, with values belonging to the same intensity value K for RVS (top row) and SVS CR (bottom row) stimulation, respectively. Similar plots but now for values belonging to the same frequency ratio (f stim /f 0 ) Á 100 are shown in Fig 8C and 8G and Fig 8D and 8H respectively.

Discussion
By systematically varying the CR stimulation frequency and intensity and comparing the stimulation outcome of the two different CR protocols, RVS and SVS CR stimulation, RVS CR proved to be more robust with respect to variations of the stimulation frequency. However, in accordance with a previous computational study, restricted to a fixed value of the stimulation frequency [84], SVS CR stimulation can induce stronger anti-kindling effects. In our study, we obtained particular parameter ranges related to particularly favorable stimulation outcome. If no closed loop adaptation for the stimulation frequency is available, RVS CR stimulation at weak intensities and with stimulation frequencies in the range of the neuronal firing rates enables to effectively and robustly achieve an anti-kindling.
To our knowledge, in our study in a plastic network the CR stimulation frequency and intensity were systematically varied for the first time to investigate the impact on acute and long-lasting stimulation outcome. Remarkably, pronounced acute desynchronization (as measured by means of the standard order parameter from Eq (12) [82,110]) does not necessarily lead to long-lasting desynchronization. On the one hand this finding might inspire future computational and pre-clinical studies aiming at specifically designing stimulation protocols for longlasting (as opposed to acute) desynchronization. On the other hand, this finding is significant for the development of clinical calibration procedures for CR stimulation, see [113].
In a previous study in networks without STDP Lysyansky and coworkers [81] considered m:n ON-OFF CR stimulation with real rather than integer m and n and varied m and n systematically. For non-integer m incomplete CR stimulation cycles are delivered, intersected by incomplete pause cycles caused by non-integer n. This type of CR stimulation has not yet been used in pre-clinical or clinical studies and is somewhat remote to the initial CR concept that builds on the periodicity of both neuronal firing and stimulus patterns [64,65]. For the majority of the CR stimulation parameters used in this work, no drastic change was observed in the firing rates. The only exception was observed for very low stimulation frequency combined with comparably high intensities (S2 Fig). Especially for the most relevant cases (weak to intermediate intensities and frequencies around the reference stimulation frequency) the firing rate of the neuron ensemble remains almost unchanged when compared with initial intrinsic firing rates before CR delivery (about up to ±3% variation of the initial intrinsic firing rate). In this study, we focused on a network of spiking Hodgkin-Huxley neurons with STDP. Compared to STDP-free networks used before [81], this is a step towards more complex and, in particular, plastic neural networks. Future studies should address yet more complex neural networks equipped with STDP to study parameter regions and stimulation protocols that are reasonably stable in different neural network models. In principle, we have to be careful about extrapolating findings obtained in one type of neural network model to network models of higher complexity. For instance, non-linear delayed feedback stimulation was introduced in globally coupled networks of limit cycle oscillators and phase oscillators [114]. It turned out to robustly cause desynchronization, nearly irrespective of the selected valued of the delay [115]. In contrast, linear delayed feedback [116] was shown to induce desynchronization only for a rather small subset of parameter pairs of delay and intensity, favoring delays close to half the intrinsic oscillation period and weak to moderate intensities [115,116].
However, in a more complex, microscopic neuronal network model consisting of a population of STN and a population of external globus pallidus (GPe) neurons [105] the parameter dependence for nonlinear delayed feedback was qualitatively different [117]. The parameter ranges of delay and intensity values associated with desynchronization were still greater for nonlinear delayed feedback as opposed to linear delayed feedback. However, in this microscopic STN-GPe network model nonlinear delayed feedback had to be properly calibrated and, in particular, the delay had to be adjusted to the intrinsic period of the neuronal oscillations, to enable desynchronization [117]. Note, the microscopic STN-GPe network model did not contain STDP [105]. Incorporating STDP to a neuronal network model substantially adds to the model's complexity (see e.g. [34]) and might, hence, further impact on the dependence Stimulation frequency, intensity and anti-kindling of the stimulation outcome on key parameters of both linear and nonlinear delayed feedback. Ultimately, we strive for using several neural networks with STDP as testbed for generating computationally based predictions and recommendations for favorable stimulus parameters and dosage protocols. However, different models may display similar or even identical spontaneous (i.e. stimulation-free) dynamics, but may have very different stimulus response properties (see e.g. [62]).
Accordingly, we cannot expect a stimulation technique to be generically effective, irrespective of the neural network model used. Nevertheless, stepwise adding further physiologically and anatomically relevant features to the neural network models employed may help to generate specific predictions and, ultimately, to further improve stimulation protocols and dosage regimes. In that sense, the finding that RVS CR stimulation at weak to moderate intensities and stimulation frequencies adapted to the neurons' intrinsic firing rates causes a desynchronization in neural network models without STDP [81] and with STDP as shown in this study, is relevant and, in fact, in accordance with pre-clinical findings [71,74]. Furthermore, the fact that SVS CR stimulation might even be more effective, but requires more careful parameter adaptation may guide future development of calibration techniques as put forward in a forthcoming study [118].
In neural networks with STDP post-stimulation transients may be complex. For instance, for stimulation dosages just reaching the level required for an anti-kindling, a rebound of excessive synchrony may occur immediately following cessation of CR stimulation, while later on a full-blown, sustained desynchronization emerges [43,67,119]. This rebound selectively relates to synchrony, rather than synaptic connectivity. This phenomenon occurs when after CR delivery the neuronal population just reaches the basin of attraction of a favorable attractor. Upon entering the basin of attraction, the synaptic connectivity is still super-critical, so that synchrony emerges in the absence of stimulation. As the neuronal network relaxes towards the favorable attractor, the initially up-regulated synaptic connectivity fades away until, finally, the synaptic connectivity remains below a critical threshold, hence, preventing the population from getting synchronized [43,67,119]. However, it remains to be shown to which extent the rebound of synchrony phenomenon might be a generic after-effect occurring for just about sufficient CR dosage or simply an epiphenomenon specific to the computational model [120] used in those studies [43,67,119], comprising networks of Morris-Lecar spike generators [120] transformed to burst mode by a slowly varying current [121,122].
We here demonstrated that over a wide range of stimulation parameters favorable acute effects do not automatically lead to favorable long-lasting, sustained after-effects. This is in agreement with a computational study in the same model, but performed in only a restricted parameter range [84], as well as with an EEG experiment performed in tinnitus patients [123]. To characterize stimulation induced effects, we here used the average synaptic weight [Eq (11)] and the average amount of neuronal synchrony [Eq (12)]. These macroscopic quantities enabled us to effectively investigate the impact of variations of stimulation parameters on the stimulation outcome. However, in Supporting Information section, we have shown that pronounced differences of the average synaptic weight do not necessarily lead to pronounced differences of the average amount of synchrony (S1 Fig). Another example in this context is the combination of weak average synaptic connectivity ( Fig 3A) combined with increased levels of average neuronal synchrony (Fig 3B) at the end of the CR-on period. To further study the relationship between connectivity pattern and synchronization processes, macroscopic quantities may not be sufficient to grasp all relevant details of the connectivity matrix and the dynamical features of the resulting synchronization processes.
In a number of previous studies, it was already shown that computational findings obtained in minimal models and, in particular, in models that are even simpler than the model studied in this manuscript, turned out to be of high clinical significance. The following computational predictions were obtained by studying minimal models, such as networks of phase oscillators with and without STDP: 1. Anti-kindling: In networks with STDP, desynchronizing stimulation may induce long-lasting, sustained desynchronizing effects that outlast cessation of stimulation by an unlearning of abnormally strong synaptic connectivity, see e.g. [33]. This fundamental prediction of long-lasting desynchronization and, in turn, therapeutic effects was verified in parkinsonian monkeys [71,74] as well as in Parkinson's patients with CR-DSB [73], in Parkinson's patients with vibrotactile CR stimulation [124], and in tinnitus patients with acoustic CR stimulation [70,77,78].
2. Cumulative effects: As shown computationally, repeated delivery of CR stimulation may have cumulative effects [67]. Cumulative effects were clinically verified in Parkinson's patients with CR-DBS [73], in Parkinson's patients with vibrotactile CR stimulation [124], and in tinnitus patients with acoustic CR stimulation [70].
3. Optimal effects of deep brain CR stimulation at weak intensities: CR-DBS has optimal desynchronizing effects at intensities that correspond to approx. a third of the intensity of standard deep brain stimulation [81]. This theoretical prediction was verified in parkinsonian monkeys [71,74].
The computational predictions obtained in minimal models were actually used to design the pre-clinical and clinical studies referred to above. From a clinical standpoint, these computationally predicted and pre-clinically and/or clinically verified findings are significant and may ultimately enable to establish superior therapies that require stimulus delivery for only a few hours, on a regular or occasional basis.
The Mexican hat connectivity was, e.g. used by [108,109] to study auditory responses, sensorineural hearing loss and tinnitus. These studies illustrate that the comparably simple Mexican hat connectivity model is able to capture relevant connectivity features. Accordingly, later on, the same connectivity profile was used in a several studies focusing on stimulus-induced desynchronization e.g. in the auditory cortex [45,84,85] In neural networks without STDP tested so far, CR stimulation works at higher intensities as well, see e.g. [81]. In that case, pronounced cluster states are induced, but coherent synchrony is reliably suppressed [81]. This is not the case in the neural network model with STDP studied here. For both RVS CR and SVS CR, for several parameters tested the long-term outcome deteriorates with increasing stimulation intensity (Figs 4, 6 and 8). Accordingly, based on our results, in pre-clinical and clinical applications stimulation at higher intensities should be avoided. Another important aspect refers to the more pronounced periodicity of SVS CR pattern. In previous papers (lacking a wider scan of the parameter space), SVS CR stimulation appeared to be superior to RVS CR stimulation [84,85]. In this paper, however, we show that SVS CR stimulation decisively depends on the appropriate choice of the stimulation frequency (Figs 6, 7 and 8). This sensitivity may significantly reduce the performance in the presence of biologically realistic variations of the neuronal firing rates and might, hence, be the very reason, why the outcome of SVS CR stimulation is significantly better for smaller numbers of sequence repetitions [preliminary results presented by Wang et al., Critical parameters determining efficacy of coordinated reset stimulation of subthalamic nucleus and related changes in primary motor cortical and subthalamic local field potentials in a parkinsonian monkey. Society for Neuroscience (2017) Poster 210.02 / I10]. There, based on first pre-clinical data, SVS CR stimulation at high numbers of sequence repetition appears to be inappropriate for an open loop application. However, its performance might be significantly improved by closed loop approaches as, e.g. computationally shown in [118].
For the development of CR stimulation, in a number of computational studies predominantly minimal models were used [33, 43-45, 61-68, 119], as opposed to biophysically realistic models [46,125]. These computational studies gave rise to qualitative non-trivial predictions, e.g. the emergence of long-lasting, sustained [33] as well as cumulative [67] effects and concerning the amplitude of the stimulation amplitude [81]. These predictions were verified in pre-clinical [71,74] and clinical studies [70,73,124]. In fact, the computational findings were used to design the study protocols and generate the underlying hypotheses. However, we have to keep in mind that the minimal-model based approach yields qualitative rather than quantitative predictions. Accordingly, we do not intend to provide "success rates" of the stimulation outcome since we do not intend to relate particular values of C av and R av with successful outcome for the following reasons. On the one hand, the mean synaptic weight C av cannot be assessed in living humans and, hence, so far, no correlation between C av and the extent of symptoms has been studied. On the other hand, the order parameter R av cannot directly be assessed either, but is related to the amplitude of macroscopic/mesoscopic signals like the local field potential. However, for instance in Parkinson's disease it is still a matter of debate whether there is a measurable quantity that reasonably represents the extent of symptoms, in a biomarker-like manner [126,127]. In fact, a number of studies provided results that are in contradiction of the biomarker notion [127][128][129][130][131][132][133]. Accordingly, so far, it is not possible to provide ranges of the amount of synchrony-reflected by R av -that correspond to physiological as opposed to abnormal, Parkinsonian states.
In recent studies, some alternative approaches have been proposed for the effective suppression of the global synchronization. In Kuramoto oscillator networks, the role of conformists, oscillators attracted to the mean field and tending to synchronize with it, and contrarians, repelled by the mean field and preferring a phase diametrically opposed to it, has been investigated, in order to suppress explosive synchronized activity. The latter refers to the transition from a non-synchronized state to a synchronized state in an abrupt/discontinuous manner (see e.g. [134][135][136][137]). Different strategies have been proposed for exploiting the local (contrarians) versus total information, the role of the negative versus positive coupling in order to achieve this goal. In our work, when implementing CR stimulation, we use inhibitory and excitatory synapses where all neurons are connected to each other. However, in [45] (Fig 12), it was shown that CR stimulation can also desynchronize effectively such networks with topologies where a fraction of neurons is excitatory and the rest inhibitory when receiving the same type of stimulus. The impact of the CR stimulation induces an increase in the inhibitory coupling weights and a decrease in the excitatory ones via the STPD. In [33], a multi-site CR stimulation, has proven to exhibit powerful long-term anti-kindling effects using a network of coupled phase oscillators with STDP. Furthermore in [66], it was shown that with CR stimulation one can achieve robust long-term curative effects, irrespectively, of the ratio between excitatory and inhibitory impact. Hence, in principle, different types of stimulation-induced modifications of plastic couplings may counteract synchronization.
Bistability and strong abnormal-explosive synchrony has also been studied recently in several physical and neuronal systems without synaptic plasticity (see e.g. [137][138][139][140] and references therein). These theoretical findings, obtained in generic networks, shed light on the underlying mechanisms on the transition from abnormal to normal activity in networks due to e.g. the interplay of local structure the internal dynamics or the critical role of the coupling strength. Moreover, experimental efforts are made in implementing such ideas, showing that conditions for explosive synchronization in the human brain could suggest a potential mechanism for rapid recovery from the lightly-anesthetized state [141].
Synaptic plasticity is another source for bistability and multistability in oscillatory networks. In fact, bistability and multistability was found in different networks with qualitatively different synaptic plasticity mechanisms. For example, in [32] a slow varying coupling matrix was used in a generalized Kuramoto model of coupled phase oscillators. In that case, the synaptic weights were governed by the time-varying phase difference of pairs of oscillators in such a way that the coupling strength increases for synchronized oscillators and weakens for nonsynchronized pairs. In [33], under spontaneous conditions two different states, a desynchronized and a synchronized state, were found in a system of coupled phase oscillators with asymmetric STDP. In addition, a similar bistability regime was found in a network of Morris-Lecar neurons with symmetric STDP in [66]. In network of phase oscillators with simplified, time-varying STDP multistability was shown to occur only for asymmetric STDP [34]. In that case, the coexistence of synchronized as well as desynchronized and cluster states depends on the distribution of the eigenfrequencies. In the network with Hodgkin-Huxley neurons with asymmetric STPD studied in this paper, a pronounced multistability was found, see also [45].
Though, not in the intended scope of our present study, the actual mechanism of action of CR stimulation deserves attention in future studies. In neural and oscillator networks without STDP, CR stimulation disrupts synchrony by causing phase resets of different subpopulations at different times [64,65,81]. The phase reset of a single subpopulation is time-locked to the stimulus affected that particular subpopulation [64,65,81]. However, in the present study we observed that CR stimulation may not just reorder the neurons' phases. Rather, for particular stimulation parameters it may even cause a significant decrease of the neuronal firing rates, intriguingly associated with a particularly pronounced anti-kindling (S2I Fig). Furthermore, in contradiction to the results obtained in networks without STDP [64,65,81], CR stimulation may cause a full-blown anti-kindling without any phase resets of the subpopulations time locked to the corresponding stimuli (S4 Fig). This is relevant for two reasons: (i) Since effective CR stimulation does not require phase resets time-locked to the individual stimuli, further computational studies should elucidate whether it makes sense to calibrate CR stimuli for preclinical and clinical applications by selecting stimulus parameters that favorably achieve phase resets. Corresponding results might be relevant for the design of calibration procedures and, in addition, challenge existing patents that are based on selecting parameters that optimally achieve phase resets of the stimuli delivered to the individual sub-populations (e.g. [142]). (ii) By the same token, our results do not only challenge current hypotheses on the mechanism of CR stimulation, but also fundamental patents in the field of invasive (e.g. [143]) as well as noninvasive (e.g. [144]) CR stimulation. Accordingly, future computational studies should focus on the mechanism of action of CR stimulation in networks with STDP in order to actually understand and possibly improve anti-kindling protocols.
Our goal is to accomplish an anti-kindling in a way as robust as possible, complying with clinically motivated constraints. For instance, striving for anti-kindling induced at minimal stimulation intensities led to the computational development of spaced CR stimulation [145] and two-stage CR stimulation with weak onset intensity [85]. The motivation behind these developments was to avoid side effects by substantially reducing stimulation intensities [85,145]. Another direction is to accomplish anti-kindling at moderate stimulation duration as computationally studied in this paper. This may be favorable from a clinical standpoint since it might help to reduce the occurrence of side effects as well as the requirement of the treatment on patients for their compliance, e.g. in terms of actually using non-invasive therapeutic devices. In this context, it might turn out to be beneficial that RVS CR stimulation causes sustained after-effects over a wide range of stimulation frequencies even at weak intensity (Fig 4). Accordingly, RVS CR stimulation might provide an appropriate stimulation protocol, in particular, if applied in an open-loop manner, without the ability to calibrate the stimulation parameters, especially the stimulation frequency by adapting it to the dominant peaks in the frequency spectrum of electrophysiological signals such as local field potentials or EEG signals. However, in a forthcoming computational study [118] we will use comparably simple closedloop control modes to significantly improve the robustness of both RVS and SVS CR stimulation and, in particular, exploit the anti-kindling potential of SVS CR stimulation.