Short-Term Dosage Regimen for Stimulation-Induced Long-Lasting Desynchronization

In this paper, we computationally generate hypotheses for dose-finding studies in the context of desynchronizing neuromodulation techniques. Abnormally strong neuronal synchronization is a hallmark of several brain disorders. Coordinated Reset (CR) stimulation is a spatio-temporally patterned stimulation technique that specifically aims at disrupting abnormal neuronal synchrony. In networks with spike-timing-dependent plasticity CR stimulation may ultimately cause an anti-kindling, i.e., an unlearning of abnormal synaptic connectivity and neuronal synchrony. This long-lasting desynchronization was theoretically predicted and verified in several pre-clinical and clinical studies. We have shown that CR stimulation with rapidly varying sequences (RVS) robustly induces an anti-kindling at low intensities e.g., if the CR stimulation frequency (i.e., stimulus pattern repetition rate) is in the range of the frequency of the neuronal oscillation. In contrast, CR stimulation with slowly varying sequences (SVS) turned out to induce an anti-kindling more strongly, but less robustly with respect to variations of the CR stimulation frequency. Motivated by clinical constraints and inspired by the spacing principle of learning theory, in this computational study we propose a short-term dosage regimen that enables a robust anti-kindling effect of both RVS and SVS CR stimulation, also for those parameter values where RVS and SVS CR stimulation previously turned out to be ineffective. Intriguingly, for the vast majority of parameter values tested, spaced multishot CR stimulation with demand-controlled variation of stimulation frequency and intensity caused a robust and pronounced anti-kindling. In contrast, spaced CR stimulation with fixed stimulation parameters as well as singleshot CR stimulation of equal integral duration failed to improve the stimulation outcome. In the model network under consideration, our short-term dosage regimen enables to robustly induce long-term desynchronization at comparably short stimulation duration and low integral stimulation duration. Currently, clinical proof of concept is available for deep brain CR stimulation for Parkinson's therapy and acoustic CR stimulation for tinnitus therapy. Promising first in human data is available for vibrotactile CR stimulation for Parkinson's treatment. For the clinical development of these treatments it is mandatory to perform dose-finding studies to reveal optimal stimulation parameters and dosage regimens. Our findings can straightforwardly be tested in human dose-finding studies.


