Computationally Developed Sham Stimulation Protocol for Multichannel Desynchronizing Stimulation

A characteristic pattern of abnormal brain activity is abnormally strong neuronal synchronization, as found in several brain disorders, such as tinnitus, Parkinson's disease, and epilepsy. As observed in several diseases, different therapeutic interventions may induce a placebo effect that may be strong and hinder reliable clinical evaluations. Hence, to distinguish between specific, neuromodulation-induced effects and unspecific, placebo effects, it is important to mimic the therapeutic procedure as precisely as possibly, thereby providing controls that actually lack specific effects. Coordinated Reset (CR) stimulation has been developed to specifically counteract abnormally strong synchronization by desynchronization. CR is a spatio-temporally patterned multichannel stimulation which reduces the extent of coincident neuronal activity and aims at an anti-kindling, i.e., an unlearning of both synaptic connectivity and neuronal synchrony. Apart from acute desynchronizing effects, CR may cause sustained, long-lasting desynchronizing effects, as already demonstrated in pre-clinical and clinical proof of concept studies. In this computational study, we set out to computationally develop a sham stimulation protocol for multichannel desynchronizing stimulation. To this end, we compare acute effects and long-lasting effects of six different spatio-temporally patterned stimulation protocols, including three variants of CR, using a no-stimulation condition as additional control. This is to provide an inventory of different stimulation algorithms with similar fundamental stimulation parameters (e.g., mean stimulation rates) but qualitatively different acute and/or long-lasting effects. Stimulation protocols sharing basic parameters, but inducing nevertheless completely different or even no acute effects and/or after-effects, might serve as controls to validate the specific effects of particular desynchronizing protocols such as CR. In particular, based on our computational findings we propose a multichannel sham (i.e., inactive) stimulation protocol as control condition for phase 2 and phase 3 studies with desynchronizing multichannel stimulation techniques.


INTRODUCTION
To establish a pharmacological therapy for clinical use, clinical trials are performed in humans that are typically classified into four phases (Friedman et al., 2015): First, in pre-clinical studies pharmacokinetic, toxicity, and efficacy are studied in non-human subjects. In phase I trials, so-called first in human-studies, safety and tolerability of a drug are investigated in healthy volunteers. Phase II trials aim to determine whether a drug can have any efficacy. More specifically, phase IIA trials typically aim at demonstrating clinical efficacy or biological activity ("proof of concept" studies), whereas phase IIB trials are dose-finding studies, performed to reveal optimum dose at which a drug has biological activity with minimal side-effects. Phase III trials investigate effectiveness and the clinical value of a new intervention in a larger patient group. In a randomized controlled trial the effect size of a new intervention is compared with state of the art treatment, if available. Finally, a phase IV trial is a postmarketing surveillance trial, performed e.g., to study whether any rare or long-term adverse effects occur within a much larger patient population and over longer time periods. Individual trials may actually comprise more than only one phase. For instance, there are combined phase I/II or phase II/III trials. Accordingly, given the different purpose of clinical trials, one may also distinguish between early phase studies and late phase trials (Friedman et al., 2015-see above). In principle, this 4-phase pattern also holds for medical technology, e.g., neuromodulation technologies.
Apart from investigating safety and tolerability, it is key to study whether a new therapeutic intervention is superior to pre-existing therapeutic options (Friedman et al., 2015). To this end, one has to take into account non-specific, placebo effects. A placebo effect is a psychobiological phenomenon that causes symptom relief after delivery of inert substances or other types of sham treatment, such as sham surgery or sham stimulation, in combination with verbal instructions suggesting clinical benefit (Price et al., 2008;Benedetti et al., 2011). Note, in clinical trials the terms placebo and sham are basically synonymous, while a placebo typically refers to an inactive substance used in pharmacological trials, whereas a sham stimulation/operation refers to a stimulation/operation without specific therapeutic effect (Price et al., 2008;Benedetti et al., 2011;Friedman et al., 2015). Real placebo effects go beyond spontaneous remission due to the natural history of a disease, regression to the mean induced by selection biases or expectation-related biases of patients and doctors (Benedetti, 2008a).
There are many different placebo effects, caused by different mechanisms and related to different types of interventions and different diseases (Benedetti, 2008b;Enck et al., 2008). Expectation, anxiety, reward, and different types of learning mechanisms may contribute to placebo effects (Benedetti et al., 2011). For instance, according to the Hawthorne effect, patients may simply improve because they are enrolled in a clinical trial (Last, 1983). Placebo treatments can decrease anxiety levels (Vase et al., 2005) and, in general, modulate emotions (Petrovic et al., 2002). Conversely, an inert substance combined with an instruction inducing negative expectations may cause a nocebo effect (Enck et al., 2008), e.g., an increase of pain (Colloca et al., 2008).
Placebo effects may actually be related to objective changes of brain action (Benedetti et al., 2011), e.g., release of endogenous dopamine (de la Fuente-Fernández et al., 2001), changes in brain glucose metabolism (Mayberg et al., 2002) or changes of the activity of specific neuronal populations (Benedetti et al., 2004).
Different types of learning mechanisms, e.g., conditioning, may play important roles in placebo mechanisms (Benedetti et al., 2011). For instance, administration of a placebo after delivery of active drugs may be more effective than placebo administration without the previous experience with the corresponding active drug (Sunshine et al., 1964;Batterman, 1966;Batterman and Lower, 1968;Laska and Sunshine, 1973;Amanzio and Benedetti, 1999;Colloca and Benedetti, 2006). Not only features related to a drug or therapeutic procedure may contribute to placebomediated clinical improvement, but also many other stimuli, related to medical environment, equipment, and personnel (Benedetti et al., 2011). From a clinical trials standpoint it is, hence, important to mimic the entire procedure of treatment delivery as well as possible, since even instructions and rituals of the treatment delivery and procedure may cause actual changes in brain activity that may be the same as those induced by the specific treatment (Benedetti et al., 2011). Accordingly, a vast majority (97%) of surveyed Parkinson's disease (PD) clinical researchers in the United States and Canada believe that even in the case of neurosurgical cell-based and gene therapies for PD double-blind, placebo-controlled trails have to be performed to assess safety and efficacy (Kim et al., 2005;Olanow, 2005). Ninety percentage of PD clinical researchers consider burr holes as justified for sham neurosurgery procedures, and a minority (<22%) even consider penetration of brain tissue to be justified for the neurosurgical sham control (Kim et al., 2005;Olanow, 2005). Hence, even in the case of clinical trials performed according to highest quality standards, e.g., in the field of deep brain stimulation (DBS) (Schuepbach et al., 2013), the comparison between qualitatively different therapeutic regimes, e.g., invasive neuromodulation plus medication vs. medication only, caused debates on whether the study design could reliably rule out placebo effects (Schüpbach et al., 2014).
There is a variety of strategies for the development of sham stimulation protocols. For instance, in the context of transcranial current stimulation a number of studies were devoted to the development of appropriate sham stimulation protocols, since current flow can elicit tingling or itching skin sensations, where different transcranial electrical stimulation methods have different cutaneous perception thresholds (Ambrus et al., 2010). In a comparative transcranial electrical stimulation study, a short-duration active protocol was used as sham, where the active stimulation was turned on only for a brief period, during which stimulation-related unwanted effects/perceptions were elicited (Inukai et al., 2016). Accordingly, the dose should be insufficient, but the patient should get the impression of receiving stimulation. Alternatively, off-target stimulation strategies were developed. In that case, the patient perceives stimuli and/or side effects thereof, but stimulation is directed to targets putatively rendering stimulation ineffective. For instance, for a sham condition for transcranial direct current stimulation a current configuration was chosen such that the current primarily traversed across the scalp, through adjacent pairs of electrodes of opposite polarity, in this way sparing cortical tissue (Richardson et al., 2014;Garnett and den Ouden, 2015). Another offtarget stimulation sham strategy is used in the field of tinnitus, where stimulation tones are delivered at sufficiently detuned pitch compared to the tonal tinnitus, putatively activating brain sites sufficiently remote from the brain regions engaged in the tinnitus-related abnormal neuronal synchrony (Tass et al., 2012a;Adamchic et al, 2017). In contrast, we here consider the situation when clinical constraints are not permitting an off-target stimulation for sham purposes. We hypothesize that an appropriate stimulation pattern may render stimulation ineffective, although the single stimuli are delivered to target sites.
This computational study is dedicated to the development of sham stimulation protocols for desynchronizing multichannel stimulation techniques, specifically coordinated reset (CR) stimulation (Tass, 2003a). The latter was computationally designed to specifically antagonize abnormal neuronal synchrony by desynchronization (Tass, 2003a,b). To this end, sequences of stimuli are administered to different neuronal sub-populations engaged in abnormal neuronal synchronization (Tass, 2003a). In computational studies it was shown that in the presence of spike-timing-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 antikindling effect (Tass and Majtanik, 2006) is caused by a CRinduced reduction of the rate of coincidences which, in turn, induces a decrease of synaptic weights, ultimately shifting the stimulated network 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).
For the clinical development, in particular, of 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 key to compare the effects of CR stimulation with an appropriate sham stimulation protocol in phase II and phase III clinical trials. The sham stimulation protocol should be reasonably similar to the CR stimulation pattern, to prevent patients from being able to distinguish between actual treatment and control. Accordingly, performing double-blind, placebo-controlled trails for non-invasive, sensory multichannel stimulation therapies requires multichannel sham stimulation protocols.
We here computationally develop a multichannel sham stimulation protocol. To this end, we investigate the anti-kindling effect of several multichannel stimulation protocols that share basic features with CR stimulation. We apply the different stimulation protocols to a one-dimensional computational network model with spiking neurons and study the stimulation effects at different levels, ranging from the macroscopic network level, via subpopulations down to the single neuron level. We obtain an inventory of qualitatively different stimulation effects elicited by the different stimulation protocols. Intriguingly, we found an inert stimulation protocol which caused only weak acute and hardly any long-lasting effects. The latter is a potential sham candidate to be tested for clinical studies in the context of desynchronizing sensory multichannel stimulation techniques.

