Multivariable closed-loop control of deep brain stimulation for Parkinson’s disease

Objective. Closed-loop deep brain stimulation (DBS) methods for Parkinson’s disease (PD) to-date modulate either stimulation amplitude or frequency to control a single biomarker. While good performance has been demonstrated for symptoms that are correlated with the chosen biomarker, suboptimal regulation can occur for uncorrelated symptoms or when the relationship between biomarker and symptom varies. Control of stimulation-induced side-effects is typically not considered. Approach. A multivariable control architecture is presented to selectively target suppression of either tremor or subthalamic nucleus beta band oscillations. DBS pulse amplitude and duration are modulated to maintain amplitude below a threshold and avoid stimulation of distal large diameter axons associated with stimulation-induced side effects. A supervisor selects between a bank of controllers which modulate DBS pulse amplitude to control rest tremor or beta activity depending on the level of muscle electromyographic (EMG) activity detected. A secondary controller limits pulse amplitude and modulates pulse duration to target smaller diameter axons lying close to the electrode. The control architecture was investigated in a computational model of the PD motor network which simulated the cortico-basal ganglia network, motoneuron pool, EMG and muscle force signals. Main results. Good control of both rest tremor and beta activity was observed with reduced power delivered when compared with conventional open loop stimulation, The supervisor avoided over- or under-stimulation which occurred when using a single controller tuned to one biomarker. When DBS amplitude was constrained, the secondary controller maintained the efficacy of stimulation by increasing pulse duration to compensate for reduced amplitude. Dual parameter control delivered effective control of the target biomarkers, with additional savings in the power delivered. Significance. Non-linear multivariable control can enable targeted suppression of motor symptoms for PD patients. Moreover, dual parameter control facilitates automatic regulation of the stimulation therapeutic dosage to prevent overstimulation, whilst providing additional power savings.


Introduction
As deep brain stimulation (DBS) enters its third decade as an approved surgical therapy for medically refractory Parkinson's disease (PD), the need to provide more effective control of patient symptoms and side-effects has focused efforts on closedloop control strategies.In contrast with conventional open-loop DBS, closed-loop DBS offers the potential to automatically adjust stimulus parameters in response to changes in patient symptoms and sideeffects, adapting the stimulation as required.Initial clinical investigations of closed-loop DBS for PD have demonstrated performance improvements and power savings with respect to open-loop DBS.On-off, dual threshold and slowly varying proportional control of subthalamic nucleus (STN) local field potential (LFP) beta activity have been demonstrated in patients for up to 8 h (Little et al 2013, 2016, Arlotti et al 2018, Velisar et al 2019), and over several days in a single case study (Nakajima et al 2021).An adaptive control scheme which adjusts stimulation based on detected sleep and medication-induced fluctuations has been demonstrated across 47 d in four patients (Gilron et al 2021).A range of different closed-loop control schemes have also been proposed and tested in silico using model simulations (Santaniello et al 2011, Grant and Lowery 2013, Pasillas-Lepine et al 2013, Popovych and Tass 2019, Su et al 2019, Fleming et al 2020a, 2020b).
To-date the focus in both simulation and experimental studies has been primarily on feedback control schemes that modify either stimulation amplitude or frequency in response to a single biomarker (Little et al 2013, Malekmohammadi et al 2016, Swann et al 2018).Although this facilitates control of symptoms or stimulation-induced sideeffects (SISEs) that are correlated with the selected biomarker, stimulation efficacy will be suboptimal when the biomarker and symptom, or side-effect, are poorly correlated (Mera et al 2011, Velisar et al 2019, Piña-Fuentes et al 2020).As the severity of different motor symptoms and side-effects can fluctuate throughout the day (Arlotti et al 2018, Zahed et al 2021) or with disease progression (Eisinger et al 2019), alternative control approaches are required to accommodate variations in multiple patient symptoms and side-effects, and to maintain long-term therapeutic efficacy.Adaptive DBS based on beta suppression has been shown to be inadequate in controlling tremor in a proportion of tremor dominant patients, suggesting that alternative adaptive control strategies may be required in patients with prominent tremor (Little andBrown 2020, Piña-Fuentes et al 2020).Adaptive control schemes enable controller parameters to adjust over time in response to changes in the system conditions (Fleming et al 2020b).However, the controller requirements may differ, depending on patient symptoms at a particular instant, or the context in which the controller is operating, for example at different times of the day or during different tasks (Tinkhauser andMoraud 2021, Fleming et al 2022).In these cases, the best control may be provided by selecting between multiple controllers depending on context or patient symptoms at that time.
For a given controller, the maximum stimulation amplitude is limited to the upper bound of the therapeutic window, defined as the range of DBS amplitudes which demonstrate therapeutic benefits without side-effects due to unintended stimulation of adjacent fiber tracts or exceeding safe charge density limits (Volkmann et al 2002, Bouthour et al 2018).Increasing the duration of the stimulation pulse while decreasing the amplitude, can help to focus stimulation on small diameter fibers located close to the electrode while avoiding larger diameter axons located further from the electrode that may be associated with SISEs (Grill andMortimer 1996, Anderson et al 2020).It follows that control approaches capable of dual modulation of stimulation pulse amplitude and duration may further help improve regulation of SISEs while maintaining efficacy in terms of symptom suppression.
Here we use these concepts to present a supervisory control structure with multi-controller architecture which uses logic-based switching to select from a bank of controllers to suppress specific symptoms.A secondary controller is employed to regulate pulse amplitude to avoid SISEs due to stimulation outside the target region by modulating the pulse duration.The controllers sit within a nested architecture, with inner and outer loops, where the primary loop targets the symptoms to be suppressed, and the slower secondary loop monitors stimulation amplitude and adjusts pulse duration to maintain the stimulation amplitude below a defined value.Within the primary control loop, a supervisor acts as a logic component that organises switching across candidate controllers, determining when to switch controllers and which controller to select.
The control strategy was implemented and tested in a computational model of the PD motor system which simulates periods of rest tremor and pathologically elevated beta activity during movement (Kehnemouyi et al 2021, Evers et al 2023).The model incorporated STN DBS, neural activity within the cortico-basal ganglia network, a motoneuron pool model and muscle activation.The results illustrate how a supervisory control scheme, used in combination with dual parameter control, can be used to target multiple symptoms and regulate SISEs during closed-loop DBS.

