Self-Tuning Deep Brain Stimulation Controller for Suppression of Beta Oscillations: Analytical Derivation and Numerical Validation

Closed-loop control strategies for deep brain stimulation (DBS) in Parkinson's disease offer the potential to provide more effective control of patient symptoms and fewer side effects than continuous stimulation, while reducing battery consumption. Most of the closed-loop methods proposed and tested to-date rely on controller parameters, such as controller gains, that remain constant over time. While the controller may operate effectively close to the operating point for which it is set, providing benefits when compared to conventional open-loop DBS, it may perform sub-optimally if the operating conditions evolve. Such changes may result from, for example, diurnal variation in symptoms, disease progression or changes in the properties of the electrode-tissue interface. In contrast, an adaptive or “self-tuning” control mechanism has the potential to accommodate slowly varying changes in system properties over a period of days, months, or years. Such an adaptive mechanism would automatically adjust the controller parameters to maintain the desired performance while limiting side effects, despite changes in the system operating point. In this paper, two neural modeling approaches are utilized to derive and test an adaptive control scheme for closed-loop DBS, whereby the gain of a feedback controller is continuously adjusted to sustain suppression of pathological beta-band oscillatory activity at a desired target level. First, the controller is derived based on a simplified firing-rate model of the reciprocally connected subthalamic nucleus (STN) and globus pallidus (GPe). Its efficacy is shown both when pathological oscillations are generated endogenously within the STN-GPe network and when they arise in response to exogenous cortical STN inputs. To account for more realistic biological features, the control scheme is then tested in a physiologically detailed model of the cortical basal ganglia network, comprised of individual conductance-based spiking neurons, and simulates the coupled DBS electric field and STN local field potential. Compared to proportional feedback methods without gain adaptation, the proposed adaptive controller was able to suppress beta-band oscillations with less power consumption, even as the properties of the controlled system evolve over time due to alterations in the target for beta suppression, beta fluctuations and variations in the electrode impedance.