MATERIALS AND METHODS
In this section, we describe the equations used to model the dynamics of our one-dimensional neuronal network, the plasticity of the synapses, and the different stimulation protocols as well as the data analysis methods.

Neuronal Network
The model we use is a one-dimensional ring composed of N spiking Hodgkin-Huxley neurons which interact via strong excitatory short-range and weak inhibitory long-range synapses (Popovych and Tass, 2012). The membrane potential V i of the i-th neuron (i = 1: N) is given by: where C denotes the membrane capacitance, the injected constant currents I i are uniformly distributed (I i ǫ [I 0 − I , I 0 + I ]) and F i represents the current induced by an external stimulation signal (see section Simulation Details for more details). The voltage-dependent rate constants α x and β x of the time-varying ion gate variables x ǫ {m, n, h} are given by α m (V) = (0.1V The coupling term S i in Equation (1) stands for the weighted ensemble average of all postsynaptic currents received by neuron i from the other neurons in the network and can be given in terms of the synaptic variable s j as: where N is the number of neurons within the network, V r,j the reversal potential of the synaptic coupling between neurons j and i and c ij is the synaptic coupling from neuron j to neuron i. M ij has the form of a Mexican hat (Wilson and Cowan, 1973;Dominguez et al., 2006;de la Rocha et al., 2008) and determines the type of the neuronal connection between neurons i and j: M ij < 0 represents an inhibitory coupling, M ij > 0 an excitatory coupling. The value of M ij determines the distance dependent maximal strength between those neurons: with σ 1 = 3.5 and σ 2 = 2.0 as in Popovych and Tass (2012), and d ij = d ·min i − j , N − i − j is the shortest distance between neurons i and j. To avoid boundary effects, the neurons form a one-dimensional ring. Therefore, the shortest distance between the neurons with indices 1 and N is 1 instead of N − 1. The lattice distance between two adjacent neurons is given by: with d 0 the length of the neuronal chain. Values used in this study are N = 200, C = 1 µF/cm 2 , maximum conductance per unit area for the sodium, potassium and leak currents g Na = 120 mS/cm 2 , g K = 36 mS/cm 2 , g l = 0.3 mS/cm 2 , sodium, potassium, and leak reversal potentials V Na = 50 mV, V K = −77 mV, and V l = −54.4 mV, reversal potential for excitatory respectively inhibitory coupling V r,j = 20 mV respectively−40 mV. For the constant injected current we used I 0 = 11.0 µA/cm 2 and I = 0.45 µA/cm 2 . The length of the neuronal chain is defined as d 0 = 10.

Spike Timing-Dependent Plasticity
The dynamical synaptic weights c ij are influenced by the precise timing of the pre-and postsynaptic spikes and are updated in an event-based manner every time a neuron spikes. This is realized by adding δ · c ij to the excitatory and −δ · c ij to the inhibitory synaptic weights c ij with learning rate δ > 0 every time neuron i or j spikes. According to the spike timing-dependent plasticity (STDP) rule (Bi and Poo, 1998) the change in synaptic weight is given by (Popovych and Tass, 2012): with t ij = t i − t j and t i is the spike time of the postsynaptic neuron i and t j the spike time of the pre-synaptic neuron j.

External Stimulation
The aim of this study is to compare the effects of different stimulation algorithms on the neuronal connectivity as well as on the synchronization of the neuronal activity. In this section, we describe six different stimulation protocols as well as a no-stimulation control protocol to investigate the influence of the different stimulation protocols on the connectivity and synchronization.

Stimulation Implemented in the Model
Each stimulation onset induces single brief excitatory postsynaptic currents, with spatial spread in the network given by a quadratic spatial decay profile: where d is the lattice distance between adjacent neurons (Equation 4), i − x k the difference in index of neuron i and index x k of the neuron at stimulation site k, and σ d the spatial decay rate of the stimulation current (Popovych and Tass, 2012). The total stimulation current induced in neuron i is given by: with the excitatory reversal potential V r = 20 mV, V i (t) the membrane potential of neuron i (Equation 1), K the stimulation intensity and D the spatial decay profile (Equation 6). G s,k is the stimulation (at stimulation site k) evoked time-dependent normalized conductance of the post-synaptic membranes defined by α-functions. Since in this study the minimal time difference between two stimulation onsets within the network is not restricted to T s /N s = 4 ms, we have adapted G s from Popovych and Tass (2012) by allowing a summation of two stimulation evoked time-dependent normalized conductances of the post-synaptic membranes if the two stimulations occur within a certain time interval. We have set this time interval to 2 · T s /N s = 8 ms, since at 8 ms the value of the αfunction is marginal (only 0.02% of the peak value). Our adapted stimulation (at stim site k) evoked time-dependent normalized conductance of the post-synaptic membranes is now given by: Here t n k is the onset of the nth activation of the kth stimulus site, τ = T s /(6N s ) represents the time-to-peak of G s,k , and min(t 1 , t 2 ) is the minimum value of t 1 , t 2 and thus represents the earliest time event of t 1 and t 2 .

Stimulation Signal Features
As a control condition we use the situation where no stimulation signal is applied and thus no stimulation current is delivered Frontiers in Physiology | www.frontiersin.org to the N neurons: F i = 0 for all i ∈ {1, .., N}. We denote the control condition as no stimulation (no-stim) protocol. The other, active stimulation protocols last 128 s and consist of stimulation ON and stimulation OFF cycles (briefly ON-and OFF-cycles), of duration T s each. N s equidistantly spaced stimulation sites are activated exactly once during one ON-cycle. After three consecutive ON-cycles, two OFF-cycles follow, before the next three ON-cycles take place. At the end of the stimulation period (briefly stim-on period) each stimulation site will have been activated exactly 4,800 times (=4,800 ON-cycles). In this study, each cycle lasts for T s = 16 ms, and for all stimulation signals the same N s = 4 equidistantly spaced stimulation sites are activated. The four stimulation sites are located at the neurons with index 25, 75, 125, and 175.
The first of the six stimulation signal approaches is the purely periodic multichannel stimulation (PPMS). The first stimulation onset of the k-th stimulation site, t 1 k , is drawn randomly (with equal distribution within an ON-cycle), the next stimulation onsets occur exactly a multiple of T s later: t n k = t 1 k + n · T s (with k ∈ {1, 2, . . . , N s } and n ∈ N). Another feature of this PPMS approach is that all four stimulation sites are activated simultaneously (see Figure 1A), which implies that t n k = t n 1 for all stimulation sites k.
The correlated multichannel noisy stimulation (CMNS) activates, like the PPMS, all four stimulation sites simultaneously, but the stimulation onsets are no longer periodic, but rather noisy: For each ON-cycle the stimulation onset is drawn randomly, with equal probability within the ON-cycle ( Figure 1B). In this case, the four stimulation sites are also active simultaneously: t n k = t n 1 for all stimulation sites k. For the uncorrelated multichannel noisy stimulation (UMNS) the stimulation onsets t n k are determined for each stimulation site k separately. To this end, for each stimulation site k and ON-cycle n the stimulation onset t n k is drawn randomly, with equal distribution within the n-th ON-cycle. The random processes that generate the stimulation onsets t n k for the different stimulation sites k ∈ {1, 2, . . . , N s } are completely uncorrelated between stimulation sites and, hence, typically do not coincide ( Figure 1C).
The other three stimulation protocols are different variants of coordinated reset (CR) stimulation (Tass, 2003a,b). For CR stimulation, within each ON-cycle the activations of the N s = 4 different stimulation sites are equidistantly spaced in time, with a time shift of T s /N s (Figures 1D-F). The different stimulation onsets are at the beginning of the ON-cycle, 4, 8, and 12 ms later. The order in which the stimulation sites are activated exactly once during an ON-cycle is called a stimulation sequence, briefly sequence. For the RVS CR the sequence randomly changes for one ON-cycle to the next. In contrast, for the fixed CR the same sequence is maintained for all ON-cycles ( Figure 1E). For the slowly varying sequences (SVS-n) CR a sequence is applied during n ON-cycles, before randomly switching to another sequence which is, in turn, used for the next n ON-cycles etc (Zeitler and Tass, 2015). Figure 1F illustrates SVS-4 CR, where each sequence is repeated four times before another sequence is drawn randomly. In the SVS CR protocols of this study we will only apply n = 100 consecutive repetitions of a sequence before The rapidly varying sequences (RVS) CR stimulation activates the stimulation sites in a random temporal order at times n · T s /4 (n ∈ {0, 1, 2, 3}) within each stimulation ON-cycles. (E) The fixed CR stimulation approach employs a fixed sequence during the entire CR-on period, and stimulation onsets times within the stimulation ON-cycles are at multiples of n · T s /4 (n ∈ {0, 1, 2, 3}). (F) In this example, for the sake of illustration, the slowly varying sequences (SVS) CR stimulation repeats the first sequence four times, before a new sequence is drawn. Different stimulation sequences are randomly drawn with equal probability. For the remainder of this study each sequence is repeated 100 times and not just four times (SVS-100 instead of SVS-4). Stimulation onsets are as for the RVS CR at n · T s /4 (n ∈ {0, 1, 2, 3}) within the stimulation ON-cycles. All stimulation protocols have four equidistantly distributed stimulation sites, stimulation ONand OFF-cycles of 16 ms, a total stimulation on period of 128 s, in which three ON-cycles are alternated by two OFF-cycles. No stimulation is applied during OFF-cycles. Stimulation ON-and OFF-cycles are separated by dashed vertical lines. A change of color indicates a change of the pattern of stimulus onsets compared to the previous ON-cycle.
we draw the next sequence. Therefore, we will use the term SVS instead of SVS-100 CR for the remainder of this work.

Simulation Details
A simulation contains four different simulation periods: first an initializing period of 2 s, then a 60 s period with STDP and without external stimulation (denoted as STDP-only period) which is followed by a stim-on period of 128 s in which the network receives an external stimulation and after withdrawal of the stimulation follows a stimulation off period (denoted as stim-off period) of 128 s. During the initializing period -the only period without STDP -a network is built by drawing random numbers for each of the N = 200 neurons from uniform distributions for the injected constant current I i ∈ [I 0 − I , I 0 + I ], for the membrane potential V i ∈ [−65, +5] ms, for the time-varying ion gate variables x i ∈ [0, 1] and the synaptic variable s i ∈ [0, 1]. The initial synaptic weights c ij are drawn from a Gaussian distribution with mean 0.5 µA/cm 2 and standard deviation 0.01 µA/cm 2 . During the initializing period the neuronal network evolves without influences of external stimulation signals or STDP.
For this selected initial distribution of the synaptic weights, during the STDP-only period the network can develop into a strongly connected network with strongly synchronized neuronal activity. The time at which the first stimulation signal is delivered to the network is defined as t = 0 s. During the stim-on period the stimulation signals are applied to the network as described in section Stimulation Implemented in the Model. After 128 s no stimulation signals are applied anymore and the evolution of the network is monitored for another 128s (stim-off period). After going through this whole process, the procedure is repeated from t = 0 s on for a different stimulation intensity, K, and/or for a different stimulation signal approach.
From previous studies Tass, 2015, 2016) we know that different initial network conditions and different stimulation signal realizations have an influence on the antikindling effects. Therefore, we draw eleven different sets of initial values from the distributions as described above at the start of the initializing period. For each of these eleven different initial network conditions we generate one realization for each stimulation signal approach. In the remainder of this study the combination of a stimulation signal realization with one set of initial network conditions is referred to as sample. All simulations were executed in Matlab R2007a. The differential equations were solved by the built-in function ODE45 with a relative tolerance of 10 −5 .

