article: Deep brain stimulation (DBS) is known to be an effective treatment for a variety of neurological disorders, including Parkinson’s disease and essential tremor (ET). At present, it involves administering a train of pulses with constant frequency via electrodes implanted into the brain. New ‘closed-loop’ approaches involve delivering stimulation according to the ongoing symptoms or brain activity and have the potential to provide improvements in terms of efficiency, efficacy and reduction of side effects. The success of closed-loop DBS depends on being able to devise a stimulation strategy that minimizes oscillations in neural activity associated with symptoms of motor disorders. A useful stepping stone towards this is to construct a mathematical model, which can describe how the brain oscillations should change when stimulation is applied at a particular state of the system. Our work focuses on the use of coupled oscillators to represent neurons in areas generating pathological oscillations. Using a reduced form of the Kuramoto model, we analyse how a patient should respond to stimulation when neural oscillations have a given phase and amplitude, provided a number of conditions are satisfied. For such patients, we predict that the best stimulation strategy should be phase specific but also that stimulation should have a greater effect if applied when the amplitude of brain oscillations is lower. We compare this surprising prediction with data obtained from ET patients. In light of our predictions, we also propose a new hybrid strategy which effectively combines two of the closed-loop strategies found in the literature, namely phase-locked and adaptive DBS.
  Symptoms of several neurological disorders, such as Parkinson’s disease (PD) and essential tremor (ET), are thought to arise from overly synchronous activity within neural populations. The severity of clinical impairment in PD is known to be correlated with an increase in the beta (13-35 Hz) oscillations in the local field potential (LFP) and in the activity of individual neurons in the basal ganglia [1,2]. The tremor symptoms associated with ET are thought to arise from synchronous activity in a network of brain areas including the thalamus [3]. In both PD and ET the muscle activity driving the tremor is coherent with local field potentials in the thalamus [4] and bursts of spikes produced by individual thalamic neurons [5]. Deep brain stimulation (DBS) is a well-established treatment option for PD and ET which involves delivering stimulation via electrodes implanted into the brain. The present generation of the technology involves manually tuning the parameters of stimulation, such as the pulse width, frequency and intensity in an attempt to achieve the best treatment. In particular, the choice of frequency is known to be crucial for efficacy, and high frequency DBS (120-180 Hz) has been found to be effective for PD and ET patients [6]. High frequency DBS is known to suppress the pathological oscillations occurring in PD [7], but despite its long history, the underlying mechanisms causing this suppression remain unclear, and several distinct theories have been proposed [8–10]. One influential theory suggests that high frequency DBS activates target neurons to such an extent that their synaptic transmission becomes saturated and they are no longer able to transmit pathological oscillations [11]. Since high frequency DBS can cause side-effects such as speech-impairments [12] and gambling tendencies [13], improvements to this treatment approach are desirable. It is thought that improvements could be achieved if future devices were to operate ‘closed-loop’, delivering stimulation only when needed and according to the ongoing symptoms of the patient [14]. A number of approaches to closed-loop DBS can be found in the literature and of these we focus on two, namely adaptive DBS [15] and phase-locked DBS [16]. In adaptive DBS high frequency stimulation is applied only when the amplitude of oscillations exceeds a certain threshold [15]. In phase-locked DBS stimulation is applied according to the instantaneous phase of the oscillations, which for ET patients corresponds to stimulation at roughly the tremor frequency (typically ∼ 5 Hz) [16]. The principles behind both of these approaches are illustrated in Fig 1. Together, these studies suggest that the effects of DBS are dependent on both the phase and amplitude of the oscillations at the time of stimulation. Modelling the effects of DBS generally poses a challenge since the brain networks involved in disorders such as ET (cortico-thalamic circuit) and PD (cortico-basal-gangla circuit) are complex and it is still debated from which parts of these circuits the pathological oscillations originate [17,18]. The task can be made more tractable by considering a simple phenomenological model which does not attempt to explicitly describe the underlying circuits, but rather focuses on general mechanisms leading to the synchronization of neurons. One example of this is the Kuramoto model, [19] where the dynamics of neurons are described using a system of homogeneously coupled oscillators, whose phases evolve according to a set of underlying differential equations. Such models are particularly attractive due to their simplicity and explicit dependence on phase, which makes them convenient for describing the effects of phase-locked stimulation. Coupled oscillator models have been used before to describe the effects of applying DBS at particular phases of ongoing oscillations [20,21]. In particular, Wilson and Moehlis [21] have described how the optimal intensity of stimulation should depend on the phase of the ongoing oscillations. The phase-dependent effects of stimulation predicted using the coupled-oscillator model [21] have been shown to generalize to other models, as they have been also observed in a biologically realistic model of a neural circuit generating oscillations related to PD [22]. However, the predictions of these models have not been directly compared with the experimental effects of closed-loop DBS. Furthermore, it has not been described how the effects of DBS should depend on another important characteristic of the ongoing oscillations—namely, the amplitude. A model describing the effects of DBS as a function of both phase and amplitude would be an important stepping stone towards designing optimal stimulation strategies. In this paper, we describe the mechanisms which give rise to the effects of phase-locked DBS using a reduced Kuramoto model. The general instantaneous response of an infinite Kuramoto system has been reported elsewhere [23,24] but has never been used in the context of DBS. Unique to this work are a number of contributions which are essential for connecting the theory with experimental data, including the relationship between the order parameter and the Hilbert transform and the theoretical interpretation of averaged response curves. We also present a framework for testing theoretical predictions on experimental data. Finally, we present a novel stimulation strategy based on our theories and existing methods found in the literature.
  We aim to describe how an underlying system of neural oscillators can give rise to oscillations, such as those found in LFP and tremor. In classic coupled oscillator models, neural oscillators correspond to individual neurons which spontaneously produce spikes [25]. In mathematical models of neurons (e. g. the Hodgkin-Huxley model) a state of a neuron is described by a set of variables, and for certain parameters, it produces spikes at regular intervals, thus its variables have periodic behaviour. A description of such a neuron can be simplified by projecting the state of the neuron on its phase, such that the state of each neuron is simply described by a one-dimensional phase variable, which in absence of external input increases with constant rate from phase 0 to 2π, corresponding to an evolution of a neuron from spike to spike [26]. The above interpretation of neural oscillators as regular-spiking neurons may be suitable for describing relatively fast beta oscillations occurring in PD, but during slower tremor oscillations thalamic neurons produce a burst of activity during a single cycle [5]. Therefore, if one interprets neural oscillators in ET as neurons, the phase would rather describe the changes in a variable governing the burst cycle (e. g. calcium level inside neuron). Alternatively, neural oscillators in ET can be interpreted as micro-circuits, which due to their internal connectivity produce oscillations in the activity of constituent neurons [27]. Since we aim here to develop a general theory, we simply consider the set of N neural oscillators with phases {θm (t) }. Nevertheless, we make two assumptions about the oscillators related to how they react to input, and how their phases determine the activity of the whole population. In the remainder of this subsection we describe these assumptions and discuss how they could be justified under different interpretations of individual oscillators mentioned above (regular-spiking neurons, bursting neurons and micro-circuits). First, we assume that if the stimulating input is provided to an oscillator when it is in phase θ, its phase will change according to a phase response function Z (θ). Biological regular-spiking neurons respond more to stimulation when they are already outside their refractory period and their membrane potential is closer to the spiking threshold [28]. Under such conditions, stimulation accelerates a neuron towards spiking. The function Z (θ) should therefore have higher values in the second part of the spiking cycle. In addition to this and under certain conditions, stimulation just after spiking and during the refractory period can slow a neuron’s spiking [28]. Therefore, Z (θ) may have negative values for small positive θ. Both of these characteristics can be captured using the phase response function Z (θ) = −sin (θ). We use this simplified phase response function in the first part of our paper and then later consider a more general form. Phase response functions for bursting neurons have also been calculated [29]. Phase response functions have been studied in a mathematical model of a micro-circuit composed of connected populations of excitatory and inhibitory neurons [30]. When such a micro-circuit receives an input, the phase of the oscillations it produces either advances or reduces depending on when within the cycle the input is provided [30]. Second, we define the average activity of neural oscillators as a superposition of cosine functions, i. e. f (t) =1N∑n=1Ncos[θn (t) ]. (1) We chose to transform the phase through a cosine function, because this periodic function has a maximum at 0, and in classic coupled oscillator model, phase 0 corresponds to the phase when neurons produce spikes [25]. When regular-spiking neurons are considered, their activity features spikes, rather than varying smoothly like the cosine function, but nevertheless, the effects of spikes on downstream neurons are prolonged in time due to non-instantaneous decay of post-synaptic potentials, thus the cosine function could be seen as a qualitative approximation of the effect of neuronal activity on downstream cells. If neural oscillators are assumed to correspond to bursting neurons or micro-circuits, the choice of the cosine function is more natural, because their rate of producing spikes varies more gradually. We would assume that the activity of bursting neurons or micro-circuits is highest at phase 0, so the function f (t) qualitatively reflects the fluctuations in firing rate. The purpose of this subsection is to relate the phases of individual oscillators to the quantities that can be measured from experimental LFP or tremor data, such as the instantaneous phase and the amplitude of the signal. The function f (t) is an abstract representation of neural activity and is not directly measurable in typical studies with patients. The experimentally measured signal fe (t) (in this case the accelerometer data from a patient’s hand) is known to be highly correlated to neural activity [5], so it is reasonable to assume fe (t) to be some transformed version of f (t). We assume this transformation to be a simple scaling and shifting, namely fe (t) =cf (t) +d. (2) It is common to subtract the mean from a signal, resulting in a signal f˜e (t) which is independent of d. This yields a simple relationship between the experimentally measured data and neural activity, f˜e (t) =cf (t). (3) The experimental signal f˜e (t) is typically filtered (around the dominant frequency) and then analysed using the Hilbert transform H^, which provides for each time point t the values of instantaneous phase ψe and the envelope amplitude ρe ρeeiψe=f˜e (t) +iH^[f˜e (t) ]. (4) This is the analytic signal for f˜e (t) —a complex function of time which has a magnitude of ρe (t) and a phase ψe (t). By inserting Eq (3) into Eq (4), we can relate the neural activity resulting from coupled oscillators to the experimental phase and amplitude ρeeiψe=c{f (t) +iH^[f (t) ]}. (5) We would like to now relate the experimental phase ψe and amplitude ρe to quantities obtained using the coupled oscillator model. The order parameter for a system of oscillators is defined to be r=1N∑n=1Neiθn. (6) Since r is a complex number, it can be written as r=ρeiψ, (7) The amplitude ρ is a measure of synchrony, with complete desynchrony and synchrony corresponding to 0 and 1, respectively. Using the Euler relation, the order parameter can be written as ρeiψ=1N∑n=1Ncos[θn (t) ]+i1N∑n=1Nsin[θn (t) ]. (8) We expect the phases {θn} to increase monotonically with time and under these conditions, H^[cos (θn) ]≃sin (θn). A justification for this can be found in the Supplementary Material S1. Using this and the expression for the time series given by Eqs (1) and (8) can be written as ρeiψ=f (t) +iH^[f (t) ]. (9) Comparing this with Eq (5), it can be seen that ρeeiψe=cρeiψ. (10) Therefore the experimental amplitude and phase is relatable to the magnitude and phase of the order parameter using ρe=cρ, ψe=ψ. (11) In summary, assuming the experimental data and neural activity are related according to Eq (2) and that the phases {θn} increase monotonically with time, we can use the Hilbert transform of the experimental data to relate the envelope amplitude and instantaneous phase to the magnitude and phase of the order parameter, respectively. The oscillation data we are concerned with arises from the correlated electrical activity of neural populations. In order to describe such systems we use a coupled oscillator model where the time-evolution for the set of N oscillators are given by the Kuramoto equations, with an additional term describing the effects of stimulation [19,20] dθldt=ωl+kN∑n=1Nsin (θn−θl) +IX (t) Z (θl). (12) The first term, ωl is the natural frequency of oscillator l, which describes the frequency in the absence of external inputs. It corresponds to the frequency with which a neuron spontaneously produces spikes or bursts (depending of the interpretation of oscillators introduced above). The second term describes the interactions between oscillators, where k is the coupling constant which controls the strength of coupling between each pair of oscillators and hence their tendency to synchronize. The third term describes the effect of stimulation. The intensity of stimulation is denoted by I and X (t) is a function which equals 1 if stimulation is applied at time t and 0 otherwise. The phase response function for a single oscillator is denoted by Z (θl). To avoid confusion, we use ‘intensity’ to refer to the magnitude of stimulation I, while the word ‘amplitude’ is used to refer to the amplitude of order parameter ρ or of the experimental signal ρe. Using the definition of the order parameter given in Eq (7), Eq (12) can be transformed to give dθldt=ωl+kρsin (ψ−θl) +IX (t) Z (θl). (13) In this form, it is clear that each oscillator has a tendency to move towards the population phase ψ and that the strength of this tendency is controlled by the coupling parameter k. To gain an intuition for this behaviour readers may wish to explore an online simulation of the model [31]. In the previous section, we outlined the conditions for which a relationship should hold between the experimental envelope amplitude and instantaneous phase and the quantities associated with the coupled oscillator model, namely the magnitude and phase of the order parameter. A patient’s response to phase-locked DBS is typically quantified using the amplitude response curve (ARC), which describes changes in the envelope amplitude of pathological oscillations (e. g. tremor) as a function of the phase at which the stimulation was delivered. Some studies also report the phase response curve (PRC), which describes changes in the phase of the pathological oscillation as a function of the stimulation phase [16]. From a theoretical perspective, the response curves arise from changes to an underlying state of oscillators which in turn gives rise to measurable changes in the ‘macroscopic’ quantities, namely the amplitude and phase of the order parameter. Therefore, in order to obtain an analytical expression for the response curves, we need to know how the order parameter evolves as a function of time. Ott and Antonsen showed that such an expression can be found under the assumption of an infinite system of oscillators and where the distribution of frequencies g (ω) is Cauchy with centre ω0 and width γ. In this section, we summarize their results but the interested reader can refer to the original paper for a full derivation. A useful review on reduced coupled oscillator models can be found in reference [32]. For an infinite system of oscillators, the order parameter can be expressed in terms of the distribution of oscillators f (ω, θ, t) r (t) =∫−∞∞∫02πf (ω, θ, t) eiθdωdθ, (14) and the time evolution of f (ω, θ, t) for the Kuramoto system given by Eq (13) is given by the continuity equation ∂f (ω, θ, t) ∂t+∂∂θ{[ω+k2i (re−iθ−r*e−iθ) +IX (t) Z (θ) ]f (ω, θ, t) }=0. (15) Central to the work of Ott and Antonsen [33] is the use of a guess, or ansatz, for the distribution of oscillators given by f (θ, ω, t) =g (ω) 2π{1+∑n=1∞[αn (ω, t) einθ+αn (ω, t) *e−inθ]}, (16) where α (ω, t) is a certain function. Inserting (16) into (15) leads to a partial differential equation for α (ω, t). Ott and Antonsen [33] considered the case of a Kuramoto system with a periodic driving term of strength I and frequency Ω, whose dynamical equations are given by dθldt=ωl+kN∑n=1Nsin (θn−θl) +Isin (Ωt−θl). (17) Using this, the ansatz and the result that r (t) = α* (ω0 − iγ, t) [23], the time evolution for the order parameter is found to satisfy (for γ = 1) drdt=12[ (kr+I) − (kr+I) *r2]−[1−i (Ω−ω0) ]r. (18) Therefore the state of the Kuramoto system for a large number of oscillators has been reduced from one being described by the set of N phases {θj} to one being described by ρ and ψ. Furthermore, Eqs (17) and (18) can be used to derive an expression for a system undergoing pulsatile stimulation, which will be discussed in the next section.
  Although the effect of phase-locked DBS may depend on the amplitude of the ongoing pathological oscillations, this dependence on amplitude has not been analysed before, due primarily to the difficulties associated with obtaining a function of two independent variables from noisy data. Instead, the averaged response curves have been reported [16,34], which are only functions of the phase and are averaged over the amplitude. Such curves are readily obtainable using standard signal processing techniques. In this subsection, we derive expressions for the ARC and PRC, where the phase response function for a single oscillator is taken to be Z (θ) = −sin (θ). By setting Ω = 0 and for general γ, Eq (17) describes a system experiencing an instantaneous impulse (i. e. X (t) = 1, during infinitesimally brief intervals at which stimulation is provided) with the phase response function Z (θ) = −sin (θ). The corresponding expression for the time evolution of the order parameter can be found by setting Ω = 0 in Eq (18), leading to drdt= (iω0−γ) r+kr2 (1−|r|2) +I2 (1−r2). (19) Inserting the expression for the order parameter Eq (7) into Eq (19) gives expressions for the time evolution of ρ and ψ dρdt=−γρ+kρ2 (1−ρ2) +P (ρ, ψ) (20) and dψdt=ω0+Ψ (ρ, ψ), (21) where P (ρ, ψ) =I2 (1−ρ2) cos (ψ) (22) and Ψ (ρ, ψ) =−I2ρ (1+ρ2) sin (ψ). (23) The functions P (ρ, ψ) and Ψ (ρ, ψ) are the instantaneous response curves for a population of oscillators with a phase response function of Z (θ) = −sin (θ). We denote the response functions by using the uppercase of the appropriate quantity. This simplified case leads to very specific qualitative predictions. For both response curves, the effects of stimulation are predicted to be more magnified when stimulation is applied at lower amplitudes. In addition to this, whether or not stimulation has an amplifying or suppressing effect on the respective quantities is dependent only on the phase. In the absence of stimulation, the first term of Eq (20) predicts the amplitude to decay with a rate proportional to the diversity of natural frequencies of individual oscillators γ whilst the second term predicts the amplitude to increase according to a term proportional to the coupling strength k. As t → ∞, ρ tends asymptotically to a non-zero finite value for k > 2γ and to zero for k < 2γ. Eq (21) predicts the phase ψ to evolve according to ω0. To demonstrate the amplitude dependent effects of stimulation predicted by Eqs (22) and (23) we simulate the Kuramoto model using a large number of oscillators (N = 3000). The stimulation amplitude I was chosen to be small at ≃ 0. 04 and numerical integration was performed using the Euler method with a time step of δt ≃ 0. 001. The parameters ω0 and γ were not expected to affect the response curves and hence were arbitrarily chosen to be ω0 = 30 and γ = 1. For such a system and in the absence of stimulation, the magnitude of the order parameter ρ will tend asymptotically to [1 − (2/k) ]0. 5 for k > 2 [33]. We can therefore fix the value of ρ by choosing an appropriate value of k. After the system has evolved to the asymptotic state, we provide pulsatile stimulation at a particular phase. The changes in ρ and ψ resulting from the perturbation (over a single time step of numerical integration) divided by δt would then approximately equal to P (ρ, ψ) and Ψ (ρ, ψ), respectively. Fig 2 shows the response functions P (ρ, ψ) and Ψ (ρ, ψ) for different amplitudes ρ and also a comparison with results from simulating a population of Kuramoto oscillators. Fig 2a and Eq (22) show the ARC is shifted with respect to Z (θ) and that, for a given ρ, the most effective reduction of oscillation amplitude is achieved when phase-locked stimulation is provided at phase π. An intuition for these effects is shown in Fig 3a. The form of the phase response function Z (θ) leads to a region of phases for which the oscillators will either slow down or speed up upon stimulation. Stimulation applied to a population of oscillators corresponds to a perturbation with a differential effect across the system of oscillators, i. e. with some oscillators responding differently to others and depending on their phase. It is this differential effect which gives rise to changes in the width of the oscillator distribution and therefore changes in amplitude. In particular, if stimulation is applied when the distribution of oscillators is centred around 180 degrees, stimulation is shown to have a desynchronising effect as half the oscillators will speed up and the other half will slow down, as illustrated in Fig 3a. An intuition for the amplitude dependent effects can also be seen in Fig 3a. If we consider the case where the system is strongly synchronised (orange curve), then stimulation can have little effect on amplitude since all oscillators would be perturbed by a similar amount. If we now increase the width of the oscillator distribution (red curve) and the amplitude reduces, then the differential effects give rise to a greater amplitude change, thus stimulation at lower amplitudes leads to magnified effects. It is also evident from Fig 2 that there exists a relationship between the ARC and slope of the PRC. In particular, the ARC is negative at those phases for which the slope of the PRC is positive. One can also see from Eqs (22) and (23) that for a given ρ, the ARC (P (ρ, ψ) ) is proportional to the negative derivative of the PRC with respect to ψ (∂Ψ (ρ, ψ) ∂ψ). We will later analyse how this relationship generalizes. In this subsection we consider the case where the phase response function Z (θ) has a general form, since the phase response curves of biological neurons may have diverse shapes [28]. We start by providing an intuition for why the qualitative effects of stimulation are expected to be different when Z (θ) contains higher harmonics. For the case of Z (θ) = sin (6θ), as shown in Fig 3b, and the simple case of oscillators distributed around 180 degrees, the effects of stimulation may actually be greater at larger amplitudes—which is in contrast to our results for the single harmonic case. In the high amplitude regime, as shown by the orange curve, the qualitative effects of stimulation are predicted to be similar to that of the single harmonic case. At lower amplitudes, some of the oscillators will be shifted away from the centre of the distribution while other towards it, reducing the overall effect of stimulation. To analyse this more formally, we now consider the case where Z (θ) takes the form of a general Fourier series Z (θ) =a02+∑m=1∞amcos (mθ) +∑m=1∞bmsin (mθ). (24) Using the results from Lai and Porter [23], an expression for the time evolution of the order parameter can be obtained, drdt= (iω0−γ) r+kr2 (1−|r|2) +iI2{a0r+∑m=1∞am[ (r*) m−1+rm+1]+i∑m=1∞bm[ (r*) m−1−rm+1]}. (25) Inserting the expression for the order parameter (Eq (7) ) into Eq (25), we find expressions for the time evolution of ρ and ψ, but now the instantaneous response curves for amplitude and phase, respectively, are given by (cf. [24]) P (ρ, ψ) =I2 (1−ρ2) ∑m=1∞ρm−1[amsin (mψ) −bmcos (mψ) ], (26) and Ψ (ρ, ψ) =I2{a0+ (1+ρ−2) ∑m=1∞ρm[amcos (mψ) +bmsin (mψ) ]}. (27) Eqs (26) and (27) describe the instantaneous response curves for a system of oscillators whose distribution satisfies Eq (16). It is worth noting that both equations are independent of the parameters of the Kuramoto model and are only dependent on the characteristics of an oscillator distribution satisfying the ansatz. It is known [23,35,36] that the presence of higher harmonic modes in the phase response function Z (θ) can cause the oscillators to cluster and lead to a breakdown of the ansatz given by Eq (16). Lai and Porter [23] investigated the effects of introducing noise through the phase response function and found good agreement between theory and simulation only when Z (θ) consisted of a dominant first harmonic mode. For Z (θ) containing higher harmonics, a deviation between theory and simulation was found [23], the extent of which was found to be proportional to the noise strength. It is important to emphasise here that the deviation between simulation and theory is only dependent on the state of the system, i. e. the distribution of oscillators, and on the phase response function Z (θ), not on whether the stimulus is deterministic. We therefore expect Eqs (26) and (27) to reasonably approximate the response when Z (θ) has a dominant first harmonic and/or the stimulation amplitude I is small [37]. Using the methodologies from before, we simulated a population of Kuramoto oscillators to demonstrate the predicted amplitude dependence of stimulation for the case of Z (θ) containing higher harmonics. Fig 4 shows an example of the ARC and PRC for the case of Z (θ) = −sin (3θ). In contrast to Fig 2, the effects of stimulation are not necessarily monotonic functions of ρ. This can be seen for the case of the ARC shown in Fig 4a, where the effects are magnified between ρ = 0. 35 and ρ = 0. 70 but are reduced between ρ = 0. 70 and ρ = 0. 90. For the PRC in Fig 4b, it is clear that the effects are monotonic but that the effects of stimulation now reduce with reducing amplitude, in contrast with the single harmonic case. A biologically realistic phase response function Z (θ) is likely to contain multiple harmonics and, as such, the qualitative effects of stimulation are expected to be different from the single harmonic case. To investigate this, we obtained the phase response function for a single spiking neuron using the Hodgkin Huxley equations [25,38]. The equations were integrated using the Euler method with a time step of δt = 0. 001. The parameters used in the equations were taken from Brown et al [25] and the phase was calculated according to the Hilbert transform. The PRC was calculated according to a single pulse of stimulation delivered at a particular phase. A general Fourier expansion of the form given by Eq (24) with 4 harmonics was fitted to the resulting PRC to obtain the Fourier coefficients. The Kuramoto equations were simulated, as before, with these coefficients. The amplitude dependent response curves are shown in Fig 5a and 5b. The phase response function Z (θ) (as fitted using 4 harmonics of a Fourier expansion) is shown in Fig 5c. The polar plane amplitude response is shown in Fig 5d. Here, the amplitude response is calculated according to Eq (26). The magnitude of a vector from the origin to a point on the plane is equal to ρ and its angle made with the horizontal axis is equal to ψ. We observe that the phase-locked stimulation suppresses oscillations most if applied at low amplitude (and at the most effective phase). We see this effect in response functions with a single first harmonic, as shown in Fig 2 and it occurs here because the phase response function shown in Fig 5c has a dominant first harmonic. The expressions for the instantaneous response curves can be used to make qualitative predictions about how a subject should respond to stimulation. Eqs (26) and (27) involve an expansion according to the harmonics of Z (θ), with each term being the product of both a phase dependent part and an amplitude dependent part. If we restrict our analysis to the simple case where Z (θ) is well-approximated to be a single harmonic mode, then it can be seen that whether stimulation has a suppressive or amplifying effect on amplitude is dependent only on the phase, while the magnitude of these effects is determined by the amplitude. Plots of the amplitude dependent part for several harmonic modes can be seen in Fig 6 for both the PRC and ARC. For the case of Z (θ) containing a single harmonic, each curve describes how the magnitude of the response curves is expected to change as a function of amplitude. Fig 6a shows the magnitude of the ARC as a function of ρ for several harmonic modes. A non-monotonic relationship is predicted only for higher modes m > 1, with the maxima occurring for each curve at ρ˜m=[m−1m+1]12. (28) In the high synchrony regime ρ>ρ˜m, the gradient in each case is found to be negative, implying that when stimulation is applied at lower ρ, the effects of stimulation increase. In the low synchrony regime ρ<ρ˜m, stimulation applied at lower amplitudes ρ is predicted to lead to smaller effects. For the case of the first harmonic mode, where ρ˜1=0 implies that the gradient is negative across the range 0 ≤ ρ ≤ 1, the predicted effects are particularly noteworthy as they are both quantitatively and qualitatively different from the other modes. Concisely, for Z (θ) consisting of a single first harmonic mode, we predict that delivering stimulation at lower ρ will result in greater effects, for all values of ρ. Fig 6b shows the magnitude of the PRC as a function of amplitude for several harmonic modes. Here we find the qualitative differences between the effects of the first harmonic mode and higher modes to be even more notable. A monotonic relationship is predicted across the harmonic modes but with differing gradients between the first and higher modes. The condition for positive (or zero) gradients is that ρ2≥− (m−2) m, which is only satisfied for m > 1, hence for higher modes, a positive (or zero) gradient can be found for all ρ in contrast to a negative gradient for all ρ for the first mode. Qualitatively, this means that for Z (θ) consisting of a dominant higher harmonic, delivering stimulation at lower ρ will result in smaller changes in phase. For Z (θ) consisting of a dominant first harmonic, the opposite is predicted, namely delivering stimulation at lower ρ will result in larger effects. These effects become particularly apparent if the oscillation amplitude is close to 0, where the ρ−1 term attached to the first harmonic begins to become very large. Finally, it is worth noting that for stimulation to have an effect on amplitude, the function Z (θ) does not need to have a region where Z (θ) < 0 (i. e. it does not need to be of type II). Although in the examples we used Z (θ) with negative regions, Z (θ) can be shifted by adding a constant a0, and this constant will not affect ARC, because it does not appear in Eq (26). A critical condition for the stimulation to have an effect on amplitude is that the function Z (θ) is not constant, which allows the stimulation to have a differential effect on oscillators in different phases. Existing experimental studies have reported the ARC and PRC averaged across amplitudes of pathological oscillations [16]. In this subsection, we study the properties of such averaged curves. In the next subsection we use this relationship to test if the response to stimulation of individual patients is well described by the Kuramoto model. Expressions for the averaged response as a function of ψ can be obtained by taking expectation values using the function h (ρ|ψ), which is the probability density function for the system being in a state ρ given a phase ψ. In the absence of stimulation, the dynamics of ρ for the Kuramoto system are phase-shift invariant. Therefore, if the effects of stimulation are small, it is reasonable to assume that h (ρ|ψ) ≃ h (ρ) and P¯ (ψ) =∫01h (ρ) P (ρ, ψ) dρ=I2∑m=1∞vm[amsin (mψ) −bmcos (mψ) ], (29) and Ψ¯ (ψ) =∫01h (ρ) Ψ (ρ, ψ) dρ=I2{a0+∑m=1∞v˜m[amcos (mψ) +bmsin (mψ) ]}, (30) where vm=∫01h (ρ) (1−ρ2) ρm−1dρ, (31) and v˜m=∫01h (ρ) (1+ρ−2) ρmdρ. (32) We first describe the relationship between the ARC and PRC for the cases where Z (θ) contains a single dominant harmonic as in these cases clear predictions can be made by the model. First, we use the derivative of Ψ¯ (ψ) giving dΨ¯dψ=−I2∑m=1∞mv˜m[amsin (mψ) −bmcos (mψ) ] (33) and also by dividing P¯ (ψ) by dΨ¯dψ leads to P¯ (ψ) =−dΨ¯dψ{∑m=1∞vm[amsin (mψ) −bmcos (mψ) ]∑m=1∞mv˜m[amsin (mψ) −bmcos (mψ) ]}. (34) Now considering the case where Z (θ) contains only the qth harmonic P¯ (ψ) =−dΨ¯dψ{vq[aqsin (qψ) −bqcos (qψ) ]qv˜q[aqsin (qψ) −bqcos (qψ) ]}=− (vqqv˜q) dΨ¯dψ, (35) which shows that in the cases where Z (θ) contains a single harmonic the averaged ARC is a scaled version of the negative gradient of the averaged PRC. Since the constants vq, vq˜ and q are all positive, it is expected that for Z (θ) containing a single or dominant harmonic, a strong positive correlation should exist between P¯ (ψ) and −Ψ¯′ (ψ). To investigate the correlation coefficient for higher harmonics, we simulated the Kuramoto model using randomly generated phase response functions. The parameters of the Kuramoto model (N, ω0, γ) were the same as in previous simulations. The number of harmonics (nh) to use in each function was chosen sequentially from 1 to 4. In each case, the harmonics of each function were chosen to be a random subset of size nh from the set {1,2, 3,4}. A random phase response function Z (θ) was generated by choosing a set of coefficients {am} and {bm} (which includes randomising a0) whose values were sampled from a standard normal distribution. For each Z (θ) with a given number of harmonics, the response curves as a function of phase were calculated for several values of synchrony ρ = {0. 4,0. 6,0. 8}. The phases were chosen from a uniformly spaced grid between 0 and 2π. The response of the system was taken to be from a single pulse of stimulation. For a given number of harmonics, 30 averaged response curves were calculated. The correlation coefficient between P¯ (ψ) and −Ψ¯′ (ψ) was then calculated by first averaging P (ρ, ψ) an Ψ (ρ, ψ) across ρ and then calculating −Ψ¯′ (ψ) by averaging the gradient in the forward and backward direction around a particular phase ψ. Fig 7 shows the value of the correlation coefficient to be only slightly affected by increases in the number of harmonics. This implies that if a system is well described by the Kuramoto model, then a strong positive correlation between P¯ (ψ) and −Ψ¯′ (ψ) should exist. The results of this subsection extend previous observations [21,39] that stimulation is most effective when applied at the phase which maximizes Z′ (θ). However, the phase response curves of individual oscillators Z (θ) are difficult to measure. Nevertheless, Eq (35) shows an analogous relationship to exist between averaged population response curves. Together with Eq (11), the above analysis predicts that there also should exist a positive correlation between experimentally measured ARC and the gradient of experimentally measured PRC. The results presented in previous sections provide a framework within which we can make some qualitative predictions about how a patient should respond to phase-locked DBS. In this subsection we compare two key predictions of the theory with experimental data: There should exist a strong correlation between averaged ARC (P¯ (ψ) ) and the negative gradient of the averaged PRC (−Ψ¯′ (ψ) ). If the phase response function Z (θ) contains a dominant first harmonic, then the effects of stimulation should be magnified if it is applied when the amplitude of oscillation ρ is low. We tested these predictions using data from the study of Cagnan et al [16]. In this study, phase-locked DBS was delivered according to the tremor measured by an accelerometer attached to the patient’s hand. Data was collected from 6 ET patients and 3 dystonic tremor patients. All patients gave their informed consent to take part in the study, which was approved by the local ethics committee in accordance with the Declaration of Helsinki. This data can be obtained through an online repository [40]. Duchet et al [41] defined a criterion for assessing significance in the averaged ARCs and PRCs from the study of Cagnan et al. In this criterion, a patient is considered to have a significant response if both the ARC and PRC are found to be significant either according to an ANOVA test or cosine model F-test. Using this, they deemed 3 out of the 6 ET patients to have a significant response curve. We restrict our analysis to these 3 patients here. A flowchart illustrating the various processes involved in calculating the response curves from experimental data is shown in Fig 8a. The tremor data was filtered using a non-causal Butterworth filter of order 2 with cut-off frequencies at ±2 Hz around the tremor frequency. Stimulation was delivered over a set of trials (typically 9), with each trial consisting of 12 blocks of 5 second phase-locked stimulation at a randomly chosen phase from a set of 12. Each block of phase-locked stimulation was also separated by a 1 second interblock of no stimulation. The envelope amplitude and instantaneous phase were calculated using the Hilbert transform. The methods for calculating the averaged ARCs and PRCs are shown in Fig 8b. The averaged amplitude response for a particular phase was calculated to be the difference between the average envelope amplitude within a 1 second window before the end of the stimulation block and the average envelope amplitude within a 1 second window prior to the onset of stimulation. The averaged phase response was calculated using a similar methodology. The unwrapped phase was calculated for the data 1 second prior to the onset of stimulation. A linear function was fitted to this phase evolution and extrapolated to the end of the stimulation block. The value for the instantaneous phase obtained using this extrapolation could be taken as the expected phase of the system in the absence of stimulation. The difference between the actual phase at the end of the stimulation block and this expected phase was taken to be the phase response. In both cases, the responses for a particular phase were averaged over all trials in the dataset. The derivative of the PRC with respect to ψ was calculated numerically by averaging the gradient in the forward and backward direction around a particular phase ψ. To determine the effects of amplitude on the response curves, we use a ‘single pulse’ method, as shown in Fig 8c. The data is binned at low, medium and high amplitudes, with each bin containing the same number of points. Within each bin we calculate the response of the system from a single pulse of stimulation. In the case of the amplitude response, we calculate the difference in the mean of the amplitude after and before the pulse. The data used for calculating the mean in each case is taken between pulses. In the case of the phase response, a straight line is fitted to the unwrapped phase evolution from before the pulse. The phase response is taken to be the difference between the actual phase and the linear extrapolation evaluated at a point after the pulse and just before the next pulse. We first tested Prediction 1 for each patient and excluded those patients which either do not have significant correlation (p < 0. 05) or have negative correlation. The lack of significant correlation between P¯ (ψ) and −Ψ¯′ (ψ) could indicate that the patient is not well described by the Kuramoto model or that the response curves have not been accurately determined. We therefore restrict subsequent testing of Prediction 2 to only those patients who exhibit significant correlation. To infer whether the phase response function Z (θ) for a given patient contains a dominant first harmonic, we used the property illustrated in Fig 6b, namely that the magnitude of the PRC should increase with reducing ρ only if Z (θ) contains a single first harmonic. For such patients, we would then expect the magnitude of the ARC to increase with reducing ρ. Table 1 gives the correlation coefficients between P¯ (ψ) and −Ψ¯′ (ψ) for each patient where DBS was found to have a significant effect on tremor. From these, significant correlation was only found for Patient 5, with the correlation coefficient being (0. 72). The response curves for this patient are shown in Fig 9. Table 1 also lists uncorrected p-values for the correlation. In summary, we did not find strong support for Prediction 1, as the correlation between the ARC and the negative gradient of the PRC was significant in only 1 out of the 3 datasets analysed. As outlined above, we then tested Prediction 2 only for Patient 5 who fulfilled Prediction 1. To determine the effects of amplitude on the response curves, we analysed the response of the system using the single pulse method. This is because the amplitude of tremor can vary substantially within the 5-second stimulation intervals. Fig 10 shows the response curves at 3 amplitude bins for Patient 5. To quantify the magnification of the response curves, we compute the standard deviation across the phase bins. The uncertainty on the standard deviation σ is calculated using the method of propagating errors [42]. We find this to be most appropriate here as it allows us to incorporate the errors of the binned response curves. For a dataset consisting of M points Y = {yj} with mean y¯ and corresponding standard errors ΔY = {Δyj}, we find the error on the standard deviation Δσ to be Δσ=1Mσ{∑j=1M[ (yj−y¯) Δyj]2}0. 5. (36) Fig 11b shows the extent of magnification across 3 amplitude bins for the PRC. For Patient 5, we find evidence for increasing magnification with reducing amplitude, which when taken together with the positive correlation between P¯ (ψ) and −Ψ¯′ (ψ), is indicative of a Kuramoto system with a phase response function Z (θ) containing a dominant first harmonic mode. For such a system, Prediction 2 states that we should also find increasing magnification of the ARC with reducing amplitude. Fig 11a shows the extent of magnification across 3 amplitude bins for the ARC. Here we can see evidence for increasing magnification with reducing amplitude for Patient 5, which agrees with our predictions and was the only patient to exhibit such effects across all the amplitude bins. For completeness, Fig 11 also shows the results of the above analysis for other patients. It is evident in Fig 11b that, for all patients, the magnification of the PRC is higher for stimulation at lower ρ. Interestingly, Fig 11a shows that the magnification of the ARC for stimulation at lower ρ is only exhibited by Patient 5, which is also the only patient who had significant correlation between the ARC and negative gradient of the PRC. A summary of our analysis is shown in Fig 12.
  We find good agreement between our theories and the data from one patient, namely a strong positive correlation between P¯ (ψ) and −Ψ¯′ (ψ) together with increasing magnification of the PRC with reducing amplitude is associated with an increasing magnification of the ARC with reducing amplitude. The lack of significant correlation between P¯ (ψ) and −Ψ¯′ (ψ) for the other patients could be due to experimental noise, which may prevent the response curves from being determined accurately. However, another possibility is that these patients are simply not well described by the models presented here. This could be due to the various assumptions, both about the nature of stimulation and the form for the distribution of oscillators, which are used to derive the theoretical response curves. Since the theoretical response curves are only valid for those distributions satisfying the ansatz given by Eq (16), the response curves for any distribution deviating from this form are expected to be different from those predicted. Distributions not represented by the ansatz include those clustered configurations which can arise through random effects or stimulation applied through phase response functions containing higher harmonic modes. Eq (13) also assumes that each oscillator responds only according to a single phase response function Z (θ), which may not be a good approximation for certain patients. Furthermore, our computational model does not capture the effects of synaptic plasticity triggered by stimulation. Presence of such plasticity may be suggested by a delayed appearance of tremor following offset of the prolonged phase-locked stimulation [16]. In addition to this, ET is known to be a heterogeneous disorder [43] and different patients are likely to have different underlying pathologies. As a result of this, the assumptions used in the model may need to differ depending on the patient. Our theories predict that an underlying Kuramoto system should be associated with a strong correlation between P¯ (ψ) and −Ψ¯′ (ψ). Since this relationship was not found in a number of patients, it is natural to wonder how the model could be modified to account for this. One assumption of the model is that each oscillator feels the same amplitude of stimulation. Since we interpret an oscillator as representing the activity of neurons or micro-circuits, it follows that these neurons should have some spatial separation in the brain [44] and hence experience stimulation differently depending on their location relative to the electrode. We investigated how providing sufficiently large stimulation to a fraction of oscillators can lead to changes in the correlation between P¯ (ψ) and −Ψ¯′ (ψ). The instantaneous response curves were calculated as before using an ensemble of 1000 oscillators and with Z (θ) = 1 − sin (θ). Fig 14a shows that when the oscillators are highly synchronized, the most effective phase to stimulate approaches 3π2, i. e. the phase at which the oscillators are most advanced by the stimulation (see Fig 14b). At this phase, half of the oscillators are maximally affected by stimulation and moved away from the rest. Fig 14c shows that the extent to which the correlation can be weakened depends on both the fraction of oscillators stimulated and the value of ρ at which the stimulation is applied. In addition to this, the amplitude dependence of the instantaneous ARC can also deviate from that predicted for a phase response function with a single first harmonic. For the case where half the oscillators are stimulated, as shown in Fig 14a, the minimum of the ARC shifts towards the maximum of the PRC with increasing synchrony. In light of our predictions, we propose a new strategy for DBS which may be effective for certain patients, namely those for whom the effects of DBS are magnified when stimulation is applied at lower amplitudes. The approach, illustrated in Fig 15, combines the aforementioned phase-locked and adaptive approaches described in Fig 1. The general idea is to only apply high frequency DBS at high amplitudes in order to drive the tremor into low synchrony regimes that are more susceptible to phase-locked DBS. With such a control strategy, the high frequency DBS would be used less often than in the current adaptive DBS approach, as it would only be used to bring the tremor into the low synchrony regime, at which point phase-locked DBS would suppress the amplitude further, thereby keeping the system in this mode. Separating the high and low synchrony regimes can be done using an amplitude threshold, in a similar way to the adaptive DBS approach [15]. When the amplitude of oscillations exceeds the threshold, then high frequency DBS would be delivered since the patient would be in the high synchrony regime. If the amplitude of oscillations falls below the threshold, then phase-locked DBS is delivered, since the patient would be in the low synchrony regime. In addition to the various parameters of stimulation, the choice of threshold is also likely to be an important factor in determining the overall efficacy of the method [45]. One possible issue concerns the phase of the system, which is known to become increasingly harder to define at lower amplitudes. At very low amplitudes we therefore expect phase-locked stimulation to be equivalent to random stimulation—the efficacy of which would need to be determined. There are a number of possible extensions to our work. On the theoretical side, investigating the effects of clustering on the response curves in addition to considering a multi-population Kuramoto model would be two possible avenues to explore. It is also worth mentioning that the instantaneous response curves given by Eqs (26) and (27) are independent of the parameters of the Kuramoto model and can be derived using only the assumption of the Ott and Antonsen ansatz. Therefore, in principle, they should be valid for other systems for which the Ott and Antonsen ansatz can be applied, such as an infinite network of theta neurons [46]. Kuramoto-like phase models arise through the phase reduction of oscillating units. By contrast, theta neurons can be in both excitable and oscillatory regimes while still being amenable to the Ott-Antonsen reduction [47]. It would be interesting to relate the phase-response curves of the population dynamics [48] to the data presented here. Furthermore, the model analysed in this paper assumes a Cauchy distribution for the natural frequencies, which has very long tails, and other distributions (e. g. Gaussian) may be a more realistic description of neuronal frequencies. It has recently been demonstrated that the Ott and Antonsen ansatz can be applied to systems where the natural frequencies are Gaussian distributed [49], and it would be interesting to extend the analysis of population response curves to this case. On the experimental side, there is the question of the hybrid DBS strategy, whose efficacy remains to be determined. To further strengthen our conclusions, we also hope to perform our analyses on more data. Given that testing our theories is conditional on finding patients for whom P¯ (ψ) is correlated to −Ψ¯′ (ψ), a study such as ours would greatly benefit from a larger dataset, both in terms of the number of patients and the length of time each patient is stimulated for. The latter would allow us to determine the response curves more accurately, which is expected to be particularly beneficial for our methods. Obtaining longer datasets poses a challenge due to the inherent difficulties associated with recording from patients, particularly the onset of fatigue. Alternatively, improvements to the accuracy of the response curves could be realised by improving the methods used to calculate them.