Supervisory and dual parameter control scheme
The proposed control scheme implemented two nested control loops, a primary loop for symptom suppression and secondary loop for control of SISEs.The primary control loop consisted of a supervisor with logic-based switching between two proportionalintegral (PI) controllers to modulate DBS amplitude to suppress pathologically elevated STN beta activity, associated with bradykinesia/rigidity during movement, or rest tremor.The secondary loop utilized a single PI controller that modulated DBS pulse duration to regulate DBS amplitude to avoid SISEs associated with activation of axons outside the STN.The supervisor selects between suppression of tremor or beta activity based on electromyographic (EMG) activity from a selected muscle.When the muscle is at rest, the supervisor selects a PI controller that is tuned to suppress rest tremor.When the muscle is active, the supervisor selects a controller tuned to suppress LFP beta activity, which will subsequently be referred to as the bradykinesia/rigidity controller.Within the secondary loop, a different PI controller simultaneously monitors DBS amplitude.If the mean DBS amplitude exceeds a given threshold value associated with induced side effects, pulse duration is increased resulting in a subsequent reduction in DBS amplitude while maintaining suppression of pathological neural oscillations.An overview of the proposed control architecture is presented in figure 1.

Supervisory control
The primary control loop consisted of a supervisor that selects between a bank of controllers, in this case two PI controllers, figure 1.At each supervisor call the average rectified value of the EMG signal, EMG ARV , over the previous 500 ms was estimated.When EMG ARV exceeded a threshold indicating voluntary movement, the supervisor selected the bradykinesia/rigidity controller; when below the threshold, the muscle was deemed to be at rest and the tremor controller was selected.During handover between controllers DBS amplitude decreased at the maximum tolerable rate before implementation of the next controller.The update intervals of both supervisor and PI controllers was 20 ms.The bradykinesia/rigidity and tremor controllers were implemented as PI controllers, as described in section 2.2, and modulated DBS amplitude to maintain their respective biomarker at or below a predefined target value.The bradykinesia/rigidity controller monitored STN LFP beta activity based on the average rectified value of the beta (22-30 Hz) filtered STN LFP during the previous 100 ms.
The biomarker for the tremor controller was similarly obtained from the average rectified value of the 4-6 Hz filtered force signal during the previous 500 ms.

Dual parameter control
To regulate SISEs, the controller in the secondary control loop modulated DBS pulse duration to constrain the mean DBS amplitude to a prespecified threshold.The SISE controller was called every 200 ms and was implemented as a PI Controller as described in section 2.2.At each controller call the mean DBS amplitude over the previous 1000 ms was estimated.A window length of 1000 ms, twice the window length used for estimating the slowest biomarker in the primary loop, was selected to minimize coupling between the dynamics of the primary and secondary control loops.

PI controller parameters
The gain parameters for each PI controller (beta, tremor and SISE) were derived using the rule-tuning method presented in (Fleming et al 2020a).The standard update law for the PI controller is defined as where u(t) is the modulated DBS parameter value at time t, u bias is a constant bias term on the output of the controller, K p is the controller proportional gain and T i is the controller integral time constant and e(t) is the normalized error between each biomarker and its respective threshold value.To prevent integral wind-up, each controller implemented conditional integration of the integral error.Integration was paused if the modulated DBS parameter reached its upper or lower parameter bounds.The constant bias term, u bias , was included if the minimum output of the controller was required to be non-zero.The integral term, T i , for each PI controller was selected based on the window length used for estimating each biomarker, l b .This window length acts as a lowpass filter on the observable dynamics in the feedback signal and so imposes a lower bound on possible T i values for each controller as (2) The proportional gain, K p , was selected based on maximum tolerable rate of change of DBS amplitude reported clinically (Little et al 2013), and the maximum error and rate of change of the error signal with DBS off.For each PI controller the proportional gain, K p , was thus selected so that the rate limit of the modulated DBS parameter was not exceeded.This was calculated by differentiation of equation (1) and setting the DBS parameter rate limit, R L , as an inequality constraint For each PI controller, the maximum proportional gain parameter was defined as where the maximum values of the absolute error, |e(t)|, and of the absolute rate of change of the error de(t)   dt of each biomarker were estimated from a 30 s simulation with no DBS applied.The clinically tolerable rate constraint, R L , was defined for each parameter as the rate required for the controller to increase the DBS parameter from its minimum to maximum value in 250 ms to avoid paraesthesia as described in (Little et al 2013).The PI controller parameters for the primary and secondary control loops are given in tables 1 and 2.

Bradykinesia/rigidity and tremor amplitude controller parameters
The PD motor network model simulated bursts of beta activity of physiological (100 ms) and pathological (>600 ms) duration in accordance with experimental observations (Tinkhauser et al 2017).T i for the bradykinesia/rigidity controller was thus selected to lie between the physiological and pathological beta burst durations at T i = 0.2 s.This corresponded to twice the window length utilized to estimate the beta biomarker and adheres to the lower parameter bound value in equation (3).For consistency, the value for the tremor controller was similarly chosen as twice the window length for estimation of the tremor biomarker, T i = 1 s.

Secondary loop control parameters
To adhere to the low T i parameter bound value, a value of 1.2 s was selected for T i for the secondary SISE controller.Three different threshold values for the maximum allowable value of the DBS mean amplitude were investigated, 2.25 mA, 1.5 mA and 1.25 mA.The parameters for each corresponding controller are given in table 2.

Motor network model
The performance of the supervisory dual parameter control scheme was investigated using a computational model of the parkinsonian motor network.The model incorporated a model of the cortico-basal ganglia network, previously presented in (Fleming et al 2020a), which projected to a model of the motoneuron pool and the first dorsal interosseous (FDI) muscle of the hand.STN DBS, simulated as monophasic, 130 Hz rectangular pulses, was applied during periods of rest tremor and elevated STN LFP beta band oscillations representative of pathological

