ORIGINAL RESEARCH article
Long-Lasting Desynchronization Effects of Coordinated Reset Stimulation Improved by Random Jitters
- Department of Neurosurgery, Stanford University, Stanford, CA, United States
Abnormally strong synchronized activity is related to several neurological disorders, including essential tremor, epilepsy, and Parkinson's disease. Chronic high-frequency deep brain stimulation (HF DBS) is an established treatment for advanced Parkinson's disease. To reduce the delivered integral electrical current, novel theory-based stimulation techniques such as coordinated reset (CR) stimulation directly counteract the abnormal synchronous firing by delivering phase-shifted stimuli through multiple stimulation sites. In computational studies in neuronal networks with spike-timing-dependent plasticity (STDP), it was shown that CR stimulation down-regulates synaptic weights and drives the network into an attractor of a stable desynchronized state. This led to desynchronization effects that outlasted the stimulation. Corresponding long-lasting therapeutic effects were observed in preclinical and clinical studies. Computational studies suggest that long-lasting effects of CR stimulation depend on the adjustment of the stimulation frequency to the dominant synchronous rhythm. This may limit clinical applicability as different pathological rhythms may coexist. To increase the robustness of the long-lasting effects, we study randomized versions of CR stimulation in networks of leaky integrate-and-fire neurons with STDP. Randomization is obtained by adding random jitters to the stimulation times and by shuffling the sequence of stimulation site activations. We study the corresponding long-lasting effects using analytical calculations and computer simulations. We show that random jitters increase the robustness of long-lasting effects with respect to changes of the number of stimulation sites and the stimulation frequency. In contrast, shuffling does not increase parameter robustness of long-lasting effects. Studying the relation between acute, acute after-, and long-lasting effects of stimulation, we find that both acute after- and long-lasting effects are strongly determined by the stimulation-induced synaptic reshaping, whereas acute effects solely depend on the statistics of administered stimuli. We find that the stimulation duration is another important parameter, as effective stimulation only entails long-lasting effects after a sufficient stimulation duration. Our results show that long-lasting therapeutic effects of CR stimulation with random jitters are more robust than those of regular CR stimulation. This might reduce the parameter adjustment time in future clinical trials and make CR with random jitters more suitable for treating brain disorders with abnormal synchronization in multiple frequency bands.
The human organism can be viewed as an integrated network where complex physiological systems continuously interact; whereby the regulatory mechanisms of one system may affect others or the organism as a whole (Bashan et al., 2012; Ivanov et al., 2016). Associations between network topology and network functionality may provide insights into how distinct pathological/physiological states emerge from nonlinear interactions between multi-compartment complex systems (Bartsch et al., 2015). On the other hand, therapeutic effects of stimulation applied to one system may spread throughout the entire network (Pfeifer et al., 2021).
In neurological disorders, such as epilepsy (Mormann et al., 2000) or Parkinson's disease (PD) (Alberts et al., 1969; Nini et al., 1995), patients suffer from pronounced motor symptoms that are caused by impaired brain activity. In PD, several strongly interconnected brain areas are involved, including the basal ganglia, the thalamus, and the sensorimotor cortex. In several of these areas, symptom-related abnormal neuronal synchrony has been observed (Nini et al., 1995; Hammond et al., 2007). An established treatment for advanced Parkinson's disease is high-frequency deep brain stimulation (HF DBS). HF DBS has been delivered to several target areas, such as the subthalamic nucleus (STN) (Krack et al., 2003) or the ventral intermediate nucleus of the thalamus (Benabid et al., 1991). HF DBS of the STN is the standard of care for treating medically refractory Parkinson's disease, however, there is no consensus on its mechanism of action (Ashkan et al., 2017; Jakobs et al., 2019; Lozano et al., 2019; Krauss et al., 2020).
Recently, DBS has been suggested as a treatment for other disorders, including obsessive-compulsive disorder (OCD) (Vicheva et al., 2020). Symptoms of OCD include uncontrollable recurring thoughts (obsessions) and repetitive behaviors (compulsions). OCD is also associated with substantial comorbidities, including substance use disorders, anxiety, and impulse-control (Ruscio et al., 2010; Vicheva et al., 2020). The exact mechanism underlying OCD and the therapeutic effect of DBS for OCD remain enigmatic (Bourne et al., 2012; Ahmari and Dougherty, 2015; Vicheva et al., 2020).
Besides several benefits, continuous HF DBS has limitations. In PD patients, HF DBS may successfully suppress symptoms while stimulation is on, however, symptoms return shortly after cessation of stimulation (Temperli et al., 2003). Permanent stimulation of the target area and surrounding tracts and nuclei as well as the corresponding medication dose adjustment may lead to DBS-induced movement disorders (Baizabal-Carvallo and Jankovic, 2016), e.g., characterized by dyskinesias, gait disorder, dysarthria, ataxia etc. (Rodriguez-Oroz et al., 2005; Temel et al., 2006; Moreau et al., 2008; Tripoliti et al., 2008; Schrader et al., 2011; Baizabal-Carvallo and Jankovic, 2016; Xie et al., 2017).
To substantially reduce the integral amount of delivered stimulation current and the risk of unwanted side effects, several studies focused on developing stimulation approaches that specifically counteract pathological synchrony. Some approaches use delayed feedback to desynchronize networks of oscillators (Rosenblum and Pikovsky, 2004a,b; Popovych et al., 2005; Hauptmann et al., 2005a,b,c; Popovych et al., 2006a,b; Pyragas et al., 2007; Popovych and Tass, 2010); clinically this could be implemented using linear or nonlinear delayed feedback envelope pulse trains (Popovych et al., 2017a,b). Other studies analyzed the nonlinear response of an ensemble of coupled oscillators to external stimuli. For instance, an ensemble of synchronized oscillators can be desynchronized by delivering a single stimulus pulse at a vulnerable phase of the collective rhythm (Mines, 1914; Winfree, 1977, 1980; Warman and Durand, 1989; Tass, 1999). During double-pulse stimulation, such a desynchronizing pulse is delivered shortly after a strong phase-resetting pulse to increase robustness (Tass, 2001, 2002; Zhai et al., 2005). Other studies suggest the delivery of periodic stimulation. It was shown that periodic stimulation with certain frequencies can desynchronize a synchronous ensemble of oscillators. This effect is known as chaotic desynchronization (Wilson et al., 2011). Phasic burst stimulation is another approach that aims at increasing the phase differences between individual oscillators by delivering stimulation bursts at certain phases of the collective rhythm. Corresponding phases were predicted using estimated phase response curves calculated for each patient. A corresponding closed-loop method was suggested in Holt et al. (2016). Other techniques deliver spatio-temporal stimulus patterns through multiple stimulation sites, one of which is coordinated reset stimulation (CR) (Tass, 2003b). During CR stimulation, desynchronization is achieved by delivering phase-shifted stimuli to individual neuronal subpopulations.
Most of these desynchronization techniques were developed for networks with fixed connectivity. In the brain, however, synaptic plasticity leads to dynamic reorganizations of neuronal networks (Doidge, 2007; Liu et al., 2015; Van Ooyen and Butz-Ostendorf, 2017; Eagleman, 2020). One plasticity mechanism is spike-timing-dependent plasticity (STDP), where the change of the strengths of synaptic connections depends on the relative timing of post- and presynaptic spikes (Markram et al., 1997; Abbott and Nelson, 2000; Caporale and Dan, 2008). In many brain regions, STDP leads to a strengthening of synapses, if the postsynaptic neuron fires shortly after the presynaptic neuron, and to a weakening in the opposite case (Markram et al., 1997; Bi and Poo, 1998). STDP may lead to the formation of strongly connected neuronal assemblies (Litwin-Kumar and Doiron, 2014) and may stabilize certain patterns of neuronal activity, e.g., synchronized activity (Karbowski and Ermentrout, 2002). Plasticity can also lead to multistability in neuronal networks and networks of oscillators. In particular, networks with coexisting stable states, such as cluster states, desynchronized states, and synchronized states have been studied (Seliger et al., 2002; Zanette and Mikhailov, 2004; Tass and Majtanik, 2006; Masuda and Kori, 2007; Maistrenko et al., 2007; Aoki and Aoyagi, 2009; Berner et al., 2020; Yanchuk et al., 2020).
Extensive theoretical and computational studies on CR stimulation of multistable plastic networks showed that CR stimulation may reshape the synaptic connectivity, and cause long-lasting desynchronization by driving the network into the attractor of a stable desynchronized state (Tass and Majtanik, 2006; Hauptmann and Tass, 2009; Popovych and Tass, 2012; Lourens et al., 2015; Manos et al., 2018; Kromer and Tass, 2020; Kromer et al., 2020). Corresponding long-lasting therapeutic effects and/or sustained reduction of neuronal synchrony were confirmed in preclinical in vitro studies (Tass et al., 2009), in preclinical in vivo studies (Tass et al., 2012b; Wang et al., 2016) as well as in clinical studies, delivering CR stimulation through implanted DBS electrodes (Adamchic et al., 2014) or noninvasively by means of vibrotactile CR fingertip stimulation (Tass, 2017; Syrkin-Nikolau et al., 2018; Pfeifer et al., 2021).
Detailed computational studies on CR stimulation of plastic neuronal networks suggested that long-lasting desynchronization effects may be sensitive to the ratio of the stimulation frequency and the dominant frequency of the pathological synchronous rhythm, frhythm (Manos et al., 2018). Presumably, this is because stimuli are delivered with fixed inter-stimulus intervals which may lead to unfavorable resonances with other time scales (Kromer and Tass, 2020). A lack of frequency robustness might restrict clinical application as individual symptoms during PD are related to pathological synchrony in different frequency bands. Specifically, synchronized basal ganglia activity in the theta band (3−10 Hz) has been associated with symptoms such as dyskinesia and tremor (Brown, 2003; Steigerwald et al., 2008; Tass et al., 2010; Contarino et al., 2012), whereas synchronized beta band activity (13−30 Hz) has been associated with symptoms such as rigidity and bradykinesia (Kühn et al., 2006; Weinberger et al., 2006). Furthermore, tremor may be associated with different central oscillators (Raethjen et al., 2000).
In order to increase the frequency robustness of long-lasting effects, spatial and temporal randomization of stimulus deliveries was suggested (Kromer and Tass, 2020). In particular, random reset (RR) stimulation was introduced (Kromer and Tass, 2020) during which stimuli are delivered at random times to randomly selected neuronal subpopulations. During RR stimulation, temporal randomization is realized by choosing stimulation times according to a Poisson spike train. Spatial randomization is achieved by delivering stimuli to randomly selected neuronal subpopulations. However, in Kromer and Tass (2020) no spatial relations between neurons were considered, this was somewhat artificial and could not be directly applied to a DBS-type setup as it implicitly assumed microscopic control, i.e., it assumed that individual neurons could be stimulated independently. Compared to CR stimulation, RR stimulation presents an extreme case with minimal temporal and spatial correlations between stimulus deliveries (Kromer and Tass, 2020; Khaledi-Nasab et al., 2021b). It was found that long-lasting desynchronization effects of RR stimulation were more robust with respect to parameter changes than those of CR stimulation, whereas acute desynchronization effects of RR stimulation were significantly weaker (Kromer and Tass, 2020; Khaledi-Nasab et al., 2021b). Also, temporal randomization during RR stimulation helped to avoid stimulus deliveries during the neurons' refractory periods. The authors reported that the resulting synaptic reshaping was more likely to drive the network into a stable desynchronized state.
It was argued that improved parameter robustness resulted from two effects: First, for CR stimulation with unfavorable stimulation frequencies, a significant portion of stimuli might be delivered during the neurons' refractory periods. This renders such stimulation protocols ineffective. In contrast, temporal variability of stimulus deliveries during RR stimulation improved the robustness of long-lasting effects of RR stimulation relative to CR stimulation. Second, RR stimulation led to a broader distribution of time lags between post- and presynaptic spike times that determined weight updates due to STDP. Thus, a bigger part of the STDP function was considered for weight updates. It was found that while CR with fine-tuned stimulation frequency led to faster weight reduction, weight reduction during RR was more robust with respect to changes of the stimulation frequency (Kromer and Tass, 2020).
In the present paper, we hypothesize that appropriately adding temporal and/or spatial randomness to CR stimulus patterns might improve the parameter robustness of long-lasting CR effects, thereby preserving its pronounced acute desynchronization effects. To this end, we analyze the effect of temporal and spatial correlations in spatio-temporal stimulus patterns on acute and long-lasting desynchronization effects of the stimulation. First, we consider a classic CR stimulation pattern (Tass, 2003a; Kromer et al., 2020). Then, we reduce temporal and spatial correlations between stimulus delivery times. Temporal correlations are reduced by adding a random jitter to the stimulus delivery times. This results in a noisy CR pattern (NCR). Spatial correlations are reduced by shuffling the sequence of stimulation site activations. This results in a shuffled CR pattern (SCR). Finally, both jitter and shuffling are applied, leading to a shuffled noisy CR pattern (SNCR). Employing theoretical analysis and computer simulations of plastic neuronal networks, we analyze and compare the effect of random jitter and spatial shuffling on the long-lasting outcome of stimulation. Our results suggest that random jitter is more suitable for increasing the frequency robustness of long-lasting desynchronization effects than spatial shuffling. Of note, random jitter does not degrade acute desynchronization effects.
This paper is organized as follows: In section 2, we introduce the model and the different stimulation patterns used throughout the manuscript. In section 2.6, we derive theoretical predictions of the stimulation-induced synaptic weight dynamics for different stimulation patterns in the limit of strong and fast stimulation. Then, in section 3, we compare theoretical predictions to results from numerical simulations of networks of leaky integrate-and-fire (LIF) neurons with STDP. Furthermore, we present results for weak stimulation. We find that random jitter is more suitable for increasing parameter robustness of long-lasting effects than shuffling. Finally, in section 4 we provide a detailed discussion of our results.
2. Model and Methods
2.1. Neuronal Network Model
Throughout the paper, we consider networks of 103 excitatory LIF neurons with STDP (Kromer et al., 2020; Khaledi-Nasab et al., 2021b). Parameters are chosen according to Kromer et al. (2020) and Khaledi-Nasab et al. (2021b) such that a stable synchronized and a stable desynchronized states coexist. See Appendix for more details. Neurons are equidistantly spaced in the interval xi ∈ [−2.5, 2.5] mm. This is motivated by the length of the short axes of an elipsoidal volume approximation of the human STN used in a detailed computational study (Ebert et al., 2014). Each neuron has Nsyn = 0.07N outgoing synapses, where the probability for a synaptic connection between neurons i and j is proportional to ∝ exp((|xj − xi|)/0.5 mm) (Ebert et al., 2014). Throughout the paper, simulation results are averaged over three different network realizations, i.e., different realizations of random network connectivity.
Initially, synaptic weights wij(t = 0) are randomly set to either one or zero, such that a mean synaptic weight of 〈w(t = 0)〉 = 0.5 is realized. Here, i and j refer to the pre- and postsynaptic neurons, respectively. Each network realization was simulated until it approached a stable synchronized state, see Kromer et al. (2020) and Khaledi-Nasab et al. (2021b).
Stimulation is applied to Ns stimulation sites. To this end, we divide the interval of possible x coordinates into Ns equal segments. Neurons with coordinates in the kth segment are assumed to receive stimuli delivered to the kth stimulation site. Throughout the paper, we will refer to these neurons as the kth subpopulation. Neurons in the same subpopulation receive stimuli simultaneously and with equal strength. We thereby neglect distance-dependent modulations of the received stimulation current and finite travel times of the electrical signal through the tissue.
2.2. Spike-Timing Dependent Plasticity
The dynamics of synaptic weights wij(t) is determined by STDP. We consider a nearest-neighbor STDP scheme in which weight updates are performed at postsynaptic spike times and presynaptic spike arrival times (Morrison et al., 2008). Corresponding weight updates wij → wij + W(tj − (ti + td)) are given by the STDP function (Kromer and Tass, 2020; Song et al., 2000)
Here, Δt = tj − (ti + td) is the time lag between the current postsynaptic spike time tj and the latest presynaptic spike arrival time ti + td (if the update is done at a postsynaptic spike time), or the time lag between the current presynaptic spike arrival time ti + td and the latest postsynaptic spike time tj (if the update is performed at a presynaptic spike arrival time). η = 0.02 scales the weight update per spike, τR = 4 yields asymmetry in STDP decay times, τ+ = 10 ms and τ− = τ+ τR, β = 1.4 scales the ratio of overall long-term depression (LTD) to long-term potentiation (LTP).
These STDP parameters lead to the coexistence of a strongly connected state with synchronized neuronal activity and a weakly connected state with asynchronous neuronal activity (Kromer and Tass, 2020; Kromer et al., 2020; Khaledi-Nasab et al., 2021a,b).
2.3. Stimulation Patterns
In the present paper, we consider four stimulation patterns: a (regular) CR pattern and three randomized CR patterns. Randomization is obtained by adding random jitters to the stimulation times (reduced temporal correlations between stimuli), shuffling of the sequence of stimulation sites (reduced spatial correlations between stimuli), and a combination of both random jitters and shuffling (reduced temporal and spatial correlations between stimuli). In the following, we introduce these four stimulation patterns in detail.
I. Regular CR stimulation (CR): CR stimulation is delivered in cycles of Ns stimuli (Tass, 2003b), where Ns is the number of stimulation sites. During each CR cycle, each site receives exactly one stimulus, and stimuli are administered at times , k = 0, 1, .., Ns − 1. Here t0 is the beginning of the CR cycle and TCR = 1/fCR is the cycle length. fCR is the stimulation frequency and corresponds to the mean frequency at which individual sites receive stimuli.
In several preclinical and clinical studies (Tass et al., 2012b; Adamchic et al., 2014; Pfeifer et al., 2021), the sequence of stimulation site activations during each cycle was chosen at random for each cycle. The resulting version of CR stimulation is referred to as CR with rapidly varying sequence in the literature (Popovych and Tass, 2012; Zeitler and Tass, 2015). For the sake of brevity, we will refer to this pattern as CR.
II. Noisy CR stimulation with random jitters (NCR): Same as CR stimulation except that there are random jitters added to the stimulation times sk. Jitters are uniformly distributed in the interval with 0 ≤ σCR ≤ 1. Thus, stimuli during a CR cycle are delivered at random times , k = 0, 1, .., Ns − 1. The case σCR = 1, refers to the limit of maximum variability of stimulus onset times. We exclude larger values of σCR as these would result in overlapping intervals of possible stimulus onset times.
A corresponding stimulation pattern with moderate jitters was used in a clinical study on vibrotactile CR stimulation of Parkinson's patients (Pfeifer et al., 2021). There, it was referred to as noisy CR (NCR). We will refer to CR with random jitter as NCR stimulation.
III. Shuffled CR stimulation (SCR): Same as CR stimulation except that stimuli are delivered to randomly selected sites. Thus, the restriction that each site receives exactly one stimulus per cycle was lifted. Sites are selected for stimulus delivery with uniform probability.
This stimulation pattern was used in a previous study on stimulation-induced desynchronization (Tass and Hauptmann, 2009). There, it was referred to as multi-site random-site stimulation (Tass and Hauptmann, 2009). We will refer to shuffled CR stimulation as SCR stimulation.
IV. Shuffled Noisy CR with random jitters (SNCR): Same as SCR but with random jitters added to the stimulus delivery times. We will refer to shuffled CR with random jitters as SNCR stimulation. In the limit of small jitter, σCR → 0, SNCR stimulation becomes equivalent to SCR stimulation.
For NCR and SNCR stimulation, the parameter σCR scales the width of the distribution for random jitters and therefore the reduction of temporal correlations between stimulus delivery times. In the limit of σCR = 0, the NCR pattern is equivalent to the CR pattern, and the SNCR pattern is equivalent to the SCR pattern.
Individual stimuli are charge-balanced and consist of excitatory and inhibitory pulses with durations of νe = 0.5 ms and νi = 3 ms, respectively. This asymmetry is motivated by preclinical and clinical studies on CR stimulation employing asymmetric pulse shapes (Wang et al., 2016; Adamchic et al., 2014). The excitatory pulse has the amplitude and the inhibitory pulse the amplitude . The two pulses are separated by a gap of 0.2 ms. Here, μ = (Vth,spike − Vreset)/〈Ci〉 and Astim is the stimulation strength. (Vth,spike − Vreset) is the voltage difference between the maximum spiking threshold Vth,spike and the voltage reset Vreset. 〈Ci〉 is the mean membrane capacitance. See Appendix for more details and see Pyragas et al. (2018) for a discussion of optimal waveforms for DBS.
Figures 1A,C shows realizations of stimulus patterns for NCR and SNCR. The shaded areas mark possible stimulus onset times. Figures 1B,D shows the distribution of inter-stimulus intervals (ISTIs) for two values of σCR. Small values of σCR lead to a narrow distribution with peaks at integer multiples of TCR/Ns, while larger σCR lead to broader ISTI distributions.
Figure 1. Illustration of stimulation patterns used throughout the manuscript and the resulting distribution of inter-stimulus intervals. (A) Possible realization of NCR stimulation for Ns = 4 stimulation sites. Colored curves indicate the stimulation currents delivered to the individual sites. The pink region marks intervals of possible stimulus onset times, for the maximum jitter σCR = 1. The limit of vanishing jitter, σCR = 0, corresponds to deterministic stimulus onsets, i.e., CR stimulation. Vertical dashed lines separate subsequent CR cycles, with cycle period TCR = 1/fCR. (B) The distribution of inter-stimulus intervals (ISTIs) for NCR stimulation for two values of stimulus jitter, σCR. (C) Possible realization of SNCR stimulation with Ns = 4 stimulation sites. The color code is the same as in (A). (D) The distribution of ISTIs for SNCR stimulation for two values of σCR.
NCR with minimal temporal correlations, σCR = 1, leads to ISTIs between zero and 2TCR with mean ISTI TCR. For SNCR, the ISTI distribution attains its global maximum at TCR/Ns and decreases for larger ISTIs. For large jitters, σCR ≈ 1 it decays approximately exponentially, see Figure 1D. Consequently, individual sites may not receive any stimulus for multiple cycles.
2.4. Quantification of Synchronization
In order to quantify the degree of neuronal synchrony, we calculate the time-averaged Kuramoto order parameter (Kuramoto, 1984)
Here, tk = 10 s is an averaging interval and N is the number of neurons. ψk(t) is a phase function associated with the spiking of neuron k. ψk(t) attains subsequent integer values at subsequent spike times and increases linearly during interspike intervals (Rosenblum et al., 2001). indicates pronounced in-phase synchronization, whereas the absence of in-phase synchronized neuronal activity leads to .
2.5. Data Evaluation
In simulations, the data for the mean synaptic weight, 〈w〉, and the time-averaged Kuramoto order parameter, , were recorded every 10 s. Acute effects of stimulation are quantified using the acute mean synaptic weight, 〈w〉ac, which is the mean weight at the end of the stimulation period. For sufficiently low 〈w〉ac the network approaches a stable desynchronized state after cessation of stimulation (see Kromer et al., 2020; Khaledi-Nasab et al., 2021b).
In order to quantify the degree of acute synchronization, we time-averaged the Kuramoto order parameter, , over the last 10 s of the stimulation duration. Accordingly, we quantified the acute after-effect of stimulation by time-averaging the Kuramoto order parameter, , over the first 10 s interval after cessation of the stimulation. Lastly, long-lasting effects of stimulation are quantified by means of the time-averaged Kuramoto order parameter, , which is evaluated over a 10 s interval 1,000 s after the stimulation ceases. indicates long-lasting desynchronization, while indicates that stimulation did not entail long-lasting desynchronization effects.
2.6. Estimated Weight Change During Strong and Fast Stimulation
We derive estimates for the stimulation-induced weight dynamics during stimulation by applying the theoretical framework presented in Kromer and Tass (2020). In the following, we present the main steps and derive specific results for NCR, SCR, and SNCR. Results for CR were previously reported in Kromer et al. (2020) and are given below as a reference for the reader's convenience.
We consider a synapse with synaptic weight wij(t), presynaptic neuron i, and postsynaptic neuron j. Its mean rate of weight change during a time interval of duration T, starting at time t, is given by Kempter et al. (1999)
W(t) is the STDP function given in Equation (1). ti and tj are the pre- and postsynaptic spike times, respectively. The sum runs over all pairs of pre- and postsynaptic spike times, that contribute to weight changes according to the given STDP scheme (Morrison et al., 2008). In the present paper, we consider a nearest-neighbor scheme in which each presynaptic spike arrival time (ti + td) is paired with the latest postsynaptic spike time (tj) and each postsynaptic spike time is paired with the latest presynaptic spike arrival time.
We are particularly interested in the expectation value , which is obtained by averaging over different realizations of the stimulation pattern. Assuming stationary dynamics and the time interval T being long compared to the interspike intervals as well as the characteristic time scale of temporal correlations in the stimulation pattern, the mean rate of weight change becomes independent of the starting point t and the length T of the time interval. Then, we can approximate by its limit for long time intervals .
Next, we restrict ourselves to the case of stimulation-controlled spiking where each spike is caused by a stimulus and each stimulus causes a spike of the stimulated neurons. In simulations of the LIF network, this can be realized for weak synaptic interaction by considering strong stimulation (Astim ≈ 1) that is fast compared to neuronal firing rates in the absence of stimulation. Furthermore, the duration of inhibitory pulses νi should be short, such that the membrane potential recovers from inhibition before the next stimulus arrives. In the case of stimulation-controlled spiking, the spike times of pre- and postsynaptic neurons can be related to stimulus delivery times by introducing the distribution of spike response times λ(ϵ), where ϵ is the time lag between stimulus delivery and triggered neuronal spiking response (Kromer and Tass, 2020). If both post- and presynaptic neurons receive stimuli at rate fCR, these assumptions allow us to rewrite Equation (3) as
is the distribution of time lags between pairs of post- and presynaptic spike times that contribute to weight changes.
We calculate for CR, NCR, SCR, and SNCR. Following, with A=CR, NCR, SCR, or SNCR will denote the distribution of time lags during ongoing stimulation with the respective stimulation protocol. To this end, we consider the statistics of time lags t′ between subsequent post- and presynaptic spikes. We set t′ = S + Δ with the inter-stimulus interval S between stimuli triggering these spikes and Δ = ϵpost − ϵpre denoting the time difference between the realizations of spike response times for the considered spikes of the post- and presynaptic neurons resulting in the time lag t′. For narrow distributions of spike response times λ(t), in Equation (4) can be approximated using the two distributions: and as
see also Kromer and Tass (2020) and Kromer et al. (2020). is the distribution of inter-stimulus intervals between stimuli delivered to the post- and presynaptic neurons for which the resulting pairs of spike times contribute to weight changes according to the STDP scheme. Note that and pA(S|Δ) are in general not normalized to one as multiple intervals per spike can contribute to weight changes. For the STDP scheme considered in the present paper, both are normalized to two, i.e., and .
Estimates of the mean rate of weight change for the four different stimulation patterns considered in the present paper can be obtained by calculating . The latter depends on the conditional probability for the postsynaptic neuron to receive a stimulus at time S after a stimulus has been delivered to the presynaptic neuron and vice versa.
For the stimulus patterns considered in the present paper, two classes of synapses exist, each characterized by a distinct statistics of stimulus delivery times to the post- and presynaptic neurons, see Kromer and Tass (2020). Intrapopulation synapses connect neurons that are affected by the same stimulation site. Hence, post- and presynaptic neurons receive stimuli simultaneously. In contrast, interpopulation synapses connect neurons that are affected by different stimulation sites. Thus, post- and presynaptic neurons receive stimuli at different times. Following, we will replace the indices i and j, referring to the connected neurons, by the terms ‘intra' and ‘inter' referring to the respective type of synapse. Specifically, we set , if neurons i and j are in the same subpopulation, and , if neurons i and j are in different subpopulations.
Results for CR stimulation have been derived in Kromer and Tass (2020) and Kromer et al. (2020) and will be presented below for the reader's convenience. In the present paper, we consider the limit of sharp distributions of spike response times λ(ϵ) = δ(ϵ), thereby focusing on the contribution of the distribution of inter-stimulus intervals to weight changes. In this limit, we derive and for A = NCR, SCR, and SNCR stimulation.
2.6.1. CR and NCR Stimulation
We calculate the distributions of inter-stimulus intervals that contribute to weight changes for CR and NCR stimulation and , respectively.
In Kromer and Tass (2020) and Kromer et al. (2020), has been derived for the case of Ns = 4 and arbitrary Ns, respectively. Note that Ns was called M in the cited papers. Furthermore, in Kromer et al. (2020) a nonlinear dependence of on stimulation parameters was found. In more detail, if S + Δ < td presynaptic spikes arrive at the postsynaptic neuron after the next spikes have been initiated by stimuli. In order to account for this effect, a correction term was suggested. In the case of Λ(Δ) = δ(Δ) the authors' approach yields
Here, and are the zeroth order term () and the first correction term accounting for the case that presynaptic spikes arrive at postsynaptic neurons after the next but before the second to next stimulus has been delivered, i.e., inter-stimulus intervals are smaller then the delay time, but second order inter-stimulus intervals are larger than the delay time. In the present paper, we neglect further correction terms. These would account for more stimulus deliveries during a single delay time td.
The results of Kromer and Tass (2020) and Kromer et al. (2020) can be expanded to NCR by considering that individual stimulation times are uniformly distributed around their mean. Consequently, while is given by a superposition of delta distributions at multiples of , shows additional variability which can be described by the distribution qσCR(S) that is centered at zero. For a given width, σCR, of the distribution of jitters, we find
Equation (8) is the distribution of the difference between two random variables that are uniformly distributed in the interval [−σCR/(2NsfCR), σCR/(2NsfCR)]. These two random variables correspond to the jitters of post- and presynaptic stimulus delivery times. The case of CR stimulation (σCR = 0) can be described by q0(S): = δ(S), which yields the results of Kromer and Tass (2020).
Accordingly, we derive the first correction term using the results for CR stimulation published in Kromer et al. (2020)
We apply the same approach to the results for interpopulation synapses. Based on the results for CR stimulation from Kromer et al. (2020), we obtain
For the first correction term, we find
Here, we introduced Buv = Ns − u + v, Cuv = 1 + u + v, and Duv = M + u + v.
2.6.2. SCR and SNCR Stimulation
Next, we consider SCR and SNCR stimulation. As stimuli are delivered to random stimulation sites independently of previous stimulus deliveries, the derivation of and is comparable to that of the results for RR stimulation presented in Khaledi-Nasab et al. (2021a). Specifically, results from
with the series
The summands contain the probabilities 1/Ns and that the neurons' subpopulations receive the kth but not the previous k − 1 stimuli.
For interpopulation synapses, we find
The first correction terms and can be derived by considering that inter-stimulus intervals of length result in negative time lags for . First, we derive the correction term for intrapopulation synapse. Both post- and presynaptic neurons receive a stimulus at time t = 0. If the neurons do not receive a stimulus at time , the presynaptic spike arrival time 0+td is paired with the postsynaptic spike at time t = 0 for a negative time lag. However, with probability 1/Ns both neurons receive a stimulus at time . Then, the presynaptic spike arrival time td is paired with the postsynaptic spike at time for a negative time lag. Proceeding accordingly for positive time lags, we find the correction term for intrapopulation synapses
and the following correction terms for interpopulation synapses
The first correction term , with A = CR, NCR, SCR, or SNCR and X = “intra” or “inter,” needs to be added if postsynaptic spikes, triggered by the next stimulus, occur before the presynaptic spike, triggered by the current stimulus, arrives at the postsynaptic neuron. This scenario changes the order of postsynaptic spike times and presynaptic spike arrival times and occurs with probability . By considering the probability for this scenario, we can generalize Equation (6) to the case of NCR and SNCR stimulation as
with σCR = 0 for CR and SCR stimulation. Using in Equation (5) yields an estimate of , which can be used in Equation (4) to estimate the mean rate of weight change for stimulation protocol A and synapses of type X.
We compare long-lasting effects of CR, NCR, SCR, and SNCRstimulation patterns using theoretical estimates of the stimulation-induced weight dynamics and simulations of networks of 103 excitatory LIF neurons with STDP.
In simulations, the network connectivity and the neurons' membrane time capacitances were initialized at random. The networks were prepared in a state of pronounced in-phase synchronization (see section 2). For each network realization, we study the acute, acute after-, and long-lasting effects of stimulation by means of the mean synaptic weight 〈w〉 and the Kuramoto order parameter , Equation (2). Stimulation was delivered for either 500 or 1,000 s, see Figure captions.
3.1. Theoretical Predictions for Stimulation-Induced Weight Dynamics
To quantify the stimulation-induced synaptic weight dynamics, we consider the mean rate of weight change during ongoing stimulation. Theoretical estimates, , of the latter have been derived in the Methods section (see Equation 4). We compare the NCR/SNCR stimulation with maximum jitter, σCR = 1 to the ones without any jitter, CR/SCR.
Results for for CR, NCR, SCR, and SNCR stimulation are shown in Figure 2 for intrapopulation synapses, Figures 2A–D, and interpopulation synapses, Figures 2A'–D', respectively. We find qualitatively different dynamics of intra- and interpopulation synapses. Intrapopulation synapses weaken () for the considered range of stimulation frequencies fCR and numbers of stimulation sites Ns, see Figures 2A–D. In contrast, interpopulation synapses weaken only during slow stimulation with a rather small number of stimulation sites, Figures 2A'–D'. For fast stimulation with larger Ns, interpopulation synapses may strengthen ().
Figure 2. Theoretical estimates of stimulation-induced weight dynamics and distributions of time lags for different CR stimulation patterns. Individual columns correspond to CR (A), NCR (B), SCR (C), and SNCR (D) stimulation patterns. Panels show the estimated mean rates of weight change , Equation (4), for intrapopulation (A–D) and interpopulation synapses (A'–D'), respectively. Here, X = “intra,” “inter” marks the considered type of synapses. White curves mark zero contour lines and indicate the boundary between strengthening () and weakening of synapses (). Corresponding estimates for the distributions of time lags that lead to weight updates (black), Equation (5), are compared to simulation results (red) in (A”–D”) for intrapopulation synapses and in (A”'–D”') for interpopulation synaspes, respectively. In (A”–D”,A”'–D”'), we set Ns = 2 and fCR = 10 Hz. Networks were simulated for 90 s of ongoing stimulation. Time lags have been recorded from 400 pairs of pre- and postsynaptic neurons. Pairs were sorted according to synapse types “intra” and “inter” and histograms were calculated using a bin size of 1 ms. Theoretical estimates for were obtained by numerical calculations of , Equation (17). To this end, the time interval [−1, 000, 1, 000] ms was discretized using a binsize of dt = 0.01 ms. Then, GX(t) was obtained using Equation (5). To compare theoretical estimates and simulation results, we plotted GX(t)dt in (A”–D”,A”'–D”') and normalized the histograms such that counts summed up to two. Parameters: td = 3 ms, η = 0.02, τ+ = 10 ms, τR = 4, β = 1.4.
We observe qualitatively different dynamics of interpopulation weights for stimulation patterns with random jitters, i.e., NCRand SNCR, and stimulation patterns with deterministic stimulation times, i.e., CR and SCR, Figures 2A'–D'. For NCR and SNCR stimulation, strengthening of interpopulation synapses () is observed for fast stimulation, whereas we find a highly nonlinear dependence of the sign of for CR and SCR stimulation.
Considering the distributions of time lags that contribute to weight updates , Equation (5), we find that stimulation protocols with random jitters, i.e., NCR and SNCR, possess broad distributions spanning a wide range of possible time lags. In contrast, for protocols with deterministic stimulation times, i.e., CR and SCR, is given by a sum of delta-like distributions at integer multiples of TCR/Ns.
We find that simulated distributions of the time lags between post- and presynaptic spikes, which lead to weight updates, are well-described by theoretical estimates for all four stimulation patterns (for interpopulation synapses, see Figures 2A”–D”, and for intrapopulation synapses see Figures 2A”'–D”'). In particular, smooth distributions of time lags were found for stimulation protocols with temporal randomization, i.e., NCRand SNCR, whereas protocols with deterministic stimulation times led to time lags of multiples of TCR/Ns.
Based on theoretical predictions for the mean rate of weight change, Figures 2A–D, A'–D', we expect long-lasting desynchronization effects of CR, NCR, SCR, and SNCRin a large portion of the parameter space, spanned by fCR and Ns.
3.2. Acute, Acute After-, and Long-Lasting Effects of Strong Stimulation
We study the acute, acute after-, and long-lasting effects of strong stimulation for the four different stimulation patterns. The limit of strong stimulation (Astim = 1) resembles the case of stimulation-controlled spiking used during the derivation of theoretical predictions. Simulation results for the acute mean synaptic weight 〈w〉ac during NCR stimulation and NSCR stimulation for different values of σCR are shown in Figures 3A–E, 4A–E, respectively.
Figure 3. Acute, acute after-, and long-lasting effects of Noisy CR (NCR) with different values of the stimulus jitter, σCR, as a function of the stimulation frequency and the number of stimulation sites for strong stimulation. (A–E) Simulation results for the acute mean synaptic weights, 〈w〉ac, at the end of the 1, 000 s stimulation duration; (A'–E') The acute Kuramoto order parameter, , time-averaged over the last 10 s of the stimulation duration; (A”–E”) The acute after-effect on synchronization as quantified by the Kuramoto order parameter, , time-averaged over a 10 s interval after cessation of the stimulation; (A”'–E”') Long-lasting desynchronization effects for respective stimulus jitters, as quantified by the Kuramoto order parameter, Equation (2), averaged over a 10 s interval 1,000 s after cessation of stimulation, . (A,A',A”,A”') show results for σCR = 0 which are similar to Kromer et al. (2020) but for longer stimulation durations. The white curves show theoretical estimates of the boundary between weakening and strengthening of interpopulation synapses, see Figure 2. Parameters: The stimulation duration was set to Tstim = 1, 000 s and Astim = 1.
Figure 4. Acute, acute after-, and long-lasting effects of Shuffled Noisy CR (SNCR) for different values of the stimulus jitter, σCR, for strong stimulation as a function of the stimulation frequency and the number of stimulation sites. (A–E) Acute mean synaptic weights, 〈w〉ac; (A'–E') The acute Kuramoto order parameter, , time-averaged over the last 10 s of the stimulation; (A”–E”) The acute after-effect of stimulation on synchronization as measured by the Kuramoto order parameter, , time-averaged over an interval of 10 s right after cessation of the stimulation; (A”'–E”') Long-lasting effects of stimulation as quantified by the Kuramoto order parameter, , time-averaged over a 10 s interval 1, 000 s after cessation of stimulation. Parameters: The stimulation duration was set as Tstim = 1, 000 s and Astim = 1.
First, we consider the extreme cases of deterministic stimulus onset times, σCR=0, and maximum variability of stimulus onset times, σCR = 1. We find that strong stimulation leads to a reduction of the mean synaptic weight for rather slow stimulation with a small number of stimulation sites. In other parameter regions, the weight dynamics strongly depends on the stimulation pattern. For patterns with deterministic stimulus onset times, i.e., CR (Figure 3A) and SCR (Figure 4A) stimulations, we find a highly nonlinear dependence of the mean synaptic weight on the stimulation frequency and the number of stimulation sites. For CR stimulation, this was previously reported in Kromer et al. (2020). In contrast, for patterns with maximum random jitters, i.e., NCR and SNCR stimulations, we find that stimulation leads to a reduction of the mean synaptic weight for a wide range of stimulation frequencies and numbers of stimulation sites, see Figures 3E, 4E.
In order to study the long-lasting outcome of stimulation, we evaluate the Kuramoto order parameter 1,000 s after cessation of stimulation. Corresponding simulation results are shown in Figures 3A”',E”', 4A”',E”'. We find that stimulation entailed long-lasting desynchronization in regions with considerable weight reduction. This indicates that the system approached the stable desynchronized state after cessation of stimulation. In regions where the mean weight was only slightly reduced, or even increased, during the stimulation, we observe long-lasting synchronization. This indicates that the system reapproached the stable synchronized state for these parameter sets.
To compare theoretical predictions to simulation results, we show the estimated boundary between stimulation-induced weakening and strengthening of interpopulation synapses () in Figures 3, 4. We find that the boundary accurately reproduces the boundary between long-lasting desynchronization and long-lasting synchronization. Deviations occur mainly for low and intermediate stimulation frequencies.
Next, we analyze the degree of stimulation-induced synchronization during stimulation (acute effects), right after cessation of stimulation (acute after-effects), and long after cessation of stimulation (long-lasting effects). Figures 3, 4 show corresponding acute effects (second row), acute after-effects (third row), and long-lasting after effects (fourth row) as quantified by the Kuramoto order parameter, Equation (2), averaged over respective time intervals.
We find that acute stimulation-induced synchrony during NCR stimulation is independent of the stimulation parameters, such as the stimulation frequency, the number of stimulation sites, and the jitter σCR. Throughout the parameter space, NCR stimulation induces acute partial synchronization with ρac ≈ 0.3. Only for very slow stimulation (fCR ≈ 1 − 2 Hz), NCR stimulation induces acute desynchronization. Remarkably, as soon as stimulation ceases, the degree of synchronization changes rapidly and becomes determined by the underlying network connectivity, i.e., strong synaptic connections lead to synchronization and weak synaptic connections to desynchronization, see the third row of Figure 3.
Similarly, for SNCR stimulation, we find acute partial synchronization. However, the degree of in-phase synchronization attains lower values for larger Ns, see Figure 4. As for NCR stimulation, the degree of synchronization changes rapidly after cessation of SNCR stimulation and attains the values determined by the underlying network connectivity. Of note, for all stimulation patterns shown in Figures 3, 4, for very low stimulation frequencies, fCR → 1 Hz, there are no long-lasting effects for most values of Ns. There, the CR stimulation frequency, fCR, is much smaller than the frequency of the synchronized rhythm. Hence even if the stimulation potentially can decouple the network, a much longer stimulation duration is needed.
3.3. Effect of Jitter on Long-Lasting Effects
Next, we consider the effect of σCR, quantifying the variability of stimulus onset times. We perform simulations for intermediate values of 0 < σCR < 1 and evaluate the theoretical prediction of the boundary between strengthening and weakening of interpopulation weights () using Equations (4) and (5), and the results for the respective stimulation patterns given in section 2. Results for NCR stimulation and SNCR are given in Figures 3, 4, respectively.
We find that for both NCR and SNCR stimulations, the nonlinear dependence of the mean synaptic weight on the stimulation frequency and the number of stimulation sites becomes less pronounced as σCR increases, see Figures 3, 4. This trend is well reproduced by our theory. We observe a similar trend for the degree of long-lasting desynchronization as quantified by the Kuramoto order parameter, see the fourth row of Figures 3, 4.
To directly compare the mean synaptic weights for different values of σCR, we fix the number of stimulation sites to Ns = 24. Results for the mean synaptic weight for different σCR as a function of the stimulation frequency are plotted in Figure 5. We find that the effect of variability of stimulus onset times on the acute mean synaptic weight depends on the stimulation frequency. In particular, we find that increasing σCR has only a minor effect on the acute mean weight for NCR stimulation with low stimulation frequencies, see Figure 5A. In contrast, it leads to a significant reduction of the acute mean weight for intermediate stimulation frequencies, (fCR ≈ 12 Hz), whereas it leads to an increase of the acute mean weight for higher stimulation frequencies, see Figure 5A. For SNCR stimulation, we find qualitatively similar results, however, these trends occur at slightly lower stimulation frequencies, see Figure 5B.
Figure 5. Acute mean synaptic weight for five values of the stimulus jitter, σCR, as a function of the stimulation frequency for Ns = 24. (A,B) show results for Noisy CR (NCR) and Shuffled Noisy CR (SNCR), respectively. These graphs represent horizontal lines at Ns = 24 in Figures 3, 4, respectively. Parameters are the same as in Figures 3, 4.
Of particular interest for clinical application is to provide parameter ranges for which adding random jitters to the stimulus onset times might improve the performance of the stimulation. In order to derive such parameter ranges, we consider the difference of the acute mean synaptic weight for deterministic stimulus onset times 〈w〉ac(CR), and the acute mean synaptic weight for maximum variability (σCR = 1) 〈w〉ac(NCR). We also calculate the corresponding difference for SCR and SNCR.
Results are shown in Figures 6A,B for NCR and SNCR stimulation, respectively. In the red regions in Figure 6, jitter improves the reduction of the mean synaptic weight during stimulation. Accordingly, we calculate the differences between the Kuramoto order parameters, , quantifying the effect of jitter on the degree of long-lasting desynchronization. Corresponding results are shown in Figures 6A',B'.
Figure 6. Differences between the outcomes of stimulation with deterministic stimulus onset times (CR/SCR) and NCR/SNCR stimulation with maximum jitter, σCR = 1. The difference maps for the acute mean weight (A,B) and for the long-lasting Kuramoto order parameter (A',B'). Parameter regions where NCR/SNCR led to smaller mean weight/values of the long-lasting Kuramoto order parameter compared to CR/SCR are marked red. Dashed vertical lines enclose the largest continuous range of stimulation frequencies where NCR/SNCR stimulation with maximum jitter led to similar or better outcome than CR/SCR stimulation. These frequency ranges are referred to as , and in the text. Data are taken from panels A,E (acute mean weight) and A”',E”' (long-lasting Kuramoto order parameters) of Figures 3, 4, respectively.
The simulation results in Figure 6 show that adding random jitters to the stimulus onset times during CR and SCR improves the performance of these stimulation patterns for intermediate frequencies. Corresponding frequency intervals (NCR) and (SNCR) are shown in Figure 6. However, for very low and high frequencies, adding random jitters tends to worsen the outcome of the stimulation.
3.4. Acute and Long-Lasting Effects of Weak Stimulation
Next, we consider weak stimulation using the four different stimulation patterns; CR, SCR, and their noisy counterparts with maximum jitter, σCR = 1. To this end, we set Astim = 0.1 and perform a similar analysis as in the previous section.
Figures 7A–D show simulation results for the acute mean weight 〈w〉ac during weak stimulation for the four different stimulation patterns. We find that stimulation leads to a reduction of the mean synaptic weight (prior to stimulation 〈 w 〉 ≈ 0.38) in the major part of the parameter space spanned by the stimulation frequency and the number of stimulation sites. In contrast to strong stimulation, weak stimulation also leads to weight reduction for fast stimulation with large numbers of stimulation sites. In the maps of Figures 7A–D, the red colors indicate that the value of the mean synaptic weight is close to that prior to stimulation in the stable synchronized state. Overall weight reduction during weak stimulation is more pronounced for sufficiently fast stimulation.
Figure 7. Acute, acute after-, and long-lasting effects of weak stimulation for the four multisite stimulation protocols. Upper (A–D) show simulation results for the acute mean weight, 〈w〉ac; (A'–D') show the acute Kuramoto order parameter, , time-averaged over the last 10 s of the stimulation duration; (A”–D”) show the acute after-effect measured by the Kuramoto order parameter, , 10 s after cessation of the stimulation; and the bottom panels (A”'–D”') show the results for the Kuramoto order parameter, , time-averaged over a 10 s time interval 1, 000 s after cessation of stimulation. Low values of the Kuramoto order parameter indicate desynchronized activity while high values refer to pronounced in-phase synchronization. The white vertical lines in (A–D) mark a stimulation frequency of fCR = 12 Hz for which we present a detailed analysis of the influence of the stimulation duration Tstim in Figure 8. The frequency of the original synchronous rhythm is approximately 3.5 Hz and it is shown by the magenta dot-dashed vertical lines. Acute mean weights are measured at the end of Tstim = 500 s stimulation period. Parameters: Astim = 0.1.
Considering the differences between the four stimulation protocols, we find that protocols with deterministic stimulation times, i.e., CR and SCR stimulations, do not lead to a substantial reduction of the mean synaptic weight for the intermediate number of stimulation sites and a wide range of stimulation frequencies. In contrast, protocols with random jitters, i.e., NCR and SNCR stimulations, lead to a reduction of the mean synaptic weight for sufficiently fast stimulation, fCR ≿ 6 Hz for NCR and fCR ≿ 10 Hz SNCR stimulation, across the range of considered numbers of stimulation sites, see Figures 7A–D.
For a vast range of stimulation frequencies, we find a nonlinear dependence of 〈w〉ac on the number of stimulation sites. While 〈w〉ac(Ns) expresses multiple extrema for constant fCR for CR and SCR stimulations, we observe only one maximum for NCR and SNCR stimulations, see Figures 7A–D.
Next, we consider the acute, acute after-, and long-lasting effects of stimulation on synchronization. Simulation results are shown in Figure 7. Similar to the case of strong stimulation, we find acute partial synchronization during weak stimulation, see the second row of Figure 7. As stimulation ceases, the network approaches the state determined by the network connectivity, i.e., synchronized activity for strong and desynchronized activity for weak connections, see the third row of Figure 7. Remarkably, in parts of the parameter space where synaptic weights have not been reduced completely during the stimulation, we find that the Kuramoto order parameter increases right after cessation of stimulation, compare rows two and three of Figure 7. As one can see from the acute after-effects, (see Figures 7A”–D”), we find pronounced in-phase synchronization in parameter regions where stimulation did not lead to a substantial reduction of the mean synaptic weight, compare Figures 7A–D, A”'–D”'. As a consequence, stimulation protocols with random jitters, i.e., NCR and SNCR stimulation, lead to long-lasting desynchronization, , in a bigger portion of the parameter space than stimulation protocols without random jitters, i.e., CR and SCR stimulations.
3.5. Required Stimulation Duration for Long-Lasting Desynchronization
We study the influence of the stimulation duration Tstim on the acute and long-lasting outcome of weak stimulation with the four different stimulation patterns. To this end, we fix the stimulation frequency to fCR = 12 Hz, see dashed vertical lines in Figures 7A–D. We deliver stimulation for a stimulation duration Tstim and record the acute mean synaptic weight 〈w〉ac. Then, we turn off the stimulation and continue the simulation for 1, 000 s to evaluate long-lasting effects using .
Figures 8A–D shows simulation results for the acute mean weight 〈w〉ac recorded after the stimulation duration Tstim. We find that the mean weight 〈w〉ac reduces rapidly for small numbers of stimulation sites, Ns < 10. For larger Ns, the evolution of 〈w〉ac(Tstim) depends on the stimulation protocol. For CR and SCR stimulation, 〈w〉ac(Tstim) reduces rather rapidly for a narrow range of large 26 < Ns < 40, whereas it reduces slowly for others Ns. In contrast, the reduction of 〈w〉ac(Tstim) only slightly depends on Ns for NCR and SNCR stimulation.
Figure 8. Acute and long-lasting effects of stimulation for different stimulation durations. (A–D) acute mean synaptic weight, 〈w〉ac, at the end of a stimulation period Tstim for different numbers of stimulation sites Ns. (A'–D') Long-lasting desynchronization effects quantified by the Kuramoto order parameter, , Equation (2), recorded 1, 000 s after a stimulation period of duration Tstim. Low values indicate desynchronized spiking activity and high values in-phase synchronization of neuronal spiking. Parameters: Astim = 0.1 and fCR = 12 Hz.
3.6. Long-Lasting Desynchronization Effects Depend on Stimulation Strength
Typically, in clinical trials, the stimulation frequency and strength can be adjusted more readily, whereas the number of stimulation sites is constrained by anatomical and physiological features of the target area. Accordingly, we study the effect of the stimulation strength, Astim, and the stimulation frequency, fCR, for fixed numbers of stimulation sites. We study the long-lasting desynchronization effects, as quantified by the Kuramoto order parameter , for CR and SCR and their noisy counterparts with maximum jitter, σCR = 1. For Ns = 4, 8, 16, and 24, we vary the stimulation frequency and strength.
Figures 9A–D shows the long-lasting desynchronization effects for CR, and Figures 9A'–D' shows the same for NCR. The third row in Figures 9A”–D” shows the difference between the long-lasting Kuramoto order parameters of CR and NCR.
Figure 9. Long-lasting effects of CR and NCR stimulations as a function of CR frequency, fCR, and simulation strength, Astim. (A–D) The long-lasting Kuramoto order parameter for CR stimulation as a function of stimulation frequency and stimulation strength for Ns = 4, 8, 16, and 24. (A'–D') Same as the top row but for NCR with σCR = 1. (A”–D”) The difference between the long-lasting Kuramoto order parameter for CR and NCR, . In the red colored parameter regions, the Kuramoto order parameter for CR was larger than that for NCR, in the blue colored regions CR led to long-lasting desynchronization whereas NCR did not. The columns correspond to different numbers of stimulation sites, Ns = 4, 8, 16, and 24, respectively.
For weak stimulation strengths (Astim → 0) the stimulation is ineffective in inducing long-lasting desynchronization. For moderate stimulation strengths, both CR and NCR lead to desynchronization over a wide range of stimulation frequencies. However, for NCR, the range of effective moderate simulation strengths is larger compared to CR, as can be seen in the second row of Figure 9. By increasing the stimulation strength (Astim → 1), we approach the limit of strong stimulation (predicted by our theory, see Figure 2). Here, we find several frequency intervals in which stimulation does not lead to long-lasting desynchronization (see black regions in Figure 9). Remarkably, we find only one of these frequency intervals for NCR, whereas several intervals occur for CR. Put differently, random jitters remove the nonlinearities in the parameter space (Compare A–D with A'–D' in Figure 9). We find that the parameter region in which stimulation leads to long-lasting desynchronization is bigger for smaller numbers of stimulation sites, i.e., Ns = 4.
Figure 10 shows the long-lasting effect for SCR and SNCR. Here, we find qualitatively similar results as in Figure 9. Comparing the results for CR and SCR, we find that SCR performs slightly better than CR at low stimulation frequencies (fCR = 1 − 3 Hz) (see Figures 9A, 10A). In the limit of strong stimulation, frequency intervals in which stimulation does not lead to long-lasting desynchronization effects are larger for SCR and SNCR than for CR and NCR. However, this reverses if the stimulation strength is reduced. Thus, shuffling increases the range of effective stimulation frequencies for moderate to strong stimulation amplitudes.
Figure 10. Long-lasting effects of SCR and SNCR stimulation as a function of CR frequency, fCR, and simulation strength, Astim. Panels show the long-lasting Kuramoto order parameter for SCR (A–D) and for SNCR (A'–D'). (A”–D”) show the difference between the long-lasting Kuramoto order parameter for SCR and SNCR, . This figure is similar to Figure 9 but for SCR and SNCR.
In the present paper, we analyze the acute, acute after-, and long-lasting effects of randomized coordinated reset (CR) stimulation patterns on plastic neuronal networks. CR stimulation has been used in preclinical in vitro studies (Tass et al., 2009), preclinical in vivo studies (Tass et al., 2012b; Wang et al., 2016; Ho et al., 2021), as well as in clinical studies (Adamchic et al., 2014; Syrkin-Nikolau et al., 2018; Pfeifer et al., 2021) to induce acute and long-lasting desynchronization effects as well as symptom relief in the context of epilepsy, Parkinson's disease, and binge alcohol drinking. We computationally study the consequences of a reduction of temporal correlations between stimulus delivery times and a reduction of spatial correlations in the stimulus pattern on the efficacy of stimulation. A reduction of temporal correlations is achieved by adding random jitters to the deterministic stimulus delivery times of the original CR pattern. We denote the resulting stimulation pattern as CR with random jitters (NCR). A corresponding stimulation pattern has been used in a recent clinical study on vibrotactile stimulation of Parkinson's patients; where a long-lasting cumulative reduction of motor symptoms was observed (Pfeifer et al., 2021). The reduction of spatial correlations is achieved by shuffling the sequence of stimulated subpopulations of the original CR pattern. The corresponding stimulation pattern is referred to as shuffled CR (SCR). Lastly, we refer to a CR pattern with both random jitters and shuffling as shuffled noisy CR (SNCR). Corresponding stimulation patterns are illustrated in Figure 1. Our detailed theoretical and computational analysis reveals a significant increase in parameter robustness of long-lasting effects due to random jitters for intermediate stimulation frequencies, whereas shuffling reduces parameter robustness.
First, we consider the limit of strong stimulation for which we accurately predict the distribution of time lags between post- and presynaptic neurons by extending a theoretical framework previously presented in Kromer and Tass (2020), see Figure 2. In particular, we find marked differences between the stimulation-induced dynamics of intrapopulation and interpopulation synapses. Here, intrapopulation synapses refer to synapses that connect neurons at the same stimulation site, while interpopulation synapses connect neurons at different stimulation sites.
For all considered stimulation patterns, we found a stimulation-induced weakening of intrapopulation synapses. This is in accordance with previous results on the stimulation-induced weight dynamics in plastic neuronal networks (Kromer and Tass, 2020; Kromer et al., 2020; Khaledi-Nasab et al., 2021a,b; Pfeifer et al., 2021). These studies found that the dynamics of intrapopulation synapses is dominated by an effect called decoupling through synchrony (Lubenov and Siapas, 2008; Knoblauch et al., 2012). For sharp distributions of spike response times, this effect leads to a reduction of synaptic weights between simultaneously stimulated neurons in networks with axonal delays. For sufficiently long axonal delays, stimulus-elicited presynaptic spikes arrive after the postsynaptic ones which leads to pronounced synaptic depression. In our setup, this effect strongly contributes to the dynamics of the mean synaptic weight for small numbers of stimulation sites Ns, i.e., when the portion of intrapopulation synapses is high. Consequently, we observe a rapid reduction of the mean synaptic weight for small Ns, see, for instance, Figure 8.
In contrast, the dynamics of interpopulation synapses is more complex and it depends on the stimulation patterns. A detailed analysis for regular CR has been provided in Kromer et al. (2020). The authors revealed a nonlinear dependence of the mean rate of weight change on the stimulation frequency, fCR, and the number of stimulation sites used for stimulus deliveries, Ns. These nonlinearities result from a delay-induced effect that leads to a change in the order of postsynaptic spikes and presynaptic spike arrival times whenever k/fCRNs < td, i.e., when spikes triggered by the next stimulus tend to occur before stimulus-triggered presynaptic spikes arrive. Here k is a natural number referring to how many stimuli may be delivered during a single delay time td. This effect results in the complex pattern of synaptic weakening () and synaptic strengthening () of interpopulation synapses, see Figure 2A'. We find the corresponding pattern for SCR stimulation, i.e., for CR with shuffled sequence of stimulation sites, see Figure 2C'. However, these nonlinearities disappear if random jitters are added, i.e., for NCR or SNCR stimulation. For these stimulation patterns, we observe one continuous region in the parameter space where interpopulation weights reduce. This happens because the condition k/NsfCR < td is only satisfied by small portions of stimuli, due to random jitters. Considering these portions in our theory led to an accurate prediction of regions with synaptic strengthening and weakening for these patterns, see Figures 2B',C'.
Using our theoretical framework, we were able to accurately predict the distribution of time lags between post- and presynaptic spikes for strong stimulation. Results for intrapopulation synapses are presented in Figures 2A”–D” and results for interpopulation synapses in Figures 2A”'–D”'. Stimulus patterns with deterministic stimulus delivery times, i.e., CR and SCR, cause distributions with several peaks at integer multiples of 1/NsfCR, i.e., multiples of the minimal inter-stimulus interval, see Figures 2A”,A”',C”,C”'. Similar distributions have been reported in Kromer and Tass (2020); Kromer et al. (2020) for CR stimulation. Kromer et al. (2020) pointed out that stimulation patterns causing such distributions of time lags require an adjustment of the inter-stimulus intervals to the STDP time scales τ+ and τ−. Random jitters, however, lead to a smoothening of these distributions, see results for NCR and SNCR in Figure 2 for the maximum jitter, σCR = 1. In Kromer and Tass (2020), a smoothed distribution of time lags was also obtained for random reset (RR) stimulation, which combines temporal and spatial randomization by delivering stimuli at random times to randomly selected subpopulations. The authors compared the performance of RR stimulation for different STDP functions, including the one used in the present paper, Equation (1), and suggested that stimulation patterns that cause smoothed distributions of time lags, may lead to weight reduction that is more robust with respect to changes of the stimulation frequency. Our results support this suggestion. In particular, we found that NCR and SNCR stimulations lead to a pronounced weakening of interpopulation synapses in a large portion of the parameter space spanned by the stimulation frequency and the number of stimulation sites used for stimulus deliveries, see Figures 2–4.
We found that stimulation that leads to a pronounced weakening of interpopulation synapses may entail long-lasting desynchronization, see Figures 3, 4. In particular, for different values of the jitter, our theory accurately predicts the boundary between regions of stimulation parameters that lead to long-lasting desynchronization and those that lead to long-lasting synchronization. We find that deviations mainly occur for slow stimulation, i.e., where stimuli are administered at a slower pace than the original synchronous rhythm (≈3.5 Hz) and right next to the boundary, where weight reduction occurs at low rates and longer stimulation durations are required to drive the network into the attractor of the stable desynchronized state. In Figures 3, 4, we also found deviations for high frequencies (fCR > 8 Hz) and large numbers of stimulation sites. These are expected as our theory only considers first-order corrections, i.e., it is restricted to the case that only postsynaptic spikes (triggered by the next stimulus) may arrive before presynaptic spikes (triggered by the current stimulus). For fast stimulation that uses a large number of stimulation sites, however, multiple stimuli may trigger postsynaptic spikes before the presynaptic spike arrives at the postsynaptic neuron.
In contrast to the long-lasting effects of stimulation, we didn't find a strong parameter dependence of the degree of acute synchronization (see the second row of Figures 3, 4). During strong stimulation, the degree of synchronization is determined by the stimulus pattern rather than the synaptic weights. In particular, in the limit of strong stimulation, the degree of acute synchronization is independent of the stimulation frequency as long as the latter is fast compared to the synchronous rhythm, see the second row of Figures 3, 4. For NCR stimulation with various jitters, we find acute partial synchronization which is in accordance with previous results on acute partial synchronization during CR stimulation, see Figure 4 in Kromer et al. (2020). For slower NCR stimulation, we find a reduction of the acute Kuramoto order parameter indicating acute desynchronization. This is due to slow but strong stimulus deliveries that perturb the original synchronous rhythm. During SNCR stimulation, we also observed partial acute synchronization, however, the Kuramoto order parameter depends on the number of stimulation sites, see the second row of Figure 4. A similar effect was observed in studies on strong RR stimulation (Khaledi-Nasab et al., 2021b). There it was argued that a low number of stimulation sites leads to synchronous spiking responses of macroscopic neuronal subpopulations. Strong RR and SNCR stimulation possess qualitatively similar correlations between stimulus delivery times compare Equations (9, 10) in Khaledi-Nasab et al. (2021b) with Equations (12, 13) in the present paper. Therefore, we expect a similar dependence on the number of stimulation sites for strong stimulation. For weak stimulation, we find a slight dependence on the stimulation frequency. In particular, partial synchronization during CR and NCR stimulation attains a maximum when fCR is close to the frequency of the synchronous rhythm (The synchronous rhythm is shown by magenta vertical lines in Figure 7). This is similar to previous computation results on CR stimulation, where weak CR stimulation was inefficient when the stimulation frequency attained multiples of the frequency of the underlying synchronous rhythm (Kromer and Tass, 2020).
In contrast to the acute synchronization effects, the acute after-effects shown in the third rows of Figures 3, 4, 7, showed a great correspondence to the long-lasting effects. Note that in the case of weak stimulation, a slight increase of the Kuramoto order parameter was observed in regions where the synaptic weights were only slightly reduced during stimulation. Thus, based on the results in our computational model, we find that acute after-effects might be more suitable for predicting long-lasting effects than acute effects. This aligns with previous studies on acoustic CR for the treatment of tinnitus; where stimulation-induced acute and acute after-effects (Adamchic et al., 2017), as well as long-lasting effects (Tass et al., 2012a) were studied. The authors suggested that significant acute after-effects might be predictive of long-lasting symptom relief (Adamchic et al., 2017).
Remarkably, we find that both the synaptic weight reduction during and the long-lasting desynchronization by weak stimulation are more robust with respect to changes of the stimulation frequency and the number of stimulation sites than those of strong stimulation. In particular, for weak stimulation, we found that all four stimulation patterns cause a pronounced synaptic weight reduction and long-lasting desynchronization for high stimulation frequencies and large numbers of stimulation sites. This is in line with the results of Kromer et al. (2020), who observed this for the CR stimulation pattern. They argued that weak stimuli delivered shortly after neuronal spiking are not strong enough to elicit spikes, which leads to longer time lags between post- and presynaptic neurons than for strong stimuli. This, on the other hand, reduces the contribution of long-term potentiation to the synaptic weight dynamics, which occurs for short positive time lags. We observe similar effects for NCR, SCR, and SNCR stimulation, compare Figures 3, 4, 7. Earlier studies also observed that weaker stimulation is more suitable for long-lasting effects. However, they considered spatial stimulation profiles and accounted for the fact that weaker stimulation affects a smaller tissue volume. For strong stimulation all neuronal subpopulations fired in response to the stimulation, and this in turn led to a weaker desynchronization and may not induce decoupling (Popovych and Tass, 2012; Lysyansky et al., 2013; Zeitler and Tass, 2015).
We point out that the stimulation duration is another relevant parameter; that is particularly important when it comes to the long-lasting effects of the stimulation. While stimulation parameters such as the number of stimulation sites, the stimulation amplitude, and the stimulation frequency determine whether the stimulation leads to a reduction of the mean synaptic weight, the time it takes to achieve a sufficient weight reduction and drive the network into the attractor of a stable desynchronized state depends on the actual rate of synaptic weight reduction. In Figure 8, stimulation is capable of inducing long-lasting desynchronization effects for all considered numbers of stimulation sites. However, the required stimulation duration to drive the network into the attractor of a desynchronized state varies by more than one order of magnitude. Thus, for too short stimulation duration, stimulation might be considered ineffective for inducing long-lasting effects even though it would be well capable of inducing such effects for a longer stimulation duration. This is particularly important as, e.g., the current parameter adjustment procedures for standard DBS focus on acute effects (Volkmann et al., 2006), which may compromise the potential long-lasting effects.
Our detailed analysis provides evidence that a reduction of spatial correlations between stimulus deliveries does not increase the robustness of long-lasting desynchronization effects with respect to changes of the stimulation frequency, see Figures 3, 4, 7. This is in line with the computational result from earlier studies where SCR stimulation was compared to CR stimulation with fixed sequence (Tass and Hauptmann, 2009). The authors reported that the weight reduction by SCR stimulation was weaker than CR stimulation. The authors further observed a strengthening of synaptic weights by strong SCR stimulation. However, stimulation parameters such as the stimulation frequency, the number of stimulation sites, as well as the stimulation duration were not varied systematically. In contrast, we provide a systematic comparison of the stimulation-induced synaptic weight dynamics, as well as the potential long-lasting desynchronization effects of CR and SCR stimulation using theoretical and computational analysis. We find that strong SCR stimulation leads to a different distribution of time lags between post- and presynaptic spikes than CR stimulation. Our theoretical analysis shows that this leads to a stronger contribution of long-term potentiation to the synaptic weight dynamics for high numbers of sites for the STDP function considered in the present paper, see Figure 2. The relative performance of CR and SCR stimulation may, however, differ for other STDP rules. In particular, Tass and Hauptmann (2009) considered a symmetric plasticity function. A detailed study of the performance of CR stimulation for different plasticity functions was provided in Kromer and Tass (2020).
Our systematic analysis provides evidence that adding random jitters to the stimulus delivery times might improve the parameter robustness of synaptic weight reduction and long-lasting desynchronization effects at intermediate stimulation frequencies. This was observed for strong stimulation, where spike times are determined by the stimulation pattern, see Figures 3, 4, and for weak stimulation see Figure 7. These results partly confirm the hypothesis that a reduction of spatio-temporal correlations in stimulus patterns increases parameter robustness of long-lasting effects, Kromer and Tass (2020); however, it also shows that this effect results mainly from a reduction of temporal correlations, while a reduction of spatial correlations by shuffling did not lead to an increase of parameter robustness. Motivated by the results of Kromer and Tass (2020), a NCR pattern of vibrotactile stimuli was studied computationally and in Parkinson's patients (Pfeifer et al., 2021). The observed results were similar as for a corresponding CR pattern, however, the number of patients included was too small to draw conclusions about the relative performance of both stimulation patterns. In Pfeifer et al. (2021) the stimulation mechanism was very different compared to the electrical stimulation used here and the jitter was only moderate. Besides this one study, the NCR pattern has not been studied before. However, other randomized versions of the classic CR pattern where considered in Tass and Hauptmann (2009) and Zeitler and Tass (2018). In Tass and Hauptmann (2009), a multi-site randomly timed reset (MRTR) pattern was considered in which activation times of individual sites were generated from a Poisson process, see Figure 8B in Tass and Hauptmann (2009). Desynchronization and synaptic weight reduction by MRTR were found to be robust with respect to variations of the stimulation amplitude, however, longer stimulation than for CR stimulation was required to achieve this effect. In Zeitler and Tass (2018) an uncorrelated multichannel noisy stimulation (UMNS) protocol was considered in which individual sites were activated at random times during each CR cycle, see Figure 1 in Zeitler and Tass (2018). The authors found that UMNS was able to reduce synaptic weights for the considered parameter combinations. In both studies, the authors did not vary the stimulation frequency and the number of stimulation sites systematically which makes a comparison to our results difficult. However, as MRTR and UMNS were obtained by temporal variation of stimulation times and shuffling, both patterns might be comparable to our SNCR stimulation, which was obtained by adding random jitters to CR with shuffled sequences. For strong and weak SNCR stimulation, we find that indeed weight reduction and long-lasting desynchronization effects are robust with respect to changes of the stimulation frequency and the number of stimulation sites, see Figures 3, 4, 7. Furthermore, the absolute mean rates of weight change in the region of slow stimulation and small numbers of stimulation sites are lower than those for CR stimulation. Consequently, although weight reduction and long-lasting effects are more robust with respect to parameter changes, longer stimulation is required to drive the network in the attractor of a stable desynchronized state.
For strong stimulation, our detailed analysis of the impact of random jitters on the performance of NCR and SNCR stimulation showed that jitters improve the parameter robustness of synaptic decoupling, as well as long-lasting desynchronization mainly in a limited range of stimulation frequencies, see Figure 6. Corresponding frequency ranges were denoted as (for NCR) and (for SNCR) and span the range of one to four times the frequency of the synchronous rhythm. In our network model, the latter is associated with a frequency of about 3.5 Hz. We find that is larger than indicating that shuffling reduces frequency robustness in the presence of random jitters. The range of these frequencies depends on the considered type of synaptic plasticity as well as the range of synaptic transmission delays. As argued above, the latter determines nonlinearities in the synaptic weight dynamics as a function of the stimulation frequency and the number of stimulation sites (Kromer et al., 2020) which, as we found in the present study, can be suppressed by inducing random jitters. We note that random jitters mainly lead to a continuous area in the parameter space where the stimulation leads to long-lasting desynchronization, see Figures 3, 4.
For weak stimulation, Astim = 0.1, our results suggest that NCR with maximum variability of stimulus onset times (σCR = 1) may improve parameter robustness of long-lasting desynchronization effects in clinical studies for intermediate stimulation frequencies (one to four times the frequency of the pathological rhythm). In this range, random jitters expand the parameter region with effective long-lasting desynchronization toward high numbers of the stimulation site. However, the precise range depends on the underlying plasticity mechanism in the target area. As we have shown, the stimulation duration is a crucial factor, and a sufficiently long duration is needed to obtain long-lasting effects. Based on our results, one might be able to utilize the NCR stimulation in existing DBS electrodes to induce long-lasting effects.
In clinical trials, it is typically easier to vary stimulation parameters such as the stimulation pattern, the frequency, and the strength, rather than varying the number of stimulation sites. The latter are given by anatomical and physiological constraints of the target area. Hence, we studied long-lasting effects as a function of stimulation frequency and strength (see Figures 9, 10). In the limit of strong stimulation, Astim → 1, the results approach the prediction of our theory (see Figure 2), i.e., we find a nonlinear dependence of long-lasting desynchronization effects on the stimulation frequency. However, for moderate stimulation strengths, 0.1 ≾ Astim ≾ 0.5, we observe the highest degree of robustness with respect to variations of the stimulation frequency for all stimulation patterns (i.e., CR, SCR, NCR, and SNCR). Comparing CR and SCR with their noisy counterparts, NCR and SNCR, we find that including random jitters expands the parameter region with long-lasting desynchronization effects toward stronger stimulation.
Our results indicate that long-lasting desynchronization effects of stimulation with a small number of stimulation sites, e.g., Ns = 4, are robust with respect to changes of the stimulation frequency, even without random jitters. For larger numbers of stimulation sites (i.e, Ns = 8, 16, 24), the most robust long-lasting effects were observed for moderate stimulation strengths for all considered stimulation patterns. We found that jitter is particularly favorable for moderate to strong stimulation with larger numbers of stimulation sites.
Our promising results on the improved parameter robustness of long-lasting effects of NCR stimulation may trigger the question, whether randomized high-frequency stimulation (HF DBS) might be suitable to induce long-lasting therapeutic effects. To date, experimental studies on temporally randomized HF DBS mostly focus on acute intra-operative effects. Furthermore, it is still a matter of debate whether temporal randomization improves the acute effects of HF DBS. In a study by Brocker et al. (2013), it was reported that irregular HF DBS led to improved performance of PD patients in a finger-tapping task (Brocker et al., 2013). However, other studies reported that, in contrast to regular HF DBS, temporally randomized HF DBS was ineffective in providing symptom alleviation (Dorval et al., 2010; Birdno et al., 2012). To the best of our knowledge, no results on the long-lasting outcome of temporally randomized HF DBS have been presented to date.
In our recent study on CR stimulation, we observed that long-lasting effects are sensitive to the number of stimulation sites (Kromer et al., 2020). This suggest that the employment of recently developed DBS electrodes comprising large numbers of stimulation contacts (Krauss et al., 2020; Steigerwald et al., 2019) might actually require complicated parameter adjustment procedures. Here, we observe a similar parameter sensitivity for SCR stimulation, see Figures 3, 4, 7. However the long-lasting effects of stimulation patterns with random jitters, i.e., NCR and SNCR stimulation, are much more robust with respect to parameter changes, e.g., the number of stimulation sites used for stimulus delivery, see Figures 3, 4, 7. This provides evidence that NCR and SNCR stimulation might be more suitable for employing multisite stimulation electrodes for inducing long-lasting therapeutic effects. However, to date, there are no studies on the long-term use of such multisite stimulation electrodes available.
In DBS, the therapeutic outcome depends on the accurate placement of electrodes (Voges et al., 2002; Saint-Cyr et al., 2002). These electrodes typically have multiple stimulation contacts that, in some designs, can be activated independently (Krauss et al., 2020). Other designs allow for simultaneous co-activation of a subset of the stimulation sites to form an arbitrary number of co-activated subpopulations (Steigerwald et al., 2019). A common target area for DBS in Parkinson's disease is the STN (Krack et al., 2003). While the spatial extension of the entire STN is of the order of several millimeters to about one centimeter (Mavridis et al., 2013), the STN possesses a non-uniform somatotopic organization (Hartmann-von Monakow et al., 1978; Nambu et al., 1996; Nambu, 2011). In addition, based on incoming cortical projections, the STN is divided into sensorimotor, associative, and limbic parts (Nambu, 2011; Tewari et al., 2016). Depending on the patient's symptoms only a small portion of the STN is considered a suitable target area. For instance, a recent paper showed that while electrodes were placed such that a total of eight stimulation contacts where in the vicinity of the STN, only activation of a few contacts led to a reduction of pathological beta oscillations (Tamir et al., 2020). Thus, only a small number of contacts may be available for delivering multisite stimulation to the target area. In this context, it is encouraging that our results indicate that even multisite stimulation using just a small number of stimulation sites is capable of inducing long-lasting desynchronization effects that are robust with respect to variations of the stimulation frequency, fCR. However, our results also indicate that stimulation protocols with spatial randomization might be more suitable as they improve frequency-robustness for moderate to strong stimulations. Given the small anatomical target dimensions of, e.g., the STN, our results indicate that qualitatively different lead topologies comprising a larger number of smaller stimulation contacts with tighter spacing are not required to induce robust long-lasting desynchronization with the stimulation patterns studied here.
To conclude, multichannel CR stimulation has shown great promise for inducing long-lasting desynchronization and therapeutic effects when delivered through DBS electrodes (Tass et al., 2012b; Adamchic et al., 2014; Wang et al., 2016) or through vibrotactile fingertip stimulators (Pfeifer et al., 2021). Our results suggest that the robustness of long-lasting effects with respect to changes of the stimulation frequency and other parameters might be increased by adding jitters to the stimulus delivery times. In contrast, shuffling the sequence of stimulus deliveries does not increase parameter robustness. In our model, acute after-effects of stimulation are strongly correlated with long-lasting effects. This suggests that acute after-effects might be suitable for predicting long-lasting effects of the stimulation. We hope that our results will lead to more clinical studies on stimulation protocols with random jitters to improve parameter adjustment procedures for brain stimulation as a treatment for neurological disorders.
Data Availability Statement
The original contributions generated for the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
AKN, JK, and PT conceived the idea, designed the study, interpreted the results, and wrote the manuscript. AKN performed the numerical simulations and analysis. JK developed the theory. All authors contributed to the article and approved the submitted version.
We gratefully acknowledge support of this study by Boston Scientific Neuromodulation (Stanford Project 127674) and the Foundation for OCD Research (New Venture Fund, 011665-2020-08-01).
Conflict of Interest
PT works as a consultant for Boston Scientific Neuromodulation.
The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
We are grateful to Stanford University and Stanford's Sherlock Computing cluster for computational resources and support that contributed to these research results.
Adamchic, I., Hauptmann, C., Barnikol, U. B., Pawelczyk, N., Popovych, O., Barnikol, T. T., et al. (2014). Coordinated reset neuromodulation for parkinson's disease: proof-of-concept study. Mov. Disord. 29, 1679–1684. doi: 10.1002/mds.25923
Adamchic, I., Toth, T., Hauptmann, C., Walger, M., Langguth, B., Klingmann, I., et al. (2017). Acute effects and after-effects of acoustic coordinated reset neuromodulation in patients with chronic subjective tinnitus. Neuroomage Clin. 15, 541–558. doi: 10.1016/j.nicl.2017.05.017
Bashan, A., Bartsch, R. P., Kantelhardt, J. W., Havlin, S., and Ivanov, P. C. (2012). Network physiology reveals relations between network topology and physiological function. Nat. Commun. 3, 1–9. doi: 10.1038/ncomms1705
Benabid, A., Pollak, P., Hoffmann, D., Gervason, C., Hommel, M., Perret, J., et al. (1991). Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus. Lancet 337, 403–406. doi: 10.1016/0140-6736(91)91175-T
Bi, G.-Q., and Poo, M.-M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. J. Neurosci. 18, 10464–10472. doi: 10.1523/JNEUROSCI.18-24-10464.1998
Birdno, M. J., Kuncel, A. M., Dorval, A. D., Turner, D. A., Gross, R. E., and Grill, W. M. (2012). Stimulus features underlying reduced tremor suppression with temporally patterned deep brain stimulation. J. Neurophysiol. 107, 364–383. doi: 10.1152/jn.00906.2010
Bourne, S. K., Eckhardt, C. A., Sheth, S. A., and Eskandar, E. N. (2012). Mechanisms of deep brain stimulation for obsessive compulsive disorder: effects upon cells and circuits. Front. Integr. Neurosci. 6:29. doi: 10.3389/fnint.2012.00029
Brocker, D. T., Swan, B. D., Turner, D. A., Gross, R. E., Tatter, S. B., Koop, M. M., et al. (2013). Improved efficacy of temporally non-regular deep brain stimulation in parkinson's disease. Exp. Neurol. 239, 60–67. doi: 10.1016/j.expneurol.2012.09.008
Contarino, M. F., Bour, L. J., Bot, M., van den Munckhof, P., Speelman, J. D., Schuurman, P. R., et al. (2012). Tremor-specific neuronal oscillation pattern in dorsal subthalamic nucleus of parkinsonian patients. Brain. Stimul. 5, 305–314. doi: 10.1016/j.brs.2011.03.011
Dorval, A. D., Kuncel, A. M., Birdno, M. J., Turner, D. A., and Grill, W. M. (2010). Deep brain stimulation alleviates parkinsonian bradykinesia by regularizing pallidal activity. J. Neurophysiol. 104, 911–921. doi: 10.1152/jn.00103.2010
Hartmann-von Monakow, K., Akert, K., and Künzle, H. (1978). Projections of the precentral motor cortex and other cortical areas of the frontal lobe to the subthalamic nucleus in the monkey. Exp. Brain Res. 33, 395–403.
Hauptmann, C., Popovych, O., and Tass, P. A. (2005a). Multisite coordinated delayed feedback for an effective desynchronization of neuronal networks. Stochastics Dyn. 5, 307–319. doi: 10.1142/S0219493705001420
Hauptmann, C., Popovych, O., and Tass, P. A. (2005c). Effectively desynchronizing deep brain stimulation based on a coordinated delayed feedback stimulation via several sites: a computational study. Biol. Cybern. 93, 463–470. doi: 10.1007/s00422-005-0020-1
Ho, A. L., Feng, A. Y., Barbosa, D. A., Wu, H., Smith, M. L., Malenka, R. C., et al. (2021). Accumbens coordinated reset stimulation in mice exhibits ameliorating aftereffects on binge alcohol drinking. Brain Stimul. 14, 330–334. doi: 10.1016/j.brs.2021.01.015
Holt, A. B., Wilson, D., Shinn, M., Moehlis, J., and Netoff, T. I. (2016). Phasic burst stimulation: a closed-loop approach to tuning deep brain stimulation parameters for Parkinson's disease. PLoS Comput. Biol. 12:e1005011. doi: 10.1371/journal.pcbi.1005011
Jakobs, M., Fomenko, A., Lozano, A. M., and Kiening, K. L. (2019). Cellular, molecular, and clinical mechanisms of action of deep brain stimulation–a systematic review on established indications and outlook on future developments. EMBO Mol. Med. 11:e9575. doi: 10.15252/emmm.201809575
Karbowski, J., and Ermentrout, G. B. (2002). Synchrony arising from a balanced synaptic plasticity in a network of heterogeneous neural oscillators. Phys. Rev. E 65:031902. doi: 10.1103/PhysRevE.65.031902
Khaledi-Nasab, A., Kromer, J., and Tass, P. A. (2021a). Nonlinear dependence of desynchronization effects of coordinated reset on the number of stimulation sites and frequency. Bull. Am. Phys. Soc. Available online at: https://meetings.aps.org/Meeting/MAR21/Session/P08.8
Khaledi-Nasab, A., Kromer, J., and Tass, P. A. (2021b). Long-lasting desynchronization of plastic neural networks by random reset stimulation. Front. Physiol. 11:622620. doi: 10.3389/fphys.2020.622620
Knoblauch, A., Hauser, F., Gewaltig, M.-O., Körner, E., and Palm, G. (2012). Does spike-timing-dependent synaptic plasticity couple or decouple neurons firing in synchrony? Front. Comput. Neurosci. 6:55. doi: 10.3389/fncom.2012.00055
Krack, P., Batir, A., Van Blercom, N., Chabardes, S., Fraix, V., Ardouin, C., et al. (2003). Five-year follow-up of bilateral stimulation of the subthalamic nucleus in advanced parkinson's disease. N. Engl. J. Med. 349, 1925–1934. doi: 10.1056/NEJMoa035275
Krauss, J. K., Lipsman, N., Aziz, T., Boutet, A., Brown, P., Chang, J. W., et al. (2020). Technology of deep brain stimulation: current status and future directions. Nat. Rev. Neurol. 17, 75–87. doi: 10.1038/s41582-020-00426-z
Kromer, J. A., Khaledi-Nasab, A., and Tass, P. A. (2020). Impact of number of stimulation sites on long-lasting desynchronization effects of coordinated reset stimulation. Chaos 30, 083134. doi: 10.1063/5.0015196
Kühn, A. A., Kupsch, A., Schneider, G.-H., and Brown, P. (2006). Reduction in subthalamic 8-35 hz oscillatory activity correlates with clinical improvement in parkinson's disease. Eur. J. Neurosci. 23, 1956–1960. doi: 10.1111/j.1460-9568.2006.04717.x
Liu, K. K., Bartsch, R. P., Lin, A., Mantegna, R. N., and Ivanov, P. C. (2015). Plasticity of brain wave network interactions and evolution across physiologic states. Front. Neural Circ. 9:62. doi: 10.3389/fncir.2015.00062
Lourens, M. A., Schwab, B. C., Nirody, J. A., Meijer, H. G., and van Gils, S. A. (2015). Exploiting pallidal plasticity for stimulation in parkinson's disease. J. Neural Eng. 12:026005. doi: 10.1088/1741-2560/12/2/026005
Lozano, A. M., Lipsman, N., Bergman, H., Brown, P., Chabardes, S., Chang, J. W., et al. (2019). Deep brain stimulation: current challenges and future directions. Nat. Rev. Neurol. 15, 148–160. doi: 10.1038/s41582-018-0128-2
Maistrenko, Y. L., Lysyansky, B., Hauptmann, C., Burylko, O., and Tass, P. A. (2007). Multistability in the kuramoto model with synaptic plasticity. Phys. Rev. E 75:066207. doi: 10.1103/PhysRevE.75.066207
Manos, T., Zeitler, M., and Tass, P. A. (2018). How stimulation frequency and intensity impact on the long-lasting effects of coordinated reset stimulation. PLoS Comput. Biol. 14:e1006113. doi: 10.1371/journal.pcbi.1006113
Moreau, C., Defebvre, L., Destee, A., Bleuse, S., Clement, F., Blatt, J., et al. (2008). Stn-dbs frequency effects on freezing of gait in advanced parkinson disease. Neurology 71, 80–84. doi: 10.1212/01.wnl.0000303972.16279.46
Mormann, F., Lehnertz, K., David, P., and Elger, C. E. (2000). Mean phase coherence as a measure for phase synchronization and its application to the eeg of epilepsy patients. Phys. D 144, 358–369. doi: 10.1016/S0167-2789(00)00087-7
Nambu, A., Takada, M., Inase, M., and Tokuno, H. (1996). Dual somatotopical representations in the primate subthalamic nucleus: evidence for ordered but reversed body-map transformations from the primary motor cortex and the supplementary motor area. J. Neurosci. 16, 2671–2683. doi: 10.1523/JNEUROSCI.16-08-02671.1996
Nini, A., Feingold, A., Slovin, H., and Bergman, H. (1995). Neurons in the globus pallidus do not show correlated activity in the normal monkey, but phase-locked oscillations appear in the mptp model of parkinsonism. J. Neurophysiol. 74, 1800–1805. doi: 10.1152/jn.19126.96.36.1990
Pfeifer, K. J., Kromer, J. A., Cook, A. J., Hornbeck, T., Lim, E. A., Mortimer, B., et al. (2021). Coordinated reset vibrotactile stimulation induces sustained cumulative benefits in parkinson's disease. Front. Physiol. 12:200. doi: 10.3389/fphys.2021.624317
Popovych, O. V., Hauptmann, C., and Tass, P. A. (2006b). Desynchronization and decoupling of interacting oscillators by nonlinear delayed feedback. Int. J. Bifurcation Chaos 16, 1977–1987. doi: 10.1142/S0218127406015830
Popovych, O. V., Lysyansky, B., Rosenblum, M., Pikovsky, A., and Tass, P. A. (2017a). Pulsatile desynchronizing delayed feedback for closed-loop deep brain stimulation. PLoS ONE 12:e0173363. doi: 10.1371/journal.pone.0173363
Popovych, O. V., Lysyansky, B., and Tass, P. A. (2017b). Closed-loop deep brain stimulation by pulsatile delayed feedback with increased gap between pulse phases. Sci. Rep. 7:1033. doi: 10.1038/s41598-017-01067-x
Pyragas, K., Fedaravičius, A. P., Pyragienė, T., and Tass, P. A. (2018). Optimal waveform for entrainment of a spiking neuron with minimum stimulating charge. Phys. Rev. E 98:042216. doi: 10.1103/PhysRevE.98.042216
Raethjen, J., Lindemann, M., Schmaljohann, H., Wenzelburger, R., Pfister, G., and Deuschl, G. (2000). Multiple oscillators are causing parkinsonian and essential tremor. Mov. Disord. 15, 84–94. doi: 10.1002/1531-8257(200001)15:1<84::AID-MDS1014>3.0.CO;2-K
Rodriguez-Oroz, M. C., Obeso, J., Lang, A., Houeto, J.-L., Pollak, P., Rehncrona, S., et al. (2005). Bilateral deep brain stimulation in parkinson's disease: a multicentre study with 4 years follow-up. Brain 128, 2240–2249. doi: 10.1093/brain/awh571
Rosenblum, M., and Pikovsky, A. (2004a). Delayed feedback control of collective synchrony: an approach to suppression of pathological brain rhythms. Phys. Rev. E 70, 041904. doi: 10.1103/PhysRevE.70.041904
Rosenblum, M., Pikovsky, A., Kurths, J., Schäfer, C., and Tass, P. A. (2001). “Phase synchronization: from theory to data analysis,” in Neuro-Informatics and Neural Modelling, Handbook of Biological Physics, Vol. 4, eds S. Gielen and F. Moss (Amsterdam: Elsevier), 279–321.
Ruscio, A. M., Stein, D. J., Chiu, W. T., and Kessler, R. C. (2010). The epidemiology of obsessive-compulsive disorder in the national comorbidity survey replication. Mol. Psychiatry 15, 53–63. doi: 10.1038/mp.2008.94
Saint-Cyr, J. A., Hoque, T., Pereira, L. C., Dostrovsky, J. O., Hutchison, W. D., Mikulis, D. J., et al. (2002). Localization of clinically effective stimulating electrodes in the human subthalamic nucleus on magnetic resonance imaging. J. Neurosurg. 97, 1152–1166. doi: 10.3171/jns.2002.97.5.1152
Schrader, C., Capelle, H.-H., Kinfe, T., Blahak, C., Bäzner, H., Lütjens, G., et al. (2011). Gpi-dbs may induce a hypokinetic gait disorder with freezing of gait in patients with dystonia. Neurology 77, 483–488. doi: 10.1212/WNL.0b013e318227b19e
Steigerwald, F., Potter, M., Herzog, J., Pinsker, M., Kopper, F., Mehdorn, H., et al. (2008). Neuronal activity of the human subthalamic nucleus in the parkinsonian and nonparkinsonian state. J. Neurophysiol. 100, 2515–2524. doi: 10.1152/jn.90574.2008
Syrkin-Nikolau, J., Neuville, R., O'Day, J., Anidi, C., Koop, M. M., Martin, T., et al. (2018). Coordinated reset vibrotactile stimulation shows prolonged improvement in parkinson's disease. Mov. Disord. 33, 179. doi: 10.1002/mds.27223
Tamir, I., Wang, D., Chen, W., Ostrem, J. L., Starr, P. A., and de Hemptinne, C. (2020). Eight cylindrical contact lead recordings in the subthalamic region localize beta oscillations source to the dorsal stn. Neurobiol. Dis. 146:105090. doi: 10.1016/j.nbd.2020.105090
Tass, P. A. (2003a). Desynchronization by means of a coordinated reset of neural sub-populationsa novel technique for demand-controlled deep brain stimulation. Progr. Theor. Phys. Suppl. 150, 281–296. doi: 10.1143/PTPS.150.281
Tass, P. A., Adamchic, I., Freund, H.-J., von Stackelberg, T., and Hauptmann, C. (2012a). Counteracting tinnitus by acoustic coordinated reset neuromodulation. Restor. Neurol Neurosci. 30, 137–159. doi: 10.3233/RNN-2012-110218
Tass, P. A., Qin, L., Hauptmann, C., Dovero, S., Bezard, E., Boraud, T., et al. (2012b). Coordinated reset has sustained aftereffects in parkinsonian monkeys. Ann. Neurol. 72, 816–820. doi: 10.1002/ana.23663
Tass, P. A., Silchenko, A. N., Hauptmann, C., Barnikol, U. B., and Speckmann, E.-J. (2009). Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation. Phys. Rev. E 80:011902. doi: 10.1103/PhysRevE.80.011902
Tass, P. A., Smirnov, D., Karavaev, A., Barnikol, U., Barnikol, T., Adamchic, I., et al. (2010). The causal relationship between subcortical local field potential oscillations and parkinsonian resting tremor. J. Neural Eng. 7:016009. doi: 10.1088/1741-2560/7/1/016009
Temel, Y., Kessels, A., Tan, S., Topdag, A., Boon, P., and Visser-Vandewalle, V. (2006). Behavioural changes after bilateral subthalamic stimulation in advanced parkinson disease: a systematic review. Parkinsonism Relat. Disord. 12, 265–272. doi: 10.1016/j.parkreldis.2006.01.004
Temperli, P., Ghika, J., Villemure, J.-G., Burkhard, P., Bogousslavsky, J., and Vingerhoets, F. (2003). How do parkinsonian signs return after discontinuation of subthalamic DBS? Neurology 60, 78–81. doi: 10.1212/wnl.60.1.78
Tewari, A., Jog, R., and Jog, M. S. (2016). The striatum and subthalamic nucleus as independent and collaborative structures in motor control. Front. Syst. Neurosci. 10:17. doi: 10.3389/fnsys.2016.00017
Tripoliti, E., Zrinzo, L., Martinez-Torres, I., Tisch, S., Frost, E., Borrell, E., et al. (2008). Effects of contact location and voltage amplitude on speech and movement in bilateral subthalamic nucleus deep brain stimulation. Mov. Disord. 23, 2377–2383. doi: 10.1002/mds.22296
Vicheva, P., Butler, M., and Shotbolt, P. (2020). Deep brain stimulation for obsessive-compulsive disorder: a systematic review of randomised controlled trials. Neurosci. Biobehav. Rev. 109, 129–138. doi: 10.1016/j.neubiorev.2020.01.007
Voges, J., Volkmann, J., Allert, N., Lehrke, R., Koulousakis, A., Freund, H.-J., et al. (2002). Bilateral high-frequency stimulation in the subthalamic nucleus for the treatment of parkinson disease: correlation of therapeutic effect with anatomical electrode position. J. Neurosurg. 96, 269–279. doi: 10.3171/jns.2002.96.2.0269
Wang, J., Nebeck, S., Muralidharan, A., Johnson, M. D., Vitek, J. L., and Baker, K. B. (2016). Coordinated reset deep brain stimulation of subthalamic nucleus produces long-lasting, dose-dependent motor improvements in the 1-methyl-4-phenyl-1, 2, 3, 6-tetrahydropyridine non-human primate model of parkinsonism. Brain Stimul. 9, 609–617. doi: 10.1016/j.brs.2016.03.014
Warman, E., and Durand, D. (1989). “Desynchronization of epileptiform activity by phase resetting,” in Images of the Twenty-First Century. Proceedings of the Annual International Engineering in Medicine and Biology Society (Seattle, WA: IEEE), 1286–1287.
Weinberger, M., Mahant, N., Hutchison, W. D., Lozano, A. M., Moro, E., Hodaie, M., et al. (2006). Beta oscillatory activity in the subthalamic nucleus and its relation to dopaminergic response in parkinson's disease. J. Neurophysiol. 96, 3248–3256. doi: 10.1152/jn.00697.2006
Xie, T., Padmanaban, M., Bloom, L., MacCracken, E., Bertacchi, B., Dachman, A., et al. (2017). Effect of low versus high frequency stimulation on freezing of gait and other axial symptoms in parkinson patients with bilateral stn dbs: a mini-review. Transl. Neurodegener. 6, 1–10. doi: 10.1186/s40035-017-0083-7
Neuronal Network Model
In our LIF model the spiking occurs when the membrane potential vi crosses a dynamic threshold potential ; where the membrane and the threshold potentials of neuron i obeys:
Here is the synaptic input current, the stimulation current, the noisy input current, Ci is the membrane capacitance, τth is the threshold time constant and the stationary threshold potential, vrest the resting potential, and gleak is the leak conductance.
Spike times are recorded at each threshold crossing. After a spike, the membrane potential is set to vspike for a time period of tspike; then, we perform an instantaneous reset: vi→vreset and .
We use the same parameters as in Kromer et al. (2020) and Khaledi-Nasab et al. (2021b): vrest = −38 mV, τth = 5 ms, mV, tspike = 1 ms, vspike = 20 mV, Vreset = −67 mV, and gleak = 0.02 mS/cm2. The capacitances Ci follow a normal distribution with mean value of 〈Ci〉 = 3 μF/cm2 and standard deviation of 0.05〈Ci〉.
Each neuron i receives noisy input, , which is obtained by feeding presynaptic Poisson spike trains with firing rate fnoise = 20 Hz into excitatory synapses (Ebert et al., 2014)
Here the noise intensity is controlled by the parameter D = 0.026 mS/cm2 scaling the strength of the Poisson input. vsyn = 0 mV is the synaptic reversal potential, τsyn = 1 ms the synaptic time scale, and the synaptic conductance.
is the excitatory synaptic input to neuron i, and it is given by
where κ = 8 mS/cm2 is the coupling strength, is the synaptic conductance, wji ∈ [0, 1] is the weight of the synapses between presynaptic neuron j and postsynaptic neuron i. The first sum runs over all presynaptic neurons, and the second sum runs over the presynaptic spikes of neuron j. We consider homogeneous synaptic delays of td = 3 ms.
Keywords: coordinated reset stimulation, spike-timing-dependent plasticity (STDP), random jitter, long-lasting desynchronization, stimulation-induced decoupling
Citation: Khaledi-Nasab A, Kromer JA and Tass PA (2021) Long-Lasting Desynchronization Effects of Coordinated Reset Stimulation Improved by Random Jitters. Front. Physiol. 12:719680. doi: 10.3389/fphys.2021.719680
Received: 02 June 2021; Accepted: 12 August 2021;
Published: 24 September 2021.
Edited by:Andras Eke, Semmelweis University, Hungary
Reviewed by:Sebastiano Stramaglia, University of Bari Aldo Moro, Italy
Ronny P. Bartsch, Bar-Ilan University, Israel
Copyright © 2021 Khaledi-Nasab, Kromer and Tass. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ali Khaledi-Nasab, email@example.com