Data Analysis
In this computational study we came up with a larger set of acute and long-lasting effects of different stimulation protocols applied to the plastic neural network as described in section Neuronal Network. We are particularly interested in whether the different stimulation protocols might induce qualitatively different antikindling effects. In this section, we will discuss the methods used to investigate the neuronal connectivity at several levels, the synchronization and a quantification of the acute stimulation and acute after-effects, as well as the phase resetting and entrainment induced by the different stimulation protocols. Matlab R2015a was used for the data analysis and for plotting the results.

Connectivity
In this study, the synaptic weights can change according to the STDP-rule (see Equation 5). On the network level the dynamics of the synaptic connectivity is monitored by the synaptic weight averaged over all synapses within the network: where N is the number of neurons within the network, sign(M ij ) is negative for inhibitory synapses and positive for excitatory synapses (with M ij defined as in Equation 3), and c ij is the synaptic coupling strength from neuron j to i. A decrease of C av over time may indicate that there is mainly a decrease in the average excitatory synaptic weights or an increase in average inhibitory weights or a combination of both. To unravel the contributions of the excitatory and the inhibitory synaptic weights, we introduce the average excitatory synaptic weight and the average inhibitory synaptic weight where [z] + stands for the half-wave rectification operation ([z] + = z if z > 0 and [z] + = 0 otherwise), N EE is the number of excitatory synapses, and N II is the number of inhibitory synapses within the whole neuronal network. On the neuronal level the connectivity matrices are analyzed at the end of the stim-on and at the end of the stim-off period, since we are interested in the acute and the long-lasting effects. Instead of considering the c ij -values, we first multiplied each c ij -value by the sign-function of the Mexican Hat, sign(M ij ) (see Equation  3 for M ij ). This allows to recognize the type of each synapse in a color plot of the connectivity matrix (negative values for inhibitory and positive values for excitatory synapses). For each synapse eleven c ij -values exist since each stimulation protocol is applied at stimulation intensity K to eleven different initial networks. By determining the median and the inter-quartilerange (IQR) of these eleven c ij -values for all i, j ∈ {1, 2, . . . , N] we obtained a larger reduction of the amount of data. Unfortunately, due to the different initial network conditions the general connectivity pattern induced by a stimulation protocol does not straightforwardly reveal at which locations neurons are coupled mainly by bidirectional weak as opposed to unidirectional strong synapses.
To get a better idea of what happens due to a stimulation signal, we first sorted for each sample the c ij (t) values such that the strongest synaptic weight between two neurons i and j is placed in the lower right triangle at location (i ′ ,j ′ ) = (min(i,j), max(i,j)) and the weaker synaptic weight between those two neurons i and j in the upper left triangle of the connectivity matrix at location (i ′ ,j ′ ) = (max(i,j), min(i,j)). This is repeated for all combinations of two neurons. After that the result is multiplied by sign(M ij ). We call this newly defined matrix the sorted connectivity matrix. After calculating the median of the eleven sorted connectivity matrices the general pattern is more evident for all stimulation protocols and the IQRs provide a clearer picture than for the unsorted connectivity matrices. By sorting the data like this, one can determine from the median if the synapses are in general strong in both directions or just in one direction or maybe even weak in both directions between two neurons. Furthermore, the IQR shows which type of synapses and at which spatial location the largest differences occur as an effect of the different samples. Note that the location (i, j) in this sorted connectivity matrix does not indicate that neuron i is the postsynaptic neuron and neuron j the pre-synaptic one as is the case for the unsorted connectivity matrices, but that if i < j (lower right triangle of the sorted matrix) than c sorted ij = max c ij , c ji with c ij , c ji elements of the unsorted connectivity matrix. For the upper left triangle of the sorted matrix where i > j this means that c sorted ij = min c ij , c ji .