Cortico-basal ganglia network
The cortico-basal ganglia model (Fleming et al 2020a) was adapted to simulate periods of centrally generated tremor and elevated beta band activity.The cortico-basal ganglia network consisted of six populations representing layer V cortical pyramidal neurons, cortical interneurons, the STN, the globus pallidus external (GPe) and internal segments and the thalamus.Each population comprised 100 neurons, represented using Hodgkin-Huxley type conductancebased models (Terman et al 2002, Otsuka et al 2004, Rubin and Terman 2004, Pospischil et al 2008, Hahn  and McIntyre 2010, Foust et al 2011, Fleming et al   2020a).Cortical pyramidal neurons projecting to the STN were simulated using multi-compartment neuron models, to allow coupling with the extracellular field.All other neurons in the network were simulated using single compartment models.Extracellular stimulation of axon collaterals projecting to the STN by DBS and the STN LFP recorded at the electrode due to aggregate synaptic activity were simulated as described in (Fleming et al 2020a).The model parameters used to simulate periods of tremor and elevated beta activity are detailed below in sections 2.3.5 and 2.3.6.Full details on the structure of the corticobasal ganglia network are given in (Fleming et al 2020a).

Motoneuron pool
Descending cortical neurons projected to a model of the FDI motoneuron pool comprised of one hundred motoneurons.Each motoneuron was randomly connected to 25 cortical pyramidal neurons and 100 independent spiking neurons.Motoneurons were simulated using a five-compartment Hodgkin-Huxley based model consisting of a cell soma and four dendritic branch compartments (Powers et al 2012, Powers and Heckman 2017, Senneff et al 2019).Motoneuron properties varied across the population from the smallest, earliest recruited motor unit (i = 1) up to the largest, latest recruited motor unit (i = 100).Conductance parameters were nonlinearly distributed between lower and upper bounds to generate a greater proportion of low threshold motor units as described in (Powers et al 2012, Powers andHeckman 2017).The parameter bounds for the motoneuron parameters are detailed in table 14 of supplemental material.
The membrane voltage of each compartment in the motoneuron model is described by where C m ir is the membrane capacitance, and v m ir the membrane potential, of the rth compartment of the ith motoneuron.I Ion ir are the ionic currents in each compartment, represented by transient and persistent sodium currents, voltage-sensitive delayed rectifier and calcium activated potassium currents and a hyperpolarization-activated mixed cation current in the soma of each motoneuron, and a slow inactivating, low-threshold L-type calcium current in the dendritic compartments.The synaptic input to the motoneuron pool was simulated as the summation of descending cortical inputs, I syn ir , and 100 independent synaptic inputs, I Ind ir , randomly selected from a population of six hundred Poissondistributed spike trains representing asynchronous cortical inputs and synaptic inputs from other pathways.Current flow into the rth motoneuron compartment from adjacent compartments was represented as I coupling ir .Membrane noise was represented using an independent zero mean gaussian signal applied to each dendritic compartment, I App ir .The current generated on the dendritic compartments, r, of motoneuron i due to synaptic input from cortical neuron k was described by where g Ctx_MU is the absolute maximum synaptic strength of the cortico-motoneuron connection, y ir k (t) represents the effective synaptic strength between the kth presynaptic cortical neuron and the rth dendritic compartment of the postsynaptic motoneuron i and E rev is the reversal potential of the excitatory synapse, i.e. 0 mV.Cortico-motoneuron synapses were implemented using a phenomenological model of short-term synaptic plasticity presented in (Tsodyks andMarkram 1997, Tsodyks et al 2000).The phenomenological synapse model evolves according to the following system of kinetic equations where x, y and z are the fractions of synaptic resources in the recovered, active, and inactive states, respectively, t sp is the spike time of presynaptic spikes, τ I is the decay constant of the postsynaptic currents and is the recovery time from synaptic depression.These equations describe the use of synaptic resources by each presynaptic spike where a fraction u of the available resources x is used by each presynaptic spike.For facilitating synapses, the effective use of the synaptic resources is described by u, where the parameter U determines the increase in u after each spike.In facilitating synapses u is increased with each presynaptic spike and returns to baseline with a time constant of τ facil .If τ facil → 0, facilitation is not exhibited, and u is identical to U for each spike.Full details regarding the synapse model can be found in its original publication (Tsodyks and Markram 1997), while the parameter values used for the synapse model in this study are detailed in table 15 of supplemental material.

Muscle force
The force produced by each motoneuron was simulated using the force model introduced by (Fuglevand et al 1993), whereby motor unit twitch force was simulated as the impulse response of a critically damped, second-order system, incorporating nonlinearity of the force-firing rate relationship.The twitch force, f i (t) , in response to an action potential generated by motor unit i is described by where T i and P i are the twitch contraction time and peak force of the ith motor unit in the pool, and nipi i (t) is the instantaneous interpulse interval of motoneuron i at time t, normalized with respect to T i .The values of the twitch contraction time and peak force were exponentially distributed across motor units with contraction time decreasing exponentially (90-30.4ms) and peak force increasing (1-95.5 mN) from the first (i = 1) to the last (i = 100) recruited motor unit.The total force produced by the muscle was calculated as the summation of the twitch forces generated by all active motor units.

Surface EMG
The surface EMG signal was simulated as the linear summation of extracellular motor unit action potentials generated in response to spiking activity in the motoneuron pool.Each motoneuron was assigned an extracellular motor unit action potential waveform simulated using an anatomical model of the FDI muscle (Pereira Botelho et al 2019).The spike train of each motoneuron was convolved with the corresponding action potential waveform for that motor unit to simulate the associated train of extracellular motor unit action potentials.Motor unit action potential trains were then summed to generate the surface EMG signal, figure 2(b).

Tremor activity
The model was configured to simulate two states, a high beta state and a rest tremor state, to demonstrate the ability of the proposed control scheme to target different symptoms (rest tremor and bradykinesia/rigidity) as required.To simulate rest tremor, a periodic excitatory input was applied to neurons in the thalamic population of the cortico-basal ganglia network to induce synchronous bursts of activity at approximately 5 Hz.The influence of synchronous thalamic oscillatory activity on the cortical pyramidal population represented a rudimentary model of the dimmer-switch hypothesis for tremor generation in PD, where oscillatory activity at the tremor frequency arises in the basal ganglia and is propagated through the cerebello-thalamo-cortical circuit (Helmich 2018).