summary: Deep brain stimulation (DBS) involves delivering electrical impulses to target sites within the brain and is a proven therapy for a variety of neurological disorders. Closed loop DBS is a promising new approach where stimulation is applied according to the state of a patient. Crucial to the success of this approach is being able to predict how a patient should respond to stimulation. Our work focusses on DBS as applied to patients with essential tremor (ET). On the basis of a theoretical model, which describes neurons as oscillators that respond to stimulation and have a certain tendency to synchronize, we provide predictions for how a patient should respond when stimulation is applied at a particular phase and amplitude of the ongoing tremor oscillations. Previous experimental studies of closed loop DBS provided stimulation either on the basis of ongoing phase or amplitude of pathological oscillations. Our study suggests how both of these measurements can be used to control stimulation. As part of this work, we also look for evidence for our theories in experimental data and find our predictions to be satisfied in one patient. The insights obtained from this work should lead to a better understanding of how to optimise closed loop DBS strategies.
keywords: medicine and health sciences
  action potentials
  myoclonus
  pathology and laboratory medicine
  nervous system
  neurodegenerative diseases
  membrane potential
  brain electrophysiology
  electrophysiology
  neuroscience
  surgical and invasive medical procedures
  deep-brain stimulation
  mathematics
  functional electrical stimulation
  signs and symptoms
  brain mapping
  bioassays and physiological analysis
  distribution curves
  research and analysis methods
  statistical distributions
  animal cells
  neural pathways
  electrophysiological techniques
  probability theory
  movement disorders
  cellular neuroscience
  diagnostic medicine
  neuroanatomy
  cell biology
  anatomy
  physiology
  neurons
  neurology
  parkinson disease
  biology and life sciences
  cellular types
  physical sciences
  neurophysiology
title: Predicting the effects of deep brain stimulation using a reduced coupled oscillator model