Synchronization
The effect of the external stimulation signals on the strongly synchronized neuronal activity is investigated on the network level by the order parameter R (Haken, 1983;Kuramoto, 1984) defined by where Φ(t) is the circular mean phase of the entire group of N neurons in the network, and (14) is a linear approximation of the phase of neuron j between its mth and (m + 1)-th spikes at spike times t j,m and t j,m+1 (Rosenblum et al., 2001). The minimum value of R is zero and indicates a complete lack of in-phase synchronization, whereas its maximum value (R = 1) indicates perfect in-phase synchronization. For our analysis, we have calculated R, Φ, and all ϕ j at each ms of the stim-on and stim-off period. In case an R-value is shown at a certain time instance t the R-values of the preceding 5 s period are averaged and denoted as R av at t.
On a mesoscopic level the amount of synchronization of the k-th subpopulation can be defined as with N k = 49 the number of neurons within subpopulation k, and ϕ j (t) the linear approximation of the phase of neuron j in subpopulation k at time t as determined by Equation (14). Using this definition, we can determine the acute stimulation effect on the k-th subpopulation by and the acute after-effect on the k-th subpopulation by A negative outcome indicates a synchronizing effect, a zero outcome indicates that there is no acute stimulation effect or after-effect, respectively, and a positive outcome means a desynchronizing effect within the k-th subpopulation due to application of that particular stimulation protocol.

Stimulus-Locked Phase Dynamics of a Subpopulation
To shed more light on the mechanisms of the different stimulation protocols and reveal, e.g., phase resetting or entrainment processes, we investigate the stimulus-locked dynamics on a mesoscopic scale. We do this by considering subpopulations of neurons as given by their proximity to the N s = 4 different stimulation sites located at neuron indices 25, 75, 125, and 175. This implies that each subpopulation k contains N k = N/N s − 1 = 49 neurons, since the subpopulations are separated by a neuron which has an equal distance to the stimulation sites in the two neighboring subpopulations, e.g., the neuron at index 50 has the same distance to the stimulation site at neuron index 25 as well as to the stimulation site at neuron index 75 and is therefore excluded from subpopulation 1 as well as from subpopulation 2. The mean phase . We focus on the distribution of the stimulus-locked phase dynamics within a time window Wof 32 ms before and up to 32 ms after each of the L stimulation onsets t n k (n = 1, .., L) for each subpopulation k separately (Tass, 2003c) k t n k + t (mod 2π) n=1,...,L .
These distributions of the stimulus-locked phase dynamics of a subpopulation will be displayed in color plots as a function of t ∈ W. To quantify the amount of stimulus locking of the phase dynamics of subpopulation k, we use the resetting index E k (t) (Tass, 2003c; see also Tallon-Baudry et al., 1996) as given by In case of a phase entrainment (a permanent stimulus-locking of the phase dynamics) E k ( t) results in a uniform distribution throughout the window W = [−32 , +32] ms. For a stimulusinduced phase reset, E k ( t) will be small in the period before stimulus onset (corresponding to a uniform phase distribution) and will increase after stimulus onset, reflecting the emergence of a unimodal phase distribution (Tallon-Baudry et al., 1996;Tass, 2003c). For a stimulus-locked disruption of a phase entrainment the phase distribution will be unimodal before stimulus onset and decrease due to stimulus onset. Since all our stimulation protocols have the same ratio of ON: OFF cycles, namely 3:2, we can split each set of the L stimulation onsets t n k (n = 1, .., L) for each subpopulation k into three subsets. The first subset contains only the L/3 stimulation onsets t n k (n = 1, 4, .., L − 2) within the first of the three consecutive ON-cycles, the second one contains the L/3 stimulation onsets t n k (n = 2, 5, .., L − 1) during the middle of each block of three consecutive ON-cycles and the last subset contains the L/3 stimulation onsets t n k (n = 3, 6, .., L/3) within the third ON-cycle of each three consecutive ON-cycles. By using these three subsets separately instead of the three subsets together, the corresponding three resetting indices can be determined in a similar way as in Equation 19. This detailed analysis shows how the phase reset and entrainment processes evolve depending on the rank order of the onset after the OFF-cycles.

Statistics
To test whether an increase of the median of the eleven C av or R av values obtained by stimulations with protocol A at stimulation intensity K A compared to stimulation protocol B applied at K B is statistically significant, we used the Matlab R2015a built-in left-sided Wilcoxon rank sum test with significance level α: [p, h] = ranksum(A,B, "alpha",0.05, "tail, " "left, " "method, " "exact"). This test is equivalent to the Mann-Whitney U-test. The built-in right-sided Wilcoxon rank sum test tests if a decrease in medians is statistically significant. In this study, we use significance level α = 0.05 and n A = n B = 11 samples unless stated otherwise.