Beta activity
Elevated beta band activity in the cortico-basal ganglia network was simulated by applying a bias current to cortical pyramidal neurons to shift the mean firing rate of the population into the beta band as described in (Fleming et al 2020a).The beta activity was amplified within the reciprocally connected STN and GPe circuit and propagated throughout the network.Cortical beta activity was modulated to generate bursts of varying duration, with pathological beta activity defined as bursts longer than 600 ms (Tinkhauser et al 2017).Beta activity was simulated during movement to emulate a bradykinetic/akinetic state, where asynchronous synaptic drive to the motoneuron pool was increased.During rest tremor, as described in section 2.3.5, desynchronization of beta activity in the cortico-basal ganglia network was simulated as observed clinically (Hirschmann et al 2013, Qasim et al 2016).Periods of elevated beta activity at rest were not simulated.

Controller performance
The performance of the DBS control scheme was quantified using the normalized mean error of the targeted biomarker signal (tremor or beta activity) and the mean power delivered by the stimulation.The normalized mean error was defined as the mean value of the half-wave rectified error of the target biomarker during the tremor (ME Tremor ) or beta (ME Beta ) period, divided by the mean error with DBS off where b est is the biomarker value estimated at time t and b target is the predefined target value for the corresponding biomarker.The aim of the controller is thus to reduce pathological tremor or beta-band oscillations to the target value.The mean power delivered during the tremor (PD Tremor ) and beta (PD Beta ) periods was estimated as where Z E is the electrode impedance, assumed to be 0.5 kΩ, T O is the duration of the observation period and I DBS is the DBS current.The mean power delivered during closed-loop DBS was normalized with respect to a reference value representing the power delivered during monophasic open-loop DBS with a pulse amplitude, duration and frequency of 2.5 mA, 60 µs and 130 Hz.

Simulation details
The behaviour of the motor network model was first investigated during continuous 130 Hz, monophasic open-loop DBS during periods of simulated tremor and beta activity.A parameter sweep was conducted to characterize the relationship between DBS amplitude and pulse duration and the magnitude of tremor and beta oscillations in the force and STN LFP signals.Ten simulations were conducted for each parameter combination with initial seeds varied between simulations.Following this, the performance of the supervisory and dual parameter control architecture was investigated during simulations which included 30 s periods of rest tremor followed by 30 s of elevated beta activity.
The ability of the primary loop supervisor to control both tremor and beta activity was investigated with the secondary SISE controller disabled.The SISE controller was then enabled, and the performance of the full supervisor and dual parameter control scheme was explored for three different values of the maximum allowable DBS amplitude.Ten simulations were conducted for each condition simulated.
The model was simulated in the NEURON simulation environment (Hines and Carnevale 1997), implemented in Python using the PyNN API package (Davison 2008) and run on a high-performance computing cluster.The model was numerically integrated using the Crank-Nicholson method with a 0.01 ms timestep for all simulations.Post-processing and signal analysis were conducted using custom scripts developed in MATLAB (MathWorks, Inc., Natick, MA).

Results
The behaviour of the model was first characterized during periods of elevated tremor and beta oscillations with open-loop DBS.The efficacy of the supervisory and dual parameter control architecture at suppressing oscillatory activity and the corresponding power were then examined.

Motor network activity during elevated tremor and beta activity
Oscillatory activity was observed across the motor network model at the tremor frequency and within the beta frequency band during tremor and beta periods, respectively, figure 3. Time-frequency analysis of the STN LFP exhibited centrally generated oscillatory activity within the corresponding frequency bands as the network transitioned from a period of tremor to beta activity, figure 3(a).
During the tremor period, the cortico-basal ganglia displayed oscillatory activity in the STN LFP at the 5 Hz tremor frequency, with peaks evident in the power spectrum at the tremor frequency and first harmonic, figures 3(a) and (b, iii).The motoneuron pool demonstrated similar oscillatory behaviour, with a mixture of doublet and single spiking observed, figure 3(b, v).The 5 Hz oscillation was present also in the muscle force and EMG signals, figure 3(a, iii-iv).Cortico-basal ganglia and motoneuron activity were coherent during the tremor period, with coherence spectra exhibiting significant coherence between the STN LFP and EMG, and between the EMG and force signals, at the tremor frequency and higher harmonics, figure 3(a, vi).
During the beta period, the cortico-basal ganglia STN LFP displayed oscillatory activity at approximately 26 Hz, figure 3(c, i-ii).Simultaneous muscle activity was also simulated during this period.A range of motoneuron firing rates were observed, with mean firing rates ranging from 5-26 Hz, figure 3(c,  v).Muscle force remained approximately constant, figure 3(c, iii), and EMG activity was predominantly asynchronous, figure 3(c, iv).Significant coherence was not observed between the STN LFP and EMG signals, figure 3(c, vi), due to the dominant contribution of independent uncorrelated synaptic inputs to the motoneuron pool during this period.Force and EMG were coherent across a range of frequencies, figure 3(c, vi).

Network oscillatory behaviour during open-loop DBS
The suppression of STN LFP beta activity with increasing DBS amplitude and pulse duration is presented in figure 4(a).The window for suppression of beta activity, defined as the DBS amplitude range which resulted in a 5%-95% reduction in the magnitude of oscillations, narrowed as pulse duration increased, from 0.75-2.75mA for a pulse duration of 60 µs to 0.25-1.5 mA at 120 µs, figure 4(a).An example of the variation in the magnitude of the beta filtered LFP during tremor and beta periods in response to open-loop DBS with 2.50 mA, 130 Hz and 60 µs amplitude, frequency and pulse duration is illustrated in figure 4(d).
The suppression of tremor in the force and EMG signals as DBS amplitude and pulse duration were increased is similarly presented in figures 4(b) and (c).With the model in the tremor state, the suppression window decreased from 0.75-2 mA to 0.5-1.25 mA as pulse duration was increased from 60 µs to 120 µs, figure 4(c).No further reduction in the width of the suppression window was observed for pulse durations greater than 120 µs.The change in EMG and muscle force due to open-loop DBS with 2.50 mA, 130 Hz and 60 µs amplitude, frequency and pulse duration is illustrated in figures 4(e) and (f).

