Algorithmic design of a noise-resistant and efficient closed-loop deep brain stimulation system: A computational approach

Advances in the field of closed-loop neuromodulation call for analysis and modeling approaches capable of confronting challenges related to the complex neuronal response to stimulation and the presence of strong internal and measurement noise in neural recordings. Here we elaborate on the algorithmic aspects of a noise-resistant closed-loop subthalamic nucleus deep brain stimulation system for advanced Parkinson’s disease and treatment-refractory obsessive-compulsive disorder, ensuring remarkable performance in terms of both efficiency and selectivity of stimulation, as well as in terms of computational speed. First, we propose an efficient method drawn from dynamical systems theory, for the reliable assessment of significant nonlinear coupling between beta and high-frequency subthalamic neuronal activity, as a biomarker for feedback control. Further, we present a model-based strategy through which optimal parameters of stimulation for minimum energy desynchronizing control of neuronal activity are being identified. The strategy integrates stochastic modeling and derivative-free optimization of neural dynamics based on quadratic modeling. On the basis of numerical simulations, we demonstrate the potential of the presented modeling approach to identify, at a relatively low computational cost, stimulation settings potentially associated with a significantly higher degree of efficiency and selectivity compared with stimulation settings determined post-operatively. Our data reinforce the hypothesis that model-based control strategies are crucial for the design of novel stimulation protocols at the backstage of clinical applications.