Acute Effects
In this section, we study the acute stimulation effects of the different spatio-temporally patterned stimulation protocols at K = 0.25 by comparing their effects on the connectivity as well as the synchronization at the end of the stim-on period (t = 128 s).
External stimulation signals can change the network's connectivity (due to STDP) and the amount of synchronized neuronal activity. The control signal (no-stim) has no influence on the average synaptic weight, C av (Figure 2A) and on the synchronization of the population activity (shown by the order parameter R av in Figure 2B). For the same initial network conditions the other stimulation signals applied at intensity K = 0.25 show an acute reduction of C av as well as of R av (see results at t = 128 s in Figures 2A,B). Only the correlated multichannel noisy stimulation (CMNS) results in an increase of C av and a small reduction of R av . The boxplots in Figures 2C,D show that these acute effects are representative for all 11 samples: compared to the control signal all stimulation protocols induce a statistically significant decrease of C av and R av , except the CMNS which induces also a statistically significant decrease of R av , but a statistically significant increase of C av (left-sided Wilcoxon rank sum test with α = 0.05). The corresponding p-values are given in Supplementary Table 1. The CMNS-induced increase of C av is the result of a statistically significant increase of the average excitatory synaptic weight c EE in combination with a statistically significant decrease of the average inhibitory synaptic weight c II (see Figures 3A,B for one set of initial network conditions and Figures 3C,D  for all 11 samples). The other stimulation protocols show an opposite behavior: the combination of a statistically significant decrease of c EE with a statistically significant increase of c II (Figure 3) explains the statistically significant decreases of C av (Figure 2C).
According to the median of the unsorted connectivity matrices induced by the control protocol (no-stim) there are many strong excitatory (and inhibitory, respectively) synapses without a spatial pattern in relation to the stimulation sites at locations (25,25), (75,75), (125,125), and (175,175) (upper panel of Figure 4A). Note, the excitatory (respectively inhibitory) synapses have positive (respectively negative) values. The corresponding IQRs show that there are large differences between the synaptic weights induced by the different samples for almost all synapses (bottom panel of Figure 4A). By first sorting for each sample the c ij (t = 128 s) values such that the strongest synaptic weight between two neurons i and j is placed in the lower right triangle and the weaker synaptic weight between those two neurons i and j in the upper left triangle (see section Connectivity for more details) and then determining the median and IQR, it becomes clear that in general the control protocol results in only one strong synapse (|c ij | ∼ = 1) between two neurons and in the reverse direction the synaptic weight is approximately zero (upper panel of Figure 4B). The corresponding IQR shows that the (sorted) synaptic strengths are rather independent of the actual sample (bottom panel of Figure 4B). So, the general pattern is that due to the control signals the synaptic weights between any two neurons are unidirectional (strong in one direction, weak in the reverse direction). This is similar for the synaptic weights at the beginning of the stim-on period (t = 0 s, result not shown). Figure 4C shows that if all N s = 4 stimulation sites are activated simultaneously and periodically (PPMS) more bidirectional strong inhibitory synapses exist at the end of the stim-on period and on the other hand some excitatory synapses are bidirectional and weak. Differences between the results induced by different samples are mainly found in the indices of the neurons which have bidirectional strong inhibitory synapses and in the indices of those neurons which have bidirectional weak excitatory synapses ( Figure 4C).
In case all N s = 4 stimulation sites are simultaneously stimulated in a noisy manner (CMNS), the median of the sorted connectivity matrices shows that besides the induction of some bidirectional strong excitatory synapses the median is similar to the median of the control signals (Figures 4B,D). However, the IQRs vary: compared to the no-stim condition (Figure 4B), different samples with CMNS will induce bidirectional strong excitatory synapses between different neuronal locations and also the induced bidirectional weak inhibitory synapses differ in locations ( Figure 4D).
If the stimulation sites are not activated simultaneously but sequentially and still noisy (UMNS) then the number of unidirectional strong excitatory neurons decreases in such a way that they exist between neighboring neurons and between neurons nearby the stimulation sites, but not between more distant neurons ( Figure 4E). Further, all inhibitory synapses are bidirectional and strong. Differences between samples are mainly found in the size of the region with strong excitatory synapses nearby the stimulation sites and between neighboring neurons and in the occurrence of bidirectional strong excitatory synapses nearby the stimulation sites ( Figure 4E).
To some extent, RVS CR is similar to UMNS, but stimulation onsets are equidistantly spaced at 0, 4, 8, or 12 ms after onset of the ON-cycle instead of at random time instances as for UMNS. The median and IQR of the sorted connectivity matrices induced by RVS CR (Figure 4F) is similar to the results induced by UMNS. Small differences are obtained for the RVS CR compared to the UMNS: the unidirectional strong excitatory neurons exist in a smaller neighborhood of each neuron as well as of the stimulation sites. These small differences are in agreement with the decrease of excitatory synaptic weights in Figure 3C and with the smaller C av -values induced by RVS CR compared to those induced by UMNS in Figure 2C.
The fixed CR activates the stimulation sites sequentially at 0, 4, 8, and 12 ms after onset of the ON-cycle, while the sequence does not change during the stim-on periods. This repetition of sequence has several effects (compare Figures 4F,G). Most prominent is the fact that due to the fixed CR ( Figure 4G) the neurons close to a stimulation site are not so strongly coupled as in case of the RVS CR ( Figure 4F). Furthermore, in general, the neurons halfway between two stimulation sites are not strongly coupled any more with their direct neighbors. Due to the repetition of the sequence strong excitatory couplings between neurons surrounding two consecutively activated stimulation sites remain present. For the eleven different samples (each with another sequence), this leads to some large IQRs for couplings between more distant excitatory synapses. By changing now and then the sequence during the stim-on period (SVS CR), these large IQRs between more distant excitatory neurons disappear ( Figure 4H).
In Figure 2D we have seen how the amount of synchronization of the complete network R av , changes under influence of the different stimulation approaches at stimulation intensity K = 0.25. On a mesoscopic level we investigate the desynchronizing effect in the four subpopulations containing N k = 49 neurons near the stimulation sites. For each stimulation protocol Figure 5A illustrates at K = 0.25 the distribution of the N s · n = 44 determined acute stimulation effects. All stimulation approaches induce a statistically significant acute stimulation effect whereby the weakest effect is induced by the CMNS (right-sided Wilcoxon rank sum test, α = 0.05, n = 44; see Supplementary Table 2 for p-values). The order of the stimulation protocols in having a stronger desynchronizing FIGURE 3 | CMNS induces reversed changes in average excitatory, respectively, inhibitory synaptic weights (c EE , respectively c II ) compared to other stimulation protocols. All stimulation protocols, except the CMNS, show an acute and long-lasting decrease in excitatory synaptic weights (A) as well as an acute and long-lasting increase in inhibitory synaptic weights (B) of the same network as studied in Figures 2A,B. CMNS induces an increase of the average excitatory synaptic weights and a decrease of the average inhibitory synaptic weights. The red horizontal bar represents the stim-on period. Boxplots confirm the results for excitatory (C) as well as for inhibitory synaptic weights (D) (11 samples). p-values of the left-sided Wilcoxon rank sum test are given in Supplementary Table 1. All stimulations are applied at K = 0.25. effect is clear for the macroscopic measure R av : only RVS CR and fixed CR have a similar effect on R av . On the mesoscopic level this order is slightly different: the PPMS turned out as good as RVS CR, while fixed CR caused a better desynchronization of the neuronal activity than RVS CR (compare Figures 2D, 5A; onesided Wilcoxon rank sum test with α = 0.05; see Supplementary  Table 1 respectively 2 for p-values corresponding by Figure 1D respectively Figure 5A). Figure 6 shows the raster plots and spike counts induced by the different stimulation protocols (K = 0.25) for the same samples as used in Figures 2A,B, 3A,B, 4 during the last 100 ms of the stim-on period. Activating the stimulation sites simultaneously does not cause a pronounced desynchronization (compare the results for no-stim with K = 0.25 in Figures 6A,B). A sequential random activation of the stimulation sites can broaden the synchronized spike-volley ( Figure 6C). RVS-CR stimulation counteracts in-phase synchronization, typically by causing cluster states, where the network forms several phaseshifted (synchronized) subpopulations (see e.g., Figure 6D). A more pronounced overall desynchronization is achieved by means of the fixed CR stimulation and the SVS-CR stimulation (e.g., Figures 6E,F).
The spike counts suggest that activating the stimulation sites simultaneously can result in stronger synchronization of the activity of the whole network, while sequential stimulation can divide the network in several synchronized, but mutually phaseshifted subpopulations, which in turn causes a pronounced overall (i.e., close to uniform) desynchronization, as reflected by R av . Since each single stimulus synchronizes the nearby neurons, while desynchronizing the entire neuronal network, on a macroscopic scale the in-phase synchronization (and hence R av ) may vanish, while the order parameters of the different subpopulations may still attain high values. Put otherwise, the acute stimulation effect will be weaker on the mesoscopic than on the macroscopic level (e.g., for UMNS and RVS CR). For the PPMS and the CMNS, which activate all stimulation sites simultaneously, the order parameter of each subpopulation represents the order parameter of the macroscopic network quite well and, hence, the acute stimulation effect is comparable on the mesoscopic and macroscopic level.
Our analysis continues on the mesoscopic level to study the mechanisms by which the different stimulation protocols influence the synchronization. For this we use the cross-trial analysis, which investigates the subpopulations' phases timelocked to the corresponding stimulus onset, averaged over all stimulus onsets of stimuli delivered to a particular subpopulation during the stim-on period. For subpopulation 2 (comprising neurons 51-99), Figures 7A-F show a clear difference between the stimulation protocols with periodic delivery pattern, PPMS, fixed CR and SVS CR, and those which have no strictly periodic stimulus delivery pattern (CMNS, UMNS, and RVS CR). The latter protocols cause phase resets of the stimulated subpopulations: Before stimulus onset (for t < 0 m) the phase distributions are close to uniform, whereas after stimulus delivery a stereotypical restart of the subpopulation phase occurs, which is reflected by the emergence of a pronounced peak of the distribution (for t > 0 ms in Figures 7B-D). In contrast, in case of the periodic stimulus patterns, PPMS, fixed CR and SVS CR, a pronounced peak of the phase distribution increasingly vanishes in the absence of stimulation (for t < 0 m) and re-occurs after stimulus delivery (for t > 0 ms) (Figures 7A,E,F). The corresponding resetting indices E 2 ( t) show that the non-periodic stimulation protocols induce a phase reset (Figures 7B-D): Following stimulus onset, the initially quite homogeneously distributed subpopulation phase turns into a unimodal phase distribution. In contrast, for the periodic stimulation protocols, the resetting indices display a completely different time course: Starting at a large value, they first decrease, then re-increase due to the first stimulus (at t = 0 ms) and further increase due to the subsequent stimulus ( t ∼ 16 ms) (Figures 7A,E,F). More precisely, at t ∼ 16 ms in only 2/3 of all stimulation onsets there is indeed a stimulation, since the third ON-cycle of each block of three consecutive ON-cycles is not directly followed by an ON-cycle but instead by two OFFcycles. This implies that after three activations of the stimulation site within the subpopulation with an inter-stimulus-interval of 16 ms, the next inter-stimulus-interval is equal to 48 ms. Accordingly, by sorting the stimulation onsets according to their  (16), acute after-effects by Equation (17). Statistically significantly stronger effects are determined by the right-sided Wilcoxon rank sum test with significance level α = 0.05 (see for p-values Supplementary Table 2). All stimulations are applied at K = 0.25. order in the blocks of three consecutive ON-cycles, reveals the effect of the different stimuli. The first of the three stimuli destroys most of the phase entrainment which is present during the OFF-cycles and then builds up the entrainment up to a lower level than the initial entrainment ( Figure 7G). The second of the three stimuli further increases the extent of phase entrainment (Figure 7H), and the third stimulus finally increases the phase entrainment to the initial level ( Figure 7I). The maximum values of the resetting indices of the periodic stimulation protocols increase within a block of three consecutive stimuli (i.e., ONcycles) from one stimulus to the subsequent one, indicating a further increase of the phase entrainment between stimulated subpopulation and corresponding stimulus train (compare the maxima in Figures 7G-I). Accordingly, the effect of the three consecutive periodic stimuli on the phase builds up within each block. In contrast, each subsequent non-periodic stimulus causes a new phase reset (Figures 7G-I).
Sequential stimulation patterns (UMNS, RVS CR, fixed CR, SVS CR) cause a more pronounced reduction of synchrony and the mean synaptic weight on the network level (Figure 2) as well as on the subpopulation level (Figures 3, 5) compared to simultaneous stimulation patterns (PPMS, CMNS). On the level of the individual neurons the sequential stimulation patterns induce strong bidirectional inhibitory synapses and subpopulations with strong unidirectional excitatory coupling of the neurons near the stimulation sites, whereas the simultaneous stimulation protocols only induce some strong bidirectional inhibitory synapses and no spatial pattern of the strong unidirectional excitatory synaptic weights (compare Figures 4E-H with Figures 4C,D).
Another important aspect refers to the sequential arrangement of stimulation sequences. Repeating a sequence many times in a row (by fixed CR or SVS CR) causes a more pronounced reduction of C av than random variations of the sequences (by UMNS or RVS CR) (Figure 2). Sufficient repetition of stimulation sequences results in a stronger reduction of c EE (Figure 3) caused by the fact that less neurons are coupled within a subpopulation around each stimulation site together with a weaker synaptic strength for those who are coupled ( Figure 4G,H compared with Figures 4E,F). Although on the network level no clear difference in the amount of desynchronization is observed (RVS CR and fixed CR induced similar C av -values; Figure 2D), on the subpopulation level the repetitive sequences induce a stronger acute stimulation effect (Figure 5). This suggests that different mechanisms of action may cause different types of macroscopic desynchronization: desynchronization between subpopulations which themselves can still be highly synchronized as opposed to a more overall, uniform desynchronization, affecting the subpopulations as well. According to the raster plots in Figure 6, a more uniform desynchronization is typically observed as a result of repetitive sequence CR (fixed CR and SVS CR). Figure 7 shows that on the subpopulation level there is an entrainment between the stimulation signal and the subpopulation activity in the networks exposed to repetitive sequence CR. In contrast, phase resets are the salient mechanism of those CR variants without sequence repetition.
Intriguingly, all stimulation protocols except the CMNS protocol induce a decrease of c EE and an increase in c II compared to the no-stim protocol (Figure 3). In contrast, the CMNS protocol induces exactly the opposite (Figure 3), resulting in an increase of C av , in contrast to all other active stimulation protocols (Figure 2C). Since also the median of the sorted connectivity matrix induced by CMNS is similar to the one of the control (no-stim) signal and not for the other stimulation protocols (Figures 4B,D), the CMNS seems to be the best candidate for a sham stimulation protocol for sensory CR stimulation. While the stimulation shares some perceptual features with CR stimulation, the acute effects of CMNS on connectivity and desynchronization are minimal.