Supervisory control
Continuous implementation of either the tremor or bradykinesia/rigidity controller alone resulted in suboptimal performance during periods in which the controller was mismatched with the dominant symptom.While the tremor controller suppressed tremor to its target value during the period of rest tremor, it was relatively insensitive to elevated beta activity, supplemental figure 1. Conversely, the bradykinesia/rigidity controller suppressed beta to its target value during the beta period, but resulted in overstimulation during the tremor period, at the cost of substantially increased power, supplemental figure 2.
When supervisory control was introduced, the supervisor successfully switched between tremor and bradykinesia/rigidity controllers depending on whether EMG ARV was above or below the selected threshold, figure 5.During periods of rest tremor, EMG lay below the threshold and the tremor controller was selected, resulting in a 86% and 57% reduction in tremor and power, ME Tremor and PD Tremor , respectively.Selection of the bradykinesia/rigidity controller during the beta period provided similar suppression of beta activity with a 83% and 28% reduction in the pathological beta and power delivered, ME Beta and PD Beta , respectively.

Supervisory and dual parameter control
Incorporation of the secondary control loop in the control architecture enabled dual parameter control of the DBS pulse amplitude and duration.Three maximum threshold values for DBS amplitude were investigated.With an amplitude threshold of 2.25 mA, the DBS amplitude did not exceed the threshold and the SISE controller did not modulate the pulse duration.For a maximum DBS amplitude of 1.50 mA, the DBS amplitude approached the amplitude threshold during both tremor and beta periods.The SISE controller then increased pulse duration to maintain efficacy while reducing pulse amplitude, figure 6(a).Pulse duration was modulated by the SISE controller between 70-80 µs (tremor period) and 90-100 µs (beta period) to maintain DBS amplitude below the 1.50 mA threshold, figure 6(a).The suppression of tremor and beta was comparable to that observed during supervisory DBS amplitude control, with modest additional improvements in power delivered.These improvements corresponded to a further 3% and 14% reduction in the PD Tremor and PD Beta relative to the power consumed during supervisory control of amplitude alone.
When the amplitude threshold was further reduced to 1.25 mA the SISE controller increased pulse duration to 90 µs during the tremor period and 120 µs during the beta period, figure 6(b).ME Tremor and ME Beta remained equivalent to the values observed during supervisory control of the DBS amplitude alone, but resulted in an additional 7% and 17% reduction in observed PD Tremor and PD Beta .
The performance of the supervisory and dual parameter control scheme in terms of the mean error and DBS power delivered is summarized in figure 7.

Discussion
Here we present a closed-loop control strategy for DBS that utilizes a supervisor to switch between two controllers tuned to suppress rest tremor and beta oscillations, along with dual parameter control of DBS pulse amplitude and duration.The control scheme targets suppression of either tremor or LFP beta activity and maintains suppression of the associated biomarker while limiting DBS amplitude to avoid stimulation-induced side effects.The performance of the control strategy was investigated in a computational model of the PD motor network which simulated periods of rest tremor and elevated beta activity associated bradykinesia.The proposed supervisory and dual parameter control strategy outperformed single biomarker and single parameter modulation strategies by delivering the level of stimulation necessary to maintain each biomarker at its target level with lower total power delivered.The use of surface EMG to identify muscle activity and switch between controllers provides an illustrative example of how peripheral or externally sensed signals can be used to provide contextual information to inform the controller selection.Depending on the requirements of individual patients, this could be substituted with alternative signals which provide information on the symptom to be targeted at a given time.

Motor network behaviour, comparison with experimental studies
Beta activity within the basal ganglia was simulated to decrease as rest tremor increased and to recover again once tremor ceased, figure 3(a), consistent with experimental observations of suppression of STN LFP beta band power with the onset of rest tremor in PD patients (Wang et al 2005, Hirschmann et al 2013).A mixture of tremor and beta oscillatory activity in microelectrode recordings from STN neurons has also been reported (Du et al 2018), indicating that oscillations at both frequencies may co-exist.Elevated coherence between the STN LFP, EMG and force signals at the tremor frequency was observed in the model, figure 3 (b, vi).Tremor in the force signal was driven by singlet and doublet firing of motoneurons at the tremor frequency, figure 3(b, v), consistent with motoneuron firing patterns observed in the FDI of PD patients (Christakos et al 2009).Tremor was simulated to occur when the muscle was at rest, while elevated beta activity was simulated during isometric muscle activity, figure 3(e).
Due to modulation of beta activity associated with movement there have been concerns that beta power may not be an appropriate biomarker for closed-loop DBS during movement.However, prolonged beta activity during movement has been associated with both decrements in movement velocity (Lofredi et al 2019) and with freezing-of-gait symptoms (Anidi et al 2018).Additionally, recent literature in both humans and parkinsonian rats has demonstrated successful use of beta activity as a biomarker for closed-loop DBS control during movement (Kehnemouyi et al 2021, Evers et al 2023).The controller thus targeted beta activity when the muscle was active, and tremor when at rest.Coherence between STN LFP and EMG signals was not observed during elevated beta periods, figure 3(c, vi).In the model, muscle activity was primarily driven by uncorrelated independent synaptic inputs to the motoneuron pool with relatively weak synchronous beta inputs from the cortex, thus limiting coherence between the STN LFP and EMG during simulated muscle activation.
Coherence between the STN LFP and EMG signals at the tremor frequency, similar to that observed in the model, has been observed in tremor dominant PD patients (Reck et al 2009, Tass et al 2010, Hirschmann et al 2013).Significant STN-LFP EMG coherence was not observed in akinesia-dominant patients, though partial directed coherence, associated with afferent activity, was observed between EMG and sites above the STN.Within the periphery, increased intermuscular coherence in the theta, alpha and beta frequency range has been observed in PD patients during isometric leg extension when compared with age-matched controls (Flood et al 2019), and at the tremor frequency in the upper limb of PD patients during rest and postural tremor (van der Stouwe et al 2015), and during dynamic motor tasks (Laine and Valero-Cuevas 2020).