INTRODUCTION
To establish a pharmacological treatment for clinical use, in humans typically a 4-phase sequence of clinical trials is performed (Friedman et al., 2010). In pre-clinical studies pharmacokinetic, toxicity and efficacy are studied in non-human subjects. In first in human-studies (phase I) safety and tolerability of a drug are studied in healthy volunteers. Proof of concept studies (phase IIA) determine whether a drug can have any efficacy, whereas dose-finding studies (phase IIB) are performed to reveal optimum dose at which a drug has biological activity with minimal side-effects. Effectiveness and the clinical value of a new intervention are studied in a randomized controlled trial (phase III), compared with state of the art treatment, if available. Finally, post-marketing surveillance trials (phase IV) are performed to detect rare or long-term adverse effects within a much larger patient population and over longer time periods. There might also be combinations of different phases.
In principle, this 4-phase pattern is also valid for medical technology, e.g., neuromodulation technologies. However, if neuromodulation technologies aim at the control of complex dynamics of e.g., neural networks, different parameters and dosage regimens may have complex, non-linear and even counterintuitive effects, (see e.g., Gao et al., 2014;Popovych et al., 2015;Gates and Rocha, 2016;Zañudo et al., 2017;Zhang et al., 2017). This computational paper illustrates how computational modeling can be used to generate hypotheses for dose-finding studies. In general, performing dose-finding studies simply by trial and error may be impossible because of the substantial parameter space to be tested, with trial durations and related costs getting out of hands.
The development of proper dosage strategies and regimens enables favorable compromises between therapeutic efficacy and detrimental factors such as side-effects or treatment duration. This is relevant, e.g., for the development of pharmaceutical (Williams, 1992;Bertau et al., 2008;Peters, 2012;Dash et al., 2014) or radiation therapy (Symonds et al., 2012). Deep brain stimulation (DBS) is the standard treatment of medically refractory movement disorders (Benabid et al., 1991;Krack et al., 2003;Deuschl et al., 2006). The clinical (Temperli et al., 2003) and electrophysiological (Kühn et al., 2008;Bronte-Stewart et al., 2009) effects of standard high-frequency (HF) DBS occur only during stimulation and cease after stimulation offset.
Coordinated reset (CR) stimulation (Tass, 2003a,b) was computationally developed to specifically counteract abnormal neuronal synchrony by desynchronization. CR stimulation uses sequences of stimuli delivered to neuronal sub-populations engaged in abnormal neuronal synchronization (Tass, 2003a,b). As shown computationally, in neuronal populations with spiketiming-dependent plasticity (STDP) (Gerstner et al., 1996;Markram et al., 1997;Bi and Poo, 1998) CR stimulation may have long-lasting, sustained effects (Tass and Majtanik, 2006;Hauptmann and Tass, 2007;Popovych and Tass, 2012). This is because in the presence of STDP, CR stimulation reduces the rate of coincidences. Accordingly, the network may be shifted from an attractor with abnormal synaptic connectivity and abnormal neuronal synchrony to an attractor with weak connectivity and synchrony (Tass and Majtanik, 2006;Hauptmann and Tass, 2007;Popovych and Tass, 2012). This process was termed anti-kindling (Tass and Majtanik, 2006). Abnormal neuronal synchronization has been shown to be associated with a number of brain diseases, for example, Parkinson's disease (PD) (Lenz et al., 1994;Nini et al., 1995;Hammond et al., 2007), tinnitus (Ochi and Eggermont, 1997;Llinás et al., 1999;Weisz et al., 2005;Eggermont and Tass, 2015), migraine (Angelini et al., 2004;Bjørk and Sand, 2008). In parkinsonian non-human primates it was shown that electrical CR stimulation of the subthalamic nucleus (STN) has sustained, long-lasting after-effects on motor function (Tass et al., 2012b;Wang et al., 2016). In contrast, long-lasting after-effects were not observed with standard HF DBS (Tass et al., 2012b;Wang et al., 2016). For instance, unilateral CR stimulation of the STN of parkinsonian MPTP monkeys, delivered for only 2 h per day during 5 consecutive days led to significant and sustained bilateral therapeutic after-effects for at least 30 days, whereas standard HF DBS had no after-effects (Tass et al., 2012b). By the same token, cumulative and lasting after-effects of electrical CR stimulation of the STN were also observed in PD patients (Adamchic et al., 2014).
HF DBS may not only cause side effects by electrical current spreading outside of the target region, but also by chronic stimulation of the target itself or by functional disconnection of the stimulated structure (Ferraye et al., 2008;Moreau et al., 2008;van Nuenen et al., 2008). Accordingly, it is key to reduce the integral stimulation current. Electrical CR stimulation of the STN employs significantly less current compared to HF DBS (Tass et al., 2012b;Adamchic et al., 2014;Wang et al., 2016). However, to further improve the CR approach, in a previous computational study the spacing principle (Ebbinghaus et al., 1913) was used to achieve an anti-kindling at subcritical intensities, i.e., particularly weak intensities rendering permanently delivered CR stimulation ineffective (Popovych et al., 2015). According to the spacing principle (Ebbinghaus et al., 1913), learning effects can be improved by repeated stimuli spaced by pauses as opposed to delivering a massed stimulus in a single long stimulation session. The spacing principle was investigated on different levels, ranging from behavioral and cognitive (Cepeda et al., 2006(Cepeda et al., , 2009Pavlik and Anderson, 2008;Xue et al., 2011;Kelley and Whatson, 2013) to neuroscientific and molecular (Itoh et al., 1995;Frey and Morris, 1997;Menzel et al., 2001;Scharf et al., 2002;Naqib et al., 2012). Computationally it was demonstrated that the spacing principle can also be applied to unlearn unwanted, upregulated synaptic connectivity at subcritical stimulation intensities (Popovych et al., 2015). In principle, the results were intriguing, but required rather long pauses and total stimulation durations (Popovych et al., 2015). Spaced CR stimulation at subcritical intensities might possibly be applied to CR DBS. However, for clinical applications, in particular, for non-invasive applications of CR stimulation (Popovych and Tass, 2012), such as acoustic CR stimulation for tinnitus (Tass et al., 2012a) or vibrotactile stimulation for PD (Tass, 2017;Syrkin-Nikolau et al., 2018), it is crucial to achieve therapeutic effects within a reasonable amount of time. Applications of non-invasive medtech devices typically rely on the patients' compliance and should favorably require short stimulation durations. Accordingly, we here set out to apply the spacing principle to CR stimulation at supercritical intensities, i.e., intensities that enable an anti-kindling for moderate stimulation duration and properly selected stimulation frequencies. The overall goal of this study is to design short-term dosage regimen that improve CR stimulation efficacy, while keeping the integral amount of stimulation as well as the overall duration of the protocols at comparably low levels.
In Manos et al. (in review) we studied the influence of the CR stimulation frequency and the intensity on the outcome of CR stimulation with Rapidly Varying Sequences (RVS) and Slowly Varying Sequences SVS (Zeitler and Tass, 2015). CR stimulation consists of sequences of stimuli delivered to each sub-population (Tass, 2003a,b). For RVS CR stimulation, the CR sequence is randomly varied from one CR stimulation period to another (Tass and Majtanik, 2006). Conversely, SVS CR stimulation is characterized by repeating a sequence for a number of times before randomly switching to the next sequence (Zeitler and Tass, 2015). In Manos et al. (in review) we demonstrated that the efficacy of singleshot CR stimulation with moderate stimulation duration depends on the stimulation parameters, in particular, on the intensity as well as the relationship between CR stimulation frequency and intrinsic firing rates. RVS CR stimulation turned out to induce pronounced long-lasting desynchronization, e.g., at weak intensities and CR stimulation frequencies in a certain range around the neurons' intrinsic firing frequencies. In contrast, SVS CR stimulation enabled even more pronounced anti-kindling, however, at the cost of a significantly stronger dependence of the stimulation outcome on the CR stimulation frequency.
Dosage regimen design is an integral part of pharmacokinetic methodology, aiming at an optimization of drug delivery and effects (Williams, 1992). By a similar token, we hypothesize that appropriate dosage regimens might further enhance the efficacy of RVS and SVS CR stimulation. To probe different dosage regimens, we here consider different stimulation singleshot and multishot CR stimulation protocols. Protocols A and B have identical integral stimulation duration, whereas Protocols C and D may require less stimulation.

Protocol A: Spaced Multishot CR Stimulation With Fixed Stimulation Parameters
Instead of one singleshot CR stimulation we deliver the identical CR shot five times, where the duration of each single pause equals the duration of each identical singleshot. Intersecting singleshot stimuli by pauses to increase stimulation efficacy, resembles the so-called spacing principle, a learning-related mechanism that is well-established in psychology (Ebbinghaus et al., 1913), education (Kelley and Whatson, 2013), and neuroscience (Naqib et al., 2012). According to the spacing principle, learning effects can be enhanced by delivering a stimulus in a spaced manner, as opposed to administering a massed stimulus in a single long stimulation session. Computationally, it was shown that subcritical CR stimulation at subcritical (ineffective) intensities may become effective if intersected by rather long pauses and delivered sufficiently often, e.g., eight times (Popovych et al., 2015). However, shorter pauses were not sufficient (Popovych et al., 2015). As yet, spaced CR stimulation at supercritical intensities was not studied. Here, we focus on comparably short stimulation protocols. Accordingly, we use CR stimulation of sufficient intensity and deliver five single CR shots intersected by pauses. We consider a symmetric dosage regimen, with identical duration of single shots and pauses.

Protocol B: Long Singleshot CR Stimulation With Fixed Stimulation Parameters
To assess the impact of the spacing principle, as a control condition we simply stimulate five times longer, without any pause and with stimulation parameters kept constant. Protocol B is shorter, but employs the same integral stimulation duration as Protocol A.

Protocol C: Spaced Multishot CR Stimulation With Demand-Controlled Variation of the CR Stimulation Frequency and Intensity
As in Protocol A, we deliver spaced CR stimulation comprising five identical CR shots, intersected by pauses, where all shots and pauses are of equal duration. However, at the end of each CR shot we monitor the stimulation effect and perform a three-stage control: (i) If no pronounced desynchronization is achieved, the CR stimulation frequency of the subsequent CR shot is mildly varied by no more than ±3%. (ii) If an intermediate desynchronization is observed, the CR stimulation frequency remains unchanged and CR stimulation is continued during the subsequent shot. (iii) If a pronounced desynchronization is achieved, no CR stimulation is delivered during the subsequent shot. Note, for stage (i) we do not adapt the CR stimulation frequency to a measured quantity. We consider two different variation types employed for stage (i): with regular and with random variation of the CR stimulation frequency. Regular variation means to increase or decrease the CR stimulation frequency in little unit steps. In contrast, random variation stands for randomly picking the CR stimulation frequency from a restricted interval.

Protocol D: Long Singleshot CR Stimulation With Demand-Controlled Variation of the Stimulation Frequency
To assess the specific pausing-related impact of the evolutionary spacing principle, as a direct control condition we perform Protocol C without pauses. To this end, we string five CR shots together, without pauses, and evaluate the stimulation effect at the end of each CR shot. If no pronounced desynchronization is achieved, the CR stimulation frequency is slightly varied by no more than ±3% for the subsequent CR shot. During each single CR shot stimulation parameters are kept constant. Only from one CR shot to the next the CR stimulation frequency can be varied.
Overall, Protocol D is shorter than Protocol C, but uses the same integral stimulation duration as in Protocols A-C.
CR stimulation and, especially, SVS CR stimulation has pronounced periodic characteristics. Accordingly, the CR stimulation frequency turned out to be a sensitive parameter, in particular, for SVS CR stimulation (see Manos et al., in review). For this reason, for stage (i) of Protocol C and D we perform a demand-controlled variation of the CR stimulation frequency to prevent from, e.g., unfavorable resonances or phase locking dynamics. Note these demand-controlled changes of the CR stimulation frequency are mild and hardly change the networks' firing rates.
In this study, we test the performance of the different Protocols A-D by selecting unfavorable stimulation parameters, which render CR stimulation ineffective according to Manos et al. (in review). By design, Protocols C and D work well for all parameter pairs (K, T s ) related to effective singleshot CR stimulation. In that case, CR stimulation actually ceases due to lack of demand. Note, in all four stimulation protocols we keep the stimulation intensity fixed. Only Protocols C and D require feedback of the stimulation outcome.
This paper is organized as follows: in the Materials and Methods section we briefly describe the computational model, the neural network (and its initialization), the synaptic plasticity rule, the CR stimulation, the analysis methods used throughout the paper as well as the summary of the CR frequency and intensity global trends which were thoroughly studied in Manos et al. (in review). In the Results section, we present all the different Protocols in detail and our main findings regarding their comparison. Finally, in the Discussion section, we discuss our findings and set this work in more general perspective, related to medical applications.

Model and Network Description
In this study we use the conductance-based Hodgkin-Huxley neuron model (Hodgkin and Huxley, 1952) for the description of an ensemble of spiking neurons. The set of equations and parameters read (see Hansel et al., 1993;Popovych and Tass, 2010): Variable The initial values of the neural synaptic weights c ij are picked from a normal distribution N(µ c = 0.5 mS/cm 2 , σ c = 0.01 mS/cm 2 ) as in Manos et al. (in review, see Popovych and Tass, 2012;Zeitler and Tass, 2015 for details). S i (t) denotes the internal synaptic input within the network to neuron i. The neurons interact via excitatory and inhibitory chemical synapses s i , by means of the post-synaptic potential (PSP) s i which is triggered by a spike of neuron i (Gerstner et al., 1996;Izhikevich, 2010) and modeled using an additional equation (see Golomb and Rinzel, 1993;Terman et al., 2002): Initially we draw s i ∈ [0, 1 ] (randomly from a uniform distribution). The coupling term S i from Equation 1a (see Popovych and Tass, 2012) contains a weighted ensemble average of all post-synaptic currents received by neuron i from the other neurons in the network: S i = N −1 N j=1 V r,j − V i c ij M ij s j , where V r,j is the reversal potential of the synaptic coupling (20 mV for excitatory and −40 mV for inhibitory coupling), and c ij is the synaptic coupling strength from neuron j to neuron i. There are no neuronal self-connections within the network (c ii = 0 mS/cm 2 ). The variable: describes the spatial profile of coupling between neurons i and j and is of a Mexican hat-type (Wilson and Cowan, 1973;Dominguez et al., 2006;de la Rocha et al., 2008) with strong short-range excitatory M ij > 0 and weak long-range inhibitory N−1 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:

Spike-Timing-Dependent Plasticity
The synaptic weights c ij are dynamical variables that depend on the time difference, t ij = t i − t j , between the onset of the spikes of the post-synaptic neuron i and the pre-synaptic neuron Frontiers in Physiology | www.frontiersin.org j, denoted by t i and t j , according to Bi and Poo (1998) and Popovych and Tass (2012): with parameters β 1 = 1, β 2 = 16, γ 1 = 0.12, γ 2 = 0.15, τ = 14 ms and with learning rate δ = 0.002, while the values of c ij are confinded to the interval [0, 1] mS/cm 2 for both excitatory and inhibitory synapses and, hence, remain bounded.

Coordinated Reset Stimulation
The term F i in Equation 1a represents the current induced in neuron i by the CR stimulation delivered at N s = 4 stimulations sites, equidistantly placed at the positions of neurons i = 25, 75, 125, 175 (Tass, 2003b). One stimulation site was active during T s /N s , while the other stimulation sites were inactive during that time window. After this, another stimulation site was active during the next T s /N s window. All N s stimulation sites were stimulated exactly once within one CR stimulation period of duration T s . The spatiotemporal activation pattern of stimulation sites is represented by the indicator functions ρ k (t) (kǫ {1, . . . , N}), taking the value 1 when the kth stimulation site is active at t and 0 else. The stimulation signals induced single brief excitatory post-synaptic currents. The evoked time-dependent normalized conductances of the post-synaptic membranes are represented by α-functions given in Popovych and Tass (2012) denotes the time-to-peak of G stim , and t k is the onset of the kth activation of the stimulation site. Note, τ stim determines the onset timing of each single signal as well as its duration. The spatial spread of the induced excitatory post-synaptic currents in the network is defined by the quadratic spatial decay profile D (i, , a function of the difference between the index of neuron i and the index x k of the neuron at stimulation site k. d is the lattice distance between two neighboring neurons, and σ d = 0.8 the spatial decay rate of the stimulation current (see Popovych and Tass, 2012 for details). Thus, the total stimulation current from Equation , where V r = 20 mV is the excitatory reverse potential, and K the stimulation intensity.

Macroscopic Measurements
We measure the strength of the coupling within the neuronal population at time t by calculating their total synaptic weight (averaged over the neuron population) , sgn is the sign-function, while C av is calculated by averaging over the last 100 · T s . The extent of in-phase synchronization within the network is assessed by the order parameter (Haken, 1983;Kuramoto, 2012) is a linear approximation of the phase of neuron j between its mth and (m + 1)th spikes at spiking times t j,m and t j,m+1 . R (t) = 1 for complete inphase synchronization, and R (t) = 0 in the absence of in-phase synchronization. Because of strong fluctuations of the order parameter, we calculate the moving average < R > over a time window of 400 · T s , to investigate the time evolution of the order parameter. Moreover, we use the quantity R av , which is the order parameter R (t) averaged over the last 100·T s of a pause following a CR shot or of the end of the post-stim epoch. For the statistical description and analysis of the non-Gaussian distributed R av data (n = 11 samples), we use boxplots (Tukey, 1977). Their Inter-Quartile Range measures the statistical dispersion around the median, which is defined as width of the middle 50% of the distribution and is represented by a box. 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.

Dependence of CR Stimulation Outcome on CR Stimulation Frequency and Intensity
This section provides a short overview of the results of a study where CR stimulation frequency and intensity were varied in detail (Manos et al., in review). That study revealed the dependence of the outcome of RVS and SVS CR on the CR stimulation frequency and intensity and, in particular, possible limitations thereof, especially for SVS CR stimulation. Based on these limitations, the present study presents an approach that enables to overcome these issues.
In the present study, for each initial network condition and its corresponding parameters (simply denoted as network), we apply RVS and SVS CR stimulation with different realizations of the CR sequence orders per network. We start the simulations with an equilibration phase without STDP, which lasts for 2 s. From this point on, the network evolves in the presence of STDP, starting with a 60 s integration with STDP only (i.e., without stimulation), where a rewiring of the connections takes place, resulting in a strongly synchronized state with intrinsic firing rate f int ≈ 71.4 Hz (corresponding to a period of T int = 14 ms). We then run four different CR stimulation protocols, resetting the starting time to t = 0 s. We use 3:2 ON-OFF CR stimulation, where three stimulation ON-cycles (with stimulation on) alternated with two OFF-cycles (without stimulation), with ON-/OFF-cycle duration of T s . 3:2 ON-OFF CR stimulation was used in a number of computational, pre-clinical and clinical studies, for details and motivation (see Manos et al., in review).
To study dosage regimens that potentially improve reliability and stimulation outcome of RVS and/or SVS stimulation, in the present study we focus on parameter ranges where RVS and/or SVS CR stimulation did not reliably elicit long-lasting desynchronization according to Manos et al. (in review). In Manos et al. (in review), we delivered CR single shots of 128 s duration followed by a 128 s CR-off period and varied the CR stimulation frequency and intensity over a wider range. In this way, we showed that RVS CR stimulation turned out to be more robust against variations of the stimulation frequency, while SVS CR stimulation can obtain stronger anti-kindling effects. This dependence on the CR stimulation intensity and frequency is summarized in Figure 1. Figures 1A,E show the boxplots for the time-averaged mean synaptic weights C av (at the end of the 128 s CR-off period) with values belonging to the same intensity value K for RVS and SVS CR stimulation, respectively, but CR stimulation frequency varying in the interval [25%f 0 , . . . , 175%f 0 ], with f 0 as defined below. The motivation for restricting the CR stimulation intensity to the interval K ǫ [0.20, . . . , 0.50] is discussed in Manos et al. (in review). In a similar manner, Figures 1B,F show boxplots for the timeaveraged order parameter R av (again at the end of the CR-off period) with values belonging to the same intensity value K for RVS and SVS CR stimulation respectively. Figures 1C,G depict the boxplots for the time-averaged mean synaptic weights C av (at the end of the CR-off period) but now with values belonging to the same frequency ratio (f stim /f 0 ) · 100 for RVS and SVS CR stimulation, respectively, but intensity value K ∈ [0.20, . . . , 0.50]. The CR stimulation frequency f stim takes values in the interval [25%f 0 , . . . , 175%f 0 ], where f 0 = 1/T 0 denotes the initial stimulation frequency. The choice of the frequency f 0 (or period T 0 ) of the CR stimulation is made with respect of the intrinsic network's firing rate frequency (or period) f int (or T int ) and is meant to have a value close to that (in this case T 0 = 16 ms with f 0 = 62.5 Hz). More details can be found in Manos et al. (in review). Figures 1D,H show the boxplots for the time-averaged order parameter R av (at the end of the CR-off period) with values belonging to the same frequency ratio (f stim /f 0 ) · 100 for RVS and SVS CR stimulation, respectively.

Simulation Description
We investigate two singleshot and two multishot, spaced CR stimulation protocols (Figure 2). The multishot Protocols A and C consist of five single CR shots of 128 s duration, each followed by a pause of 128 s, respectively (Figures 2A,C). The CR singleshot Protocol B consists of a long singleshot of 5 × 128 s followed by a pause of 5 × 128 s (Figures 2B). The CR singleshot Protocol D consists of a long singleshot consisting of five single shots of 128 s duration, strung together without pauses in between, followed by a pause of 5 × 128 s (Figures 2D).
The integral stimulation duration is identical for Protocols A and B. In Protocols A and B all stimulation parameters are kept constant. In contrast, in Protocol C at the end of each pause the amount of synchrony is evaluated in a time window of 100 stimulation periods length ( Figure 2) and a three-stage control scheme is put in place: (i) If the amount of synchrony does not fall below a pre-defined threshold, the CR stimulation frequency is mildly varied. (ii) If the desynchronization effect is moderate, the CR stimulation frequency remains unchanged. (iii) If desynchronization is achieved, the stimulation intensity is set to zero for the subsequent shot. Analogously, in Protocol D at the end of each single shot the amount of synchrony is evaluated in a time window of 100 stimulation periods length (Figure 2) and the three-stage control scheme is executed. The difference between Protocol C and D is that the evaluation for the control intervention is performed in a pause subsequent to a single shot (Protocol C) as opposed to during a single shot (Protocol D).
For the stage (i) control, the variation of the CR stimulation frequency is not adapted to frequency characteristics of the neuronal network. Rather a minor variation of the CR stimulation frequency is performed to make a fresh start with the subsequent single CR shot. These minor changes of the CR stimulation frequency do not lead to changes of the neurons' intrinsic firing rates of more than ±3%.
Due to the stage (iii) control, the demand-controlled shutdown of CR stimulation, the maximum integral stimulation duration of Protocol C and D can reach the level of Protocols A and B, but may well fall below. We use the order parameter to assess the amount of synchronization (see Materials and Methods).

Protocol A: Spaced Multishot CR Stimulation With Fixed Stimulation Parameters
For this stimulation protocol all stimulation parameters are kept constant (Figure 2). Accordingly, the CR stimulation period T s remains constant, too. We study the stimulation outcome of only five symmetrically spaced consecutive single CR shots. To this end, for both RVS CR and SVS CR stimulation we consider two unfavorable parameter pairs of fixed CR stimulation period and intensity, respectively. One example refers to cases where CR stimulation induces acute effects, but no long-lasting desynchronizing effects (Cases I and IV). The other example concerns the case where CR stimulation causes neither acute nor long-lasting desynchronizing effects in a reliable manner (Cases II and III). At a stimulation duration of 128 s these parameters caused only an acute, but no long-lasting desynchronization in the majority of networks studied [ Figure 4B of Manos et al. (in review), where T s = 11 ms corresponds to ∼127% of the intrinsic firing rate (or ∼91 Hz)]. Case II: (K, T s ) = (0.20, 28). In the majority of networks tested, these parameters did neither lead to acute nor long-lasting desynchronization after administration of a single CR shot [ Figure 5B of Manos et al. (in review), where T s = 28 ms corresponds to ∼50% of the intrinsic firing rate (or ∼36 Hz)]. For both cases, we investigate the order parameter < R > averaged over a sliding window for 11 different networks (marked with different color/line types) (Figures 2A,C). Boxplots of the order parameter R av averaged over a window of length 100 · T s at the end of each pause demonstrate the overall stimulation outcome for all tested 11 networks (Figures 2B,D).

Case I:
RVS CR stimulation induces a desynchronization during the CR shots ( Figure 3A), but no reliable, long-lasting desynchronization in the subsequent pauses ( Figure 3B). The spacing protocol with five identical RVS CR shots does not significantly improve the desynchronizing outcome of a single RVS CR shot. In fact, in the boxplots the large dispersion around the median value remains almost unchanged in the course of this protocol ( Figure 3B). Case II: Neither during the RVS CR shots nor during the subsequent pauses a sufficient desynchronization is observed (Figures 2C,D). The spacing protocol does not cause an improvement of the stimulation outcome in this case, too ( Figure 3D). SVS CR stimulation: Case III: (K, T s ) = (0.20, 9) Single shot SVS CR stimulation with these parameters caused neither pronounced acute nor long-lasting desynchronization [ Figure 6D of Manos et al. (in review), where T s = 9 ms corresponds to ∼156% of the intrinsic firing rate (or ∼111 Hz)]. Case IV: (K, T s ) = (0.20, 14). Single shot SVS CR stimulation with these parameters led to acute, but no longlasting desynchronization in the majority of networks tested [ Figure 7B of Manos et al. (in review), where T s = 28 ms corresponds to ∼50% of the intrinsic firing rate (or ∼36 Hz)]. For both cases we performed the same analysis as shown in Figure 3.

Case III:
SVS CR stimulation neither induces a pronounced and reliable desynchronization during the CR shots ( Figure 4A) nor during the pauses (Figure 4B). In fact, the dispersion around the median value is increased during the fourth and fifth pause ( Figure 4B).

Case IV:
During the SVS CR shots a desynchronization occurs ( Figure 4C). However, no reliable and pronounced desynchronization is observed during the pauses ( Figure 4D).
In summary, for both RVS CR and SVS CR stimulation the spacing protocol with five consecutive CR shots does not cause an improvement of the long-lasting desynchronization (assessed after cessation of stimulation). We performed the same analysis  for a larger set of (K, T s ) pairs in the parameter plane analyzed in Manos et al. (in review), with K ranging from 0.2 to 0.3 (weak intensities) and T s ranging from 9 to 28 ms (around the intrinsic period). For all parameter pairs tested, the spacing Protocol A did not improve the long-term desynchronization effect.

Protocol B: Long Singleshot CR Stimulation With Fixed Stimulation Parameters
For this stimulation protocol all stimulation parameters are kept constant, too (Figure 2). Instead of five single CR shots of 128 s duration each (Figures 1-3), we deliver one fivefold longer singleshot of 5 × 128 s duration (Figure 2). For both RVS CR and SVS CR stimulation we consider the corresponding two unfavorable parameter pairs of fixed CR stimulation period and intensity already studied above (Cases I-IV). This is to study whether a fivefold prolongation of the stimulation duration leads to an improvement of the stimulation outcome.

RVS CR stimulation:
For comparison, we consider the cases studied above. Case I: (K, T s ) = (0.30, 11). Case II: (K, T s ) = (0.20, 28). We study the order parameter < R > averaged over a sliding window for 11 different networks (marked with different color/line types in Figure 5A). The overall stimulation outcome for all tested 11 networks is illustrated with boxplots of the order parameter R av averaged over a window of length 100 · T s at the end of the post-stimulus epoch ( Figure 5B).

Case I:
RVS CR stimulation induces a desynchronization during the long RVS CR singleshot ( Figure 5A), but no reliable, longlasting desynchronization in the subsequent pauses ( Figure 5B). The median of the order parameter of the long-term outcome hardly changes, but the dispersion around the median value is greater for the post-stim order parameter (as the IQR of the boxplots in Figure 5B show). Note, the overall longterm desynchronization for the long singleshot ( Figure 5B) is more pronounced compared to the spaced RVS CR stimulation Protocol A (Figure 3B).

Case II:
Neither during the long RVS CR singleshot nor during the subsequent stimulation-free epoch a reliable and pronounced desynchronization is observed (Figures 4C,D). Interestingly, the network that undergoes an acute desynchronization during the singleshot relaxes back to a synchronized state (Figure 5B, green curve). Conversely, the only network that displays a long-term desynchronization does not undergo a pronounced desynchronization during the singleshot (Figure 5B, magenta curve).

SVS CR stimulation:
For comparison, we consider the time course of the timeaveraged order parameter < R > (Figures 5A,C) and the corresponding boxplots of the order parameter R av averaged over

Case III:
In the majority of networks SVS CR stimulation does not induce a pronounced and reliable desynchronization during the fivefold longer SVS CR singleshot as well as in the post-stim epoch ( Figure 6A). Three out of 11 networks display a pronounced acute and long-lasting desynchronization ( Figure 6A). Accordingly, the dispersion around the median is large during and after the singleshot (Figure 6B).

Case IV:
During the long SVS CR singleshot a pronounced desynchronization occurs (Figure 6C), as reflected by the small dispersion around the small median in the corresponding boxplot ( Figure 6D). However, in the post-stimulation epoch most of the networks relax to a synchronized state, with only a few networks remaining in a long-term desynchronized state ( Figure 6D). Accordingly, there is a large dispersion around a large median in the boxplot (Figure 6D).
In summary, for both RVS CR and SVS CR stimulation the fivefold increase of the stimulation duration does not lead to a reliable and pronounced long-lasting desynchronization. Again, we performed the same analysis for a larger set of (K, T s ) pairs in the parameter plane analyzed in Manos et al. (in review), with K ranging from 0.2 to 0.3 (weak intensities) and T s ranging from 9 to 28 ms (around the intrinsic period value). For all parameter pairs tested, the spacing Protocol B did not lead to a reliable and pronounced long-term desynchronization.

Protocol C: Spaced Multishot CR Stimulation With Demand-Controlled Variation of Stimulation Period T s and Intensity
We study the stimulation outcome of only five symmetrically spaced consecutive single CR shots with stimulation period T s and intensity varied according to a three-stage control scheme. To this end, for both RVS CR and SVS CR stimulation we consider two unfavorable parameter pairs of fixed CR stimulation period and intensity, respectively. One example refers to cases where CR stimulation induces acute effects, but no long-lasting desynchronizing effects (Cases I and IV). The other example concerns the case where CR stimulation causes neither acute nor long-lasting desynchronizing effects in a reliable manner (Cases II and III). We consider a regular and a random type of demandcontrolled variation of the CR stimulation period T s . Note, in both cases the CR stimulation period is not adapted to frequency characteristics of the network. We consider the time courses of the time-averaged order parameter < R > and R av , the order parameter averaged over a window of length 100 · T s at the end of pause.

Demand-controlled regular variation of the CR stimulation period and demand-controlled variation of the intensity
At the end of each pause we calculate the order parameter R av averaged over a window of length 100 · T s . We vary the CR stimulation period and intensity according to the amount of synchrony, based on a three-stage control scheme: (i) Insufficient desynchronization: If R av > 0.4, we decrease the CR stimulation period of the subsequent RVS shot by T s j + 1 = T s j −1 ms, where the index j stands for the j-th CR shot. As lower bound we set T s = 9 ms (corresponding to ∼156% of the intrinsic firing rate), in order to avoid undesirably high CR stimulation frequencies.
In a previous computational study the latter turned out to be unfavorable for desynchronization (see Manos et al., in review). As soon as T s reaches its lower bound of 9 ms, it is reset to T s (1) + 1 ms. (ii) Moderate desynchronization: If 0.2 ≤ R av ≤ 0.4, we preserve the CR stimulation period for the subsequent CR shot: T s j + 1 = T s j , where the index j denotes the j-th CR shot. 0.2 ≤ R av ≤ 0.4 is considered to be indicative of a desynchronization effect. (iii) Sufficient desynchronization: If R av < 0.2, the CR stimulation is suspended for the subsequent shot by setting K = 0 for the next shot and until 0.2 ≤ R av . R av < 0.2 is considered a sufficient desynchronization.

Spaced Multishot RVS CR Stimulation With Demand-Controlled Regular Variation of the Stimulation Period T s and Demand-Controlled Variation of the Intensity
In both Cases (I and II) this protocol reliably induces a desynchronization for all networks tested (Figures 7A,C). After the second RVS CR shot the median of the time-averaged order parameter R av at the end of the corresponding pauses falls below 0.4, with moderate dispersion (Figures 7B,D). Note, already after the first mild variation of the CR stimulation period T s the amount of synchrony is strongly reduced. In several networks and pauses, the desynchronization criterion, R av < 0.2, is fulfilled, so that during the subsequent CR shots no stimulation is delivered (Figures 7A,C). Accordingly, Protocol C enables to reduce the integral amount of stimulation.

Spaced Multishot SVS CR Stimulation With Demand-Controlled Regular Variation of the Stimulation Period T s and Demand-Controlled Variation of the Intensity
This protocol causes a desynchronization for all networks tested in both Cases (III and IV) (Figures 8A,C). After the third (Case III, Figure 8B) or the second SVS CR shot (Case IV, Figure 8D) a pronounced desynchronization is achieved, as reflected by a median of R av close to 0.2 (Figures 8B,D). Accordingly, about half of the networks fulfilled the desynchronization criterion R av < 0.2 after the third SVS CR shot and, hence, did not require further CR stimulation. The mean firing rate, measured at the end of each pause did not deviate from the baseline firing rates by more than ±3%, irrespective of the extent of protocol-induced variation of the stimulation period T s (data not shown). We do not adapt the CR stimulation period T s to frequency characteristics of the stimulated network. To further illustrate this aspect, we replace a regular, increasing or decreasing variation of the stimulation period by a random variation.

Demand-Controlled Random Variation of the CR Stimulation Period and Demand-Controlled Variation of the Intensity
Again, at the end of each pause we calculate the order parameter R av averaged over a window of length 100·T s . A random variation of the CR stimulation period is performed, depending on the amount of synchrony detected. To this end, we select the interval [T s (1) − 4 ms, T s (1) + 4 ms], where T s (1) denotes the CR stimulation period of the first shot. By design, this interval has a lower bound at 9 ms. The three-stage control scheme is governed by: (i) Insufficient desynchronization: If R av > 0.4 at the end of the pause of the j-th CR shot, we randomly pick T s j + 1 and skip inefficient values used before.
(ii) Moderate desynchronization: If 0.2 ≤ R av ≤ 0.4, we preserve the CR stimulation period for the subsequent CR shot: T s j + 1 = T s j , where the index j denotes the j-th CR shot. (iii) Sufficient desynchronization: If R av < 0.2, the CR stimulation is suspended for the subsequent shot by setting K = 0 for the next shot and until 0.2 ≤ R av .
The feasibility of this protocol is demonstrated by considering one example for RVS CR stimulation (Case II, Figure 9) and one for SVS CR stimulation (Case IV, Figure 10). For both cases we additionally provide the mean firing rate of the networks at the end of each shot and at the end of each subsequent pause to demonstrate that deviations do not exceed ±3% (Figures 9C, 10C). In the RVS case (Figure 9), the time course of the order parameter < R > ( Figure 9A) and the corresponding boxplots of R av (Figure 9B) display a similar pattern of reliable desynchronization as obtained by Protocol C with regular variation of the CR stimulation duration T s (Figure 7). In principle, the SVS case (Figure 10) provides similar findings as with a regular variation of the CR stimulation duration T s (Figure 8). However, one network relaxes back to a strongly synchronized state (Figure 10A, dashed blue line). Delivering a sixth SVS CR shot with randomly varied T s to that network caused a desynchronization (data not shown). This example illustrates that a sequence of five SVS CR shots might not be sufficient to induce desynchronization in all possible networks. In summary, for the five-shot RVS CR as well as SVS CR stimulation Protocol C with regular as well as random variation of the CR stimulation duration T s we observed a pronounced desynchronization, with the exception of one network ( Figure 10A, dashed blue line). Our analysis was performed for a larger set of (K, T s ) pairs in the parameter plane analyzed in Manos et al. (in review), with K ranging from 0.2 to 0.3 (weak intensities) and T s ranging from 9 to 28 ms (around the intrinsic period). For all parameter pairs tested, the spacing Protocol C with regular and random variation of T s led to a reliable and pronounced long-term desynchronization in the vast majority of networks tested.

Protocol D: Long Singleshot CR Stimulation With Demand-Controlled Variation of the Stimulation Frequency
Protocol D consists of five consecutive shots. Unlike in Protocol C, there are no pauses between the five consecutive shots, so that they form one long singleshot.

Demand-controlled regular variation of the CR stimulation period and demand-controlled variation of the intensity
At the end of each shot we calculate the order parameter R av averaged over a window of length 100 · T s . We vary the CR stimulation period and intensity according to the amount of synchrony, based on the three-stage control scheme as used for Protocol C (see above).

Long singleshot RVS CR stimulation with demand-controlled variation of the stimulation frequency
In Case I this protocol seems to perform similarly well (Figures 11A,B) as Protocol C (Figures 7A,B) and Protocol B (Figures 5A,B) which is also an alternative long singleshot but with fixed T s . After the second RVS CR shot almost all networks reach a moderate or sufficient desynchronization which is maintained fairly well after the RVS CR is ceased. Nonetheless, this particular protocol does not perform equally well for Case II (Figures 11C,D). Even in the cases, where the variation of T s leads to some improvement, the overall long-lasting effect is worse than with Protocol C (Figures 7C,D).

Long singleshot SVS CR stimulation with demand-controlled variation of the stimulation frequency
In Case III this protocol does not show any systematical improvement (Figures 12A,B). In fact, this is partly due to the fact that in some cases the network gets trapped in an unfavorable parameter variation loop, bouncing between T s = 9 ms and T s = 10 ms. In Case IV (Figures 12C,D) the global evolution is quite similar to the one found for Protocol B (Figures 6C,D), i.e., pronounced desynchronization during a single shot, with a FIGURE 9 | Protocol C: Spaced multishot RVS CR stimulation with demand-controlled random variation of the stimulation period T s and with demand-controlled variation of the intensity. (A) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots. The CR stimulation period is randomly varied depending on R av , by randomly picking a value from a narrow interval around the start period (see text). The horizontal solid red lines indicate the CR shots, while the horizontal dashed gray lines highlight the two control thresholds (see text). (B) Boxplots for the time-averaged order parameter R av at the end of each pause, illustrate the overall outcome for all tested 11 networks. (C) For each (vertically aligned) network the table presents CR intensity (K = 0.2 if CR is ON or K = 0 if CR is OFF during a CR shot) and stimulation period T s used for each CR shot (indicated by red bars) together with the mean firing rate of the network at the end of each CR shot ("ON") and at the end of the subsequent pause ("OFF"). The mean firing rate was strongly fluctuating and, hence, calculated in a window of length 200 · T s . Case II stimulation parameters are unfavorable for anti-kindling: (K, T s ) = (0.20, 28) (see text). tendency to relapse back to the synchronized state while some of the networks remain desynchronized. However, the overall final outcome is rather poor as the corresponding boxplot (blue color) at the end of the CR-off period indicates.
In summary, for the singleshot RVS CR as well as SVS CR stimulation Protocol D with regular variation of the CR stimulation duration T s (without pauses between two consecutive shots) did not lead to a reliable and systematic long-lasting desynchronization.

DISCUSSION
By comparing spaced CR stimulation with fixed stimulation parameters (Protocol A) and massed, continuous CR stimulation with equal integral duration (Protocol B) with a flexible spaced CR stimulation with demand-controlled variation of CR stimulation frequency and intensity (Protocol C), and with a flexible non-spaced CR stimulation with demand-controlled variation of CR stimulation frequency and intensity (Protocol D), we demonstrated that Protocol C enables to significantly improve the long-term desynchronization outcome of both RVS and SVS CR stimulation, even at comparatively short integral stimulation duration. Remarkably, spacing alone (Protocol A) is not sufficient to provide an efficient short-term dosage regimen (Figures 2, 3). In fact, in particular cases fivefold longer stimulation duration might even be more efficient than five consecutive single CR shots with identical integral stimulation duration, at least for RVS CR stimulation (Figure 2B vs. Figure 4B). The low performance of pure spacing (Protocol A) might be due to the low number of single CR shots, here five, as opposed to slightly larger numbers of CR shots, say eight, tested for the case of subcritical CR stimulation before (Popovych et al., 2015). However, more important might be the approximately fiftyfold longer stimulation and pause duration used for the spaced subcritical CR stimulation protocol (Popovych et al., 2015). The long spaced subcritical CR stimulation protocol might be FIGURE 10 | Protocol C: Spaced multishot SVS CR stimulation with demand-controlled random variation of the stimulation period T s and with demand-controlled variation of the intensity. (A) Time evolution of the order parameter < R > averaged over a sliding window during five consecutive RVS CR shots with demand-controlled random variation of T s (see text). (B) Boxplots for the time-averaged order parameter R av at the end of each pause, illustrate the overall outcome for all tested 11 networks. (C) For each (vertically aligned) network the table presents CR intensity (K = 0.2 if CR is ON or K = 0 if CR is OFF during a CR shot) and stimulation period T s used for each CR shot (indicated by red bars) and the mean firing rate of the network at the end of each CR shot ("ON") and at the end of the subsequent pause ("OFF"). Case IV stimulation parameters are unfavorable for anti-kindling: (K, T s ) = (0.20, 14) (see text). Same format as in Figure 9.
beneficial for invasive application, such as DBS, and help reduce side-effects by substantially reducing stimulation current intake of the issue.
However, computationally we show that a spacing with rigid five-shot timing structure, but flexible, demand-controlled variation of stimulation frequency and intensity (Protocol C) provides a short-term dosage regimen that significantly improves the long-term desynchronization outcome of RVS and SVS CR stimulation (Figures 6-9). At the end of each pause between CR shots, the stimulus after-effect is assessed. If the desynchronization is considered to be insufficient, a mild variation of the CR stimulation frequency is performed to possibly provide a better fit between network and CR stimulation frequency, without actually adapting the stimulation frequency to frequency characteristics of the network stimulated. If desynchronization is considered to be moderate, the subsequent CR shot is delivered with parameters unchanged. If desynchronization is sufficient, CR stimulation is suspended during the subsequent shot. Intriguingly, in the vast majority of parameters and networks tested, this short-term dosage regimen induces a robust and reliable long-lasting desynchronization (Figures 6-9). This protocol might be a candidate especially for non-invasive, e.g., acoustic (Tass et al., 2012a) or vibrotactile (Tass, 2017;Syrkin-Nikolau et al., 2018), applications of CR stimulation to increase desynchronization efficacy, while keeping the stimulation duration at moderate levels.
Demand-controlled variation of CR stimulation frequency and intensity (Protocol D) alone (i.e., without inserting pauses) is not sufficient to significantly improve the outcome of RVS and SVS stimulation (Figures 11, 12). Hence, introducing pauses significantly improves the effect of the demand-controlled variation of CR stimulation frequency and intensity.
In principle, stimulation parameters other than the CR stimulation frequency might be varied depending on the stimulation outcome. However, in this study we have chosen to vary the CR stimulation frequency, since the latter turned out to be a sensitive parameter, especially for SVS CR stimulation (see Manos et al., in review). In fact, the short-term dosage regimen FIGURE 11 | Long singleshot CR stimulation with demand-controlled variation of the stimulation frequency. (A,C) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots. If R av , the order parameter averaged over a window of length 100 · T s at the end of a shot, exceeds 0.4, the CR stimulation period of the subsequent RVS shot is decreased by T s → T s −1 ms (see text). Stimulation parameters are unfavorable for anti-kindling in Case I (A,B) and Case II (C,D) (see text). (B,D) Boxplots for the time-averaged order parameter R av at the end of each shot (red color) and at the end of CR-off period (blue color), illustrate the overall outcome for all tested 11 networks. Case I: (K, T s ) = (0.30, 11). Case II: (K, T s ) = (0.20, 28). Same format as in Figure 7.
with demand-controlled variation of stimulation parameters (Protocol C) might help to turn SVS CR stimulation in a method that causes a particularly strong anti-kindling in a robust and reliable manner.
Protocol C does not require a direct adaption of the CR stimulation frequency to measured quantities reflecting frequency characteristics of the network. We have chosen this design, since it might be an advantage not to rely on specific biomarker-type of information. For instance, in the case of PD a number of relevant studies were devoted to closed-loop DBS (Graupe et al., 2010;Rosin et al., 2011;Carron et al., 2013;Little et al., 2013;Priori et al., 2013;Yamamoto et al., 2013;Hosain et al., 2014;Rosa et al., 2015). A relevant issue in this context is the availability of a biomarker adequately reflecting the individual patient's extent of symptoms (Beudel and Brown, 2016;Kühn and Volkmann, 2017). In fact, it is not clear whether low or high frequency beta band oscillations might be more appropriate as biomarker-type of feedback signal (Beudel and Brown, 2016). For several reasons, beta band oscillations might possibly not be an optimal feedback signal Kühn and Volkmann, 2017). Enhanced beta band oscillations are not consistently found in all PD patients (Kühn et al., 2008;Kühn and Volkmann, 2017). The clinical score of PD patients might more appropriately be reflected by the power ratio of two distinct bands of high frequency oscillations around 250 and 350 Hz (Özkurt et al., 2011). Appropriate biomarkers might depend on the patient phenotype (Quinn et al., 2015): In tremor dominant (compared to akinetic rigid) PD patients resting state beta power may decrease during tremor epochs (Bronte-Stewart et al., 2009;Quinn et al., 2015). By a similar token, theta and beta oscillations interact with high-frequency oscillations under physiological (Yanagisawa et al., 2012) as well as pathological (Yang et al., 2014) conditions. Also, quantities assessing the interaction of brain oscillation, e.g., phase amplitude coupling (PAC) might be used as biomarker to represent the amount of symptoms (Beudel and Brown, 2016). Also, activity in the beta band might be relevant for compensatory purposes, as recently shown in a parkinsonian monkey study with sensorimotor rhythm neurofeedback (Philippens et al., 2017).
It might be another potential advantage for clinical applications that the three-stage control of the proposed short-term dosage regimen (Protocol C) could possibly be approximated by scores reflecting the patient's state or the amount symptoms. A simple three-stage rating of the patient's state (bad, medium, and good) might replace the feedback signal-based stages (i), (ii), and (iii). Assessments of the patient's state might be performed in a pause after a CR shot. Depending on the rating, the CR stimulation frequency or intensity of the subsequent CR shot may be varied. In particular, for non-invasive application of CR stimulation a non-invasive assessment of the stimulation effect might straightforwardly be realized.
In realistic biological systems intrinsic (model) parameters typically vary over time. These variations may be of complex dynamical nature (see e.g., Timmer et al., 2000;Yulmetyev et al., 2006). To obtain some indication as to whether Protocol C is robust against low-amplitude intrinsic variations of the neuronal firing rates, we added a low-amplitude term I var = A·sin(2π ·f ·t) to the right-hand side of Equation 1a. In the stimulation-free case, I var causes variations of the neurons' firing rates in the order of ±3% and no qualitative changes of the network dynamics (data not shown). For different frequencies f this type of variation does not significantly affect the long-term desynchronization outcome of Protocol C (f = 0.004, 4, and 20 Hz in Supplementary  Figure 1). By the same token, the neuronal firing rates are not significantly altered by the additional periodic force (data not shown).
Note, this is not intended to be a comprehensive study of the impact of periodic forcing of arbitrary frequency on the spontaneous or stimulation-induced dynamics of the model network under consideration. Rather, the slow oscillatory forcing is meant to model slow physiological modulatory processes in an illustrative manner. In the extreme case of f = 0.004 Hz the slow oscillatory modulation acts on the same time scale as a cycle comprising shot and pause and, hence smoothly emulates the step-wise modulation of the CR stimulation frequency in Protocol C.
Conversely, intrinsic variations of sufficient size might naturally mimic variations of the relationship between CR stimulation frequency and intrinsic neuronal firing rates as introduced on purpose in Protocol C. Accordingly, already the purely spaced stimulation without demand-controlled variability (Protocol A) might display some variability of the relationships between intrinsic firing rates and CR stimulation frequency simply due to the intrinsic variability. However, at least with the frequencies 0.004, 4, and 20 Hz in the low-amplitude term I var = A · sin(2π · f · t) added to the right-hand side of Equation 1a, we were not able to observe any substantial improvement of the desynchronizing outcome of Protocol A (Supplementary Figure 2). However, more physiological patterns of firing rate modulations might have a more significant impact on the stimulation outcome of Protocol A. In future studies typical variations of the signals relevant to a particular preclinical or clinical application might be taken into account to further improve desynchronizing short-term dosage regimen. The additional periodic forcing considered here was meant to illustrate the stability of the suggested control approach. However, future studies could also provide a detailed analysis of the interplay of one or more periodic inputs and noise, thereby focusing on stochastic resonance and related phenomena (e.g., Pikovsky and Kurths, 1997;Gammaitoni et al., 1998;Manjarrez et al., 2002;Torres et al., 2011;Bordet et al., 2015;Yu et al., 2016;Guo et al., 2017;Uzuntarla et al., 2017 and references therein). The number of stimulation sites and CR stimulation spatial decay was based on Lysyansky et al. (2011). In accordance to that study, adding more stimulation sites does neither lead to qualitatively different results nor does it improve the stimulation outcome (see Supplementary Figure 3).
The short-term dosage regimen proposed here provides a closed-loop CR stimulation concept that enables to significantly increase the robustness and reliability of the stimulation outcome. Our results motivate to further improve the CR approach by closed loop or feedback-based dosage regimen. Compared to the computationally developed initial concept of demand-controlled CR-induced desynchronization of networks with fixed coupling constants (Tass, 2003a,b), the focus will now be on a feedback-adjusted modulation of synaptic patterns to induce long-lasting therapeutic effects. Currently, clinical proof of concept (phase IIa) is available for deep brain CR stimulation for the therapy of Parkinson's disease (Adamchic et al., 2014) and acoustic CR stimulation for the treatment of chronic subjective tinnitus (Tass et al., 2012a). In addition, promising first in human (phase I) data are available for vibrotactile CR stimulation for the treatment of Parkinson's disease showing pronounced and highly significant sustained therapeutic effects (Syrkin-Nikolau et al., 2018). For the clinical development of these treatments it is mandatory to perform dose-finding studies (phase IIb) to reveal optimal stimulation parameters and dosage regimens, for comparison (see Friedman et al., 2010). The latter are required to get properly prepared for large efficacy (phase III) trials (Friedman et al., 2010). Since CR stimulation modulates complex neuronal dynamics, dosefinding studies are sophisticated, since stimulation parameters as well as dosage patterns have to be chosen appropriately. Selecting appropriate stimulation parameters and dosage regimens by trial and error may neither be effective nor affordable, since it would require a huge number of patients. In contrast, our manuscript illustrates the important role of computational medicine in generating hypotheses for dose-finding studies. Specifically, we show that spacing (i.e., adding pauses in between stimulation epochs) as well as moderate and unspecific parameter variations adapted in the course of the therapy are not sufficient to overcome limitations of CR stimulation. Intriguingly, the combination of both, spacing plus adaptive moderate parameter variation increases the robustness of the stimulation outcome in a significant manner. This computational prediction can immediately be tested in dose-finding studies and, hence, help to optimize the CR therapy, shorten the development time and reduce related costs.

AUTHOR CONTRIBUTIONS
TM, MZ, and PT: Conceived and designed the experiments; TM: Performed the experiments; TM, MZ, and PT: Analyzed the data; TM and MZ: Contributed reagents, materials, and analysis tools; TM, MZ, and PT: Wrote the paper.