Long-Term Effects
It is a key goal for non-invasive neuromodulation techniques to cause long-lasting, sustained effects that persist after cessation of stimulation. Sufficiently pronounced long-lasting effects may open up the possibility to deliver stimulation only regularly or occasionally e.g., for a few hours only, to maintain substantial relief. Therefore, we are particularly interested in the effects of the different stimulation protocols after withdrawal of stimulation. In this section, we study these long-term effects at the end of the stim-off period (t = 256 s) unless stated otherwise. Again, the stimulations during the stim-on period were applied at intensity K = 0.25.
For the same initial network configuration, a stimulation epoch of duration t = 128 s can have different acute and longterm effects on the average synaptic weight C av , depending on the stimulation protocol selected (Figure 2A). To focus on the longterm effects, the distributions of the C av (t = 256 s) at the end of the stim-off period are shown in Figure 8A for the different stimulation protocols. The SVS CR protocol induces the greatest reduction of C av (t = 256 s), while the CMNS induces even a small increase. These long-term effects ( Figure 8A) are qualitatively comparable with the acute effects (Figure 2A) and statistically significant. None of the stimulation protocols induces C av effects statistically equivalent to the no-stim protocol. However, CMNS induces R av (t = 256 s) values similar to the control signal (no-stim) (Figure 8B), and there are no statistically significant differences. In contrast, the SVS CR stimulation induces the smallest R av values at t = 128 s (acute effect; Figure 2B) as well as at t = 256 s (long-term effect; Figure 8B). Therefore, the CMNS protocol turns out to be the best candidate for sham stimulation. The p-values obtained by the left-sided Wilcoxon rank sum tests for the long-term effects as shown in Figures 8A-D are given in Supplementary Table 3.
The strongest C av is induced by the CMNS protocol ( Figure 8A) and can be disentangled in the strongest average excitatory synaptic weight c EE of the tested protocols ( Figure 8C) in combination with the weakest inhibitory synaptic weight c II of the tested protocols ( Figure 8D). The weakest C av as shown in   (Figure 8C) and the strongest c II (Figure 8D) of the tested protocols. In general, a weaker C av value obtained by a particular stimulation protocol compared to another protocol can be explained by a significantly reduced c EE and increased c II of the first stimulation protocol except for the comparison between fixed CR and SVS CR. These two protocols induce a similar c II and thus only the stronger reduced c EE induced by the SVS CR stimulation contributes to the smaller C av value compared to the fixed CR stimulation (Figures 8A,C,D).
Comparing the medians and IQRs of the unsorted connectivity matrices at t = 256 s (long-term effect) for CMNS ( Figure 9A) with those for the no-stim protocol at t = 128 s ( Figure 4A) show similar patterns. Note that the patterns do not change for the no-stim protocol at different times e.g., at t = 0 or 256 s (results not shown). At the end of the stim-off period the medians and IQRs of the sorted connectivity matrices of CMNS ( Figure 9D) and of the no-stim protocol ( Figure 9B) look also similar and are comparable to those of the no-stim protocol at the end of the stim-on period ( Figure 4B): strong unidirectional synaptic connections, with only a very small difference between the sorted connectivity matrices of different samples. Compared to the acute effect induced by CMNS ( Figure 4D) the network has lost its strong bidirectional excitatory synapses during the stim-off period ( Figure 9D) for all samples.
In case the stimulation onsets of the simultaneous stimulation patterns are not random as for CMNS but periodic (PPMS), FIGURE 8 | Long-lasting anti-kindling effects for different stimulation protocols. (A) Boxplots of C av (t = 256 s) for different stimulation protocols show a statistically significant decrease compared to the control condition (no-stim) for the stimulation protocols with sequential activation of the four stimulation sites. In contrast, CMNS causes no decrease of C av. (B) Boxplots of R av (t = 256 s) show a statistically significant desynchronization induced by the different stimulation protocols compared to the control signal (no-stim) except for the inert CMNS. All stimulation protocols, except CMNS, induce a long-lasting decrease in the average excitatory synaptic weight c EE (C) and a long-lasting increase in the average inhibitory synaptic weight c II (D). Each boxplot represents the results of 11 samples. The horizontal line within the box represents the median, the length of the box the IQR (middle 50%), whereas the whiskers below and above the box indicate the first and last 25%. Outliers are defined as 1.5 times the length of the box below or above the box and are represented by open circles. p-values of the left-sided Wilcoxon rank sum test are given in Supplementary Table 3. K = 0.25 for all panels.
bidirectional strong inhibitory synapses exist at the end of the stim-off period and most excitatory synapses are unidirectional and strong, while some of them are bidirectional and weak ( Figure 9C). Indices of the neurons which have bidirectional strong inhibitory synapses or bidirectional weak excitatory synapses can differ for other samples (Figure 9C).
For the CR protocols, the medians of the sorted connectivity matrices have lost their connectivity patterns of the shortrange excitatory synapses during the stim-off period (compare Figures 9E-H with Figures 4E-H). Mainly the CR protocols without repetition (Figures 9E,F) show a loss of some bidirectional strong inhibitory synapses, which were formed during the stim-on period (Figures 4E,F). For the CR protocols the IQRs of the sorted connectivity matrices are different at the end of the stim-off period compared with those at the end of the stim-on period (compare the bottom panels of Figures 9E-H with those of Figures 4E-H): at the end of the stim-off period large IQR are not only found for the strongest synapse between neighboring neurons (diagonals) or nearby stimulation sites, but without a clear spatial pattern in the lower right triangle large IQRs can be found for the excitatory synapses and in the upper left triangle for the inhibitory synapses influenced mainly by UMNS or RVS CR (bottom panels of Figures 9E-H) .
Although on the macroscopic neural network level there is no significant difference between the desynchronization induced by the CMNS and by the no-stim protocol (Figure 8B), on the subpopulation level, there is a statistically significant acute aftereffect ( Figure 5B): An increase in the amount of synchronization is induced by CMNS. All other stimulation protocols have a stronger desynchronizing effect than the no-stim protocol, whereby the acute after-effects of the RVS, fixed and SVS CR are increased compared to their acute stimulation effects (Figures 5A,B). These increases are statistically significant (see Supplementary Table 2 for the p-values of the right-sided Wilcoxon rank sum test with significance level α = 0.05).
Accordingly, we conclude that sequential activation of the stimulation sites also induces more pronounced long-term antikindling effects except for UMNS, which tends to induce better long-term anti-kindling effects than PPMS, but the improved medians are not statistically significant (network level see Figures 8A,B; subpopulation level see Figures 5B, 8C,D). At the end of the stim-off period the results induced by consecutively repeating each sequence many times (fixed CR and SVS CR) are still better than without repetition (UMNS and RVS CR) on the network level (Figures 8A,B), as well as on the subpopulation level (Figures 8C,D for the connectivity; Figure 5B for the desynchronizing effect).
Even at the end of the stim-off period the CMNS protocol results in the strongest c EE and the weakest c II , and therefore gives rise to the largest C av value of all stimulation protocols ( Figures 8A,C,D). The long-term R av shows no difference with the control signals (no-stim; Figure 8B), but it does show a synchronizing effect on the subpopulation level, which implies that within the subpopulations the amount of synchronization has increased compared to the beginning of the stim-on period ( Figure 5B). Despite these small differences the median and IQR of the sorted connectivity matrices induced by the no-stim and CMNS protocols appear to be similar. From all the investigated stimulation protocols at K = 0.25 the CMNS protocol turns out to be the one which induces the most similar results as the control signal (no-stim), despite small but statistically significant differences.