Influence of pulse duration on stimulation therapeutic window
The application of continuous high-frequency openloop DBS led to a reduction in tremor, figures 4(b) and (c), and beta activity in the STN LFP figure 4(a).This is consistent with clinical investigations of openloop DBS for PD, where STN DBS reduces tremor (Sturman et al 2004, Blahak et al 2007), and suppresses STN LFP beta band power in parallel with improvements in bradykinesia and akinesia (Kühn et al 2006, Kehnemouyi et al 2021).Exploration of the window for suppression of oscillatory activity revealed a unique window for tremor and beta oscillations.The suppression window ranged from 0.75-2 mA for tremor and from 0.75-2.75mA for beta oscillations for a 60 µs pulse duration, figures 4(a)-(c).A similar differential response of motor symptoms to DBS has been reported clinically, with higher stimulation amplitudes required for full suppression of bradykinetic symptoms than for tremor (Mera et al 2011).In the presented model, reducing pulse duration increased the width of the suppression window for both tremor and beta oscillations, figures 4(a)-(c).Therapeutic window expansion has been observed in Computational investigations have observed that longer pulse durations (>60 µs) result in more focal activation of small fibers responsible for the therapeutic efficacy of DBS at lower DBS amplitudes (Anderson et al 2020).This is in line with classic biophysical studies of electrical stimulation where short pulse durations selectivity for large diameter, side-effect inducing fibers (Gorman and Mortimer 1983).Moreover, the selectivity for these fibers is increased with increasing distance from the point of stimulation (Grill and Mortimer 1996).It has been suggested that the observed clinical benefits of shorter pulse durations may be related to the charge or energy equivalence metrics utilized to maintain the therapeutic stimulation dosage during open-loop DBS parameter tuning (Anderson et al 2020).These metrics may not maintain the intended stimulation dosage when pulse duration is varied as they result in a greater or reduced, but never equivalent, volume of tissue activated (Anderson et al 2020).
Here we have made use of the ability of long pulse durations to focus stimulation on small diameter fibers lying close to the electrode and away from larger diameter fibers lying further away.The controller acts to identify the minimum current required to suppress beta activity via stimulation of STN afferent fibers lying close to the electrode.Large diameter fibers lying adjacent to this region may also be stimulated, resulting in potential side-effects.Increasing the pulse duration, will temporarily result in over stimulation of the adjacent small diameter fibers and in beta activity falling below the target for therapeutic efficacy.To compensate, the amplitude controller will reduce DBS amplitude allowing the beta activity to return to the target level, figures 6(e) and (j).The updated amplitude will lie below the threshold for stimulation of the more distant large diameter axons thereby reducing the risk of side-effects while maintaining beta suppression and therapeutic efficacy.This implementation of two control loops for symptom and side-effect suppression by separately controlling both pulse amplitude and duration, takes advantage of the nonlinear relationship observed between these parameters in the stimulation strength-duration curves where pulses of equivalent charge per phase but of different amplitude and durations do not necessarily activate the same neurons.

Controller performance
The model was configured to simulate two states, a high beta state during movement and a rest tremor state, associated with beta desynchronization (Hirschmann et al 2013, Qasim et al 2016, Asch et al 2020), to demonstrate the ability of the control scheme to target different symptoms as required.Both tremor and bradykinesia/rigidity controllers were effective at suppressing the associated biomarker to its target value.The implementation of a supervisor to select between the tremor or bradykinesia/rigidity controller, based here on detected muscle activity, elicited superior performance when compared with the use of a single controller, figure 7(a).The supervisor selected the controller which provided the necessary therapeutic dosage to maintain the relevant biomarker at the target level, while avoiding over or under stimulation, figure 7(a).The controllers, however, were less efficient when the biomarkers did not correlate with the dominant pathological oscillation, resulting in either over-or under-stimulation, with an associated increase in either power consumption or error, figures 7(a) and (b).This error would be observed as a failure to control symptoms, similar to the breakthrough tremor which has been reported in a proportion of tremor dominant patients during betabased adaptive DBS (Little andBrown 2020, Piña-Fuentes et al 2020).
As beta desynchronization is observed during rest tremor (Hirschmann et al 2013, Qasim et al 2016, Asch et al 2020), co-occurrence of pathological beta activity and rest tremor is unlikely to occur.Action tremor may occur during movement but may be suppressed by DBS acting to suppress beta activity.More generally, co-occurrence of multiple symptoms could lead to suboptimal performance of the control strategy.In this scenario, the supervisor would prioritize suppression of the symptom it is programmed to target in each state, while providing inefficient or suboptimal control of the other symptom.
A state corresponding to elevated beta activity at rest, in the absence of tremor, was not considered here.Beta activity is present during rest with longer and higher amplitude beta bursts, and higher averaged beta power, observed during rest than during movement.Suppression of beta activity during rest when tremor is not present could be incorporated by setting the supervisor to default to the beta controller when tremor is not detected, and the tremor controller is inactive, while at rest.Switching the controllers based on tremor alone could lead to undesirable, rapid switching between controllers as tremor is repeatedly suppressed and reemerges.
Although the characteristic desynchronization and subsequent resynchronization of beta activity prior to and immediately following movement was not simulated, we anticipate that this would have minimal influence on the performance of the proposed control strategy.In (Lofredi et al 2019) the authors demonstrated that while beta bursts occurred less frequently and were shorter in duration during movement than at rest, they still surpassed the threshold defined at rest (Lofredi et al 2019).
Additionally, as the duration of beta suppression prior to movement is relatively short, the lowpass filter used to extract the envelop of beta activity in the STN LFP in conjunction with the integral component of the PI controller in the proposed control strategy would both act to minimize the influence of these fast fluctuations in beta activity on the output produced by the controller.
When a maximum tolerable stimulation amplitude was applied to the controller, dual parameter control adapted the stimulation pulse duration to constrain the DBS amplitude without loss of therapeutic efficacy, figures 6 and 7(b).Moreover, due to power being proportional to the square of current, increasing pulse duration while lowering pulse amplitude resulted in further power consumption reductions during each oscillatory period.When the amplitude was limited to 1.25 mA an additional 7% and 17% reduction was observed in the PD Tremor and PD Beta values, figure 7(b).More generally, the results suggest that simultaneous modulation of DBS amplitude and pulse duration may be beneficial for closed-loop DBS even when induced side-effects are not a primary concern.