Introduction
The use of electrical deep brain stimulation (DBS), during approximately the last 30 years, has been proven to provide striking benefits for patients with advanced Parkinson's disease (PD), essential tremor and dystonia [1][2][3][4] who have failed conventional therapies. In the interim, promising applications of this technique for the treatment of neuropsychiatric disorders have emerged, including treatment-refractory obsessive-compulsive disorder (OCD), Tourette's syndrome, major depressive disorder, drug addiction and anorexia nervosa [5][6][7][8][9]. Challenges however exist and are principally related to the optimization of the efficiency of stimulation through refined strategies at multiple peri-operative levels. Particularly, in addition to the appropriate patient selection and anatomical target determination [10], the outcome of DBS may be strongly influenced by the quality of post-operative clinical management, i.e., the adjustment of stimulation parameters and the selection of the optimal contact, usually over periods of weeks to months [11]. Apart from being considerably time consuming, this trialand-error procedure may not necessarily yield the optimal trade-off between maximal therapeutic benefit and minimal stimulation-induced side-effects [12]. Moreover, it fails to keep pace with the fact that movement and neuropsychiatric disorder symptoms may fluctuate over significantly shorter time-scales of seconds to days. Chronically, the open-loop nature and monomorph pattern of conventional high-frequency stimulation appears to favor tolerance/ habituation phenomena, while being associated with a significant rate of power consumption [13].
Against this background, closed-loop DBS is emerging as a more robust alternative and one of the most promising breakthroughs in the field of neuromodulation [14,15]. In an optimal closed-loop-stimulation scenario, delivery of maximally efficient stimulation protocols is adjusted to the fast dynamics of movement and neuropsychiatric disorder symptoms through utilization of specific biomarkers that capture the patient's clinical state in real time [16]. Hence, any algorithm designed for a maximally efficient closed-loop DBS system shall conceptually satisfy two core specifications [17]: the reliable assessment of optimal neurophysiological biomarkers for feedback control and the robust identification of alternative stimulation protocols that may be more therapeutically-and energy-efficient compared with the conventional pattern of stimulation [18,19].
Nonlinear coupling across multiple frequency bands in the basal ganglia and in cortical structures is being increasingly highlighted as a potentially predictive biomarker of PD and OCD pathophysiology [20][21][22][23][24][25]. To date, assessment of this biomarker has largely relied on evaluation of phase-amplitude coupling by means of the Hilbert transform combined with linear band-pass filtering [20,26]. Remarkably however, the respective phase reconstruction method may be characterized by a high level of susceptibility to measurement noise and a high rate of artificial phase slips [27,28], thereby discarding possibly rich information that would be revealed by employing noise-resistant or phase reconstruction-free methods [29,30]. Meanwhile, model-based control policies for the determination of temporally alternative stimulation protocols [31][32][33][34][35][36][37][38][39][40][41], though still limited, are most commonly oriented towards the minimum energy desynchronizing control of neuronal activity. The rationale behind this objective lies in indications that temporally alternative DBS waveforms, including stochastic waveforms, hold the potential to drive the neuronal dynamics within the basal ganglia back to the normal desynchronized state-namely to more irregular and less burst-like firing patterns [19,42]-thereby outperforming the action of standard DBS waveforms, the mechanism of which has been principally attributed to the reinforcement-driven regularization of neural firing patterns in the vicinity of the stimulated nucleus [43][44][45].
Considering these observations collectively and employing concepts drawn from dynamical systems and control theory, in this study we elaborate on the algorithmic aspects of a noiseresistant and efficient closed-loop neuromodulation system for advanced PD and treatmentrefractory OCD (Fig 1A and 1B). Specifically, we first state a series of methods robust to the presence of internal and measurement noise that are employed in order to reliably assess significant nonlinear coupling between beta and high-frequency subthalamic activity, as a biomarker for feedback control in the closed-loop neuromodulation scheme. We then suggest a model-based control strategy through which optimal parameters of stimulation for minimum energy desynchronization of neuronal activity are being identified.
In our analysis, we opt for a phase-reduced bursting neuron model [40,46,47]. Our motivation for this particular selection is twofold. First, phase reduction theory constitutes a powerful mathematical framework for the analysis of the synchronization and response properties of nonlinear oscillatory activity based on a single phase variable [48]. Second, since bursting activity is a prominent characteristic of subthalamic neuronal activity in PD and OCD (Fig 1C) [ [49][50][51][52], a qualitative model of neuronal bursting, like the well-established Hindmarsh-Rose model for bursting, may be a highly appropriate point of reference for capturing the respective neuronal dynamics [31,53]. It should be noted that, depending on parameter selection, the Hindmarsh-Rose model may capture a wide range of neuronal dynamics: from regular spiking to bursting to chaotic regimes and fixed-point behavior. However, in this study, we focus on a computational model able to qualitatively capture pathological neuronal dynamics, i.e., bursting behavior. Accordingly, a major part of the phase-response dynamics of the reduced model has been determined based on the Hindmarsh-Rose model for bursting [46,47]. Importantly, the employed phase-reduced model, which simulates the effect of stimulation on pathological neuronal activity, is data-driven, i.e., microelectrode recordings (MERs) acquired during subthalamic nucleus (STN) DBS surgical interventions for PD and OCD are used to estimate the unknown model parameters off-line. Thereby, the ability of the model to simulate realistic neuronal dynamics is further enhanced. A data-driven phase-reduced model of subthalamic neuronal activity was employed in [40], where we adopted a measure of the invariant density (steady-state phase distribution) of the simulated dynamical system, as a quantity inversely related to the desynchronizing effect of temporally alternative patterns of stimulation, and further provided evidence for a possible correlation of this measure with clinical effectiveness of stimulation in PD. Determination of the precise optimal parameters of stimulation is accomplished through the application of a derivative-free optimization algorithm, in particular a model-based pattern search method guided by simplex derivatives [54,55]. This approach is motivated by the fact that the neural response to DBS parameters is expected to be a complex, non-differentiable function [19,56].
Finally, extending the results of our previous work, we attempt to provide indications for a possible correlation of the invariant density measure with clinical effectiveness of stimulation in treatment-refractory OCD. Overall, the results of this study corroborate the ability of the presented modeling approach to identify stimulation settings potentially associated with a significantly higher degree of efficiency and selectivity compared with stimulation settings determined post-operatively, while guaranteeing a relatively low computational cost.

Data description
We used MER data acquired during 8 STN-DBS surgical interventions for advanced PD at Evangelismos General Hospital of Athens and 8 STN-DBS surgical interventions for treatment-refractory OCD at Grenoble University Hospital. Surgery was performed after provision Applicability of cross-frequency coupling as a biomarker for feedback control in case of treatmentrefractory OCD may be subject-specific. Accordingly, the presence of bursting neuronal activity was included as an alternative biomarker for feedback control in the respective closed-loop neuromodulation scheme; *short interburst interval and high of written informed consent by patients with advanced PD and after consideration of strict ethical guidelines and inclusion criteria for patients with treatment-refractory OCD (the consent procedure was approved by the National Consultative Ethics Committee on Health and Life Sciences 2002). Patients' clinical characteristics have been described in detail elsewhere (cases 1-5 and 7-9 [28]; cases O1-O6 and O9 [49]; case P5 [57]). In case of PD, medication was withdrawn at least 12 hours before surgery. All recordings were obtained, while patients were awake and at rest, using five tungsten microelectrodes (2mm apart, tip diameter < 25 μm, Medtronic Inc., Minneapolis, MN (Evangelismos General Hospital of Athens) / FHC Inc., Bowdoinham, ME, USA (Grenoble University Hospital)). Signals acquired during STN-DBS for advanced PD were pre-amplified, band-pass filtered between 0.1 Hz and 10 kHz, sampled at 24 kHz (Leadpoint TM Neural Activity Monitoring System, Medtronic Inc., Minneapolis, MN) and off-line band-pass filtered at 1-300Hz and 300Hz-6 kHz, applying four-pole Butterworth filters (Matlab, Mathworks, Natick, MA). Signals acquired during STN-DBS for treatment-refractory OCD were pre-amplified, band-pass filtered at two frequency bands (1-300 Hz and 300 Hz-6 kHz) and sampled at 3 kHz and 48 kHz, respectively (Neurotrek System, Alpha-Omega Engineering, Nazareth, Israel). A total of 31 acceptable MER trajectories (i.e., trajectories traversing the broadest extent of the nucleus) acquired during STN-DBS for PD and 12 acceptable MER trajectories acquired during STN-DBS for OCD were selected for offline analysis. Preprocessing of each MER included its subdivision into three distinct neuronal populations: spiking activity acquired by employing a five-point spike template [28,58], background unit activity reconstructed according to [59], and local field potential (LFP) activity . In order to keep pace with the employed phase reduced bursting neuron model, following the assessment of the biomarker for feedback control, optimal parameters of stimulation were determined only for sites at which bursting activity was recorded. A bursting or burst-like firing pattern of neuronal activity was identified according to [60] (Fig 1C).
Noise-resistant assessment of cross-frequency coupling as a biomarker for feedback control The first objective of this study was to propose a method for the assessment of nonlinear coupling between beta and high-frequency activity characterized by reduced sensitivity to internal and measurement noise compared with the combined application of linear band-pass filtering and the Hilbert transform [20-22, 24, 26]. We therefore employed a two-part technique for the designed scheme (Fig 1). First, the beta-band-frequency (13-30Hz) envelope of the highpass filtered (200-300Hz) LFPs (or, alternatively, of the high-frequency signal component (300-500Hz)) was assessed. Particularly, the high-pass filtered signal was full-wave rectified, mean subtracted and downsampled to 1kHz. The derived signal was band-pass filtered at 13-30Hz by applying the complex-valued filter proposed by [29]. This filter is designed based upon minimization of the relative variance where z is the filter output. In particular, z = f Ã x, where f is the filter's impulse response, x is the given signal and '' Ã " denotes convolution. The robustness of this filter lies at its property to increase the signal to measurement noise ratio with respect to the complete dynamics of its impulse response. Moreover, this filter has been proven to cope with strong internal noise that constitutes a prominent characteristic of the recorded subthalamic neuronal activity. Accordingly, the presented method accounts not only for the presence of measurement noise, but also for the presence of intrinsic noise which is relevant to PD and OCD pathophysiology. Following the employment of the complex-valued filter, we applied the 0-1 test for chaos [30] to a logarithmic transformation of the complex magnitude of the filter output in order to assess the presence of significant nonlinear coupling between beta and high-frequency activity in the STN of patients with PD or OCD. Nonlinear coupling corresponds to regular, non-chaotic dynamics indicated by a test outcome approximately equal to 0. In addition to being a phase reconstruction-free method for the determination of regular or chaotic dynamics in a deterministic dynamical system, the 0-1 test for chaos retains the advantage, over the traditional methods for detecting chaos (using the maximal Lyapunov exponent), of displaying reduced sensitivity to measurement noise [61]. Briefly, for the first n = 1, . . ., n max = 1000 samples of the input signal and N c = 100 values of c chosen randomly in the interval (0, π), we evaluated the translation variables VðjÞcosðjcÞ and q c ðnÞ ¼ where V is the amplitude of the input signal. Secondly, considering the presence of measurement noise, we quantified for n n max 10 ¼ n cut the damped mean squared displacement of the translation variables, as follows where M c (n) is the mean squared displacement of the translation variables, defined as EV is the expectation of V, while parameter h was defined based upon optimization of the outcome of the test across a subset of 12 MER trajectories in PD and 12 MER trajectories in OCD. We next computed the strength of correlation of D Ã c ðnÞ with linear growth as The outcome of the test, K t , was given by A test outcome, K t < 0.1 indicated the presence of regular dynamics [62], i.e., the presence of significant cross-frequency coupling.
With respect to parameter determination, we used 1s (i.e., 1000 samples) of the input signal, since this value yielded the best trade-off between low computational cost and optimal outcome of the test. In addition, N c = 100 different values of c have been proven to constitute an appropriate variable selection in [30].
The performance of the 0-1 test outcome was compared with the performance of an alternative measure of cross-frequency coupling, the modulation index. This index is based on the Kullback-Leibler (KL) distance between two distributions and its calculation involves application of the Hilbert transform combined with linear band-pass filtering [26,63].

Model-based control strategy for the identification of the optimal stimulation protocol
The phase-reduced model. We used a previously published stochastic phase-reduced model [40], inclusively allowing for the phase-response dynamics of a bursting neuron in both weak and strong perturbation regimes [46,47]. The phase-reduced model is described by the following Ito stochastic differential equation (derivation of Eq (6) is provided in S1 Text): where Considering a rectangular stimulation waveform, parameter β may be expressed as: β = wI 0 / C [64], where w represents the stimulus pulse width (expressed in μs), I 0 is the stimulus current amplitude (expressed in A) and C = 1μF/cm 2 . Variable ϕ 2 [0,1) denotes the oscillator's phase, ω is its natural frequency, while K, r and ψ symbolize the coupling strength, the degree of synchrony and the mean phase of the oscillators, respectively, in the surrounding neural population. These parameters were evaluated based on the processing of the MERs, as described in [28,40]. Parameter α represents the phase shift inherent to nonlinear coupling. This parameter was considered equal to ¼, so as to better capture the partially synchronous quasiperiodic dynamics (0 < r < 1) of the subthalamic neuronal activity in the pathological state [28,65]. ξ(t) is the zero mean Gaussian white noise, added independently to the oscillator, and η(t) is the colored (common) noise with zero mean, unitary variance and autocorrelation function C(t). Parameters σ I and σ C denote the intensity of independent and common noise, and were determined based on the spiking activity and the LFP signal component, respectively [40]. Phase sensitivity functions R I (ϕ) and R C (ϕ) were evaluated according to [40,46,66]. Function Δ(ϕ, β) represents the phase response curve (PRC) to a single (DBS) impulse and was evaluated according to [47] (S1 Fig). Values of β were appropriately scaled according to the size of perturbations upon which the PRC was constructed. Variable τ k denotes the input times (0 k < 1). We considered that the inter-impulse interval (IPI) Δτ n = τ n+1 − τ n obeys a Poisson distribution with parameter λ and that stimulus trains have a mean frequency f (expressed in Hz).
We emphasize that the MERs used for the estimation of model parameters are exclusively those for which cross-frequency coupling and/or bursting activity was identified, as described in the previous section. Through this assignment, phase model (6) is being appropriately elaborated to capture the regular dynamics of pathological neuronal activity, in addition to simulating the effect of stimulation on this activity.
The recorded neuronal activity had to follow a bursting-like pattern in order to be compatible with the model, which is a bursting neuron model. Essentially, the first part of the closedloop scheme is crucial for model parameter estimation in the second part.
After solving phase Eq (6), we employed the Perron-Frobenius operator, P, in order to extract the stochastic phase map from one stimulus cycle to the next: where p n (ϕ) and p n+1 (ϕ) are the densities of the phases at the time of the nth and (n + 1)th impulses, respectively [40]. By discretizing the phase into m = 500 bins, operator P was approximated using a m × m transition matrix (or stochastic kernel), A = [α ij ]. The iterative process (7) converges to the steady-state phase distribution, i.e., the invariant density, p st (ϕ). The invariant density vector, p st 2 R m !0 , is the eigenvector corresponding to the dominant (unit) eigenvalue of the transition matrix. In accordance with [40], we assessed the variance of the elements of the invariant density vector, s 2 p st , as a quantity inversely related to the desynchronizing effect, but potentially also to clinical effectiveness of stimulation. We may express this variance in terms of the Euclidean norm of p st , as where m p st is the mean of the elements of the invariant density vector and 1 is the m by 1 vector of ones. It should be noted that we did not assess the variance of the phase variable, ϕ, which would be analogous to the desynchronizing effect of stimulation. Rather, we wanted to place emphasis on the properties of the derived largest eigenvector of the transition matrix by assessing the variance of its elements. Identification of the optimal stimulation protocol was based on minimization of the cost function wfÁR is the stimulation power [67]. Parameter g is a penalizing scalar (we set g = 0.25), while R represents the electrode impendance (we considered R = 1000O).
Model-based derivative-free optimization. We considered the determination of the d = 4 optimal stimulation parameters for minimum energy desynchronizing control of neuronal activity as a constrained optimization problem defined as: subject to 30 x 1 210, 0:001 Determination of the pulse-width constraints was based on evidence that pulse durations lower than 60 μs may lead to increased selectivity of stimulation, i.e., activation of the targeted neural elements without activation of distant pyramidal tract fibers, and therefore possibly also to an increased therapeutic window of DBS [68].
Cost function (9) is expected to exhibit non-smooth behavior. We therefore employed a clever derivative-free optimization algorithm, in particular a model-based pattern search method guided by simplex derivatives (SID-PSM) [54,55]. This algorithm belongs to the general class of direct search methods of directional type. Its distinguishing feature is the use of past objective function evaluations to improve the computational algorithmic efficiency.
Direct search methods of directional type are iterative algorithms, where the process of finding a new iterate (x k+1 ) can be organized in a search step and a poll step. The search step is optional, not required to ensure the convergence of the algorithm, and typically used to improve the numerical efficiency. In SID-PSM, past function evaluations are used to compute a local quadratic model of the objective function, which is minimized in a region of interest (a ball around the current iterate x k ). The corresponding minimizer is then evaluated for the original objective function and, if it decreases the value of the current iterate F(x k ), it is accepted as the new iterate x k+1 . In this case, the iteration is declared as successful and the poll step is skipped.
If the search step fails in obtaining a better point, the algorithm will obligatorily perform the poll step, where a local search around the current best point x k is considered by evaluating the feasible points belonging to the poll set P k = {x k + α k d k : d k 2 D k }. In this case α k represents a step size parameter and D k a set of directions with good geometrical properties, typically corresponding to a positive generating set [69]. This evaluation process is opportunistic, meaning that a new point is accepted as a new iterate once it decreases the value of the objective function, without evaluating the remaining poll points. Thus, the order by which the directions are tested is relevant for the algorithmic performance. Using previous evaluated points, SID-PSM computes a descent indicator, based on simplex derivatives, and directions are tested according to the angle that they make with this descent indicator. If a better point is found during this testing procedure, the new point is accepted as a new iterate x k+1 , and the iteration is declared as successful. Otherwise, x k+1 = x k and the iteration will be unsuccessful.
At the end of each iteration, the step size is updated: decreased at unsuccessful iterations and maintained or increased for successful ones. Points evaluated during both search and poll steps at iteration k are stored in a list, X k , allowing the computation of the quadratic models and of the descent indicators at no further expense in terms of function evaluations.
Incorporating quadratic models. Given a sample set Y 0 k ¼ fy 00 k ; y 01 k ; . . . ; y 0pk k g X k (where y 00 k ¼ x k ), a quadratic polynomial basis e = {e 0 (x), e 1 (x), . . ., e q (x)} and a quadratic polynomial model mðy 0 k Þ ¼ a T k eðy 0 k Þ, the condition for quadratic polynomial interpolation can be expressed as where : System (11) is determined, if p k = q = (d + 1)(d + 2) / 2 − 1, overdetermined, if p k > q, and underdetermined, if p k < q. The quadratic polynomial model may also be written as where g k and H k represent the gradient and the Hessian of the model, respectively. The quality of the quadratic model as approximation to the original function is strongly dependent on the norm of H k . Thus, when a reasonable number of points is already available, but system (11) is still underdetermined (d + 1 p k < q), a minimum Frobenius norm (MFN) solution is computed by minimizing the Frobenius norm of the Hessian H k subject to the interpolation conditions: . . . ; p k : If q < p k (d + 1)(d + 2) − 1, a regression quadratic model is considered by solving the problem: When there are more points available than the ones required to build the overdetermined model, a subset of X k is selected, with 80% of the points chosen near to the current iterate, x k , and the remaining 20% chosen far from it. The search step will be defined by minimizing the quadratic model in kdk. Parameter σ k equals to 1, if the previous iteration was unsuccessful, or 2, otherwise.
Incorporating descent indicators. At the beginning of each poll step, a sample set with some desirable geometric properties may be identified. This sample set should be part of a ball of the same or larger radius of the smallest ball enclosing the poll set P k . The simplex gradient of F at x k , g k ¼ r S k Fðx k Þ, is the solution of the following system where S k ¼ ½y 1 k À x k Á Á Á y r k k À x k and δ k ðF; S k Þ ¼ ½Fðy 1 k Þ À Fðx k Þ Á Á Á Fðy r k k Þ À Fðx k Þ T . Again, this system allows determined (r k = d), underdetermined (r k < d) or overdetermined (r k > d) solutions. Underdetermined and overdetermined forms of simplex gradients are computed by where U k Σ k V k T is the singular value decomposition of the scaled matrix S T k =D k and Δ k is defined as in the previous section.
Importantly, the quality of simplex gradients as approximations to some form of real function derivatives has been established even in the non-smooth case [70] and depends on the geometrical properties of the sample set. SID-PSM uses the geometrical notion of Λ-poisedness to determine the quality of the geometry of the sample set and considers that a sample set Y k is Λ-poised, if kΣ À 1 k k L, for some positive constant Λ. Thus, the negative simplex gradient À r S k Fðx k Þ may be considered as a direction of potential descent, namely this gradient constitutes a descent indicator, which is used for ordering the poll vectors. In particular, the polling procedure will start by first testing the poll vectors that make the smallest angle with the negative simplex gradient. S2 Fig displays an

Assessment of a possible correlation of the invariant density measure with clinical effectiveness of stimulation in OCD
In order to assess a possible correlation of the invariant density measure with clinical effectiveness of stimulation in treatment-refractory OCD, we evaluated s 2 p st simulating the application of regular 130 Hz stimulation and based upon the model parameters estimated for two subsets of recordings: a total of 39 MERs of subthalamic neuronal activity acquired during DBS for OCD and characterized by a high mean discharge rate (39.7 ± 14.71 Hz), a high intraburst frequency and a short interburst interval (μ ISI = 0.0289 ± 0.0114s, Var ISI = 0.0038 ± 0.0056) vs. a total of 39 MERs of subthalamic neuronal activity characterized by a low mean discharge rate (13.53 ± 7.13Hz), a low intraburst frequency and a long interburst interval (μ ISI = 0.1072 ± 0.093s, Var ISI = 0.0265 ± 0.0542). This specific approach was based on indications correlating the efficacy of standard STN-DBS for OCD with locations of neuronal activity characterized by a high discharge rate and intraburst frequency, and a short interburst interval [50]. Statistical significance was determined by means of the Mann-Whitney U test. For assessment of the modulation index, after linear band-pass filtering the LFP signals between 13 and 30 Hz, the respective phase series is estimated by means of the Hilbert transform. As illustrated in Fig 2B, the derivative of the phase so estimated, i.e. the instantaneous angular frequency, is characterized by a high rate of singularities reflecting a high rate of artificial phase slips. At a singularity (or phase slip), the instantaneous angular frequency reaches values that exceed the limits imposed by the band-pass filter [71] (inset; Fig 2B). The high rate of artificial phase slips is probable to render the calculation of cross-frequency coupling unreliable, particularly in the presence of increased noise levels. On the contrary, a phase-reconstruction-free method does inherently not suffer from the ambiguity associated with phase singularities. This fact is demonstrated through four representative cases in Fig 2C. While results obtained by means of the modulation index and the test outcome are in good agreement in case of a relatively high signal to noise ratio (SNR) (cases P5, P6; Fig 2C), the modulation index fails to discriminate sites with significant non-linear coupling from sites without, in case of a low SNR (cases P7, P8; Fig 2C).

Noise-resistant assessment of nonlinear coupling
In the two representative cases of Fig 3A, employment of the 0-1 test for chaos following the application of the complex valued-filter, singled out sites with significant nonlinear coupling between beta and high-frequency activity (indicated by a test outcome smaller than 0.1). Conversely, following the application of a conventional Butterworth band-pass filter, the 0-1 test for chaos did not discriminate sites with significant non-linear coupling from sites without. This result was corroborated across the total of the MER trajectories examined and highlighted the robustness of the two-part technique, i.e., the combined application of the complex-valued filter and the 0-1 test for chaos. In case of PD, reduced sensitivity to measurement noise in the 0-1 test for chaos was warranted by assigning a positive value to parameter h in Eq (3) (h = 1, Fig 3B). On the contrary, this assignment did not prove to be a prerequisite in case of OCD (h = 0, Fig 3B).
Nonlinear coupling may be a reliable biomarker for feedback control in case of STN-DBS for PD. Cross-frequency coupling was identified at a total of 18 MERs-sites within the STN of 8 patients with PD (case P1: 2 sites; case P2: 1 site; case P3: 3 sites; case P4: 2 sites; case P5: 3 sites; case P6: 2 sites; case P7: 1 site; case P8: 4 sites). Approximately 67% of these sites (n = 12) was located at the dorsal border of the STN (Fig 3C). These results are rather predictable given that beta-HFO coupling is closely correlated with the pathophysiology of PD and strongest at the dorsal border of the STN [20][21][22][23][24]. They further corroborate the potential appropriateness of nonlinear coupling between beta and high-frequency neuronal activity as a biomarker for feedback control in PD. Neuronal activity at 13 out of the 18 sites with cross-frequency coupling followed a bursting or burst-like firing pattern (case P1: 1 site; case P2: 1 site; case P3: 3 sites; case P4: 1 site; case P5: 3 sites; case P6: 2 sites; case P7: 1 site; Comparative assessment of the modulation index and the 0-1 test outcome in four representative cases of PD. Results obtained by means of both measures are in good agreement in case of a relatively high signal to noise ratio (SNR) (cases P5 and P6: identification of significant cross-frequency coupling at +0.5 mm). On the contrary, in the presence of increased noise levels, the high rate of artificial phase slips renders the calculation of cross-frequency coupling by means of the modulation index unreliable. In particular, the index fails to discriminate sites with significant non-linear coupling from sites without, in case of a low SNR (cases P7 and P8: significant cross-frequency coupling at -3 mm and at -1.5 mm/+1 mm, respectively, identified only by means of the 0-1 test for chaos).
doi:10.1371/journal.pone.0171458.g002  (Eq (3)), based upon optimization of the outcome of the 0-1 test for chaos across a subset of 12 MER trajectories in PD and 12 MER trajectories in OCD. According to the results, in case of PD, sensitivity to measurement noise had to be further decreased by assigning a unitary value to parameter h. (C) Cross-frequency coupling was identified at at least 1 site within the STN of each patient with PD (total:18 MERs). Approximately 67% of these sites was located at the dorsal border of the STN, while at 72.2% of these sites neuronal activity followed a bursting or burst-like firing pattern and was considered for further processing in the phase-reduced bursting neuron model. Contrary to the case of PD, cross-frequency coupling was identified within the STN of only 2 patients with treatment-refractory OCD (total:4 MERs).
doi:10.1371/journal.pone.0171458.g003 case P8: 1 site). These sites were considered for further processing in the bursting neuron model. At the remaining 5 sites a rather irregular firing pattern was observed, and therefore these sites were excluded from subsequent analysis.
Nonlinear coupling may display subject-specific applicability as a biomarker for feedback control in case of STN-DBS for OCD. Contrary to the case of PD, cross-frequency coupling was identified at only 4 MERs-sites within the STN of 2 patients with OCD (case O2: 2 sites; case O3: 2 sites) (Fig 3C). The latter fact may be attributable to the lower number of acceptable MER trajectories in case of STN-DBS for OCD. Otherwise, it implies that nonlinear coupling between beta and high-frequency activity may not consistently be an appropriate biomarker for feedback control in closed-loop STN-DBS for treatment-refractory OCD [25] and that an alternative biomarker should, therefore, additionally be considered (Fig 1B). For this reason, on the basis of evidence pointing to a correlation of subthalamic bursting neuronal activity, characterized by certain features, with symptom severity and stimulation efficacy in OCD [50], we assessed, for the remaining of the cases wherein no cross-frequency coupling was identified, the presence of bursting neuronal activity with specific characteristics, i.e., a short interburst interval and a high intraburst frequency (μ ISI = 0.0242 ± 0.0113s, Var ISI = 0.0059 ± 0.0083). Specifically, we considered for further processing a total of 12 MERs (case O1: 1 site; case O2: 2 sites; case O3: 2 sites; case O4: 1 site; case O5: 1 site; case O6: 1 site; case O7: 2 sites; case O8: 2 sites).

Performance of the model-based control strategy in terms of efficiency, selectivity of stimulation and computational cost
The performance of the model-based direct search method (SID-PSM) in the determination of the optimal parameters of stimulation in cases of PD and OCD was compared with the performance of a non-model-based generalized pattern search method ( [72]; Matlab, Mathworks, Natick, MA), in terms of the resulting values of the invariant density, stimulation power and total computation time. At each site, we acquired five evaluations of optimal stimulation parameters by means of each distinct solver and assessed the respective mean values illustrated in Fig 4A-4C. In both optimization procedures, current amplitude was consistently maintained at its minimal value (I 0 = 0.001A). The optimal pulse width determined by means of the model-based direct search method in case of PD proved to be equal to 33.36 ± 1.06 μs (mean ± standard error mean) and in case of OCD, equal to 33.75 ± 1.29 μs (mean ± standard error mean) (Fig 4D). Given that pulse durations lower than 60 μs have been associated with increased selectivity of stimulation [68], this result indicates a potentially outstanding performance of the model-derived stimulation parameters in terms of selectivity of stimulation. We should further comment on the fact that, following the application of the model-based direct search method, the mean optimal stimulation frequency proved to be significantly higher in case of OCD compared with PD (p<0.05, Mann-Whitney U test), while a similar outcome was obtained with respect to the mean optimal pulse width and the mean optimal Poisson parameter ( Fig 4D). We suggest that differences in the underlying pathophysiology [49,50] may have led to the observed differences in optimal stimulation frequency in case of PD vs. treatment-refractory OCD. The derived mean optimal stimulation frequency of 39 ± 3.43 Hz in case of PD is in partial accordance with Class III evidence that STN-DBS at patient-specific frequencies in the range 31-100 Hz may improve motor Unified Parkinson's Disease Rating Scale (mUPDRS) scores equally effectively with high-frequency stimulation in patients with PD [73].
Statistical analysis corroborated a significantly higher performance of the model-based direct search method, in terms of both stimulation power and computation time corresponding to the optimal stimulation settings, compared with the non-model-based generalized pattern search method (p<0.0001, Mann-Whitney U test), while an almost equivalent effect was observed on the invariant density measure (p PD = 0.3299, p OCD = 0.4705, Mann -Whitney U test) (Fig 5).
The results corresponding to the stimulation settings determined by means of the modelbased control strategy (combined application of the stochastic phase-reduced model and the model-based direct search method) were further compared with the results obtained by simulating application (through employment of the stochastic phase-reduced model) of the stimulation settings determined post-operatively, during the last follow-up of patients having undergone STN-DBS for PD or OCD (Tables 1 and 2). The comparison was performed in terms of the values of the invariant density measure and stimulation power. Statistical analysis corroborated the ability of the model-based control strategy to identify stimulation settings that yield significantly lower values of the invariant density measure and stimulation power compared with the respective values acquired by simulating application of the stimulation settings determined post-operatively (p<0.0001, Mann-Whitney U test) (Fig 6A and 6B). This result combined with the reported possible correlation of the invariant density measure with  clinical effectiveness of stimulation in PD [40], but probably also in OCD (see next section), points to a potentially superior performance of the model-based stimulation parameters in terms of therapeutic and energy efficiency of stimulation. The differential desynchronizing effect on neuronal activity exerted by the model-based stimulation settings vs. the stimulation settings determined post-operatively is qualitatively reflected in the distinct form of the respective stochastic kernels (Fig 6C and 6D).  Possible correlation between the invariant density measure and clinical effectiveness of stimulation in OCD Fig 7 displays the results obtained by assessing the invariant density measure based on a total of 39 MERs of subthalamic neuronal activity acquired during DBS for OCD and characterized by a high discharge rate, a high intraburst frequency and a short interburst interval vs. a total of 39 MERs of subthalamic neuronal activity characterized by a low discharge rate, a low intraburst frequency and a long interburst interval. Remarkably, the desynchronizing effect of standard 130Hz stimulation proved to be significantly stronger in the former case compared with the latter (p<0.01, Mann-Whitney U test). This result points to a possible correlation of the invariant density measure with clinical effectiveness of stimulation in OCD, since values of this measure are proven to be lower at locations of neuronal activity that have been correlated with the best clinical outcome of STN-DBS for OCD [50]. Bikson et al. (2015) [74] remark: "Approaches using closed-loop stimulation are inherently state dependent and require computational neurostimulation." Elaborating on this concept and considering the implications of the current approach, we make the following two key observations: first, though evidence about the pathophysiology of medically refractory movement and neuropsychiatric disorders remains to date to a large extent inconclusive, a growing body of basic and clinical work supports the important role of nonlinear coupling between beta and high-frequency activity in the pathophysiology of PD [75], thereby pointing to a possible utility of this measure as a state biomarker in closed-loop neuromodulation approaches for PD. Nevertheless, any attempt to reliably assess this biomarker should be made by carefully considering the presence of strong internal and measurement noise in the recorded neural activity. In this study, we presented an innovative technique drawn from dynamical systems theory  guaranteeing low sensitivity to noise, and corroborated the presence of cross-frequency coupling in each case with advanced PD. Thereby, we provided indications for a possible appropriateness of this approach in optimization of closed loop neuromodulation systems. Second, throughout this paper, we attempted to provide compelling evidence for the critical role of computational neurostimulation in closed-loop identification of novel stimulation protocols [56,76]. The computational model employed operates on the principles of phase reduction and phase-resetting that are inherently characterized by simplicity and analytical tractability [48,[77][78][79], and further incorporates the dynamics of neuronal bursting activity that constitutes a hallmark of PD and OCD pathophysiology. In addition to the employment of the phase-reduced bursting neuron model, employment of direct search optimization based on quadratic modeling has significantly contributed to the performance of the presented approach. Crucially, since the quality of simplex gradients as approximations to the real derivatives of the objective function has been demonstrated in the non-smooth case [70], SID-PSM comprised an ideal choice for the problem under consideration, i.e., the non-smooth dynamics of the neuronal response to stimulation.

Discussion
In previous work, we provided important indications for the realistic substructure of the stochastic phase-reduced model and further highlighted a possible correlation of the invariant density measure with clinical effectiveness of stimulation in PD [40]. By extending the latter result to the case of treatment-refractory OCD, we here prove that the proposed model-based control strategy holds the potential to exhibit remarkable performance in terms of therapeutic and energy efficiency of stimulation for both pathologic conditions. By yielding a mean optimal pulse width equal to~33μs, the model-based control strategy may further achieve outstanding performance in terms of selectivity of stimulation. Importantly, application of modelbased direct search has been associated with a significantly higher computational speed compared with non-model-based derivative-free optimization.
There are several methodological considerations inherent to this work. Even though this study may provide a rigorous theoretical foundation for the design of a therapeutically-and energy-efficient closed-loop neuromodulation system, it should be pointed out that the stochastic model employed here does not capture large-scale neuronal interactions within key circuits involved in the pathophysiology of PD or OCD. Nonetheless, similar approaches have led to valid interpretations of the neuronal response to stimulation [80]. Another consideration is related to the potential effect of dopaminergic treatment on the performance of the proposed methodology. This assessment was not feasible in the current study, since the available data were acquired after a long period of withdrawal of medication. Third, clinical validation of the presented predictions should comprise a basic priority in future studies. Clinical validation may be possible once novel neural probes or systems offering the capability of concomitant DBS and microelectrode recording are introduced in clinical practice [81]. We also emphasize the necessity to test whether the findings of this study are replicated in macroelectrode recordings, and whether the model-based predictions are validated clinically based on current neuromodulation technologies [82]. Furthermore, modeling approaches similar to those proposed in this study may display greater fidelity in the framework of constant current stimulation. The transition from the use of constant-voltage to constant-current DBS devices is being motivated by the rationale that constant-current stimulation will accommodate for interpatient and temporal fluctuations in the impedance of the tissue and electrode-tissue interface [83,84].
Last, special mention should be made of further important factors for efficient closed-loop neuromodulation. These include the optimization of the stimulation waveform [85], of the circuit topology [86,87] and of contact selection [88], as well as the incorporation of neurochemical control [89,90] and the adjustment of closed-loop DBS systems to the phenotypical heterogeneity of movement and neuropsychiatric disorders [91,92]. Cumulative research towards these directions may ultimately favor the optimization of less invasive, groundbreaking treatment options including closed-loop optogenetic control [93,94].