Robustness Against Stimulation Intensity
In this section, we investigate acute and long-term effects elicited by stimulation intensities weaker than K = 0.25. From the previous sections, we can conclude that a small difference in stimulation protocols might result in completely different acute as well as long-term effects. For instance, the simultaneous and noisy stimulation (CMNS) does not decrease C av as the simultaneous and period stimulation (PPMS) does. By the same token, at weaker intensities (K = 0.10, 0.15, and 0.20) PPMS induces statistically significantly smaller values of C av and R av values than CMNS (Figures 10A,B). See Supplementary  Tables 1, 3, 5-7 for the corresponding p-values of the left-sided Wilcoxon rank sum test with significance level α = 0.05. The basic difference in the stimulation protocol feature between UMNS and fixed CR is the same as between PPMS and CMNS, namely noise or periodic stimulation. By comparing the pvalues for the UMNS and the fixed CR it follows that for K = 0.10 UMNS induces smaller C av and R av values than fixed CR (Supplementary Table 5), but for K = 0.25 it is just the other way around (Supplementary Tables 1,3). Since for the investigated K-range PPMS always induces smaller C av and R av values than the CMNS, we can only conclude that the effect of noise stimulation also depends on how the stimulation sites are activated: simultaneously or sequentially. In case of sequential stimulation it depends strongly on K.
Fixed CR and SVS CR differ in the number of different sequences applied to the network. For fixed CR one sequence is applied 4,800 time during the on-period, whereas for SVS CR each sequence is applied 100 times before the next sequence is applied. For weak stimulation intensities up to K = 0.15 there is no difference in the anti-kindling effects. In contrast, at higher intensities changing the sequence from time to time, as for SVS CR, decreases the C av and R av values even more (Figures 10E,F) than the fixed CR does.
Cross-trial analysis shows that already for K = 0.10 a weak phase reset is observed for CMNS, UMNS, and RVS CR and a weak entrainment followed by a stimulus-locked disruption of the weak phase entrainment for PPMS, fixed CR, and SVS CR (Figure 11). Increasing the intensity also increases the amount of entrainment as well as the strength of the phase reset (e.g., Figures 11A,D,G). During the first ON-cycle, PPMS, fixed CR as well as SVS CR stimulation destroy the entrainment of the subpopulation phase-dynamics with the stimulation signal in the entire investigated K-range. During the second ON-cycle, this entrainment is partly recovered by the stimulation except for stimulations at K = 0.15. In the latter case the entrainment is destroyed even further. During the third ON-cycle entrainment is restored in the entire K-range (K ∈ {0.10; 0.15; 0.20; 0.25}), where the amount of entrainment increases with increasing K. The phase rest seems to increase slightly from the first to the second and, finally, to the third ON-cycle, but clearly increases with increasing K (K ∈ {0.10; 0.15; 0.20; 0.25}).