Control strategy
The control strategy presented is a simple example of a context-dependent control approach for PD motor symptom suppression during voluntary movement.The supervisor selected between tremor and bradykinesia/rigidity controllers based on EMG activity, figures 5 and 7(a).In practice, there are a wealth of alternative signals which may be used in contextdependent strategies to select controllers based on the time of day (Toth et al 2020), accelerometry during specific motions and actions (Opri et al 2020, Zamora et al 2022, Borzì et al 2021) or to target the suppression of other axial disease symptoms such as action tremor or freezing of gait (Horak and Mancini 2013).Due to the heterogenous expression of symptoms between patients, future implementations of the proposed context-dependent approach may however require patient-specific tailoring of the controllers included in the controller bank to target each patient's individual dominant symptoms.In addition, selection of the most contextually appropriate controller from the bank of potential controllers may also require more sophisticated classification methods than the simple threshold approach applied (Haddock et al 2019).Moreover, more nuanced approaches for handling controller transitions such as sample-and-hold approaches or interpolation, may be preferable.In practice, ramping down from the last controller output prior to handover may result in an undesirable re-emergence of symptoms during the handover period.Lastly, alternative controllers such as those based on phase-based (Holt et al 2019), adaptive (Fleming et al 2020b) or frequency modulation methods (Xie et al 2015, Jia et al 2017) may be readily incorporated, in addition to varying electrode contacts or utilizing current steering approaches.Similarly, utilization of biomarkers correlated with side-effects such as speech impairment (Khan et al 2014) or hyperkinesia narrowband gamma (Swann et al 2016) may provide more appropriate measures of side-effect severity than the tolerable amplitude limit implemented here.
Other important considerations when developing control strategies for closed-loop DBS is the definition or quantification of stimulation efficacy.Open-loop DBS provided the greatest reduction in beta-band and tremor oscillations in the investigated model, figure 7(a).However, it also resulted in the highest power consumption, figure 7(b).The generalizability of the improvements in efficacy and power observed with the proposed control strategy will vary depending on the biomarkers and controllers implemented.As most DBS devices still rely on primary cell technology, quantification of the tradeoff between clinical efficacy and power consumption is important for effective control strategy design.Algorithms which stimulate with higher power than necessary may result in a more frequent need for battery replacement surgery than open-loop DBS, while algorithms which consume less power can potentially extend battery life, subsequently reducing the frequency of replacement surgery when compared to open-loop DBS.In addition, the therapeutic effects of open-loop stimulation may vary over time due to system variations (Fleming et al 2020b, van Rheede et al 2022), which can be accommodated by the proposed control strategy.Thus, metrics that suitably quantify the trade-off between symptom suppression and power consumption over time are required.
The PI controller parameters were derived based on the dynamics of the tremor and beta biomarkers and the clinical rate constraints of the modulated DBS parameter.This resulted in controllers which modulated DBS parameters with dynamics similar to the associated biomarker.Modulation of DBS amplitude was investigated by Arlotti et al (2018) where the DBS amplitude was linearly modulated in response to slow biomarker variations due to levodopa-induced modulations of beta oscillations (Arlotti et al 2018).Previous clinical investigations of closed-loop DBS using on-off or dual threshold control have increased DBS amplitude at the maximum clinically tolerable rate limit (Little et al 2013, 2016, Velisar et al 2019), without consideration of the biomarker dynamics.A mismatch between biomarker and parameter modulation dynamics may, however, lead to stimulation parameters being adjusted faster than the resulting changes are reflected in the biomarker signal potentially inducing oscillations within the system.To minimize coupling between the dynamics of the two controllers here, the pulse duration controller was called at every tenth call of the amplitude controller.Mean DBS amplitude was estimated over a 1000 ms window, greater than twice the duration of the lowest motor symptom oscillation, i.e. the 5 Hz tremor.This ensured that the biomarker used for control of pulse duration displayed slower dynamics than the slowest motor symptom biomarker.Minimizing the coupling between control loops in the control strategy ensured that simultaneous modulation of the DBS amplitude and pulse duration did not unintentionally introduce additional oscillations into the system due to interactions between the control loops.In practice, successful implementation of chronic closed-loop DBS would require that the emergence of side-effects due to stimulation is rare.

Model considerations
The cortico-basal ganglia network model presented simulated both tremor and beta oscillatory activity.Previous computational models have simulated either tremor (Terman et al 2002, Rubin and Terman 2004, Santaniello et al 2011) or elevated beta activity (Hahn and McIntyre 2010, Nevado-Holgado et al 2010, 2014, Pavlides et al 2012, Grant and Lowery 2013, Davidson et al 2016, Kumaravelu et al 2016), however, few studies have explored the oscillatory activity associated with both motor symptoms in the same model (Kang and Lowery 2013).Extension of the model to include spinal motoneurons, EMG and force in the motor network model facilitated the investigation of tremor-based closedloop DBS approaches which better approximate clinical investigations where tremor is monitored using accelerometry signals (Malekmohammadi et al 2016, Haddock et al 2018).It also opens the possibility of exploring EMG as a potential biomarker for closedloop control.In this way, the motor network model can provide a surrogate testbed for the development of beta-or tremor-based closed-loop DBS approaches prior to clinical investigation in patients.
A number of model limitations should be considered when interpreting the results.Afferent feedback from the muscle to the spinal cord or central nervous system was not incorporated, though this can also contribute to tremor and potentially beta activity within the motor network.The integration of synaptic inputs to the motoneuron pool from multiple sources, including cortical and reticulospinal drives, spinal interneurons and afferent inputs, has been simplified.In addition, the generation of tremor and beta activity within the cortico-basal ganglia network was simplified by introducing oscillatory activity at each frequency to the network.A more physiological representation of the generation of beta and tremor activity could be included by extending the model of the cortex and thalamus to capture changes in coupling within and between both structures (Farokhniaee and Lowery 2021).Moreover, the model at present does not simulate action tremor, pathologically elevated beta-band oscillations during rest (Kühn et al 2006) or dynamic variation in beta-band oscillations at movement onset (Kühn et al 2004).Lastly, for the purposes of illustrating the dual parameter control mechanism, here we have assumed that SISEs occur when the stimulus amplitude exceeds a given threshold associated with activation of large-diameter axons that anatomically project to different regions of the cortex or neighboring networks.However, as these large diameter axons are assumed to not project directly onto the STN they have not been explicitly included in the model.In practice, the threshold for induced side-effects will likely vary with pulse width (Steigerwald et al 2018) thus modelling an alternative biomarker, directly related to side-effects, would provide a better biomarker for constraining stimulus amplitude.

