Original Research ARTICLE
Multisite Delayed Feedback for Electrical Brain Stimulation
- 1Institute of Neuroscience and Medicine, Brain and Behaviour (INM-7), Research Centre Jülich, Jülich, Germany
- 2Department of Neurosurgery, Stanford University, Stanford, CA, United States
Demand-controlled deep brain stimulation (DBS) appears to be a promising approach for the treatment of Parkinson's disease (PD) as revealed by computational, pre-clinical and clinical studies. Stimulation delivery is adapted to brain activity, for example, to the amount of neuronal activity considered to be abnormal. Such a closed-loop stimulation setup might help to reduce the amount of stimulation current, thereby maintaining therapeutic efficacy. In the context of the development of stimulation techniques that aim to restore desynchronized neuronal activity on a long-term basis, specific closed-loop stimulation protocols were designed computationally. These feedback techniques, e.g., pulsatile linear delayed feedback (LDF) or pulsatile nonlinear delayed feedback (NDF), were computationally developed to counteract abnormal neuronal synchronization characteristic for PD and other neurological disorders. By design, these techniques are intrinsically demand-controlled methods, where the amplitude of the stimulation signal is reduced when the desired desynchronized regime is reached. We here introduce a novel demand-controlled stimulation method, pulsatile multisite linear delayed feedback (MLDF), by employing MLDF to modulate the pulse amplitude of high-frequency (HF) DBS, in this way aiming at a specific, MLDF-related desynchronizing impact, while maintaining safety requirements with the charge-balanced HF DBS. Previously, MLDF was computationally developed for the control of spatio-temporal synchronized patterns and cluster states in neuronal populations. Here, in a physiologically motivated model network comprising neurons from subthalamic nucleus (STN) and external globus pallidus (GPe), we compare pulsatile MLDF to pulsatile LDF for the case where the smooth feedback signals are used to modulate the amplitude of charge-balanced HF DBS and suggest a modification of pulsatile MLDF which enables a pronounced desynchronizing impact. Our results may contribute to further clinical development of closed-loop DBS techniques.
High-frequency (HF) deep brain stimulation (DBS) is the standard therapy for medically refractory Parkinson's disease (PD), where electrical pulse trains are permanently delivered via depth electrodes at high frequencies (> 100 Hz) (Benabid et al., 1991, 2009; Kuncel and Grill, 2004; Johnson et al., 2008). The mechanism of action of HF DBS is still debated (Johnson et al., 2008; Gradinaru et al., 2009; Deniau et al., 2010). HF DBS may cause side effects by stimulation of the target as well as surrounding structures (Ferraye et al., 2008; Moreau et al., 2008; van Nuenen et al., 2008; Xie et al., 2012). It is, hence, desirable to reduce the integral current delivered. Accordingly, different types of closed-loop, demand-controlled and adaptive DBS (aDBS) have been developed in computational and engineering studies (Tass, 2001, 2003; Rosenblum and Pikovsky, 2004a,b; Hauptmann et al., 2005a,b; Popovych et al., 2005, 2006; Kiss et al., 2007; Pyragas et al., 2007; Tukhlina et al., 2007; Luo et al., 2009; Popovych and Tass, 2010; Montaseri et al., 2013). Closed-loop aDBS approach received recent development in pre-clinical and clinical studies (Graupe et al., 2010; Rosin et al., 2011; Carron et al., 2013; Little et al., 2013; Priori et al., 2013; Yamamoto et al., 2013; Grahn et al., 2014; Hosain et al., 2014; Rosa et al., 2015).
Closed-loop aDBS was successfully tested in parkinsonian monkeys under acute conditions (Rosin et al., 2011). In the considered setup the globus pallidus internal (GPi) was stimulated by a short pulse train delayed by 80 ms following an action potential recorded from the primary motor cortex. Under such conditions aDBS was shown to be more effective in reducing pallidal discharge rate and pathological oscillatory neuronal activity as well as in alleviation of akinesia than the conventional continuous HF DBS (cDBS). In a proof of principle study in PD patients (Little et al., 2013), aDBS was switched on and off depending on whether the amplitude of the subthalamic nucleus (STN) local field potential (LFP) in the beta band increased above or decreased below a certain threshold. During aDBS the clinical motor scores strongly improved by about 30% better than during cDBS, whereas aDBS was switched on during 44% of the time (reduced by 56%) as compared to 100% of cDBS (Little et al., 2013). Demand-controlled aDBS was applied for suppression of essential tremor (Graupe et al., 2010). The onset of the tremor was predicted from the measured electromyographic (EMG) signal, which was used to initiate aDBS stimulation epochs. In patients with intention tremor aDBS was switched on and off based on the threshold crossing by EMG power (Yamamoto et al., 2013), where the accurately triggered switching of aDBS resulted in a complete control of the tremor.
Instead of the on-off strategy of the papers (Graupe et al., 2010; Rosin et al., 2011; Little et al., 2013; Yamamoto et al., 2013) mentioned above, the stimulation intensity can also be adapted in real time to the amplitude of the synchronized neuronal activity. Such an approach was used in a clinical study (Rosa et al., 2015), where the voltage of the stimulation was adapted to the beta-band power of the LFP each second (Rosa et al., 2015). The latter approach resembles closed-loop methods that have been developed in the past for the specifically desynchronizing control of abnormal neuronal synchronization that is characteristic for several neurological disorders including PD (Nini et al., 1995; Hammond et al., 2007), essential tremor (Schnitzler et al., 2009), epilepsy (Wong et al., 1986), and tinnitus (Llinas et al., 1999; Weisz et al., 2005; Eggermont and Tass, 2015). These techniques are feedback approaches utilizing the mean field of the synchronized population, which is measured and processed (e.g., filtered, delayed, amplified, etc.), and then fed back as stimulation signal to desynchronize neuronal populations (Rosenblum and Pikovsky, 2004a,b; Hauptmann et al., 2005a,b; Popovych et al., 2005, 2006; Kiss et al., 2007; Pyragas et al., 2007; Tukhlina et al., 2007; Luo et al., 2009; Popovych and Tass, 2010; Montaseri et al., 2013). Direct electrical stimulation of the neuronal tissue with smooth and slowly oscillating feedback signal may however cause an irreversible charge deposit in the neuronal tissue that can exceed safety limits (Harnack et al., 2004; Kuncel and Grill, 2004; Merrill et al., 2005). Two desynchronizing delayed feedback methods, linear delayed feedback (LDF) and nonlinear delayed feedback (NDF) were recently adapted and computationally tested for electrical closed-loop DBS (Popovych et al., 2017a,b). In both cases, the amplitude of charge-balanced short pulses composing the stimulation signal of the standard HF DBS was modulated by the slow feedback signal. The feedback method with such a pulsatile stimulation signal is referred to as a pulsatile feedback stimulation that can be used for electrical DBS.
In principle, abnormal neuronal synchronization can be counteracted in different ways. For instance, LDF (Rosenblum and Pikovsky, 2004a,b) and NDF (Popovych et al., 2005, 2006; Popovych and Tass, 2010) aim at restoring incoherent neuronal activity. In contrast, in this study we consider a multisite linear delayed feedback (MLDF) which has been designed for the control of excessive neuronal synchronization (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008). In previous modeling studies it was observed that MLDF stimulation can counteract the synchronized dynamics by inducing clustering states, which may lead to a variety of spatio-temporal patterns of neuronal activity. Such patterns of the neuronal activity are important, for example, in the context of central pattern generators (CPG) (Marder and Calabrese, 1996; Yuste et al., 2005), where synchronized neuronal discharges have to be well-coordinated both in space and time. The MLDF stimulation approach was suggested for inducing and control of such a patterned activity, for example, in the case when physiological CPG dynamics needs to be restored (Hauptmann et al., 2007a,b; Omel'chenko et al., 2008).
We introduce a pulsatile MLDF for electrical brain stimulation and test it on a physiologically motivated model of interacting neuronal populations of STN and external globus pallidus (GPe) neurons (Terman et al., 2002; Rubin and Terman, 2004). We show that for four-site stimulation setup of MLDF with smooth stimulation signals, a weak clustering, mostly two-cluster states can be observed in a limited parameter range. For pulsatile MLDF stimulation the stimulation-induced clustering becomes even less pronounced such that the main impact of the pulsatile MLDF stimulation consists in a desynchronization, i.e., a suppression of in-phase synchronization in the stimulated population. The pulsatile MLDF stimulation is however less effective in inducing desynchronization than the pulsatile LDF. For effective desynchronization, we here introduce differential pulsatile MLDF stimulation and show this stimulation method can effectively and robustly desynchronize the model STN neurons.
We consider a model network of STN and GPe neuronal populations suggested by Terman et al. (2002), where the dynamics of individual neurons is described by the following system:
In Equations (1–3), v is a membrane potential of the neuron, the currents IL, IK, INa, IT, ICa, IAHP, Isyn, and Iapp are the corresponding leak, potassium, sodium, low threshold calcium, high threshold calcium, afterhyperpolarization potassium, synaptic, and external current, respectively. [Ca] is the intracellular concentration of Ca2+ ions, and X = n, h, r are the gating variables.
The following currents are given by the same expressions for STN and GPe neurons:
On the other hand, current IT is different for the excitatory STN neurons for the inhibitory GPe neurons:
where b∞(r) = 1/(1 + exp[(r − θb)/σb]) − 1/(1 + exp[−θb/σb]). The functions X∞(v) and τX(v) used in the above equations read
For GPe neurons τr(v) = τr is a constant parameter.
The STN and GPe neuronal populations contain N = 200 neurons each and arranged on 1Dim lattices with periodic boundary conditions. Each GPe neuron receives an exitatory input from a single neighboring STN neuron and inhibits three neighboring STN neurons. The considered model was introduced to study pathological neuronal dynamics in PD and was investigated in several papers (Terman et al., 2002; Rubin and Terman, 2004; Park et al., 2011; Popovych et al., 2017a,b). The neurons are interacting via chemical synapses, where the synaptic currents are
Summation is taken over all corresponding presynaptic neurons, where j is the index of neurons. The coupling strength between neurons is defined by parameters of synaptic weights gS→G = 0.4 nS/μm2 (from STN to GPe) and gG→S (from GPe to STN, will be specified below). The reversal potential for the excitatory coupling from STN to GPe vS→G = 0 mV, and vG→S = −100 mV for the inhibitory coupling from GPe to STN. Synaptic variables sj are governed by
The neurons are nonidentical such that the applied currents Iapp = Iapp, j for STN neurons are Gaussian distributed around the mean 10 pA/μm2 and with the standard deviation 0.015 pA/μm2. For GPe neurons the parameter ε = εj are also Gaussian distributed around the mean 0.0055 ms−1 and with the standard deviation 2 · 10−5 ms−1. The other parameters for the STN and GPe neurons and their values are listed in Table S1.
2.2. Synchronized Dynamics of STN Neurons
In this study we investigate how the synchronized dynamics of the STN-GPe network can be controlled by an external stimulation. We estimate the extent of synchronization by the order parameters (Haken, 1983; Kuramoto, 1984; Tass, 1999).
where ψj(t) is the phase of neuron j, which attains the values ψj(tn) = 2πn, n = 0, 1, … at the burst onset time moments tn of the jth neuron. The phase linearly increases between two consecutive bursts ψj(t) = 2π(t − tn)/(tn+1 − tn) + 2πn for t ∈ (tn, tn+1), n = 0, 1, … (Pikovsky et al., 2001). The order parameters Rk(t) range from 0 to 1, where the values of the first order parameter R = R1 correspond to the extent of in-phase synchronization in the population. Large values of the k-th order parameter Rk together with small values of the order parameters Rn of smaller degree n < k are characteristic for a k-cluster state, where the oscillators are in-phase synchronized within the clusters, but the clusters are time (and phase) shifted with respect to each other equidistantly in the oscillation period.
Examples of the collective dynamics of STN neurons without stimulation (Istim = 0 in Equation 1) are illustrated in Figure 1. Depending on the coupling strength as given by the values of parameter gG→S, STN neurons can exhibit synchronization of different extents and forms. For weak coupling, e.g., gG→S = 1.28 nS/μm2, the neurons are weakly and intermittently synchronized, and the order parameter fluctuates around small values 〈R〉 ≈ 0.2, see Figure 1A (black curve). STN neurons exhibit desynchronized bursting dynamics (Figure 1B, black dots), where the individual bursting frequencies (number of bursts per second) of STN neurons are relatively broadly distributed in the range 9.91 ± 0.017 Hz (mean ± standard deviation), see Figure S1. The firing patters exhibit strong variation as time evolves (Figure S1). The inter-bursts intervals (IBI, time intervals between the first spikes of two consecutive bursts) vary irregularly from one burst to the next, which is illustrated in Figure 1E (black circles), where the next IBIn+1 are plotted vs. the previous IBIn as a scatter plot of the first return map.
Figure 1. Collective dynamics of STN-GPe neurons (Equations 1–4) without stimulation. (A) Time courses of the order parameter R = R1 of the STN neurons obtained for the same initial conditions, but for different coupling parameter gG→S as indicated in the legend. (B–D) The corresponding spike raster plots of N = 200 STN neurons, where the spike onsets are indicated by black dots for (B) gG→S = 1.28 nS/μm2, (C) gG→S = 1.38 nS/μm2 and (D) gG→S = 1.7 nS/μm2. The raw LFP (ensemble-averaged synaptic variables sj) and filtered LFP (variable of the linear oscillator, Equation 6) scaled by the factor 1,000 are depicted by blue and red curves, respectively. (E) Scatter plot of the first return map of the inter-burst intervals (IBI, time intervals between the first spikes of two consecutive bursts) of the STN neuron j = 100 for the values of the coupling parameter indicated in the legend.
An increase of the coupling leads to synchronized dynamics of bursting STN neurons as illustrated by the spike raster plot in Figure 1C (black dots) for gG→S = 1.38 nS/μm2. The order parameter fluctuates around a larger value 〈R〉 ≈ 0.8 (Figure 1A, blue curve). The individual bursting frequencies become much narrowly distributed in the range 9.75±0.003 Hz indicating that most neurons become frequency synchronized, see Figure S1. Although the firing patters still vary in time, they clearly demonstrate an in-phase synchronized dynamics of STN neurons (Figure 1C and Figure S1). The behavior of the IBIs remains however irregular (Figure 1E, blue diamonds).
Only stronger coupling can regularize the dynamics of IBIs. For example, for gG→S = 1.7 nS/μm2 the IBIs can attain only two values as illustrated in Figure 1E (red triangles) such that IBIs periodically alternate between them in a period-2 manner. STN neurons synchronize at the same frequency in the very narrow range 9.06 ± 2 · 10−6 Hz, and the firing pattern demonstrates periodic dynamics, see Figure S1. In such a way the STN neurons become periodically synchronized for strong coupling (Figure 1D, black dots), and the order parameter is nearly constant with 〈R〉 ≈ 0.7 (Figure 1A, red curve).
The extent of synchronization is also reflected by the amplitude dynamics of the local field potential (LFP). The latter can be modeled as an ensemble-averaged synaptic activity of the neuronal population (Buzsaki, 2004), where sj(t) are the synaptic variables (Equation 4). For a more sophisticated approach see the papers (Lindén et al., 2011; Parasuram et al., 2016). The measured raw LFP(t) can be on-line filtered by means of a linear damped oscillator
where ω approximates the mean frequency of the LFP oscillations, ω = 2π/T, and T is the LFP mean period. Variable of Equation (6) has a zero phase shift with respect to the input raw LFP signal (Tukhlina et al., 2007), and we consider it as the filtered LFP. The other parameters of Equation (6) were chosen αd = kf = 0.008, which approximately preserves the amplitude of the input LFP signal (Popovych et al., 2017a,b).
The dynamics of raw and filtered LFP is illustrated in the raster plots in Figures 1B–D (blue and red curves) for the above three considered values of the coupling parameter, where large- and small-amplitude oscillations of LFP are in correspondence with strong and weak neuronal synchronization, respectively.
2.3. Delayed Feedback Stimulation
We stimulate the considered STN-GPe model neuronal network by multisite linear delayed feedback (MLDF), where the stimulation is administered to the STN neurons only. This stimulation techniques has been suggested and investigated in the papers (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008). We assume that, for example, 4 stimulation sites are implanted in the STN population of N = 200 neurons at the equidistant lattice coordinates (index of neurons) j = 25, 75, 125, and 175 as schematically illustrated in Figure 2A (upper plot). The feedback stimulation signal Si(t) administered via the i-th stimulation site is calculated as (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008)
The delays τi are considered in such a form in order to achieve a time shift by T/4 between the feedback signals for τ = T. Indeed, if the measured signal x(t) is periodically oscillating with period T, the corresponding oscillating feedback signals (Equation 7) with neighboring indices will be time shifted with respect to each other by T/4 (e.g., τ1 − τ2 = T/4) including S1 with respect to S4, where the latter signal is considered over the next oscillation period, i.e., T + τ4 − τ1 = T/4.
Figure 2. Stimulation setup of MLDF. (A) Four stimulation sites (color dots) equidistantly located within the 1-Dim lattice of STN neurons (upper plot) and example of the corresponding feedback signals Si for delay τ = 102 ms and K = 1 (lower plot). The filtered LFP is depicted by the dashed red curve. (B) Pulsatile stimulation signal, where the amplitude of charge-balanced biphasic pulses composing a high-frequency stimulation train (solid lines) is modulated by the feedback signal (dashed curve). For details see Popovych et al. (2017a,b). The shape of a single pulse is shown in the insert, where the pulse can contain an interphase gap between the negative and positive phases of the pulse.
To calculate the feedback signals (Equation 7), the LFP of synchronized STN neurons is measured, filtered by the linear oscillator (Equation 6), i.e., the signal from Equation (6) is delayed with delays τi from Equation (7), and amplified by the factor K that is a dimensionless feedback gain and will be referred to as parameter of the stimulation intensity.
Instead of four delays, we use only two of them, e.g., τ1 and τ2, and calculate the feedback signals S1 and S2 according to Equation (7). The other two feedback signals are defined by reversing polarity S3 = −S1 and S4 = −S2 (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008). We thus have two stimulation parameters that can be varied: The stimulation intensity (feedback gain) K and the stimulation delay τ. We also assume that neurons within the same sub-population assigned to the corresponding stimulation site receive approximately the same signal administered via that stimulation site, which is called a segmental stimulation, see Hauptmann et al. (2005a,b), Popovych et al. (2006), Hauptmann et al. (2007a,b), and Omel'chenko et al. (2008). More precisely, the neurons j = 1, 2, …, 50 are stimulated with the same signal S1, neurons j = 51, 52, …, 100 are stimulated with the same signal S2, and so on.
The feedback signals Si are illustrated in Figure 2A (lower plot) together with filtered LFP (dashed red curve). For the considered delay τ ≈ T, where T ≈ 102 ms is the mean period of LFP oscillations, the neighboring feedback signals are time shifted by approximately T/4 with respect to each other if the filtered LFP signal is periodic or close to that. Delayed feedback stimulation with such smooth signals may be referred to as smooth feedback stimulation. Direct electrical stimulation of the brain with such signals might violate safety requirements and cause an irreversible charge deposit to the neuronal tissue and could lead to its damage (Harnack et al., 2004; Kuncel and Grill, 2004; Merrill et al., 2005). This problem was studied recently for single-site delayed feedback stimulation (Popovych et al., 2017a,b). By a similar token, we here use a high-frequency stimulation pulse train of the standard HF DBS consisting of biphasic charge-balanced pulses (Volkmann et al., 2002; Kuncel and Grill, 2004; Butson and McIntyre, 2007), whose amplitude is modulated by the slowly oscillating feedback signals Si(t) as schematically illustrated in Figure 2B, where an example of the pulsatile stimulation current Istim in Equation (1) is shown. The cathodic and anodic phases of the pulses administer the same charge of opposite polarity, and a charge-balanced stimulation is realized in this way. The resulting zero net charge injection after each short biphasic pulse can prevent from damaging nervous tissue (Lilly et al., 1955; Harnack et al., 2004; Kuncel and Grill, 2004; Merrill et al., 2005). Each pulse can contain an interphase time gap between its cathodic and anodic phases (Figure 2B, insert). We refer to the stimulation with such pulsatile signal whose amplitude is modulated by the smooth MLDF signals Si(t) as pulsatile MLDF stimulation.
3.1. Smooth MLDF
The impact of the smooth MLDF is illustrated in Figure 3, where synchronized STN neurons are directly stimulated by the smooth feedback signals (Equation 7), i.e., the stimulation currents Istim = Si in Equation (1) for STN neurons. Depending on the parameter of the stimulation delay τ, the stimulation can induce several qualitatively different dynamical regimes as reflected by the values of the order parameters R1, R2, and R4 (Figure 3A). The stimulation can desynchronize the STN neurons, where all order parameters are small as illustrated in Figure 3B for τ = 72 ms. Another stimulation-induced regime is a two-cluster state characterized by small values of R1 and large values of R2, see Figure 3A (red circles and green squares). One of such regimes is illustrated in Figure 3C for τ = 10 ms, where the stimulation splits the stimulated neuronal population into two groups of nearly in-phase synchronized neurons which are shifted by approximately 48 ms with respect to each other. These two clusters are thus nearly in anti-phase to each other (LFP oscillation period T ≈ 103), which results in a small first order parameter R1 and large second order parameter R2 (Figure 3C). Further regimes mimic a four-cluster state, where the fourth order parameter R4 attains greater values as compared to R1 and R2, see Figures 3A,D for τ = 114 ms. Although for such parameters the stimulation clearly splits the stimulated population into four distinct groups phase shifted with respect to each other (with the borders at the lattice coordinates i = 1, 50, 100, and 150, see Figure 3D, right panel), the neurons within these clusters can be far from in-phase synchronization. This leads to relatively small values of R4 and little pronounced four-cluster states.
Figure 3. Impact of smooth MLDF on the STN-GPe neurons (Equations 1–4). (A) Time-averaged order parameters 〈R1〉, 〈R2〉, and 〈R4〉 (as indicated in the legend) vs. parameter of the feedback delay τ for fixed stimulation intensity K = 10. The horizontal dashed line indicates the value of the time-averaged order parameter 〈R1〉 of the STN neurons without stimulation (K = 0). (B–D) Corresponding time courses of the order parameters (left plots) and spike raster plots of STN neurons (right plots) for fixed delay (B) τ = 72 ms, (C) τ = 10 ms, and (D) τ = 114 ms. The stimulation starts at t = 20 s. In plot (B) the time course of the order parameter R1 of the stimulation-free STN population is also shown (black curve). Coupling gG→S = 1.38 nS/μm2.
3.2. Pulsatile MLDF
As mentioned above, the safety requirements for the electrical stimulation of neuronal tissue may be violated for direct electrical stimulation with smooth and slowly oscillating feedback signals. We thus utilize a pulsatile stimulation protocol for MLDF, see section Methods and Figure 2B. Such a pulsatile MLDF stimulation is administered to synchronized STN neurons (Figure 1) for different values of the stimulation parameters τ and K. The time-averaged first order parameter of the obtained stimulation-induced regimes is depicted in Figure 4 in color vs. parameters (τ, K) for width of the interphase gap GW = 0 ms (Figure 4A) and 5 ms (Figure 4B). As compared to the case of smooth MLDF stimulation (Figure 3) the first order parameter R1 exhibits much more pronounced alterations when the delay parameter τ is varied such that several desynchronization regions emerge in the parameter space characterized by small values of R1.
Figure 4. Desynchronization of STN-GPe neurons (Equations 1–4) by pulsatile MLDF. The values of the time-averaged order parameter 〈R1〉 are encoded in color vs. the stimulation intensity K and the stimulation delay τ for the width of the interphase gap (A) GW = 0 ms and (B) 5 ms. Coupling gG→S = 1.38 nS/μm2.
When stimulation pulses include an interphase gap (Figure 2B) the desynchronizing effect of the pulsatile MLDF stimulation can significantly be improved as illustrated in Figure 4B for the interphase gap GW = 5 ms. In fact, for vanishing interphase gap the desynchronizing effect is only moderate (Figure 4A). The structure of the parameter space with desynchronization regions is preserved, but the extent of the stimulation-induced desynchronization is enhanced, as reflected by the values of the first order parameter, i.e., the values of R1 get smaller, compare Figure 4A and Figure 4B. This indicates a favorable effect of the interphase gap on the desynchronization outcome of the pulsatile MLDF.
A detailed consideration by calculating all order parameters R1, R2, and R4 reveals that pulsatile MLDF stimulation with zero gap does not induce any kind of clustering. This is implied by relatively large values of the first order parameter R1 and small values of the other order parameters, see Figure 5A. As mentioned above (Figure 4), an increase of the interphase gap results in a decrease of the first order parameter, which can also be observed in Figures 5A–C (red circles). Simultaneously the other order parameters increase. For example, for the stimulation intensity K = 20, delay τ = 30 ms, and interphase gaps GW = 0 ms, 2 ms, and 5 ms the order parameters R1 ≈ 0.54, 0.4, 0.3, R2 ≈ 0.23, 0.41, 0.5, and R4 ≈ 0.14, 0.16, 0.26, respectively. A relatively large R2 and small R1 are indicative for a two-cluster state. The pattern of the neuronal firing of STN neurons for such parameters is illustrated in Figure 5D, where two groups of neurons (clusters) with different patterns of activity can be distinguished for neuron indices i < 100 and i > 100. Another example of a two-cluster state is illustrated in Figure 5E for the gap width GW = 5 ms and delay τ = 212 ms from another parameter region of a two-cluster regime of large R2 and small R1, see Figure 5C.
Figure 5. Effect of the interphase gap on the dynamics of STN-GPe neurons (Equations 1–4) induced by pulsatile MLDF. (A–C) Time-averaged order parameters 〈R1〉, 〈R2〉, and 〈R4〉 (as indicated in the legends) vs. the stimulation time delay τ for fixed stimulation intensity K = 20 and width of the interphase gap (A) GW = 0 ms, (B) 2 ms, and (C) 5 ms. The first order parameter 〈R1〉 for K = 10 is also shown for comparison (black solid curves). (D,E) Raster plots of the stimulation-induced cluster dynamics for K = 20, GW = 5 ms, and (D) τ = 30 ms and (E) τ = 212 ms. Coupling gG→S = 1.38 nS/μm2.
We compare the desynchronizing impact of the pulsatile MLDF to that of pulsatile linear delayed feedback (LDF). The smooth and pulsatile LDF administered to synchronized STN neurons has been investigated in Popovych et al. (2017a,b) together with smooth and pulsatile nonlinear delayed feedback (NDF). The feedback signal S(t) of the differential LDF can be obtained as (Rosenblum and Pikovsky, 2004a,b):
where, as before, the variable x(t) is a filtered LFP and calculated by means of Equation (6), i.e., , and K and τ are the parameters of the stimulation intensity and delay, respectively. The smooth feedback signal S(t) of LDF is then used to modulate the amplitude of the stimulation pulses as discussed above, where we assume that all STN neurons receive the same stimulation current Istim depicted in Figure 2B.
STN neurons can be desynchronized by the pulsatile LDF stimulation as illustrated in Figure 6. The parameter space of the pulsatile LDF contains large desynchronization regions, where the first order parameter R1 exhibits pronounced minima (Figures 6A,B, red circles). An interphase gap of a finite width can enhance the desynchronizing effect of pulsatile LDF, where the desynchronization regions are enlarged and deepened such that the stimulation induces stronger desynchronization as reflected by smaller values of R1 for larger interphase gap, compare Figure 6A and Figure 6B (red circles). The other order parameters of higher degree are however not significantly affected by the values of the interphase gap. For example, the values of the second order parameter R2 are nearly preserved within the desynchronization regions when the width of the interphase gap increases from GW = 0 to 5 ms, compare Figure 6A and Figure 6B (green squares). Albeit the suppressed second order parameter R2 can still be somewhat larger than the first order parameter R1 (Figure 6B), the stimulated neurons do not exhibit any consistent clustering as illustrated in Figure 6C. This indicates that the pulsatile LDF stimulation does not induce clustering among stimulated STN neurons for any widths of the interphase gap, which is different to the impact of the pulsatile MLDF, see Figure 5.
Figure 6. Desynchronization of STN-GPe neurons (Equations 1–4) by pulsatile LDF. (A,B) Time-averaged order parameters 〈R1〉 and 〈R2〉 (as indicated in the legend) of the STN neurons stimulated by pulsatile LDF vs. the stimulation delay τ for fixed stimulation intensity K = 20 and interphase gap (A) GW = 0 ms and (B) 5 ms. The first order parameter 〈R1〉 for pulsatile MLDF stimulation is also shown for comparison (black solid curves, see Figures 5A,C, red circles). (C) Example of the time courses of the order parameters R1 and R2 (left plot) and the corresponding spike raster plot of STN neurons (right plot) for pulsatile LDF stimulation for fixed K = 20, τ = 60 ms, and GW = 5 ms (from plot B). Coupling gG→S = 1.38 nS/μm2.
The pulsatile LDF is more efficient in inducing desynchronization than the pulsatile MLDF. For the same stimulation intensity K and width of the interphase gap, the pulsatile LDF can induce much stronger desynchronization than the pulsatile MLDF as given by the values of the first order parameter R1 within desynchronization regions, see Figures 6A,B and compare black solid curves (R1 for MLDF) to red circles (R1 for LDF). For a more detailed comparison, we fix optimal stimulation delay τ = 90 ms for pulsatile MLDF and τ = 60 ms for pulsatile LDF, where the stimulation induces strongest desynchronization, see Figures 4–6, and increase the stimulation intensity K. We find that both pulsatile MLDF and LDF stimulations with larger intensity can induce stronger desynchronization, and the first order parameter R1 decreases when K increases, see Figure 7A for MLDF and Figure 7B for LDF (empty symbols). We also observe that R1 decreases much faster with increasing K and can reach much smaller values for pulsatile LDF than for MLDF for the same range on the stimulation intensity.
Figure 7. Desynchronizing outcome of pulsatile MLDF and LDF stimulations administered to STN-GPe neurons (Equations 1–4). (A,B) Time-averaged first order parameter 〈R1〉 and the absolute value 〈|S|〉 of the feedback signal (A) Equation (7) of MLDF and (B) Equation (8) of LDF vs. parameter K of the stimulation intensity for two widths GW = 0 ms and 5 ms of the interphase gap as indicated in the legends. The scale of 〈|S|〉 is indicated on the right vertical axis. (C,D) Administered amount of the stimulation as given by the values of 〈|S|〉 vs. the reached extent of the stimulation-induced desynchronization as given by values of 〈R1〉 for (C) pulsatile MLDF and (D) pulsatile LDF for different width of the interphase gap as indicated in the legends. Stimulation delay (A,C) τ = 90 ms for MLDF and (B,D) τ = 60 ms for LDF. Coupling gG→S = 1.38 nS/μm2.
An important characteristics of the stimulation of the neuronal tissue is the amount of the administered stimulation. We thus estimate it for the considered feedback stimulations by the time-averaged absolute value 〈|S|〉 of the smooth feedback signals (Equation 7) for MLDF and (Equation 8) and LDF. The amount of the administered stimulation is depicted in Figure 7A for MLDF and Figure 7B for LDF (filled symbols) vs. parameter K of the stimulation intensity. When K increases, the amount of the administered stimulation also increases such that stronger desynchronization can be obtained at stronger stimulation. Since the amplitude of the feedback signals depends on the amplitude of the LFP, it also inversely relates to the extent of the stimulation-induced desynchronization. For a given value of the stimulation intensity K, the amount of the administered stimulation 〈|S|〉 will be smaller if the stimulation-induced desynchronization is stronger, i.e., when the order parameter R1 and amplitude of the LFP are smaller. Therefore, for the considered range of parameter K, 〈|S|〉 increases more slowly for pulsatile LDF than for pulsatile MLDF with increasing stimulation intensity, see Figures 7A,B (filled symbols). For both stimulation methods larger interphase gap leads to better desynchronization and smaller amount of the administered stimulation.
The latter claim is also supported by Figures 7C,D, where the amount of the administered stimulation 〈|S|〉 is depicted vs. the extent of the stimulation-induced desynchronization as reflected by the first order parameter R1. As follows, the same extent of desynchronization can be obtained at smaller amount of the administered stimulation for larger interphase gap. In such a way we can compare the efficacy in inducing desynchronization of MLDF and LDF methods by comparing the amount of the administered stimulation necessary to achieve the same extent of the stimulation-induced desynchronization. Comparing the depicted data for pulsatile MLDF (Figure 7C) and pulsatile LDF (Figure 7D) we conclude that the pulsatile LDF is more effective in inducing desynchronization than pulsatile MLDF, where desynchronization can be obtained for much smaller amount of the administered stimulation.
The above results obtained for the synchronized regime of STN-GPe neurons with irregular interburst intervals (Figure 1 for coupling gG→S = 1.7 nS/μm2) are preserved for other parameters and synchronized regimes in the considered populations of the STN-GPe neurons. We consider, for instance, stronger coupling gG→S = 1.7 nS/μm2, where the STN-GPe neurons are periodically synchronized, see Figure 1A (red curve) and Figure 1D. The desynchronizing effect of the pulsatile MLDF and LDF on the periodically synchronized STN-GPe neurons is illustrated in Figure 8A. The structure of the parameter space is preserved for both stimulation methods except for that the pulsatile MLDF induces a somewhat weaker desynchronization as compared to the above case of a weaker coupling. The clustering induced by the pulsatile MLDF for the considered stronger coupling also becomes less pronounced. The efficacy of the pulsatile MLDF and LDF in inducing desynchronization is compared in Figure 8B, where the amount of the administered stimulation 〈|S|〉 is depicted vs. the time-averaged first order parameter 〈R1〉. As for the case of the irregular synchronization for weaker coupling (Figure 7), the interphase gap has the same favorable impact on the stimulation outcome, and the pulsatile LDF is apparently superior to MLDF and can induce stronger desynchronization for smaller amount of administered stimulation.
Figure 8. Desynchronization of periodically synchronized STN-GPe neurons (Equations 1–4) by pulsatile MLDF and LDF stimulations. (A) Time-averaged order parameter 〈R1〉 of the STN neurons stimulated by pulsatile MLDF with K = 30 and pulsatile LDF with K = 10 vs. the stimulation delay τ for two widths of the interphase gap GW = 0 ms and 5 ms as indicated in the legend. (B) Administered amount of the stimulation as given by the values of 〈|S|〉 vs. the reached extent of the stimulation-induced desynchronization as given by values of 〈R1〉 for pulsatile MLDF (empty symbols) with τ = 50 ms and LDF (filled symbols) with τ = 70 ms. The width of the interphase gap is indicated in the legend. Coupling gG→S = 1.7 nS/μm2.
3.3. Differential MLDF
To overcome the limitations of pulsatile MLDF revealed above, we suggest to use a differential MLDF. The feedback signals are constructed by analogy with the differential LDF (Equation 8) and read.
where the signal x(t) is the filtered LFP from Equation (6), and the delays are as in Equation (7). As before, we reverse the polarity of the two feedback signals, such that S3 = K · (−x(t − τ1) − x(t)) and S4 = K · (−x(t − τ2) − x(t)). An example of the feedback signals Si is illustrated in Figure 9A. As for the case of direct MLDF (Equation 7), the feedback signals of the differential MLDF are time shifted with respect to each other, but, in contrast to the direct MLDF, they may however have very different amplitude.
Figure 9. Impact of the differential smooth MLDF on collective dynamics of STN-GPe neurons (Equations 1–4). (A) Example of the feedback signals Si of the differential MLDF (Equation 9) for delay τ = 90 ms. The filtered LFP is depicted by dashed red curve. (B) Time-averaged order parameters 〈R1〉, 〈R2〉, and 〈R4〉 (as indicated in the legend) vs. parameter of the feedback delay τ for fixed stimulation intensity K = 10. The first order parameter 〈R1〉 of the STN neurons for direct smooth MLDF stimulation (Equation 7) is also shown by the black solid curve for comparison (copied from Figure 3A, red circles). The horizontal dashed line indicates the value of 〈R1〉 without stimulation (K = 0). (C,D) Spike raster plots of STN neurons for fixed stimulation intensity K = 10 and delays (C) τ = 10 ms and (D) τ = 90 ms. Coupling gG→S = 1.38 nS/μm2.
Stimulation by differential smooth MLDF can perturb the neuronal synchronization of the stimulated STN neurons as illustrated in Figures 9B–D, where the time-averaged order parameters 〈R1〉, 〈R2〉, and 〈R4〉 are plotted vs. the stimulation delay τ. Based on the values of the first order parameter R1, smooth differential MLDF is less effective in inducing desynchronization than smooth direct MLDF, compare values of R1 in Figure 9B of differential MLDF (red circles) to those of direct MLDF (black curve). For differential MLDF large values of the other order parameters are accompanied by large values of R1. This indicates that clusters of a possible clustered state are not equidistantly spaced over the oscillation period. Indeed, this is illustrated in Figure 9C for delay τ = 10 ms, where the fourth order parameter R4 attains maximal values, see Figure 9B. In this example the stimulation divides the neurons into two groups, where the burst onsets of the two clusters are time shifted with respect to each other by ΔT ≈ 26 ms (Figure 9C). Such firing patterns result in large values of the order parameters R1 and R4, but in small values of R2 as depicted in Figure 9B for the mentioned value of the stimulation delay. The best desynchronization can be achieved, for example, at τ = 90 ms as reflected by small values of R1 (Figure 9B). The corresponding firing pattern of STN neurons is illustrated in Figure 9D.
In contrast to the stimulation with smooth signals, the differential pulsatile MLDF can induce stronger desynchronization as compared to the direct pulsatile MLDF. Two-parameter diagrams for differential pulsatile MLDF, where the order parameters are depicted in color in the (τ, K)-parameter plane, are shown in Figure 10. For the same range of the stimulation parameters the first order parameter R1 (Figures 10A,B) exhibits smaller values as compared to those depicted in Figure 4 for direct pulsatile MLDF. Interphase gap enhances the extent of the stimulation-induced desynchronization, compare Figure 10A and Figure 10B. For large interphase gap the (τ, K)-parameter plane contains a narrow region around τ ≈ 15 ms, where the second order parameter R2 is relatively large (Figure 10C). For such parameters of the best two-cluster regime, the pattern of the neuronal firing looks similar to that in Figure 5D, i.e., there are no pronounced two- and four-cluster states.
Figure 10. Desynchronization of STN-GPe neurons (Equations 1–4) by the differential pulsatile MLDF (Equation 9). The values of the time-averaged order parameters (A,B) 〈R1〉, (C) 〈R2〉, and (D) 〈R4〉 are depicted in color vs. the stimulation intensity K and delay τ for the width of the interphase gap (A) GW = 0 ms and (B–D) 5 ms. Coupling gG→S = 1.38 nS/μm2.
To evaluate the stimulation outcome of the differential pulsatile MLDF we calculate the amount of the administered stimulation 〈|S|〉 that is approximated by the average of the absolute values of the MLDF feedback signals (Equation 9). In Figure 11A 〈|S|〉 is depicted together with the extent of the stimulation-induced desynchronization as given by the values of the first order parameter 〈R1〉 vs. parameter of the stimulation intensity K. For the differential pulsatile MLDF, the order parameter 〈R1〉 decreases much faster than for the direct pulsatile MLDF, and, as a results, 〈|S|〉 increases more slowly, compare Figure 11A to Figure 7A. The efficacy of the differential pulsatile MLDF in inducing desynchronization is comparable with that of the pulsatile LDF, see Figure 7B. Indeed, this conclusion apparently follows from Figure 11B, where a given extent of the stimulation-induced desynchronization as given by values of R1 can be obtained at approximately the same amount of the administered stimulation for differential pulsatile MLDF (Figure 11B, empty symbols) and for pulsatile LDF (Figure 11B, filled symbols).
Figure 11. Stimulation outcome of the differential pulsatile MLDF stimulation administered to STN-GPe neurons (Equations 1–4). (A) Time-averaged first order parameter 〈R1〉 and the averaged absolute value 〈|S|〉 of the feedback signals (Equation 9) of differential MLDF vs. parameter K of the stimulation intensity for two widths GW = 0 ms and 5 ms of the interphase gap as indicated in the legend. The scale of 〈|S|〉 is indicated on the right vertical axis. (B) Administered amount of the stimulation as given by the values of 〈|S|〉 vs. the reached extent of the stimulation-induced desynchronization as given by values of 〈R1〉 for differential pulsatile MLDF (empty symbols) for delay τ = 90 ms and pulsatile LDF (filled symbols, copied from Figure 7D for comparison) for delay τ = 60 ms for different width of the interphase gap as indicated in the legend. Coupling gG→S = 1.38 nS/μm2.
We also verified that differential pulsatile MLDF stimulation is robust with respect to the extent of the initial synchronization in the stimulated neuronal population. If the stimulation is administered to weakly synchronized neurons as, for example, for the coupling parameter gG→S = 1.28 nS/μm2, see Figure 1, the synchronization can further be suppressed by differential pulsatile MLDF practically irrespective of the values of the delay parameter τ, see Figure 12. Such a reduction of an already weak and intermittent neuronal synchronization is comparable with or even slightly better than that of pulsatile LDF, see also Popovych et al. (2017a,b).
Figure 12. Impact of the differential pulsatile MLDF on collective dynamics of weakly coupled and weakly synchronized STN-GPe neurons (Equations 1–4). Time-averaged order parameters 〈Ri〉, i = 1, 2, 4, are plotted (as indicated in the legend) vs. parameter of the feedback delay τ for fixed stimulation intensity K = 10. The first order parameter 〈R1〉 of the STN neurons stimulated by pulsatile LDF is also shown by the black solid curve for comparison. The horizontal dashed line indicates the value of 〈R1〉 without stimulation (K = 0), see Figure 1A (black curve). Coupling gG→S = 1.28 nS/μm2 and interphase gap GW = 5 ms.
Multisite linear delayed feedback (MLDF) has been suggested for control of neuronal synchronization patterns (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008). As shown computationally, stimulation by MLDF can suppress synchronization in a stimulated neuronal population and was, hence, suggested for counteracting abnormal neuronal synchronization characteristic for several neurological disorders (Hauptmann et al., 2005a,b, 2007b; Popovych et al., 2006). The desynchronization induced by MLDF stimulation was found to be accompanied by the emergence of several interacting clusters of neurons equidistantly distributed over the oscillation period and space. Depending on the stimulation setup and parameters the stimulation-induced spatio-temporal patters can consist of two or four clusters, for example, for four-site MLDF stimulation (Hauptmann et al., 2005a,b, 2007a,b; Popovych et al., 2006; Omel'chenko et al., 2008). These properties of MLDF made the method appropriate for the control of spatio-temporal patterns of neuronal activity, for example, for regulating activity of central pattern generator (CPG) in case of its malfunction as suggested in Hauptmann et al. (2007a,b) and Omel'chenko et al. (2008).
In this study we adapted the MLDF technique for electrical stimulation of the neuronal tissue. As mentioned above, direct electrical stimulation with smooth feedback signals may violate safety aspects like charge density limits (Harnack et al., 2004; Kuncel and Grill, 2004; Merrill et al., 2005). The feedback signals are slow, so that an irreversible charge can be deposited into the neuronal tissue during the comparably long feedback stimulation periods, which can exceed safety limits. We resolved this problem and showed that the demand-controlled character and desynchronizing impact of the MLDF feedback technique can be preserved together with gaining the advantages of the pulsatile HF DBS signal with the charge-balanced property. For this, the slow feedback signal is used to modulate the amplitude of the HF train of charge-balanced pulses of HF DBS. We computationally illustrated the desynchronizing properties of smooth and pulsatile MLDF in a network of STN-GPe neurons suggested to model parkinsonian neuronal dynamics (Terman et al., 2002; Rubin and Terman, 2004).
We showed that both smooth and pulsatile MLDF stimulation can suppress neuronal synchronization in the stimulated STN neurons. While smooth MLDF could induce relatively well pronounced two-cluster states in some parameter ranges, the expected four-cluster states were only weakly expressed and could be observed in a limited parameter range. For pulsatile MLDF we found that interphase gap in the stimulation pulses could significantly enhance the desynchronizing impact of the stimulation. The clustering state was observed for large interphase gap only, where some two-cluster states could be induced by the stimulation for some selected parameter values. We thus conclude that the pulsatile MLDF is mostly a desynchronizing stimulation rather than inducing coordinated spatio-temporal clustering patterns. We however showed that the efficacy of the pulsatile MLDF in inducing desynchronization was much lower than that of the pulsatile LDF. Therefore, in the standard realization the pulsatile MLDF can be suggested as effective method neither for desynchronization nor for the control of spatio-temporal clustering patterns.
We therefore proposed to use a differential pulsatile MLDF. Such a modified pulsatile MLDF turned out to hardly induce clustering for any width of the interphase gap. Nevertheless, we showed that differential pulsatile MLDF can effectively and robustly desynchronize the stimulated neurons. We verified that the efficacy of the differential pulsatile MLDF in inducing desynchronization is comparable with that of the pulsatile LDF, and we suggest this technique for closed-loop desynchronizing DBS together with pulsatile LDF and NDF investigated in recent papers (Popovych et al., 2017a,b). The differential pulsatile MLDF is robust with respect to variations of the stimulation parameters, in particular, if initially weakly synchronized neuronal populations need to be further desynchronized. However, MLDF requires several stimulation sites to be placed in the neuronal target population. In the case of small targets this might be difficult, so that single-site stimulation techniques such as pulsatile LDF and NDF may be more appropriate in those cases. It would be interesting to investigate how the effectiveness and efficacy of pulsatile NDF stimulation are affected if a multisite stimulation protocol is adapted for this method. It would also be interesting to use a realistic 3-Dim reconstruction of the target neuronal structures, e.g., STN and GPe, as well as localization of the stimulation sites within these structures to explore the spatio-temporal patterns induced by a multisite stimulation (Ebert et al., 2014), which however essentially requires the usage of a supercomputer.
However, as a word of caution, it should be noted that the approach presented in this study relies on the assumption that abnormal neuronal synchrony is recordable and represents the patient's individual symptoms in a sufficient manner, like a biomarker (Beudel and Brown, 2016; Kühn and Volkmann, 2017). For instance, it is doubtful that beta band oscillations might be a biomarker-like feedback signal (Özkurt et al., 2011; Johnson et al., 2016; Kühn and Volkmann, 2017). Beta band oscillations are no stand-alone oscillations, but interact with brain oscillations in other frequency bands under both physiological (Yanagisawa et al., 2012) and pathological (Yang et al., 2014; Beudel and Brown, 2016) conditions. Different PD phenotypes might require different biomarkers, since the amplitude of beta band oscillations may decrease during tremor epochs in tremor dominant PD patients (Quinn et al., 2015). Changes of the amplitude of abnormal brain oscillations in the course of physiological processes have to be taken into account, too. For instance, in an MPTP monkey study HF DBS and closed-loop DBS (CL-DBS) reduced rigidity to a comparable extent, where CL-DBS reduced the DBS ON time by approx. 50% (Johnson et al., 2016). However, only HF DBS improved bradykinesia during a cued reaching task, likely because the amplitude of beta band oscillations was reduced related to the reaching process, in this way reducing the extent of the presumed biomarker (Johnson et al., 2016). Also, beta band oscillations need not be entirely pathological. Rather activity in the beta frequency range might be key for compensatory purposes, as demonstrated in an MPTP monkey study with sensorimotor rhythm neurofeedback (Philippens et al., 2017).
Apart from merely reducing the stimulation current, differential pulsatile MLDF may be beneficial because of its specifically desynchronizing effect. As shown computationally in the context of Coordinated Reset (CR) stimulation (Tass, 2003), desynchronizing stimulation may cause an anti-kindling, where abnormal synaptic connectivity and neuronal synchrony can be unlearned, ultimately leading to sustained desynchronizing effects (Tass and Majtanik, 2006; Popovych and Tass, 2012). In accordance with these theoretical predictions, long-lasting therapeutic effects were demonstrated in pre-clinical studies in MPTP monkeys with CR-DBS (Tass et al., 2012b; Wang et al., 2016) as well as in a clinical proof of concept study with CR-DBS in Parkinson's patients (Adamchic et al., 2014). Analogously, long-lasting therapeutic effects were observed in a proof of concept study with acoustic CR stimulation in tinnitus patients (Tass et al., 2012a) as well as in a first in man study with vibrotactile CR stimulation in patients with Parkinson's disease (Syrkin-Nikolau et al., 2018).
The requirements for CR-DBS and MLDF are quite different. CR stimulation can be administered in open loop as well as closed loop, e.g., demand-controlled manner (Tass, 2003). In particular, CR stimulation does not require a feedback signal. So far, pre-clinical (Tass et al., 2012b; Wang et al., 2016) and clinical (Adamchic et al., 2014) proof of concept of CR-DBS were obtained with open loop CR-DBS. In contrast, MLDF requires a reliably measurable clean biomarker signal sufficiently representing the amount of abnormal synchronization. Despite first positive results (Little et al., 2013; Rosa et al., 2015), several findings indicate that beta-band STN LFP does not provide a reliable biomarker (Özkurt et al., 2011; Quinn et al., 2015; Johnson et al., 2016; Kühn and Volkmann, 2017; Philippens et al., 2017), see above. In addition, MLDF requires a challenging registration-stimulation setup: An LFP signal, representative for the entire neuronal target population, has to be measured at one site, while stimuli have to be delivered to different sites of the target population. Because of stimulation artifacts this might be difficult. However, to overcome this limitation, alternatively, one might try to measure a representative LFP and stimulate different parts of fibers projecting on the target population. One could also separate stimulation and recording in time (Ratas and Pyragas, 2014).
PT conceived HFS amplitude modulation by feedback. OP performed the experiments, analyzed the data, and prepared the initial draft of the manuscript. PT contributed to the numerical analysis and extended the manuscript. All authors reviewed the manuscript.
The study was funded by the Helmholtz Society and by the John A. Blume Foundation.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00046/full#supplementary-material
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
Benabid, A. L., Chabardes, S., Mitrofanis, J., and Pollak, P. (2009). Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson's diseaseegge. Lancet Neurol. 8, 67–81. doi: 10.1016/S1474-4422(08)70291-6
Benabid, A. L., Pollak, P., Gervason, C., Hoffmann, D., Gao, D. M., Hommel, M., et al. (1991). Longterm suppression of tremor by chronic stimulation of ventral intermediate thalamic nucleus. Lancet 337, 403–406.
Butson, C. R., and McIntyre, C. C. (2007). Differences among implanted pulse generator waveforms cause variations in the neural response to deep brain stimulation. Clin. Neurophysiol. 118, 1889–1894. doi: 10.1016/j.clinph.2007.05.061
Deniau, J.-M., Degos, B., Bosch, C., and Maurice, N. (2010). Deep brain stimulation mechanisms: beyond the concept of local functional inhibition. Eur. J. Neurosci. 32, 1080–1091. doi: 10.1111/j.1460-9568.2010.07413.x
Ferraye, M. U., Debû, B., Fraix, V., Xie-Brustolin, J., Chabardès, S., Krack, P., et al. (2008). Effects of subthalamic nucleus stimulation and levodopa on freezing of gait in parkinson disease. Neurology 70, 1431–1437. doi: 10.1212/01.wnl.0000310416.90757.85
Grahn, P. J., Mallory, G. W., Khurram, O. U., Berry, B. M., Hachmann, J. T., Bieber, A. J., et al. (2014). A neurochemical closed-loop controller for deep brain stimulation: toward individualized smart neuromodulation therapies. Front. Neurosci. 8:169. doi: 10.3389/fnins.2014.00169
Graupe, D., Basu, I., Tuninetti, D., Vannemreddy, P., and Slavin, K. V. (2010). Adaptively controlling deep brain stimulation in essential tremor patient via surface electromyography. Neurol. Res. 32, 899–904. doi: 10.1179/016164110X12767786356354
Harnack, D., Winter, C., Meissner, W., Reum, T., Kupsch, A., and Morgenstern, R. (2004). The effects of electrode material, charge density and stimulation duration on the safety of high-frequency stimulation of the subthalamic nucleus in rats. J. Neurosci. Methods 138, 207–216. doi: 10.1016/j.jneumeth.2004.04.019
Hauptmann, C., Omelchenko, O., Popovych, O. V., Maistrenko, Y., and Tass, P. A. (2007a). Control of spatially patterned synchrony with multisite delayed feedback. Phys. Rev. E 76:066209. doi: 10.1103/PhysRevE.76.066209
Hauptmann, C., Popovych, O., and Tass, P. A. (2005a). Delayed feedback control of synchronization in locally coupled neuronal networks. Neurocomputing 65–66, 759–767. doi: 10.1016/j.neucom.2004.10.072
Hauptmann, C., Popovych, O., and Tass, P. A. (2005b). 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
Hauptmann, C., Popovych, O., and Tass, P. A. (2007b). Demand-controlled desynchronization of oscillatory networks by means of a multisite delayed feedback stimulation. Comput. Visual. Sci. 10, 71–78. doi: 10.1007/s00791-006-0034-9
Johnson, L. A., Nebeck, S. D., Muralidharan, A., Johnson, M. D., Baker, K. B., and Vitek, J. L. (2016). Closed-loop deep brain stimulation effects on Parkinsonian motor symptoms in a non-human primate - is beta enough? Brain Stimul. 9, 892–896. doi: 10.1016/j.brs.2016.06.051
Johnson, M. D., Miocinovic, S., McIntyre, C. C., and Vitek, J. L. (2008). Mechanisms and targets of deep brain stimulation in movement disorders. Neurotherapeutics 5, 294–308. doi: 10.1016/j.nurt.2008.01.010
Kiss, I. Z., Rusin, C. G., Kori, H., and Hudson, J. L. (2007). Engineering complex dynamical structures: sequential patterns and desynchronization. Science 316, 1886–1889. doi: 10.1126/science.1140858
Llinas, R. R., Ribary, U., Jeanmonod, D., Kronberg, E., and Mitra, P. P. (1999). Thalamocortical dysrhythmia: a neurological and neuropsychiatric syndrome characterized by magnetoencephalography. Proc. Natl. Acad. Sci. U.S.A. 96, 15222–15227.
Luo, M., Wu, Y. J., and Peng, J. H. (2009). Washout filter aided mean field feedback desynchronization in an ensemble of globally coupled neural oscillators. Biol. Cybern. 101, 241–246. doi: 10.1007/s00422-009-0334-5
Merrill, D. R., Bikson, M., and Jefferys, J. G. R. (2005). Electrical stimulation of excitable tissue: design of efficacious and safe protocols. J. Neurosci. Methods 141, 171–198. doi: 10.1016/j.jneumeth.2004.10.020
Montaseri, G., Yazdanpanah, M. J., Pikovsky, A., and Rosenblum, M. (2013). Synchrony suppression in ensembles of coupled oscillators via adaptive vanishing feedback. Chaos 23:033122. doi: 10.1063/1.4817393
Moreau, C., Defebvre, L., Destée, A., Bleuse, S., Clement, F., Blatt, J. L., 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
Nini, A., Feingold, A., Slovin, H., and Bergmann, 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.
Omel'chenko, O. E., Hauptmann, C., Maistrenko, Y. L., and Tass, P. A. (2008). Collective dynamics of globally coupled phase oscillators under multisite delayed feedback stimulation. Physica D 237, 365–384. doi: 10.1016/j.physd.2007.09.019
Özkurt, T. E., Butz, M., Homburger, M., Elben, S., Vesper, J., Wojtecki, L., et al. (2011). High frequency oscillations in the subthalamic nucleus: a neurophysiological marker of the motor state in Parkinson's disease. Exp. Neurol. 229, 324–331. doi: 10.1016/j.expneurol.2011.02.015
Parasuram, H., Nair, B., D'Angelo, E., Hines, M., Naldi, G., and Diwakar, S. (2016). Computational modeling of single neuron extracellular electric potentials and network local field potentials using lfpsim. Front. Comput. Neurosci. 10:65. doi: 10.3389/fncom.2016.00065
Park, C., Worth, R. M., and Rubchinsky, L. L. (2011). Neural dynamics in parkinsonian brain: the boundary between synchronized and nonsynchronized dynamics. Phys. Rev. E 83:042901. doi: 10.1103/PhysRevE.83.042901
Philippens, I., Wubben, J., Vanwersch, R., Estevao, D., and Tass, P. (2017). Sensorimotor rhythm neurofeedback as adjunct therapy for Parkinson's disease. Ann. Clin. Transl. Neurol. 4, 585–590. doi: 10.1002/acn3.434
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
Priori, A., Foffani, G., Rossi, L., and Marceglia, S. (2013). Adaptive deep brain stimulation (aDBS) controlled by local field potential oscillations. Exp. Neurol. 245, 77–86. doi: 10.1016/j.expneurol.2012.09.013
Pyragas, K., Popovych, O. V., and Tass, P. A. (2007). Controlling synchrony in oscillatory networks with a separate stimulation-registration setup. Europhys. Lett. 80:40002. doi: 10.1209/0295-5075/80/40002
Quinn, E. J., Blumenfeld, Z., Velisar, A., Koop, M. M., Shreve, L. A., Trager, M. H., et al. (2015). Beta oscillations in freely moving Parkinson's subjects are attenuated during deep brain stimulation. Mov. Disord. 30, 1750–1758. doi: 10.1002/mds.26376
Rosa, M., Arlotti, M., Ardolino, G., Cogiamanian, F., Marceglia, S., Di Fonzo, A., et al. (2015). Adaptive deep brain stimulation in a freely moving parkinsonian patient. Mov. Disord. 30, 1003–1005. doi: 10.1002/mds.26241
Rosenblum, M. G., and Pikovsky, A. S. (2004b). 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
Rosin, B., Slovik, M., Mitelman, R., Rivlin-Etzion, M., Haber, S. N., Israel, Z., et al. (2011). Closed-loop deep brain stimulation is superior in ameliorating parkinsonism. Neuron 72, 370–384. doi: 10.1016/j.neuron.2011.08.023
Rubin, J. E., and Terman, D. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. J. Comput. Neurosci. 16, 211–235. doi: 10.1023/B:JCNS.0000025686.47117.67
Schnitzler, A., Munks, C., Butz, M., Timmermann, L., and Gross, J. (2009). Synchronized brain network associated with essential tremor as revealed by magnetoencephalography. Mov. Disord. 24, 1629–1635. doi: 10.1002/mds.22633
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–180. doi: 10.1002/mds.27223
Tass, P. A., Adamchic, I., Freund, H.-J., von Stackelberg, T., and Hauptmann, C. (2012a). Counteracting tinnitus by acoustic coordinated reset neuromodulation. Rest. Neurol. Neurosci. 30, 137–159. doi: 10.3233/RNN-2012-110218
Tass, P. A., Qin, L., Hauptmann, C., Doveros, 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
van Nuenen, B. F. L., Esselink, R. A. J., Munneke, M., Speelman, J. D., van Laar, T., and Bloem, B. R. (2008). Postoperative gait deterioration after bilateral subthalamic nucleus stimulation in parkinson's disease. Mov. Disord. 23, 2404–2406. doi: 10.1002/mds.21986
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
Weisz, N., Moratti, S., Meinzer, M., Dohrmann, K., and Elbert, T. (2005). Tinnitus perception and distress is related to abnormal spontaneous brain activity as measured by magnetoencephalography. PLoS Med. 2:e153. doi: 10.1371/journal.pmed.0020153
Xie, T., Kang, U. J., and Warnke, P. (2012). Effect of stimulation frequency on immediate freezing of gait in newly activated stn dbs in parkinson's disease. J. Neurol. Neurosurg. Psychiatry 83, 1015–1017. doi: 10.1136/jnnp-2011-302091
Yamamoto, T., Katayama, Y., Ushiba, J., Yoshino, H., Obuchi, T., Kobayashi, K., et al. (2013). On-demand control system for deep brain stimulation for treatment of intention tremor. Neuromodulation 16, 230–235. doi: 10.1111/j.1525-1403.2012.00521.x
Yanagisawa, T., Yamashita, O., Hirata, M., Kishima, H., Saitoh, Y., Goto, T., et al. (2012). Regulation of motor representation by phase-amplitude coupling in the sensorimotor cortex. J. Neurosci. 32, 15467–15475. doi: 10.1523/JNEUROSCI.2929-12.2012
Yang, A. I., Vanegas, N., Lungu, C., and Zaghloul, K. A. (2014). Beta-coupled high-frequency activity and beta-locked neuronal spiking in the subthalamic nucleus of parkinson's disease. J. Neurosci. 34, 12816–12827. doi: 10.1523/JNEUROSCI.1895-14.2014
Keywords: neuronal synchronization, delayed feedback, deep brain stimulation, desynchronization, electrical pulse stimulation, closed-loop stimulation
Citation: Popovych OV and Tass PA (2018) Multisite Delayed Feedback for Electrical Brain Stimulation. Front. Physiol. 9:46. doi: 10.3389/fphys.2018.00046
Received: 19 October 2017; Accepted: 15 January 2018;
Published: 01 February 2018.
Edited by:Alexey Zaikin, University College London, United Kingdom
Reviewed by:Anuj Agarwal, Signal Solutions LLC, United States
Mikhail Ivanchenko, N. I. Lobachevsky State University of Nizhny Novgorod, Russia
Copyright © 2018 Popovych 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 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: Oleksandr V. Popovych, email@example.com