DISCUSSION
We studied acute and long-lasting effects of six different stimulation protocols and compared the observed effects on a plastic neural network with a no-stimulation control condition. While sharing the same average rate of stimuli per channel, the tested stimulation protocols differ with respect to their amount of periodicity as opposed to randomness, both between ON-cycles and between stimulation sites. One of the tested stimulation protocols, CMNS, turned out not to induce desynchronization. In fact, comparing CMNS with the no-stimulation control condition showed that CMNS is nearly inert. More precisely, during stimulation (Figure 2), CMNS caused a significant increase of the strength of the mean synaptic connectivity by 49% compared to the no-stimulation control condition. Intriguingly, during CMNS the overall synchrony nevertheless was 8% smaller than for the no-stimulation control. Remarkably, there was hardly any long-term post-stimulation effect of CMNS compared to the no-stimulation control condition (Figure 8): While the mean synaptic weight increased significantly, but only slightly by 7%, the spontaneous (i.e., stimulation-free) synchronization did not significantly differ between the CMNS and the nostimulation condition. Whether a slight, but significant increase of the mean synaptic weight might be relevant in a nonspontaneous context, where the network is, e.g., subjected to other types of stimuli, remains to be tested. Ultimately, clinical studies will provide the necessary tests.
In this study, we have shown that stimulation protocols, that differ by just one, putatively minor feature, may cause massively different anti-kindling effects, robustly over a range of stimulation intensities as well as for different samples at the end of the stim-on period and also at the end of the stimoff period (see e.g., Figure 10). Stimulating simultaneously or sequentially has a big influence on the synchronization and connectivity of the network at all levels (e.g., compare CMNS vs. UMNS in Figures 2C,D, 3C,D, 4D,E, 5, 6B,C, 8, 9D,E, 10). Another influential feature of some stimulation protocols is repetition (compare e.g., PPMS vs. CMNS or RVS vs. fixed CR), which goes along with a different dynamical stimulation mechanism: entrainment caused by repetition vs. phase reset otherwise (Figures 7, 10). In this study, we restricted ourselves to cycle durations of 16 ms, which is slightly above the intrinsic firing period of the individual neurons (14 ms). For a more pronounced mismatch of ON-cycle duration and firing periods, stimulus-locked entrainment of subpopulations may probably become more difficult if not completely impossible. Another drawback of the SVS and fixed CR is that they have to be applied at slightly stronger intensities than UMNS and RVS to be effective. However, SVS and fixed CR may, in principle, induce more pronounced long-term anti-kindling effects than UMNS and RVS CR (see e.g., Figures 8A,B).
Compared to all stimulation protocols tested in this study, the anti-kindling effects induced by the CMNS protocol are most similar to the control (i.e., no-stim) protocol (Figures 2C,D,  3C,D, 4D,E, 5, 6B,C, 8, 9D,E, 10), although the similarity is not always statistically significant. The CMNS protocol is also the only protocol which induced opposite effects on the connectivity, both on the network as well as subpopulation level, compared to all other active protocols (Figures 3, 8). Although the implementation of the stimulation in the model used here is for sensory stimulation, we expect similar results for electrical stimulation. For comparison between electrical and sensory stimulation effects and their qualitative similarities in our model, see Popovych and Tass (2012).
In addition, we also performed simulations which extended the finished trials with the CMNS protocol (K = 0.20) by adding a second, additional 128 s lasting stim-on period at t = 256 s with the RVS CR protocol as well as with the SVS CR protocol, each at K = 0.20. Anti-kindling effects induced at t = 384 s (i.e., at the end of the second stim-on period) as well as at t = 512 s (the end of the second stim-off period) were statistically significantly similar with the anti-kindling effects at t = 128 s FIGURE 10 | Boxplots of C av and R av show the effect of one difference between stimulus signal features. (A,B) PPMS and CMNS differ in the time between two consecutive stimulation onsets: PPMS has a constant period of 16 ms in consecutive ON-cycles, whereas CMNS has changing periods between stimulation onsets. The induced anti-kindling effects induced by PPMS are better than those induced by CMNS at the end of the stim-on period (A) as well as at the end of the stim-off period (B). (C,D) For UMNS stimulation onsets are random and uncorrelated between different sites, whereas for RVS CR they are restricted to only four equidistant moments within the ON-cycle. RVS CR induces a stronger reduction of C av and R av than UMNS for most K-values at t = 128 s (C) as well as at t = 2 56 s (D). (E,F) Fixed CR applies the same sequence during the stim-on period, whereas SVS CR randomly draws a new sequence after 100 consecutive repetitions of a sequence. For K = 0.10 and K = 0.15 there are no statistically significant differences in the distributions of C av and R av at t = 128 s (E) and t = 256 s (F). However for K = 0.15 and K = 0.20 SVS CR reduces the medians and IQRs of C av and R av more than the fixed CR does. Each boxplots represents eleven samples. p-values of the left-sided Wilcoxon rank sum test are given in Supplementary Tables 1, 3, 5-7.
(the end of the first and only stim-on period) and t = 256 s (the end of the first and only stim-off period) for the RVS and SVS CR only protocols applied during only one stim-on period as shown in Figures 10C-F (data not shown). Accordingly, from a computational standpoint there is no reason to assume that application of CMNS sham stimulation might render subsequent delivery of RVS or SVS CR stimulation ineffective.
Apart from the stimulation-related aspects, our findings are relevant with respect to the assessment of changes of synchronization and synaptic connectivity patterns and their mutual interrelation. Although the macroscopic measures C av and R av were often strongly correlated, an increase in C av did not necessarily imply an increase in R av . For example, in Supplementary Table 4, we show that for the CMNS protocol applied at K = 0.25 during the stim-off period a statistically significant decrease in C av is combined with a statistically significant increase in R av . Accordingly, it is not sufficient to exclusively monitor the connectivity in our model. Rather, we also have to determine the amount of synchronization. By a similar token, comparable amounts of macroscopic synchronization/desynchronization, as assessed by the order parameter, may differ on the mesoscopic, i.e., subpopulation level (see e.g., Figure 7). Accordingly, there is not just one type of (e.g., uniform) desynchronization. In summary, the analysis of both synchrony and synaptic connectivity on macroscopic as well as mesoscopic network levels may further the development of both active and inactive (sham) stimulation protocols.
In the field of drug development, a placebo treatment is realized by delivering a substance with no active therapeutic effect (Friedman et al., 2015). Accordingly, placebo effects can, for instance, be assessed with a parallel group design by comparing a placebo group with a natural history (i.e., no-treatment) group (Wager and Atlas, 2015). According to our computational results, CMNS is a promising candidate for a sham (i.e., inactive) stimulation protocol in the field of desynchronizing multichannel stimulation. In the present study, we did not aim at developing a biophysical, microscopic model. Rather, we used a minimal model, equipped with robust spontaneous multistable dynamics, comprising synchronized and desynchronized states, that served as testbed for different stimulation protocols to generate first and experimentally testable hypotheses. By a similar token, several previous computational studies in the field of CR stimulation were carried out in minimal models and led to a number of clinically significant predictions, e.g., concerning cumulative effects and stimulation intensity (Hauptmann and Tass, 2009;Lysyansky et al., 2011), that were verified in preclinical and clinical studies (Tass et al., 2012b;Adamchic et al., 2014;Wang et al., 2016). In the same manner, the present computational study yields the testable hypothesis that CMNS might serve as sham stimulation protocol. Accordingly, the effects of CMNS should be tested in pre-clinical and, in particular, in clinical studies in order to disentangle CMNS effects from placebo effects in the best possible way. For instance, in a phase 1 (first in man) study feasibility and tolerability could be tested. In addition, in a phase 2 study an assessment of effects should be performed, possibly in comparison to a no-stimulation control and/or an active control group. Obviously, this adds to the complexity of the clinical development of neuromodulation treatments.
Furthermore, CMNS may be useful in clinical studies focusing on revealing predictive EEG markers for optimizing stimulation parameters of desynchronizing neuromodulation interventions, see Adamchic et al (2017). CMNS does neither cause substantial acute nor long-lasting effects. Nevertheless, CMNS might have unspecific EEG effects, e.g., on brain rhythms and/or brain areas less important to the disease related network dynamics. In this way, CMNS might help to separate EEG responses related to core stimulation effects from concomitant EEG responses. In addition, CMNS may also be helpful in pre-clinical studies to elucidate mechanisms by which stimulation protocols cause a desynchronization. Furthermore, CMNS might also help to reveal mechanisms that might be related to, but go beyond desynchronization, such as therapeutic rewiring (Tass and Majtanik, 2006) or neuroprotective effects (Musacchio et al., 2017).
Wherever appropriate and possible, sham procedures and double-blind protocols should be developed to scrutinize specific effects of neuromodulation interventions by adequate clinical trials and rule out placebo effects. However, placebo effects should not just be considered as a nuisance, requiring cumbersome clinical study protocols. Rather, the mechanisms underlying different placebo effects (Benedetti et al., 2011) could actually be specifically exploited to better neuromodulation techniques. For instance, similar to conditioned immunomodulation (Metal'nikov and Chorine, 1926;Ader and Cohen, 1975;Ader, 2003), one could condition specific desynchronization stimulation delivered invasively and causing long-lasting effects with unspecific non-invasive, sensory stimuli (Tass, 2011).
For the clinical implementation of multichannel sham stimulation protocols, one should take into account possible side effects related to the technical generation of the stimuli. For instance, in the field of transcranial magnetic stimulation (TMS) the current gold standard for sham TMS appears to be the use of a shielded coil, generating characteristic stimulusrelated auditory stimuli, but no magnetic brain stimuli, together with surface electrodes for skin stimulation, mimicking magnetic skin stimulation (Duecker and Sack, 2015). Intriguingly, sham TMS may have specific side effects Sack, 2013, 2015). Obviously, TMS is not just a purely magnetic stimulation modality, but may constitute a compound stimulation approach which may cause specific effects caused by stimuli of different modality Sack, 2013, 2015). For the clinical development of multichannel sham stimulation protocols such aspects might be relevant, for instance, when vibrotactile mulitchannel stimulation causes auditory and possibly other sensory side effects.
This computational study is a first step for the development of a sham stimulation protocol for multichannel desynchronizing stimulation techniques. For comparison, for the development of CR stimulation, in computational studies predominantly minimal models were used (Tass, 2003a,b;Tass and Majtanik, 2006;Tass, 2007, 2009;Lysyansky et al., 2011;Popovych and Tass, 2012;Tass, 2015, 2016), as opposed to biophysically realistic models (Ebert et al., 2014). These computational studies revealed non-trivial predictions, e.g., concerning the emergence of long-lasting, sustained (Tass and Majtanik, 2006) as well as cumulative (Hauptmann and Tass, 2009) effects and concerning the amplitude of the stimulation amplitude (Lysyansky et al., 2011). These predictions were verified in pre-clinical (Tass et al., 2012b;Wang et al., 2016) and clinical studies (Tass et al., 2012a;Adamchic et al., 2014;Syrkin-Nikolau et al., 2018). The computational predictions were used as hypotheses for the design of the corresponding study protocols. Analogously, the development of sham stimulation protocols requires a combined effort, comprising computational, pre-clinical, and clinical studies. The effect of a sham stimulation protocol may depend on the type of the stimulated neurons, the target area of the stimulation and the mechanism of the stimulation. Accordingly, future studies should also use other neuronal network models, e.g., network of FitzHugh-Rinzel bursting neurons (Rinzel, 1987;Izhikevich, 2001). By the same token, sham stimulation should ideally be inactive also in the presence of additional features and mechanisms, such as synaptic noise (Destexhe et al., 2003), propagation delays (Madadi Asl et al., 2017) as well for different stimulation mechanisms, e.g., excitatory vs. inhibitory stimulation (Popovych and Tass, 2012). Hence, future computational studies should take into account these refinements, too.

AUTHOR CONTRIBUTIONS
PT came up with the initial ideas for this work; MZ designed and performed the simulations and analyzed the data; MZ and PT discussed findings and interpretations, and wrote the manuscript.

ACKNOWLEDGMENTS
PT gratefully acknowledges the generous research funding support (donation) from the John A. Blume Foundation.