Conclusion
A multivariable closed-loop control strategy to control motor symptoms and stimulation-induced side effects is presented.The control strategy simultaneously modulated DBS pulse amplitude and duration to suppress either tremor or beta activity, while maintaining stimulation amplitude below the threshold at which side effects are induced.The control strategy was investigated in a model of the PD motor network which simulated periods of centrally generated tremor and beta activity.The proposed strategy outperformed single biomarker-based strategies by maintaining effective biomarker suppression as the dominant pathological oscillation varied.Performance was maintained when DBS amplitude was limited by modulating pulse duration, leading to further savings in power consumption than during amplitude modulation alone.These results suggest that the proposed control strategy may provide improved performance over current single biomarker and single parameter modulation closedloop DBS methods.

Figure 1 .
Figure1.Supervisory and dual parameter control architecture.The control architecture is comprised of a primary and secondary control loop.The supervisor in the primary control loop implements logic-based switching to select between two DBS amplitude controllers tuned to suppress either tremor, θ(t), or beta-band, β(t), oscillatory activity depending on EMG activity.The secondary SISE controller monitored the mean stimulation amplitude, u1 (t), and modulated pulse duration, u2(t), to maintain the DBS amplitude below a defined threshold, AL0, to avoid side-effects.

Figure 2 .
Figure 2. Motor network model overview.(a) The cortico-basal ganglia network model was coupled to (b) a model of the motoneuron pool, muscle force and EMG of the FDI muscle.Muscle force and EMG were estimated as the summation of the force twitch and motor unit action potentials across the population.Tremor and beta activity in the motor network were centrally generated and LFP, EMG and force signals were simulated using the model.

Figure 3 .
Figure 3. Motor network activity.Simulated periods of centrally generated tremor and beta oscillatory activity.(a) Continuous wavelet transform of the STN LFP during periods of elevated tremor (10-40 s) and beta activity (40-70 s).During the tremor period (b) the LFP, muscle force and EMG signals displayed coherent activity at the 5 Hz tremor frequency.Tremor oscillations in the muscle force and EMG are reflected in single and doublet firing at the tremor frequency in the motoneuron pool.During the beta period which coincided with isometric muscle contraction (c) the LFP displayed synchronous beta-band activity at 26 Hz.Coherence was observed between the muscle force and EMG, but not the LFP and EMG signals, during the elevated beta period.

Figure 4 .
Figure 4. Motor network model oscillatory behaviour during open-loop DBS.The magnitude of (a) STN LFP beta oscillations (b) EMG tremor and (c) force tremor in response to 130 Hz monophasic DBS as pulse amplitude and duration were varied.The magnitude of the beta oscillations was estimated during the beta period (45-50 s), while the magnitude of the tremor oscillations was estimated during the tremor period (40-45 s).(d) Beta filtered STN LFP EMG and (f) muscle force signals during the transition between tremor and beta periods in response to 2.50 mA, 130 Hz and 60 µs open-loop stimulation (black) and with DBS off (orange).Prior to 40 s the motor network model was simulated to exhibit steady state tremor oscillatory activity as illustrated in the 40-45 s tremor period.

Figure 5 .
Figure 5. Supervisory logic-based switching of DBS controllers based on muscle activity.Panels (a)-(d) present the average rectified value (ARV) of the STN LFP beta activity, EMGARV, force tremor and DBS amplitude with supervisory switching between tremor and bradykinesia/rigidity amplitude controllers during the tremor (15-45 s) and beta periods (45-75 s).

Figure 6 .
Figure 6.Supervisory dual parameter DBS control for varying SISE controller amplitude thresholds.Panels (a)-(e) demonstrate the behaviour of the beta ARV, EMG ARV, muscle force tremor ARV, DBS pulse amplitude and duration for a 1.50 mA amplitude threshold.Similarly, panels (f)-(j) illustrate the behaviour of the same signals in response to a 1.25 mA amplitude threshold.

Figure 7 .
Figure 7. Summary of supervisory and dual parameter control architecture performance.(a) Normalized mean error, ME Tremor/Beta , and (b) normalized mean power delivered, PD Tremor/Beta , during the tremor and beta periods due to amplitude modulation using the tremor controller alone, the bradykinesia/rigidity controller alone and supervisory switching between the tremor and bradykinesia/rigidity controllers.(c) ME Tremor/Beta and (d) PD Tremor/Beta during the tremor and beta periods due to supervisory and dual parameter control of DBS pulse amplitude and duration.Values are presented for three different DBS amplitude thresholds.The mean error and mean power delivered have been normalized with respect to the mean error and mean power delivered during DBS off and open-loop DBS, respectively.

Table 1 .
PI controller parameters for the primary control loop which modulated DBS amplitude by switching between controllers designed to suppress oscillatory activity associated with bradykinesia/rigidity, based on LFP beta activity, or rest tremor symptoms.

Table 2 .
PI controller parameters for the secondary control loop which modulated pulse duration to constrain the DBS mean amplitude at a given threshold.