INTRODUCTION
Deep brain stimulation (DBS) is a clinically effective treatment used in patients with advanced Parkinson's disease (PD) to supplement or replace pharmacological treatment of symptoms. It consists of high-frequency stimulation of neurons within the basal ganglia, with the subthalamic nucleus and globus pallidus being the most common targets, using a chronically implanted electrode and a subcutaneous pulse generator. At present, DBS is delivered clinically in an open-loop fashion where stimulation parameters remain fixed over time. This approach, however, may lead to overstimulation, inducing side effects and shortening battery life. Regular retuning of parameters is also required, involving a time-consuming trial and error process (Lozano et al., 2019).
The search for new approaches to improve the attenuation of patient symptoms and minimize side effects while increasing battery life has motivated a growing interest in closed-loop DBS over the last several years (Santaniello et al., 2010;Rosin et al., 2011;Santos et al., 2011;Carron et al., 2013;Beuter et al., 2014;Shah et al., 2018;Eitan et al., 2019). In a closedloop paradigm, the stimulation waveform is modified based on specific biomarkers which are used as a surrogate measure of symptom severity. To date, beta-band (13-30 Hz) activity in the STN local field potential (LFP) has been one of the most widely investigated biomarkers for closed-loop DBS during Parkinson's disease (Parastarfeizabadi and Kouzani, 2017). Increased betaband power in the STN LFP is correlated with motor impairment symptoms in Parkinson's disease, and its suppression, due to medication or high frequency DBS, with improved motor performance (Kühn et al., 2006(Kühn et al., , 2008Hammond et al., 2007;Eusebio et al., 2011).
Closed-loop stimulation utilizing LFP beta-band power has been experimentally validated in patients for short periods (Rosin et al., 2011;Little et al., 2013Little et al., , 2016Rosa et al., 2015;Arlotti et al., 2018;Velisar et al., 2019), however, longer term investigations have yet to be conducted. The improvement in symptoms has been comparable to that with conventional openloop DBS while the energy required has been substantially reduced. In a number of studies, the stimulation signal has been delivered in an "on-off " fashion, switching the stimulation on or off depending on whether or not the biomarker exceeded a specified threshold (Little et al., 2013(Little et al., , 2016. Using this approach, suitable stimulation parameters must be first identified, as in the case of open-loop DBS. Using two thresholds, Velisar et al. (2019) increased or decreased stimulation amplitude to maintain the biomarker within a given range. Proportional feedback approaches, where the stimulation amplitude is proportional to the measured biomarker, have also shown potential in both clinical (Rosa et al., 2015;Arlotti et al., 2018) and computational studies (Tukhlina et al., 2007;Chaillet et al., 2017;Popovych and Tass, 2019). Clinical studies in patients have confirmed that proportional stimulation responding to slowly changing beta-band LFP activity is not only effective and well-tolerated by patients, but might also help avoid stimulation-induced dyskinesia when patients are on medication (Rosa et al., 2015;Arlotti et al., 2018). Similar to on-off strategies, which utilize fixed stimulation parameters, proportional feedback requires a fixed controller gain parameter to be identified where the stimulation amplitude at a given time is controlled by this gain and an estimated biomarker value. Identification of the controller gain parameter is a potentially complicated and time-consuming postoperative process. Furthermore, although the selected gain may be suitable for the operating point at which it was initially set, it may provide suboptimal performance and require retuning when properties of the system change, for example in response to disease progression or changes in the properties of the electrode-tissue interface.
To address this problem, we propose a self-tuning control strategy inspired by adaptive control theory, where the value of the controller gain evolves based on the measured pathological activity. The controller gain automatically increases until the detected pathological oscillations are sufficiently suppressed and then begins to dissipate when the oscillation amplitude is low enough. In this manner, the proportional controller self-adapts its gain to the lowest value that guarantees suppression of the oscillations to the desired level. Using tools from control theory, namely Lyapunov-Krasovskii analysis, we mathematically show in a firing-rate model that the proposed controller guarantees disruption of the pathological oscillations provided that the internal coupling within the GPe is sufficiently weak. The controller is derived first under the assumption of endogenously generated oscillations, arising from increased coupling within the STN-GPe network, and is then extended to the case where oscillatory activity arises in the STN-GPe loop due to exogenous cortical inputs.
While firing-rate models provide a means to represent the activity of the network in a mathematically tractable manner, they lack the physiological detail that enables them to be easily related to the underlying processes at the cellular level. To overcome this limitation, we computationally test its performance in a more physiologically relevant context, by implementing it in a network of conductance-based neuron models. We first demonstrate through simulations that the firingrate and conductance-based network models studied here have similar characteristics in terms of the emergence of oscillations as a function of connection strength between STN and GPe, and that both models have similar qualitative frequency responses. Firing-rate model parameters are then identified based on data obtained from simulations of the conductance-based model to link the two models and demonstrate that the latter fulfills the theoretical criterion of low internal GPe connectivity. Finally, the ability of the adaptive controller to "self-tune" to maintain the suppression of pathological oscillations is assessed numerically in three different scenarios: changing the target suppression level, modifying the background beta activity in the network, and varying the electrode impedance.

Firing-Rate Model
The firing-rate model, which we derive the self-tuning DBS strategy and mathematically prove its efficacy, is inspired by the STN-GPe loop model originally proposed in Nevado-Holgado et al. (2010) to study the emergence of pathological beta oscillations observed in the parkinsonian basal ganglia. The model is defined as follows: The instantaneous firing activities (in pulses per second) of the STN and GPe, are respectively represented by x 1 (t) and x 2 (t). τ 1 , τ 2 > 0 are time constants. For each i, j ∈ {1, 2}, the constant c ij ≥ 0 represents the synaptic connection strength from population j to population i and δ ij ≥ 0 is a time delay that occurs due to finite velocity of axonal and synaptic transmission from population j to population i. u 1 and u 2 represent the influence of cortical (into STN) and striatal (into GPe) inputs to the system, respectively, modulated by the synaptic weights b 1 ≥ 0 and b 2 ≥ 0. All the coupling constants c ij and b j being non-negative, the sign represents whether the neurons in the presynaptic population have excitatory (STN and cortex) or inhibitory (GPe and striatum) effect on the postsynaptic one. The activation functions S 1 and S 2 encode the response of the neuronal populations to their input. Although they were considered in Nevado-Holgado et al. (2010) as sigmoids, the theoretical analysis provided below allows them to be any increasing function with bounded derivative.

Self-Tuning Controller Derivation
It was shown in Nevado-Holgado et al. (2010) and Pasillas-Lépine (2013) that if the coupling weights c 12 and c 21 between STN and GPe are sufficiently high, the model exhibits sustained endogenous oscillations, which fall in the beta frequency band for appropriate values of the model parameters. The controller is first derived and assessed under that assumption and is then further explored in the case where oscillations within the loop arise from exogenous inputs to the STN-GPe network. Under the assumption that the cortical and striatal inputs are constant (u(t) = (u 1 (t), u 2 (t)) T =ū), consider the change of variables u ← u −ū, and x ← x −x, wherex is an equilibrium value of x = (x 1 , x 2 ) T for the inputū, whose existence is guaranteed by Pasillas-Lépine (2013). Adding a feedback µ(t) representing the influence of artificial stimulation (DBS) on STN and shifting the activation functions S j such that S j (0) = 0 for j ∈ {1, 2}, we obtain the following dynamics: The model described by (2) has been studied by Pasillas-Lépine et al. (2013) and Haidar et al. (2016), who showed that a proportional feedback acting only on STN is capable of disrupting exaggerated oscillations and stabilizing the system. More precisely, assuming that the inner GPe interconnections are sufficiently weak, namely: where ℓ 2 denotes the maximum slope of the GPe activation function S 2 , it was demonstrated in those papers that the system (2) is asymptotically stable under proportional feedback from the STN to itself, namely: where the proportional gain θ > 0 should be chosen to be sufficiently large. The original result was obtained using linearization techniques. Since then, it has been extended to take full account of the nonlinear effects induced by the activation functions S j (Chaillet et al., , 2019. In particular, it was shown in Chaillet et al. (2019), using Lyapunov-Krasovskii methodology, that the network described in (2) in closed loop with (4) is globally exponentially stable, provided that θ is above some minimum value θ * > 0 and condition (3) is satisfied. Global exponential stability of the origin means that there exist η, γ > 0 such that the solutions of (2) satisfy for any initial state x 0 . The model being a delay differential equation, the state here is not a point in R 2 , but rather a history function defined as x t (s) = x(t + s) for all s ∈ [−δ, 0], where δ denotes the maximum delay δ ij involved in the dynamics. This state (hence, the initial state) belongs to the set of all continuous functions from [−δ, 0] to R 2 and we employ the norm x t = max s∈[−δ,0] |x(t + s)| on this functional set. Global exponential stability is thus a very strong stability property, as it imposes an exponential convergence to the equilibrium and a transient overshoot proportional to the magnitude of the initial state, no matter where the system initially lies. Imposing such a property on the firing-rate model (2) impedes the existence of steadystate pathological oscillations. Global exponential stability is also known to induce robustness properties with respect to exogenous inputs for a wide class of systems (Yeganefar et al., 2008), which may prove useful for the problem considered here due to the inherent imprecision and variability of biological models. One of the major limitations of the stimulation strategy proposed above is that we do not know a priori the value of the minimum effective gain θ * . In Chaillet et al. (2019), the following estimate of θ * was proposed: but it is a conservative approximation. Moreover, this value depends on the connection parameters c ij that would be very hard to estimate accurately in practical applications, due to the high level of abstraction of the considered model. The alternative to pure proportional control that we propose here is inspired from adaptive control theory and involves updating the gain parameter θ based on the measured state of the system, namely: where σ , τ θ > 0 are control parameters governing how fast the control gain reacts to changes in the system. First note that, similar to (4), this control law requires recording of the STN activity only. In the same way, it involves only stimulation of the STN. In practice, these constitute important features in terms of limited insertion of measurement and stimulation electrodes.
The idea behind the adaptive law (6) is simple. As long as x 1 is not at zero (meaning that STN activity has not reached its equilibrium), the term |x 1 (t)| ≥ 0 increases the proportional gain θ . With that, θ (t) will eventually overpass θ * (no matter what its precise value is) and cause the state to converge to zero exponentially. On the other hand, the dissipation (or leakage) term −σ θ (t) decreases the value of θ whenever σ θ (t) ≥ |x 1 (t)|, due to either a too high value reached by the gain θ or because x 1 has reached a sufficiently low value (as desired). These balanced effects are designed in such a way that the control law automatically adjusts its gain around the (unknown) value θ * .
The dissipation term −σ θ (t) present in (6b) is in the spirit of what is known in the control theory literature as the "σ -modification" (Ioannou and Kokotovic, 1984). It was introduced to increase robustness of adaptive control to external disturbances and unmodeled dynamics. It has been shown to guarantee that all the closed-loop signals are bounded and that their mean values converge to a residual set, whose size can be made arbitrarily small with an appropriate choice of σ (Ioannou and Fidan, 2006), even for certain classes of nonlinear systems (Fradkov et al., 1999). Until recently, this methodology was confined to delay-free systems but, for the purpose of the present study, we have extended it to nonlinear time-delay systems (Orłowski, 2019). More precisely, we have the following result (the interested reader is referred to Orłowski, 2019 for more details on the mathematical aspects of this result).
Proposition. Let ℓ i denote the maximum slope of the activation functions S i and letθ 0 = θ 0 − θ * . Under the condition that ℓ 2 c 22 < 1, there exists q > 0 such that, for any σ ≥ 0 small enough and any initial conditions x 0 and θ 0 , the solution of system (2) in closed loop with (6) is bounded and satisfies the following property for all t, T ≥ 0: This statement ensures that solutions are bounded and the system is "stable in the mean." This latter property guarantees that the mean value of the solution, taken over a sufficiently long time window T, converges to a neighborhood proportional to σ , regardless of the initial state. Since σ is a tunable parameter in our controller, this means we can arbitrarily decrease the average amplitude of steady-state oscillations by picking a sufficiently small σ (picking σ as zero would annihilate steadystate oscillations, but would impede the ability to decrease the proportional gain θ whenever possible). The key assumption under which this stabilization is made possible is that ℓ 2 c 22 < 1, meaning that the GPe self-coupling should be reasonably low. More discussion on this assumption and its biological meaning is provided in section 4.2.
In order to selectively attenuate pathological oscillations, with moderate effect on other frequency bands, we propose the following frequency-sensitive version of (6): where β(x 1t ) ≥ 0 is a biomarker detection function. β is implemented as the peak-to-peak amplitude of a bandpassfiltered signal (Butterworth filter, order 5, 15-30 Hz) from t −500 ms to t. This version is similar in spirit to the self-tuning controller (6), but the increase of the gain θ depends only on the STN activity within the targeted frequency band (beta). In a situation when no beta activity is present in the STN, the leakage term −σ θ (t) would cause the gain θ (t) to converge to zero, in which case no DBS signal would be delivered.

Endogenous and Exogenous Generation of Beta Oscillations
The mathematical derivations in section 2.1.1 assume that external inputs into the STN-GPe loop are constant and pathological oscillations arise in an endogenous manner due to too strong synaptic coupling between STN and GPe (Plenz and Kital, 1999;Nevado-Holgado et al., 2010;Pavlides et al., 2012). However, the origin of beta oscillations is still a subject of much debate and there is increasing evidence supporting a role for an exogenous mechanism in which oscillations originating in the cortex or striatum are transmitted to the STN-GPe network and amplified within the network (Magill et al., 2001;Sharott et al., 2005;Mallet et al., 2008;Tachibana et al., 2011;Corbit et al., 2016).
To simulate oscillations arising due to exogenous inputs to the STN, the synaptic coupling between STN and GPe were taken as sufficiently low such that the firing-rate model is not only stable but also incrementally stable  which implies that its steady-state solutions in response to any T-periodic inputs are themselves T-periodic. Thus, any oscillatory input (from cortex or striatum) can entrain the STN-GPe network and cause it to oscillate at the same frequency. This observation enables us to identify the frequency characteristics of the STN-GPe network (despite its nonlinear nature) and which frequency bands, if any, are preferably amplified (see section 3.1.2).
In the firing-rate model, the exogenous and endogenous hypotheses of beta oscillations generation thus correspond to two distinct dynamical behaviors: instability for the former and incremental stability (entrainment) for the latter. Both mechanisms are studied in the simulations presented in this paper (Figures 2-5), using the parameter values presented in Table 1. Additionally, the time constants of the populations were set to τ 1 = 6 ms, τ 2 = 14 ms, and the delays in the system were set to δ 12 = δ 21 = 6 ms and δ 22 = 4 ms. The functions S i were implemented as sigmoids with slope 1 Frontiers in Neuroscience | www.frontiersin.org  where the constants were set to M 1 = 300, M 2 = 400, B 1 = 17, B 2 = 75.

Model and Controller Implementation
In view of its self-tuning capacities, induced robustness and limited requirements in terms of recording and stimulation electrodes, the proposed adaptive closed-loop DBS strategy is promising. Nevertheless, the strong abstraction of the model considered, which summarizes the activity of an entire neuronal population by a unique variable (its firing rate), makes it difficult to assess whether this strategy would be effective clinically. To address this, the performance of the proposed controller was assessed in a biophysically realistic network model of the cortical basal ganglia network comprised of conductance-based neurons. The model is presented in Fleming et al. (2020) and consists of a population of multicompartment cortical pyramidal neurons and single compartment models of cortical interneurons as well as STN, GPe, GPi, and thalamic neurons which have been previously validated and used in other modeling studies (Terman et al., 2002;Otsuka et al., 2004;Rubin and Terman, 2004;Pospischil et al., 2008;Hahn and McIntyre, 2010;Foust et al., 2011;Lowery, 2013, 2014;Kumaravelu et al., 2016). Each population was comprised of 100 neurons, where synaptic connections between neurons were modeled by spike detectors in presynaptic neurons coupled to synapses in postsynaptic neurons by a time delay. Synaptic connections in the model were either excitatory (AMPAergic) or inhibitory (GABAergic) depending on the synapse type (Destexhe et al., 1994). Striatal input to the network was represented as Poisson-distributed spike trains to GPe neurons with a mean firing rate of 3 Hz. An overview of the model structure is presented in Figure 1.
The model captures key features of the cortical basal ganglia network required for simulating clinical implementations of closed-loop DBS including: (i) the extracellular DBS electric field, which is required to accurately model changes in the DBS amplitude, (ii) antidromic and orthodromic activation of STN afferent fibers, and (iii) the STN LFP detected at nonstimulating contacts of the DBS electrode. In Fleming et al. (2020), the model parameters were tuned to match key features observed during experimental investigations of DBS including cortical desynchronization (Li et al., 2012), GPe entrainment (McConnell et al., 2012), and a gradual suppression of betaband power detected in the STN LFP for increasing stimulation amplitude (Davidson et al., 2016). The model is described in full detail in its original publication (Fleming et al., 2020) and available to download from ModelDB (https://senselab.med.yale. edu/modeldb/) at ascension number 262046.
The oscillatory properties of the model's STN-GPe network were first examined to explore the duality between its oscillatory behavior and that of the firing-rate model. Beta-band activity was then configured to remain fixed, or was varied according to the three numerical scenarios detailed below which may require controller adaptation in vivo.
In line with (7), the proposed self-tuning controller was implemented in the conductance-based model to adapt the amplitude of the stimulation waveform as follows: where µ(t) is the controller output and represents the instantaneous stimulation amplitude, θ (t) is the controller gain, x 1 (t) is the biomarker measurement [i.e., the average rectified value (ARV) of a 100 ms epoch from the beta-band filtered STN LFP, which was filtered using a fourth order Chebyshev band-pass filter with an 8 Hz bandwidth, centered on 25 Hz as described in Fleming et al., 2020], |e(t)| is the half-wave rectified error signal calculated as the difference between the measured biomarker and the desired target suppression level, and τ θ and σ represent tuning parameters which were fixed at 100 ms and 0.00875, respectively, for all controller simulations in the conductance-based model.
A key difference with (7) is that the DBS in (9) is delivered as a positive feedback on the biomarker [as indicated by the positive sign in µ(t)]. This difference is due to the implementation of DBS in the conductance-based model, whereby increasing DBS amplitude results in a stronger suppression of pathological oscillations.
The model was simulated in the NEURON simulation environment (Hines and Carnevale, 1997) and implemented in Python using the PyNN API package (Davison et al., 2009). The model was numerically integrated using the Crank-Nicholson method with a 0.01 ms timestep for all simulations. Simulations were run on the UCD Sonic high-performance computing cluster.

Numerical Scenarios
The self-tuning controller (9) was tested in three independent scenarios to simulate practical situations in which adaptation of controller parameters in vivo may be required to maintain the biomarker (the ARV of LFP beta activity) at a target value. All scenarios were simulated for a 130 s duration.
In the first scenario, background beta-band activity in the model was set to its maximum value for the duration of the simulation while the target level for beta suppression was varied. The beta-band activity (prior to DBS activation) was set by fixing the firing rates of cortical neurons to 26 pulses per second, which resulted in a peak in the LFP power spectrum at 26 Hz. Target values of 10, 0.2, 0.05, 0.15, and 10 µV were then considered over the time intervals 0-10, 10-40, 40-70, 70-100, and 100-130 s, respectively. In the second scenario, the target value for the LFP biomarker was fixed at 0.1 µV for the duration of the simulation while background beta-band activity in the model was modulated to be low during the 0-10, 40-70, and 100-130 s time intervals and high during the 10-40 and 70-100 s time intervals. The intracellular cortical neuron bias current was varied to shift the mean cortical neuron firing rate between 14 and 26 pulses per second during the low and high beta-band activity periods, respectively. Thus, beta-band activity in the model was modulated to display a 14 Hz peak in the LFP power spectrum during the low beta-band activity periods, which shifted to 26 Hz during the high beta-band activity periods. As the bandwidth of the biomarker filter was centered at 25 Hz, modulation of the background beta-band activity in this manner led to lower beta ARV measurements during the low beta-band activity periods.
The third scenario considered a linear variation of the electrode impedance over the simulation period, while background beta-band activity in the model was fixed at its maximum value and the beta ARV target value remained constant at 0.1 µV. In the simulation, the electrode impedance remained constant at 0.5 k up to t = 30 s, after which it was linearly increased to a maximum value of 2.5 k at t = 130 s.

Performance Measures
Controller performance was quantified using two measures: the error while tracking the target value and the mean power consumption. The mean squared error (MSE) was utilized to measure the controller's ability to track the target level. It is defined as where T sim is the simulation duration (T sim = 130 s) and e(t) is the normalized error signal between the measured LFP beta ARV and the target value, as used in (9). For simplicity, the MSE value for the controllers in each scenario are reported as a percentage of the MSE value that was measured in each respective scenario when DBS was off. Power consumption (PC) was measured as where Z E is the electrode impedance, assumed to be 0.5 k in the non-varying electrode impedance scenarios, and I DBS is the delivered DBS current.

Parameter Identification
One of the necessary conditions for practical stabilization using the adaptive controller, as recalled in section 2.1.1, is that the internal connections within GPe are weak, as expressed by c 22 ℓ 2 < 1 in the firing-rate model. It was, therefore, first checked whether this condition was fulfilled in the conductance-based model. The value of c 22 from the firing-rate model is related to the maximum conductance of the GPe-GPe synapsesḡ GPeGPe , in the conductance-based model, which was set to 0.015 µS. ℓ 2 is the maximum slope of the GPe activation function S 2 . Thus, in order to verify the stabilizability criterion, we need to identify the slope of the GPe activation function in the conductance-based model, as well as the value of c 22 corresponding to the value ofḡ GPeGPe used in the conductance-based model.
To that aim, the GPe neurons of the conductance-based model were disconnected from the rest of the network, the only remaining connections being the excitatory projections from STN to GPe. In the firing-rate model, this translates in the following dynamics: where u(t) represents synaptic inputs to GPe from STN. Settinḡ g GPeGPe = 0 in the conductance-based model leads to c 22 = 0 in the firing-rate model, thus yielding We conducted a series of simulations for different values of cortical input to STN and we estimated the firing rate u(t) of STN and the firing rate x 2 (t) of GPe. We have rescaled the data to obtain a slope 1, by fitting a linear function and dividing the values by the obtained slope a = 1.29. We then fitted a sigmoid of the form (8) to the obtained data and rescaled it by the same factor a, to obtain an estimate of the activation function S 2 . Next, we obtained a similar set of data forḡ GPeGPe = 0.015 µS. Using the activation function S 2 determined in the previous step, we found the equilibrium of (12) for different values of c 22 and constant STN input u using a numerical solver. We compared the curves obtained for the different values of c 22 with the steadystate data from the conductance-based model, to identify the value that minimized the normalized square error: where STN[i] and GPe[i] represent the firing rate of STN and GPe taken from simulation i, and f c 22 (STN[i]) is the solution of for a given c 22 with u = STN[i], meaning the steady-state solution of (12) for these generated inputs from STN. The best fit was reached for c 22 = 0.35. Since the GPe activation S 2 identified in (8) has maximum slope ℓ 2 = 1.29, the stabilizability criterion (3) is satisfied.

Endogenous Oscillations
For constant striatal and cortical inputs, beta-band oscillations emerged in the firing-rate model when the STN-GPe and GPe-STN connectivity strength was sufficiently increased (Figure 2A). In the conductance-based model, the intensity of beta power in the spectra of the cumulative STN and GPe population spike trains similarly increased with increasing STN-GPe and GPe-STN connectivity strengths ( Figure 2B). The increase in beta power within the STN and GPe was accompanied by an increase in synchronization of the two populations within the beta band. The beta-band coherence of the STN and GPe neural spike trains increased similarly with increased connectivity strengths (Supplementary Figure 1).
Both models exhibited beta-band oscillations as STN-GPe and GPe-STN connectivity increased, indicating that the firingrate abstraction well captures the oscillatory dynamics of the conductance-based model. Differences between the two models are, however, observed. In particular, in the firingrate model, the transition from non-oscillatory to oscillatory behavior occurs more abruptly than in the conductance-based model and oscillatory conditions are associated with a relatively stronger coupling from STN to GPe than from GPe to STN. Also, in the conductance-based model, fluctuations in the frequency and amplitude of the beta oscillations are apparent as GPe-STN connectivity increases. Nevertheless, the overall behavior of the models is qualitatively similar with oscillations emerging in both models as STN-GPe connectivity and GPe-STN connectivity increase.

Exogenous Oscillations
When synaptic coupling between STN and GPe is sufficiently low, the firing-rate model of the STN-GPe loop is not only stable but also entrainable, meaning that any T-periodic input (whether cortical or striatal) generates T-periodic steady-state solutions . While this feature is guaranteed for stable linear systems, the nonlinear nature of the firingrate model makes it less straightforward. This entrainability is a fundamental requirement for constructing frequency profiles of the STN-GPe network. By considering a sinusoidal input at a given frequency, it is indeed possible to measure the magnitude of the resulting steady-state oscillations, and thus the amplification of the network at this specific frequency (Pavlov et al., 2007). Repeating this procedure across a range of input frequencies, we obtained the nonlinear Bode plots depicted as solid lines in Figure 3A.
When DBS is off (solid curves), a clear resonance is observed in the beta frequency band, thus indicating that the network preferably amplifies beta components of the cortical input. This resonance can therefore be interpreted as the beta generation mechanism in the exogenous hypothesis. Due to the nonlinear nature of the firing-rate model, this resonance strength depends on the amplitude of the applied cortical input with more pronounced resonance occurring for stronger mean cortical input.
Akin to the resonance behavior in the firing-rate model, the influence of an external cortical oscillatory input to the STN-GPe network was investigated in the conductance-based model. The connectivity strengths between the STN and GPe populations were selected to lie below the threshold for which the network generated endogenous beta-band oscillations (at 0.11 µS for both), analogous to the entrainable state in the firing-rate model. The frequency of cortical inputs to the STN were varied from 3 to 100 pulses per second, while striatal input to the network remained fixed at 3 pulses per second. The frequency response of the STN-GPe network due to synchronous cortical inputs through the hyperdirect pathway was examined by estimating the power of the cumulative population spike trains and the spike train coherence between STN and GPe populations within a 4 Hz window centered on the mean frequency of the cortical input. Additionally, resonant network activity was calculated by estimating the power in the simulated STN LFP at the cortical input frequency. Resonance effects were examined as the strength of cortical connectivity to the STN was systematically increased from 0.03 to 0.12 µS ( Figure 3B).
As the strength of the hyperdirect pathway was increased, a beta-band resonance emerged in the STN and GPe populations and in the power spectrum of the STN LFP. This was accompanied by synchronization across the STN and GPe populations, as evidenced by a peak in the coherence between STN and GPe spike trains for cortical inputs in the beta frequency range (Figure 3B). Further strengthening of the hyperdirect pathway led to a broadening of the frequency band at which resonance occurred in the cumulative spike trains of the STN and GPe populations and the STN LFP, extending beyond the beta-band. Synchronous cortical inputs to the STN at low connectivity strengths in the beta band and at frequencies outside this range resulted in coherent activity in the subnetwork, but with relatively low power ( Figure 3B). Consistent with the firing-rate model (Figure 3A), the beta-band resonance was more pronounced for stronger cortical inputs. FIGURE 3 | Resonance behavior of the firing rate and conductance-based models in response to synchronous cortical drives. (A) Bode plot illustrating amplification of cortical input signals of varying frequency by the (i) STN and (ii) GPe populations in the firing rate model. The amplitude ratio is defined as the amplitude of the input oscillation (the cortical input signal) to the amplitude of the steady state oscillation in the (i) STN and (ii) GPe populations at each cortical input frequency. Solid lines represent the frequency response of each population when DBS is off. Dashed lines represent the frequency response of each population when the self-tuning DBS (7) is implemented with σ = 0.1 and τ θ = 50 ms. Two mean values of the cortical input signal are represented. (B) Resonance plot of the power centered at the cortical input frequency in the (i) STN and (ii) GPe cumulative spike trains, (iii) the STN LFP power spectrum and (iv) the coherence between the STN and GPe populations for varying cortical input frequencies. Power at the cortical input frequency was estimated in the cumulative spike trains and LFP power spectra by integrating the power in a 4 Hz window centered on the input frequency in the respective power spectra. The coherence between the STN and GPe populations at the cortical input frequency was estimated from pairs of composite spike trains randomly chosen from the STN and GPe populations (Farina et al., 2014;McManus et al., 2019). The spike trains in the STN and GPe were summed to obtain two composite spike trains. The magnitude squared coherence between the two composite spike trains was then calculated with 1 s overlapping Hamming windows. This was repeated for 200 randomly chosen combinations of spike trains from the STN and GPe populations, as each combination will generate a slightly different coherence estimate. The coherence between the populations was then estimated as the median coherence spectrum over all 200 combinations. The coherence at the cortical input frequency was determined by integrating the resulting coherence spectrum in the 4 Hz window centered on the input frequency. Four cortical to STN connectivity strengths are represented, where the resonance responses are normalized between 0 and 1 in each panel separately.

Endogenous oscillations with beta-band activity variation
Suppression of endogenous beta oscillations during DBS was assessed first in the firing-rate model, with synaptic weights c 12 and c 21 between STN and GPe increased to generate beta oscillations within the network. A performance comparison between self-tuning DBS (6) and proportional DBS (4) is presented in Figure 4. Both strategies successfully disrupt pathological oscillations between 200 and 750 ms. At t = 750 ms, an additional increase of the cortical input to STN was artificially FIGURE 4 | Performance comparison of self-tuning and proportional controllers on endogenous oscillations in the firing-rate model. During simulation, DBS was initially off and then switched on at t = 200 ms, where either the (A) self-tuning controller or (B) proportional controller was implemented. At t = 750 ms, the magnitude of the constant cortical input to STN was artificially increased by from 27 to 42, amplifying the endogenous oscillations present in the STN and GPe. The cortical input is presented in green at the top of the figure, aligned with respect to the simulation time with the plots underneath. Left panels illustrate the firing rate behavior of the STN and GPe populations and the right panels correspond to the instantaneous peak-to-peak oscillation amplitude of each population (measured with a sliding window of 100 ms) due to either (A) the self-tuning controller (7) with τ θ = 75 ms and σ = 0.19 or (B) the proportional controller (4) with θ = 2.
introduced to simulate an increase in beta oscillations. While the self-tuning DBS automatically adapts the proportional gain θ to maintain the attenuation of beta oscillations, pure proportional DBS is unable to do so, resulting in strong beta oscillations that cannot be counteracted without manual tuning of the proportional gain. The self-tuning controller thus outperforms the proportional controller in terms of robustness to disease evolution in the simulated endogenous mechanism scenario.

Exogenous oscillations with beta-band activity variation
Selective disruption of pathological beta oscillations generated through exogenous inputs to the STN using the self-tuning controller (7) was then confirmed. With the self-tuning DBS (dashed curves) ( Figure 3A) the beta-band resonance was eliminated in the firing-rate model, while the frequency profile of the STN-GPe network remained essentially unaltered in other frequency bands, as DBS remains off when beta activity is not detected within the STN. Figure 5 illustrates that, similar to the endogenous case, self-tuning DBS (6) outperforms proportional DBS (4) when faced with changes in the exogenous oscillations. After the stimulation is turned on at t = 200 ms, both controllers successfully decrease the amplitude of the pathological oscillations. After the mean level of the cortical input as well as the amplitude of oscillations is increased at t = 750 ms, the self-tuning controller achieves higher damping of the oscillations than the proportional controller.

Scenario 1: beta-band target variation
The self-tuning controller maintained the LFP beta ARV around the desired target values between 10 and 100 s, At t = 750 ms, the amplitude of the cortical oscillations was increased from 10 to 60 and the mean level of the oscillations is raised from 50 to 60, amplifying the exogenous oscillations present in the STN and GPe. The cortical input is presented in green at the top of the figure, aligned with respect to the simulation time with the plots underneath. Left panels illustrate the firing rate behavior of the STN and GPe populations and the right panels correspond to the instantaneous peak-to-peak oscillation amplitude of each population (measured with a sliding window of 100 ms) due to either (A) the self-tuning controller (7) with τ θ = 5 ms and σ = 0.01 or (B) the proportional controller (4) with θ = 25. and turned off from 0 to 10 and 100 to 130 s, where the target value was high and stimulation was not required (Figure 6A). The self-tuning controller resulted in a mean power consumption of 11.0 µW and a 82.1% reduction in the MSE. Proportional controllers with fixed gain were able to maintain beta ARV at a single target value, however they were unsuitable for target values other than the one for which they were tuned and resulted in either under or over stimulation when attempting to maintain beta ARV (Figures 6B-D). Proportional controllers with fixed gains at 7, 10 and 26 resulted in mean power consumption values of 6.9, 9.3, and 28.2 µW and reductions in the MSE of 79.4, 85.1, and 81.1%, respectively.

Scenario 2: beta-band activity variation
The LFP beta ARV was maintained at the target value by the self-tuning controller as the background beta activity was varied between low and high activity periods ( Figure 7A). The selftuning controller consumed 17.9 µW and resulted in a 91.4 % reduction of the MSE. The proportional controller with a fixed gain value of 10 was able to maintain the LFP beta ARV at the target value during low background beta activity periods, but did not provide sufficient stimulation to suppress to the target value during periods of high background beta activity ( Figure 7B). The proportional controller with a fixed gain value of 40 maintained the LFP beta ARV at the target value during both low and high background beta activity periods ( Figure 7C). The proportional controllers with fixed gain values of 10 and 40 resulted in MSE reductions of 73.4 and 95.0 % and power consumption values of 6.6 and 34.8 µW, respectively. When maintaining the LFP beta activity at the target value, the proportional controller with a fixed gain value of 40 resulted in over stimulation during the low background activity periods, consuming more power than necessary to maintain the LFP beta ARV. FIGURE 7 | Performance comparison of self-tuning and proportional controllers in response to background beta-band activity variations in the conductance-based model. The beta ARV target value was fixed to 0.1µV for the duration of the simulation while the background beta-band activity was modulated to be at its maximum value during the 10-40 and 70-100 s time periods and at its minimum value during the 0-10, 40-70, and 100-130 s time periods, respectively. The mean ARV of the background beta-band activity measured in the LFP when DBS was off is highlighted by the gray segmented bar at the top of the figure. Left panels illustrate the target value (dashed red), the beta ARV measured from the STN LFP when DBS was off (gray) or on (black), where the DBS amplitude was modulated by the corresponding controller. Right panels illustrate the time evolution of the DBS amplitude (black) and controller gain θ (blue) during simulation. (A) Self-tuning DBS controller with τ θ = 100 ms and σ = 0.00875. (B) Proportional controller with a fixed gain value of θ = 10. (C) Proportional controller with a fixed gain value of θ = 40.

Scenario 3: electrode impedance variation
The self-tuning controller maintained the LFP beta ARV at the target level as the electrode impedance was linearly increased over the course of the simulation to five times its initial impedance value, i.e., from an initial value of 0.5-2.5 k ( Figure 8A). Background beta activity and the target value of 0.1 µV remained fixed over the course of the simulation. The self-tuning controller resulted in a MSE reduction of 93.5% and a power consumption of 11.2 µW. Proportional controllers with fixed gain values of 10 and 40 lead to MSE reductions of 77.0 and 96.0 % while consuming 6.2 and 27.7 µW, respectively. Similar to scenarios 1 and 2, the self-tuning controller was able to tune its gain to the required level to maintain beta ARV at the target level. Proportional control with fixed gain of 10 became less effective over the course of the simulation, while proportional control with fixed gain of 40 was effective throughout the simulation, but consumed more power than necessary (Figures 8B,C). A summary of the controller performance under the different scenarios considered is presented in Table 2.
FIGURE 8 | Performance comparison of self-tuning and proportional controllers in response to electrode impedance variations in the conductance-based model. The beta ARV target value was fixed to 0.1µV and the background beta-band activity was fixed to its maximum value for the duration of the simulation. The electrode impedance was fixed at 0.5 K during the 0-30 s period and was then linearly increased from 0.5 to 2.5 K over the 30-130 simulation period. The electrode impedance variation is illustrated by the purple bar at the top of the figure. Left panels illustrate the target value (dashed red), the beta ARV measured from the STN LFP when DBS was off (gray) or on (black), where the DBS amplitude was modulated by the corresponding controller. Right panels illustrate the time evolution of the DBS amplitude (black) and controller gain θ (blue) during simulation. (A) Self-tuning DBS controller with τ θ = 100 ms and σ = 0.00875. (B) Proportional controller with a fixed gain value of θ = 10. (C) Proportional controller with a fixed gain value of θ = 40.

Firing-Rate and Conductance-Based Models
The proposed firing-rate model facilitated the derivation of a robust control law capable of disrupting pathological beta-band oscillations in the STN-GPe network in Parkinson's disease and the analytic establishment of its efficiency. Firing-rate models capture the average behavior of neural populations and facilitate tractable mathematical analysis of network behavior (Destexhe and Sejnowski, 2009). These models, however, summarize the neuronal population to a single variable: the number of spikes it emits per unit time. They are thus unable to capture cellular-level features such as sub-threshold activity, LFP activity, specific responses induced by bursting, antidromic activation of STN afferent inputs or the interaction between neurons and DBS Off 100 0 Self-tuning controller 7.9 11.0 Proportional controller (θ = 7) 20.6 6.9 Proportional controller (θ = 10) 14.9 9.3 Proportional controller (θ = 26) 18.9 28.2 Beta variation: MSE (%) Power consumed (µW) DBS Off 100 0 Self-tuning controller 8.6 17.9 Proportional controller (θ = 10) 26.6 6.6 Proportional controller (θ = 40) 5.0 34.8 Electrode impedance variation: MSE (%) Power consumed (µW) DBS Off 100 0 Self-tuning controller 6.5 11.2 Proportional controller (θ = 10) 23.0 6.2 Proportional controller (θ = 40) 4.0 27.7 In each scenario, the MSE value is reported as a percentage of the DBS Off case, where the DBS Off case corresponds to a 100% MSE.
the induced extracellular potential. In contrast, networks of conductance-based neuron models offer the ability to capture more complex network dynamics and interactions, but are usually too complex to analyze mathematically. While recent approaches have attempted to bridge the gap between these two types of models (di Volo et al., 2019), we have empirically assessed the similarities and discrepancies between them and shown that key features needed to analyze and control pathological oscillations in Parkinson's disease are indeed captured by both models (Figures 2, 3). We also demonstrated in section 2.3 the possibility to check the theoretical stabilization condition (3), derived on the firing-rate model, using observations from the conductance-based model. By utilizing both modeling approaches, the limitations of each model are complemented by the strengths of the other. The joint analysis of these two models also supports the relevance of abstract firing-rate models for the derivation of advanced DBS strategies aiming at counteracting a targeted brain oscillation, which can then be validated in computationally detailed models before preclinical investigations.

Physiological Interpretation of the Stabilizability Condition
The theoretical condition obtained on the firing-rate model to ensure stabilizability by the self-tuning DBS signal reads ℓ 2 c 22 < 1 (see Proposition). In other words, the synaptic weights from GPe to itself should be sufficiently low. This condition ensures that the GPe does not act as a pacemaker on its own, as low internal coupling is a standard sufficient condition for stability of a neuronal population (see for instance Faye and Faugeras, 2010). The necessity to impose that GPe does not generate pathological oscillations on its own is quite reasonable: considering the extreme case when STN is not connected to the GPe, it would be impossible to attenuate self-generated GPe oscillations by stimulating STN only.
It is worth noting that the proposed condition precludes GPe self-oscillatory activity no matter the value of its internal delay. This constitutes a noteworthy feature of our mathematical result as this delay does not need to be estimated. Nevertheless, for realistic values of internal GPe delays (of the order of few milliseconds), no such self-oscillatory GPe activity is observed even if the condition is violated, at least for reasonable values of the striatal input. This was confirmed in numerical simulations of the conductance-based model in which the internal GPe synaptic weights were artificially increased by two orders of magnitude. Even in that case, GPe was unable to autonomously generate oscillations and the proposed self-tuning DBS successfully disrupted network beta oscillations (data not shown).

Endogenous and Exogenous Generation of Oscillations
Beyond employing two modeling approaches of the structures involved, the paper also investigated two possible mechanisms of pathological oscillations generation: the emergence of endogenous oscillations, in which the STN-GPe network acts as a pacemaker, and the generation of oscillations through the interaction of the network with inputs originating from other structures such as the cortex.
The ability of the STN-GPe network to endogenously generate beta-band oscillations was consistent across both the firingrate and conductance-based models (Figure 2). The STN-GPe network has been proposed as a potential source of pathologically increased beta-band oscillations in Parkinson's disease, where connectivity changes in the reciprocally connected network leads to the endogenous generation of beta-band oscillations. This behavior has been previously explored in modeling studies utilizing firing-rate models where the progression of Parkinson's disease is represented as an increase in the synaptic coupling strengths between the STN and GPe neuron populations (Nevado-Holgado et al., 2010;Pavlides et al., 2012;Pasillas-Lépine, 2013). The firing-rate model presented here is consistent with these previous studies in which a Hopf bifurcation occurs and leads to beta-band oscillations in the network when the synaptic coupling strengths are sufficiently increased (Figure 2A). This behavior is well-matched in the conductance model, where increases in the synaptic coupling strengths between the two populations also leads to the endogenous emergence of beta-band oscillations in the network (Figure 2B).
Although modeling studies support the hypothesis that the reciprocally connected STN-GPe network is capable of generating the beta-band oscillatory activity, investigations of isolated STN-GPe cell cultures in vitro have observed the emergence of endogenous oscillations at much lower frequencies (Plenz and Kital, 1999). Furthermore, increasing evidence from experimental studies in patients and animal models suggest that external inputs to the STN-GPe loop may play a key role in the generation of elevated beta-band oscillations in the parkinsonian cortex and basal ganglia. The STN-GPe network occupies a crucial location in the cortical basal ganglia network receiving inputs from both cortical and striatal structures through the hyperdirect and indirect pathways, respectively. Due to connectivity changes in the STN-GPe loop during Parkinson's disease, it is thus hypothesized that exogenous beta oscillations are locally amplified by the loop, with subsequent oscillations then propagating throughout the full cortical basalganglia network (Magill et al., 2001;Mallet et al., 2008;Corbit et al., 2016;West et al., 2018). Consistent with this, Sharott et al. (2005) observed that the power and coherence of beta oscillations in the cortex and STN were elevated during dopamine depletion in a parkinsonian rat model, while Litvak et al. (2011) observed cortical beta activity leading to STN beta activity in parkinsonian patients. In primate studies, oscillatory activity in the GPe was observed to be mainly due to excitatory activity from the STN, while oscillatory activity in the STN was primarily due to excitatory cortical input (Tachibana et al., 2008(Tachibana et al., , 2011. Both the firing-rate and conductance-based models in the present study demonstrated the STN-GPe network's ability to entrain to an external cortical rhythm (Figure 3). A Bode plot of the input-output relationship in the firing-rate model showed a marked peak in the beta frequency band, with this peak increasing in magnitude and broadening with increasing strength of the inputs to the STN ( Figure 3A). This behavior was consistent with the conductance-based model behavior, where the STN-GPe network displayed a resonant peak in the beta frequency band for low connectivity strength between the cortex and STN (Figure 3B). Increasing the strength of the hyperdirect pathway led to increased resonance in the model and also to a widening of the frequency bands where resonance was observed ( Figure 3B). This observed resonance is consistent with other modeling investigations of the STN-GPe network, where firing-rate models (Nevado-Holgado et al., 2014;Detorakis and Chaillet, 2017;Liu et al., 2017Liu et al., , 2020 and conductancebased neuron models (Ahn et al., 2016;Shouno et al., 2017;Koelman and Lowery, 2019) were used to investigate the behavior of STN-GPe network in response to external drives. Although those studies have illustrated the resonant capabilities of the STN-GPe network of both firing-rate and conductance-based models separately, this study is the first to demonstrate that both modeling approaches lead to comparable results (Figure 3).

Self-Tuning DBS Controller
Having established qualitative consistency in the behavior of the firing-rate and conductance-based models, the performance of the proposed self-tuning DBS controller was first proven mathematically and validated in the firing-rate model where the controller was capable of disrupting both exogenously and endogenously generated beta-band activity in the STN-GPe network (Figures 3A, 4, 5). The simulations confirmed that the self-tuning controller autonomously adapts its gain value to the minimal value required to counteract pathological oscillations, thus avoiding over-stimulation and allowing for adaptation to possible changes in the system properties associated with disease progression (Figures 4, 5).
The performance of the controller to maintain network beta-band oscillations at a target level was then assessed in the conductance-based model in three example conditions, which emulated practical situations in which gain adaptation may be required in vivo: modification of the beta-level target, variation of the beta oscillations intensity, and alteration of the electrode impedance.

Adaptation to Target-Level Changes
The self-tuning controller adapted the controller gain in response to changes in the target value ( Figure 6A). For each target value, the controller identified the necessary gain to maintain the biomarker at the target level. In contrast, proportional controllers with fixed gain were unable to track target changes and resulted in less reduction in the MSE than the self-tuning controller (Figures 6B-D, Table 2). The self-tuning controller consumed more power than the controllers with low fixed gain and less power than the controller with high fixed gain, but was able to maintain low error as the target changed ( Table 2).
Similar issues were identified by Su et al. (2019) who investigated the ability of a proportional-integral controller to modulate DBS frequency to track dynamic changes in a target signal during closed-loop DBS using a conductance-based model. While the proportional-integral controller considered in that study was able to successfully track dynamic changes in the target beta signal, it required different controller gain values for each beta-band target level considered. The adaptive controller proposed here overcomes this issue as the controller self-tunes its gain value to find the gain necessary for maintaining the beta ARV at the target values ( Figure 6A).

Adaptation to Beta Oscillation Fluctuations
The self-tuning controller was able to maintain the biomarker at the target level while beta activity in the network varied ( Figure 7A). Between t = 0 and t = 10 s, the self-tuning controller adapted its gain value to the low beta activity. Once beta activity in the network increased, during the 10-40 s period, the controller increased its gain to suppress the beta activity to the target level. The proportional controller with low fixed gain value (θ = 10) was able to maintain the biomarker at the target level during low beta activity periods, but was unable to during high beta periods ( Figure 7B). In contrast, the proportional controller with a high fixed gain value (θ = 40) was able to maintain the biomarker close to the target during both low and high beta activity periods (Figure 7C), but at the cost of increased power consumption ( Table 2).
In line with the power consumption in scenario 1, the controllers with low and high fixed gain resulted in the lowest and highest power consumption, respectively, and the power consumption of the self-tuning controller lay between these values (Table 2). Essentially, the self-tuning controller identified the gain required to maintain the beta ARV at the target level in both the high and low beta activity periods and consumed the necessary power to maintain beta at the target.
Clinical investigations of proportional control strategies in patients with Parkinson's disease (Rosa et al., 2015;Arlotti et al., 2018), have observed attenuation of DBS when subjects were on and off medication, motivating the need for self-tuning in response to varying background beta activity. Arlotti et al. (2018) investigated proportional control over an 8 h period and showed a 30-45 % improvement in patient UPDRS III motor scores. Although these clinical studies showed promising results, we hypothesize that additional benefits may result from using the self-tuning controller proposed here. The proportional control scheme implemented in Rosa et al. (2015) and Arlotti et al. (2018) utilized a fixed gain value, and thus may lead to under or over stimulation if a suitable gain is not selected.

Adaptation to Electrode Impedance Variations
Clinical measurements of electrode impedance usually vary between 0.5 and 1.5 k (Obeso et al., 2001;Volkmann et al., 2002). However, many factors contribute to the electrode impedance, and its variability between subjects. These factors include the surface properties of the electrodes, electrical double layer, conductivity of the bulk tissue medium and thickness of the surrounding encapsulation layer (Butson et al., 2006). The electrode impedance plays a crucial role in determining the current delivered to tissue during voltage-controlled stimulation. Clinical investigations of closed-loop DBS have not yet investigated the impact of electrode impedance variations on DBS efficiency as the timescales at which these variations take place are much longer than the current timescales at which closed-loop DBS has been investigated clinically. For closedloop DBS to remain effective when chronically implemented, it is necessary to utilize controllers which can accommodate such changes and avoid a need for clinical retuning of the controller parameters. This third scenario, therefore, aimed to assess this robustness to electrode impedance variations. In response to the simulated changes, the self-tuning controller gradually increased the controller gain to maintain the biomarker at the target level ( Figure 8A).
In contrast, the proportional controller with low fixed gain was able to maintain the beta ARV at the target up until t = 60 s, but became less effective thereafter (Figure 8B). With a higher fixed gain value, the proportional controller was able to accommodate changes in the electrode impedance and remained effective at suppressing the beta ARV to the target for whole the duration of the simulation ( Figure 8C).
In line with the other examples, the power consumption of the self-tuning controller was greater than that of the proportional controller with low fixed gain and less than the proportional controller with the high fixed gain ( Table 2). Although both the proportional controller with high fixed gain and the self-tuning controller successfully maintained beta activity at its target level, more power consumption was needed for the former than for the latter.

Clinical Implementation of the Self-Tuning DBS
The self-tuning controller presented here offers several advantages for experimental implementation. First, it relies only on data from the STN, no additional recording electrodes in other brain structures are required. Moreover, although the firing-rate model utilized the number of STN pulses per second, the simulations conducted in the conductancebased model showed the efficiency of the approach when only STN LFP is accessible. More importantly, the self-tuning controller relies on very limited information regarding the system parameters. The theoretical result of Proposition simply requires that the GPe internal synaptic weights are sufficiently low. This condition, also present in theoretical investigations on fixed-gain proportional DBS , does not require precise knowledge on the exact shape of activation functions, or the values of the time constants, synaptic weights, or delays.
Nonetheless, experimental validation of the self-tuning controller also comes with challenges. First, both the stimulation and the recording are assumed to be within the STN. This would lead to stimulation artifacts that should be removed. Several techniques are available to address this issue (Rossi et al., 2007;Stanslaski et al., 2012;Basir-Kazeruni, 2017) The proposed controller also requires some embedded computational capabilities in order to filter the STN LFP and implement the control law (7). While more demanding than classical open-loop DBS, the required computational power is comparable to the algorithms employed for proportional or on-off stimulation. In addition, the upper and lower bounds on the applied current or voltage would be the same as those used during conventional stimulus parameter setting.
Despite its self-tuning nature, the control law (7) requires two parameters to be chosen: the time constant τ θ and the dissipation parameter σ . Both these parameters have a strong impact on the beta attenuation performance. For the sake of illustration, these parameters have been chosen here to demonstrate quick adaptation of the gain over the simulation duration. In an experimental setting, τ θ should be chosen on the order of the timescale at which the targeted variations occur. Generally speaking, a large value of τ θ is more convenient for adaptation to slow variations. σ should be chosen as a compromise between rapid decrease of the gain when the beta level is low (high σ ) and limited gain overshoot in response to an increase of the beta level (low σ ). In practice, the controller parameters could also be obtained through an optimization process, such as the dual-loop framework outlined in Grado et al. (2018), though an optimization approach would itself also require the identification of appropriate objective functions which may be non-trivial. In summary, the proposed self-tuning DBS controller offers increased robustness to variations in system properties including changes in the strength of beta oscillations and or the electrode-tissue interface, and can accommodate alterations in the desired beta level. Compared to fixedgain proportional control, it provides a better compromise between power consumption and efficient attenuation of beta oscillations, with little additional computational complexity or technological requirements.

DATA AVAILABILITY STATEMENT
The models presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: The firing-rate model is available to download from GitHub (https://github.com/ silvanx/self-tuning-control). The conductance-based model is available to download from ModelDB (https://senselab.med.yale. edu/modeldb/) at ascension number 262046.