In silico development and validation of Bayesian methods for optimizing deep brain stimulation to enhance cognitive control

Testing a wide range of optimization strategies directly in patients undergoing VCVS DBS is impractical. It would require thousands of task trials across hundreds of runs, which is unlikely to be tolerated by most patients. Thus, we developed a simulation testbed for our cognitive control paradigm. The simulation is meant to model a patient continuously performing the MSIT (or a similar cognitive control assay), while stimulation is adjusted to improve task performance (to minimize RTs without inducing errors). The overall system (figure 1(A)) includes a data generator, a sensor, and an optimizer. The generator models a patient performing the MSIT and emitting behavioral data, in the form of RT values that are simultaneously influenced by stimulation, task factors, and noise processes. The sensor models the real-time estimation we would then need to do with a patient, analogous to the processing in [37]. It attempts to filter out the noise and task processes to identify stimulation-driven RT changes, without explicit information about the current stimulation settings. This smoothed derivative of RT (sometimes referred to as a 'cognitive state') is then used by the optimization algorithm to update its distribution estimate. Based on the update, the optimizer then modifies the simulated stimulation applied to the generator. This process continues until the algorithm converges.

Zoom In Zoom Out Reset image size

In the current work, we focused on the problem described in [37], optimizing the stimulation location (choosing among available DBS contacts or combinations of contacts), assuming that other parameters such as amplitude, pulse width, and frequency are already reasonable. We had empirical data to constrain this problem (see below), making it a logical starting point.

The dataset includes six participants with long-standing pharmaco-resistant epilepsy. They were reported as participants 8–13 in [37]. We relabel them here as S1–S6. Participants voluntarily enrolled with fully informed consent (obtained by a member of the study staff who was not the participant's primary clinician) in accordance with guidelines and procedures approved by the local institutional review boards at Partners Healthcare (Massachusetts General Hospital), with secondary review from the US Army Human Research Protections Office.

Participants performed the MSIT with simultaneous behavior and LFP recordings. Stimuli for the MSIT were presented on a computer screen with either Presentation software (neurobehavioral systems) or Psychophysics Toolbox. Each participant performed 1–3 sessions, and each session consisted of multiple blocks of 32 or 64 trials, (ranging between 5–10 blocks) with brief rest periods (median break time of 5.25 min) in between the blocks. They were instructed to be as fast and as accurate as possible. Median success rates of 100% ± 2.47% and 97.1 ± 5.52% were reported during congruent and incongruent conditions, respectively [37]. The task contained a roughly equal number of congruent and incongruent trials. Stimuli were presented for 1.75 s, with an inter-trial interval randomly jittered within 2–4 s. Simultaneously, LFPs were recorded using implanted depth electrodes; these data were not considered in the current study. The electrode locations were determined by a team of caregivers solely on clinical grounds and were not in any way modified for the research.

Only one site was stimulated during each MSIT block. Stimulation was a 600 ms long train of symmetric biphasic (charge balanced) 2–4 mA, with 90 μs square pulses at a frequency of 130 Hz, delivered using a Cerestim 96 (Blackrock Instruments). Stimulation was delivered through a neighboring pair of contacts on a single depth electrode (bipolar) with parameters set manually by the experimenter. Depth electrodes (Ad-Tech Medical or PMT) had a diameter of 0.8–1.0 mm and had 8–16 contacts (platinum/iridium), each 1–2.4 mm long. The stimulations were triggered by a separate PC that was either delivering or monitoring the task or behavioral stimuli. All stimulations were delivered at the image onset to influence a decision-making process that begins with that onset. Participants S1–S6 completed 343, 378, 320, 383, and 440 trials, respectively, with 233, 254, 192, 255, 224, and 249 trials of stimulation and 110, 124, 128, 128, 159, and 191 trials of non-stimulation. S1-S6 performed 6, 6, 5, 8, 9, and 10 blocks respectively. All stimulations were given to varying sites within the internal capsule as described in [37] where sites were chosen to have one site each in the ventral and the dorsal IC. Up to four sites were tested per patient, with 1–3 blocks per tested site. All stimulation was administered at the image's onset in order to impact a decision-making process that began at that time. Prior to task-linked stimulation, all participants were tested at 1, 2, and 4 mA for 1 s, repeated five times, with 5–10 s between each 130 Hz pulse train. Participants were told that stimulation was active and asked to describe any acute perceptions. The absence of epileptiform after-discharges was confirmed, ensuring that the individuals could not detect stimulation, such as odd feelings. If participants reported any sensation (for example, tactile sensations in the limbs or head), task-linked stimulation was reduced to the next lowest intensity. Note that although this experiment used intermittent IC stimulation, we have previously shown that the same effects hold with continuous IC stimulation as would be used in clinical DBS [36].

The MSIT is described extensively in [36, 37]. In brief, it is a standard cognitive control/conflict paradigm where roughly 50% of trials contain distractors that induce pre-potent (incorrect) responses. These are referred to as 'incongruent', 'interference', or 'high conflict' trials. Successful task performance requires suppression of those responses, both through reactive processes (e.g. stopping an incorrect action once started) and proactive processes (e.g. engaging attentional systems to detect impending errors). In this context, task RTs measure a participant's ability to apply cognitive control—faster RTs without errors imply more efficient control systems.

Task RTs are positive random variables with right-skewed distributions [49, 50]. In particular, RTs are often well described by Gamma distributions, with parameters that vary depending on the trial type (e.g. the high vs. low conflict trials in MSIT). Thus, we model the RT of a participant performing MSIT as described in [50],

$$\begin{align}{Z_{{\text{RT}}}}\sim \Gamma \left( {\mu ,v} \right) + \alpha \end{align} \tag{ 1.1 }$$

where ZRT is the observed RT, α is the offset term representing the lower bound of the RT ($\alpha \geqslant 0$), 1/v is the dispersion observed in the RT, and μ is the mean RT. Here, the rate parameter for the Gamma distribution will be, $\beta = v/\mu$ and its shape parameter is defined by, $\alpha = v$.

That mean RT varies both systematically and stochastically. Systematically, it is higher on high conflict trials. Stochastically, it has both random variation (captured in $v$) and slow change over time. The latter may be due to learning processes (participants get slightly better at the task over the first 100–200 trials) or, more commonly, due to the effect of brain stimulation that alters the participant's 'baseline' RT. As in [37], we express this in terms of two unobservable 'cognitive states', ${x_{{\text{base}}}}$ and ${x_{{\text{conflict}}}}$, connected to the RT by a log link function [50, 51],

$$\begin{align}{\text{ln(}}\mu ) = {b_0} + {b_1}{x_{{\text{base}}}} + {b_2}{I_{{\text{conflict}}}}{x_{{\text{conflict}}}}^{}\end{align} \tag{ 1.2 }$$

${x_{{\text{base}}}}$ is the baseline cognitive state, which represents the mean RT when no conflict is present. ${x_{{\text{conflict}}}}$ is the mean amount of increase in RT during high conflict. ${I_{{\text{conflict}}}}$ is an indicator variable representing presence or absence of conflict in a given trial of MSIT. b₀, b₁, b₂ are then weights (effectively, regression coefficients) that relate cognitive states to the distribution mean μ. b₀ provides a bound on RT, whereas b₁, b₂ are scaling factors. We have demonstrated this model as well suited to MSIT datasets in multiple prior papers [37, 48, 51].

The temporal evolution of cognitive states at trial k is defined by a first order autoregressive model, AR(1), as in [50, 51].

$$\begin{align}{x_{{\text{base}},k + 1}} = {a_1}{x_{{\text{base}},k}} + {b_{in1}}{u_m} + {w_{1,k}}\end{align} \tag{ 2.1 }$$

$$\begin{align}{x_{{\text{conflict}},k + 1}} = {a_2}{x_{{\text{conflict}},k}} + {b_{in2}}{u_m} + {w_{2,k}}\end{align} \tag{ 2.2 }$$

$$\begin{equation*}{w_{1,k}}\sim N\left( {0,{\sigma _1}} \right)\end{equation*}$$

$$\begin{equation*}{w_{2,k}}\sim N\left( {0,{\sigma _2}} \right).\end{equation*}$$

The next states ${x_{{\text{base}},k + 1}}$, ${x_{{\text{conflict}},k + 1}}$ are dependent on the current states ${x_{{\text{base}},k}}$, ${x_{{\text{conflict}},k}}$, the input ${u_m}_{}$ (which represents the presence of applied DBS at one of the m stimulation sites), and the Gaussian state noise processes $w$1, $k$, $w$2, $w$ with zero mean and variance σ1, σ2 respectively.

Here, we simulated a case analogous to our empirical experiments in [37], where we tested different contacts along a cylindrical lead while holding other stimulation parameters (amplitude, pulse width) constant. In this simulation, the input vector ${u_m}$ has length equal to the number of available contacts, all of which are assumed to be stimulated in monopolar mode (as in a clinical monopolar survey). The elements of u are 1 when a given contact is active, and 0 when it is not. The weighting vectors ${b_{in1}},{b_{in2}}$ then describe how stimulation alters the state variables. Note that in this formulation, effective DBS changes these underlying generating states, but does not act directly on the mean RT $\mu$. Thus, when a new stimulation regime is applied, the states ${x_{{\text{base}}}}$ and ${x_{{\text{conflict}}}}$ will not immediately change but will drift to a new steady state. This is consistent with our observations in [37]. The rate of drift will depend on a and ${b_{{\text{in}}}}$. As further described below, this influences the design of the optimization approach.

The state space model of (1) and (2) is applied to generate new samples of behavior (MSIT RTs) that model those we would expect to observe from a patient undergoing DBS. The evolution of behavior is represented by stepping (2.1) and (2.2) forward in time. The modeling requires identifying reasonable values for the parameters in these equations.

We perform parameter and state estimation using the COMPASS toolkit [48], which uses a combination of an expectation maximization (EM) algorithm and a Kalman-type filtering/smoothing approach to simultaneously infer model parameters and cognitive states given the observed behavior in a block of trials.

The parameter estimation problem for MSIT as described in (1) and (2) would be

$$\begin{align}\phi = \left( {\alpha ,v,{b_0},{b_1},{b_2},{a_1},{a_2},{w_1},{w_2},{b_{in1}},{b_{in2}}} \right).\end{align} \tag{ 3 }$$

This forms a complex estimation problem with a large number of unknown parameters. This increases the possibility of converging to a local maximum, resulting in sub optimal estimates. To prevent this, we reduce the parameter estimation search space. First, we temporarily ignore the input weightings ${b_{in1}}$ and ${b_{in2}}$, by setting them to 0. As described below in section 'Optimizer test with simulated stimulation response surfaces', we will assign new values to these parameters when performing the simulation, making it irrelevant to estimate them from patients' data. Note that we can still use data from patients receiving stimulation (as in [37]) to estimate the model—any drift in ${x_{{\text{base}}}}$ and ${x_{{\text{conflict}}}}$ caused by stimulation is captured by the dispersion 1/v and state Gaussian noise W. For convenience, here the generator models were estimated only from non-stimulation blocks within the empirical dataset.

As a further simplification, we apply no scaling to the cognitive states, which is done by setting scaling parameters ${b_1}$ and ${b_2}$to 1. Finally, we recognize that the state noise process $W$ is the most challenging free parameter to estimate, as a wide range of values can be consistent with observed data. We therefore assign values to ${\sigma _1}$ and ${\sigma _2}$ that we identify as working well across patients; these are described further in 'Grid search for model noise processes' below. This simplifies the parameter search space to

$$\begin{align}{\phi _{{\text{generator}}}} = \left( {\alpha ,v,{b_0}} \right).\end{align} \tag{ 4 }$$

The problem (4) is then estimated using the COMPASS EM algorithm. The convergence of the EM algorithm is determined by assessing the maximum likelihood (ML) of the data at each iteration, given the current parameters ${\phi _{generator}}$. The convergence criterion for EM algorithm is a plateau in the ML change given by (4.1)

$$\begin{align}\left| {\frac{{M{L_t} - M{L_{t - 1}}}}{{M{L_{t - 1}}}}} \right| < 0.001.\end{align} \tag{ 4.1 }$$

To demonstrate the models' ability to simulate new participants that are very similar to the empirical data, 1000 new models ('generator pool') were estimated based on each participant's data (S1–S6). For each model fit, the initial conditions ${x_0}$ in (2) were selected randomly from a normal distribution N(0,1), autoregressive parameters ${a_{}}$ were selected close to 1 [${a_1}$ = 0.9999, ${a_2}$ = 0.9999] and for $W$ values for ${\sigma _1}$ and ${\sigma _2}$ were randomly selected from the scale range [7, 16] and [13, 27] respectively identified in 'Model noise process determination via grid search'.

To validate the generator model, we computed the Kolmogorov–Smirnov (KS) distance between the empirical distribution of those simulated RTs and the empirical probability distribution function (pdf) of a gamma distribution fit to the simulated RTs (MATLAB fitdist function). A KS distance/test at p > 0.05 would suggest that the generator is producing RT-like, gamma-distributed data. We performed two versions of the simulation, one with model tuning and other without model tuning. Model tuning was performed to demonstrate that, if necessary, the generator model could closely replicate the actual distribution obtained from any given empirical participant. Model tuning was a trial and error method of modifying model parameters A, v, w and $\alpha$ by hand. We simulated both the tuned model and untuned model each for 1000 trials as described above. These simulations randomly generated interference trials by setting ${I_{{\text{conflict}}}}$ to 0 or 1 on each trial. No stimulation effects were included. This was repeated 1000 times for each participant, followed by the KS test as just noted. We report the fraction of simulations that had p > 0.05. Further, a KS distance/test was performed between the simulated RT from tuned models and empirical data to verify if the models can be tuned to generate distribution similar to the empirical data. Again we reported with a fraction of simulation that had p > 0.05.

In actual patients, we do not have direct access to their underlying cognitive states but infer these by fitting the model of (1) and (2) to the observed RTs. In our simulations, however, the ground truth is known, and we can confirm that our sensor model accurately recovers it. The sensor model fitting problem, analogous to the generator problem, is given by

$$\begin{align}{\phi _{{\text{sensor}}}} = \left( {\alpha ,v,{b_0}} \right).\end{align} \tag{ 5 }$$

The difference between the sensor and the generator model is the data used for modeling and evaluation. For the generator model, we fit it and tested its deviance against the empirical data from S1–S6. For the sensor model, we used data simulated from a pool of randomly initialized generators. We used the same COMPASS EM algorithm, with convergence rule (4.1), and ${\sigma _1}$, ${\sigma _2}$ drawn from the scale range [7, 12] and [13, 27]; the origin of this choice is described below in 'Grid search for model noise processes'. Further, we added 0.0625 (equal to scale parameter 8) to both ${\sigma _1}$ and ${\sigma _2}$. This slight increase in noise covariance is helpful in later modeling steps, when we attempt to capture the effect of stimulation. Because the sensor model does not explicitly represent stimulation, a slightly higher state noise leads to more rapid adaptation when stimulation is applied and the state variable is forced away from its default AR(1) behavior.

For each sensor model, the initial state ${x_0}$ in (2) was randomly selected from a normal distribution N(0,1). To further accommodate a wide range of generators, we set the lower bound $\alpha$ of RT to 0.1. We simulated 1000 fits of a sensor model to a generator model. The generator for each simulation was chosen from the 'generator pool' described above under 'Generator model'. We then evaluated sensor tracking ability on this randomly selected generator model using normalized root mean squared error (NRMSE), as described in 'Sensor tracking validation' below.

For the problem described in (3), a large number of unknown parameters increases the probability of incorrect convergence of the EM algorithm. In initial pilots, we noticed that this was particularly problematic when attempting to estimate the covariance of the state noise process W. (Here, W is a diagonal matrix containing parameters ${\sigma _1},\,{\sigma _2}$). This estimation often converged to different values, with diagonal elements between [${10^{ - 1}}$, ${10^{ - 7}}$]. All other parameters in $\phi$ would be altered in turn. This wide range of W values is not reasonable on face value, suggesting a need for a constraint. We theorized that there might be acceptable values of W that may work well across all patients, simplifying the problem of (3)–(5).

We tested this theory by conducting a two-dimensional grid search using a logarithmically spaced grid. Grid search was selected above other approaches such as random, gradient descent, or Bayesian search because it performs an exhaustive search for the underlying hyperparameter W. The advantages of the other search methods is faster convergence, but for an offline problem with no time constraints, grid search is the better choice. Further grid search has the advantage of revealing information on the underlying unknown structure of the response surface as a function of W. At each grid point, we set the diagonal elements of W to the specified values, estimated ${\phi _{{\text{generator}}}}$ using COMPASS, then calculated the deviance between observed RTs and those predicted by the model (μ from equation (1)). The logarithmic grid is given by

$$\begin{align}{\text{W}} = \left| {\begin{array}{*{20}{c}} {\frac{1}{{{{\sqrt 2 }^{{\text{scal}}{{\text{e}}_1}}}}}}&0 \\ 0&{\frac{1}{{{{\sqrt 2 }^{{\text{scal}}{{\text{e}}_2}}}}}} \end{array}} \right|\end{align} \tag{ 6 }$$

where ${\text{scal}}{{\text{e}}_1}\epsilon$ [1, 40] and ${\text{scal}}{{\text{e}}_2}\epsilon$ [1, 40]. The range was selected empirically based on values observed in previous studies [37]. Numerically, ${\sigma _1},{\sigma _2}$ takes values from [9.5367 × 10⁻⁰⁷, 0.707] for logarithmic search. For each point in the grid search, the problem (4) was estimated by fitting to data including both non-stimulation and stimulation trials from the empirical dataset using the COMPASS EM algorithm. We then evaluated the deviance between observed RTs and the one predicted by the model. This fit was performed separately for each dataset for participants S1–S6. It was repeated five times for each participant, using a random ${x_0}$ selected from N(0,1), plus an additional run with initial condition ${x_0}$ set to 0. For this model inference we set ${A_k}$ close to 1 [${a_1}$ = 0.9999, ${a_2}$ = 0.9999], $W$ was set based on the current grid point, b was 0 and the overall distribution was assumed to be gamma.

In previous reports, we validated sensor models by showing that the residuals of predicted RTs were consistent with a white noise process [37, 51]. Here, because we had access to the ground truth generator internal states, we directly measured whether the sensor model accurately tracked those states. We calculated the NRMSE between the ${x_{{\text{base}}}}$ and ${x_{{\text{conflict}}}}$ values in the generator and those inferred by the sensor. NRMSE is the root mean square error between the generator and sensor states after both variables have been re-scaled multiplicatively to the [−1, 1] interval. Because our optimization is concerned primarily with changes in a state variable relative to its own baseline, and not with the specific numeric value of that variable, NRMSE is a more appropriate metric. It tracks whether the sensor's estimate changes in the same direction and scale as the generator ground truth.

NRMSE was calculated between the sensor model (predicted RT and cognitive states) and the corresponding ground truth from the generator. We selected generators randomly from the generator pool (1000 models from 6 participants as noted above). To ensure that we demonstrated that the sensors tracked well over all assumed noise process values, we enforced a further constraint that ${\sigma _1}$ in the generator be drawn uniformly from the selected scale range of [7, 12]. The generator-sensor models were simulated for 1000 trials each, for 1000 randomly selected generators and sensors. During these trials, interference variable ${I_{{\text{conflict}}}}$ was set randomly on each trial, and the value of this indicator was visible to both sensor and generator.

To illustrate the sensor's timescale for settling to a steady state when stimulation is changed, we also simulated the two-site, high contrast stimulation ${b_{{\text{SS}}}}$ as described in 'Optimizer test with simulated stimulation response surfaces'. This simulation also used 1000 trials, but stimulation was changed every 25 trials between the two elements of ${b_{{\text{SS}}}}$.

Identifying the optimal stimulation contact on a cylindrical DBS lead, assuming that other parameters are held constant during this survey, is a multi arm bandit (MAB) problem. MABs consist of competing discrete options (here, k independent stimulation sites from available stimulation sites in u), of which one of the options produces the maximum effect (here, reduction in the RT). Common algorithms for solving MABs include a greedy algorithm [52, 53], $\epsilon$-greedy algorithm [52, 53], UCB1 [52], Bayesian upper confidence bound (Bayes-UCB) [54]and Thompson sampling (TS) [53].

In general, the goal of these algorithms is to minimize the approximate cost function $Q\left( u \right)$ by performing a trade off between exploring (selecting arms with uncertain payoff) and exploiting (selecting the arm currently estimated to have the best payoff). Algorithms differ in the cost functions and the steps they take to perform the trade off between exploration-exploitation, as described below.

2.9.1.  Greedy and $\epsilon$ -greedy algorithm

The pure greedy algorithm solely exploits and is a degenerate case of $\epsilon$-greedy with $\epsilon$ = 0. The $\epsilon$-greedy algorithm is an improvement to the greedy algorithm which primarily exploits, but also explores a new stimulation site ${u_k}$ randomly with a small probability $\epsilon$ (here selected as 0.1). The cost function $Q\left( u \right)$ is adapted from [52]

$$\begin{align} Q\left( u \right) & = \frac{1}{{{N_t}\left( u \right)}}\mathop \sum \limits_{t = 1}^T RT_{}^t \nonumber\\[8pt] {u_k}^t & = \frac{{argmin}}{{{u_k} \in u}}\left( {Q\left( u \right)} \right)\,\,\,p = 1 - \epsilon \nonumber\\[8pt] {u_k}^t & = random\left( u \right)p \geqslant \epsilon \end{align} \tag{ 7 }$$

where u is the input vector to (2), with 1 indicating the selected site ${u_k}$, ${N_t}\left( u \right)$ is the number of times a site is selected, T is the overall number of trials done so far, p is a random probability drawn uniformly from [0,1], and ${u_k}^t$ is the stimulation site k selected at time t based on cost function.

2.9.2. UCB1

UCB1 explores and exploits based on the Hoeffding inequality [52] that accounts for the number of times a particular stimulation site ${u_k}$ has been selected. It maintains an estimate of each site's worst-case performance, with a confidence interval around that estimate that shrinks as the site is selected more. Over time, as all sites have been repeatedly selected, UCB1 becomes less exploratory and more exploitative. The cost function [52] for UCB1 is

$$\begin{align} Q\left( u \right) & = {\mu _t}\left( u \right) - \sqrt {\frac{{2ln\left( T \right)}}{{{N_t}\left( u \right)}}} \nonumber\\[8pt] {u_k}^t & = \frac{{{\text{argmin}}}}{{{u_k} \in u}}\left( {Q\left( u \right)} \right) \end{align} \tag{ 8 }$$

where ${\mu _t}\left( u \right)$ is the vector of the mean RT of stimulation sites in u.

2.9.3. Bayesian inference and conjugate prior

Bayesian algorithms differ from UCB1 in that UCB1 directly maintains an estimate of the distribution of observed RTs, whereas Bayesian algorithms maintain an estimate of the parameters defining that distribution.

The exploration-exploitation of stimulation sites is based on these belief dynamics about the quality of each available site. This belief is expressed in terms of prior ${\Pi^o}$ and posterior distributions ${\Pi^t}$. The ${\Pi^o}$ is the probability distribution before the RT is observed at t, whereas the ${\Pi^t}$ is the probability distribution after observation at t. On each trial t, a random hypothetical RT ${h^t}$ is drawn for each stimulation site ${u_k}$in u, using the prior estimate of the distribution parameters ${\theta _k}$, and the site with the best draw is selected by (9)

$$\begin{align} {h^t}_{} & = X\left( {{\theta _k}} \right){\text{for}}\,{\text{each}}\,k \nonumber\\[4pt] u & = {\text{min}}\left( {{h^t}} \right) .\end{align} \tag{ 9 }$$

The Bayesian algorithms assume the RT observed on each trial to be an independent identically distributed draw. With each observation, they perform a posterior update of the parameter belief to obtain a posterior belief ${\Pi^t}$, following Bayes' rule. This ${\Pi^t}$ then is the prior for the subsequent trial at t + 1.

Prior and posterior distributions are selected from conjugate pairs. For conjugate distributions the prior (conjugate prior) and the posterior distribution fall in the same probability distribution family [53]. Here, the posterior update is performed in conjunction with a rule described in (10), which transforms the observed RT into a reward/payoff metric for the action of selecting stimulation site ${u_k}$.

For any arbitrary distribution parameters ${\theta _1},$ ${\theta _2} \in {\theta _n}$, if $R{T_k}^t \geqslant R{T_{val}}$: update ${\theta _1}$; else: update ${\theta _2}$

$$\begin{align} {\text{Where }}R{T_{{\text{val}}}} & = R{T^{t - 1}}\,{\text{if }}{u_k}^t \ne {u_k}^{t - 1} \nonumber\\[4pt] & = {\text{ min}}\,{\text{(}}\mu \left( {{u_k}^{^{\prime}}} \right))\,{\text{if }}{u_k}{u_k}^t = {u_k}^{t - 1} \end{align} \tag{ 10 }$$

where ${u_k}^{^{\prime}} \subseteq u|{u_{k}} \ne {u_t}$. (10) states that if the current and previously selected stimulation site are different, then the threshold $R{T_{{\text{val}}}}$ is the previous $R{T^{t - 1}}$ Otherwise, $R{T_{{\text{val}}}}$ is the minimum RT among the mean RTs of the rest of the stimulation sites. The parameters ${\theta _1}$ and ${\theta _2}$ depend on the specific choice of conjugate prior. The relevant distributions and their update rules are expressed in table 1.

Table 1. Conjugate prior and posterior distributions and update rules for Bayesian multi-arm bandit algorithms.

Algorithms	Likelihood	Conjugate prior ${\Pi^0}$	Posterior hyperparameters update ${\Pi^t}$
			RT $\geqslant$ X	RT $<$ X
Bayes-UCB	Normal RT ∼ N($\mu ,{\sigma ^2}$)	Gamma $\left( {\alpha ,\beta } \right)$	$\alpha = \alpha$ + $\frac{n}{2}$	$\beta = \beta + \frac{{\mathop \sum \nolimits _{t = 1}^T{{(R{T^t} - \mu )}^2}}}{2}$
TS-Bernoulli	Bernoulli ${r_t}\sim Ber\left( p \right)$	Beta ($\alpha ,\beta )$	$\alpha = \alpha + {r_t}$ ${r_{t}} \subset \left[ {0,1} \right]$	$\beta = \beta + 1 - {r_t}$ ${r_{t}} \subset \left[ {0,1} \right]$
TS-Poisson	Poisson ${r_{t}}\sim Pois\left( \lambda \right)$	Gamma $\left( {\alpha ,\beta } \right)$	$\alpha = \alpha + R{T^t}$	$\beta = \beta + 1 - {r_t}$ ${r_{t}} \subset \left[ {0,1} \right]$
TS-Normal	Normal RT ∼ N($\mu$,$\frac{1}{\tau }$)	Gamma $\left( {\alpha ,\beta } \right)$	$\alpha = \alpha$ + $\frac{n}{2}$	$\beta = \beta + \frac{{\mathop \sum \nolimits _{t = 1}^T{{(R{T^t} - \mu )}^2}}}{2}$
C-TS	Normal RT ∼ N($\mu ,{\sigma ^2}$)	Normal (${\mu _z}$, $\frac{1}{{1 + {k_z}}}$)	${\mu _z} = \frac{{{\mu _z}{k_z} + R{T^t}}}{{{k_{z}} + 1}}$	${k_z} = {k_z} + 1$

2.9.4. Bayesian UCB

Bayes-UCB assumes that the process generating the observed RT is normally distributed, where the parameters mean $\mu$ and variance $\sigma$ for that distribution are unknown. A gamma distribution $\Gamma \left( {\alpha ,\beta } \right)$ is an appropriate choice of conjugate prior (${\Pi^o}$) for such a normally distributed process [54]. Prior to each trial t, Bayes-UCB estimates $\mu$ and $\sigma$ using the prior shape ($\alpha$) and rate ($\beta$) hyperparameters (11). It then calculates a cost function/hypothetical RT ${h^t}$ to select the best stimulation according to the current belief. Once the actual RT is observed, as per Bayes theorem the prior belief is updated using posterior update ${\Pi^t} = \Gamma \left( {\alpha + \frac{n}{2},\beta + \frac{{\mathop \sum \nolimits _{t = 1}^T{{(R{T_t} - \mu )}^2}}}{2}} \right)$ in combination with (10),

$$\begin{align} \mu & = \frac{\alpha }{\beta }\,,\sigma \, = \,\frac{\alpha }{{{\beta ^2}}}, \nonumber\\[6pt] {h^t} & = {\mu _t}^*\left( u \right) - c*\frac{\sigma }{{\sqrt {{N_t}\left( u \right)} }} \nonumber\\[6pt] {u_k} & = \frac{{argmin}}{{{u_k}\epsilon u}}\left( {{h^t}} \right) \end{align} \tag{ 11 }$$

where c is a confidence level value (here selected as 1.96 corresponding to 95% confidence level) which adds to the exploration-exploitation. The margin of error $\frac{\sigma }{{\sqrt {{N_t}\left( u \right)} }}$ decreases as the sample size increases and algorithm is more certain about the mean ${\mu _t}\left( u \right)$.

2.9.5. TS

2.9.6. Improper priors

2.9.7. Brute force (BF)

The performance of a given optimizer depends on the problem structure. For instance, if there is a single obvious global best setting, the greedy and BF algorithms will suffice. If the response surface is shallow, with local minima very close to the global minimum, the more complex optimizers are more likely necessary. Greedy and BF algorithms also are expected to scale less well as the number of options (stimulation sites) expands, e.g. as DBS leads move from legacy quadripolar designs to modern multi-contact segmented designs [57, 58].

We thus simulated a range of increasingly complex optimization problems, from 2 to 8 stimulation sites (table 2). We modeled the effects of individual sites after the empirical effects observed in [37]. The source data for changes in the RT due to the stimulation are shown in supplementary table 3(a) of [37]. We assumed that all sites either decrease ${x_{{\text{base}}}}$(the desired effect) or have no effect. Table 2 gives the values for ${b_{{\text{in}}}}$ 's effect on ${x_{{\text{base}}}}$; the effect on ${x_{{\text{conflict}}}}$ was assumed to be 10-fold smaller (consistent with the empirical findings that VCVS stimulation mainly affects ${x_{{\text{base}}}}$). For each re-simulation of the problem (see below), the order of these sites was randomly shuffled.

Another key performance factor for an optimization algorithm is its ability to distinguish between two closely spaced minima. We therefore also simulated the two-site case, but with increasingly small separation between the two options, all the way to a degenerate case with no difference between options. These values b1–b7 are given in table 3.

Table 3. Stimulation effects for two stimulation site locations for the optimization problems, with decreasing differences between the two sites. ${b_{{\text{SS}}}}$ is a special case used only for demonstrating steady-state response/settling of the sensor model. It is specifically chosen to have very large and visibly different effects.

Stimulation sites	Stimulation effect		Stimulation sites	Stimulation effect
	1	2		1	2
b1	0.05	−0.05	b5	0	−0.02
b2	0	−0.10	b6	0	−0.01
b3	0	−0.05	b7	0	0
b4	0	−0.03	${b_{ss}}$	0.25	−0.25

The primary metric for optimizer evaluation is accuracy, defined as the fraction of times that a given algorithm identifies the known best site.

Optimizer accuracy may be sensitive to the steady-state settling behavior of the specific problem, the convergence rules, and aspects of the problem. We tested the above algorithms across a range of assumptions (table 4), including:

2.11.1. Effect of block-wise vs. trial-wise stimulation

2.11.2. Effect of convergence criterion

2.11.3. Effect of number of stimulation sites

2.11.4. Effect of state noise/signal to noise ratio (SNR)

The SNR influences how well bandit algorithms perform. A noisy system would produce noisy RTs, and it would be difficult for the algorithm to identify any changes caused by the stimulation effect. We evaluated how different noise floors affect the algorithm accuracy. For this simulation, rather than using the values from the grid search, we selected $W$ (${\sigma _1}$ and ${\sigma _2}$) as described in table 5. These values were perturbed by N(0, 0.0001) in order to use a slightly different generator model in each simulation. These values correspond to different SNRs, expressed as the ratio of the largest stimulation effect in b to the size of the noise variance. We simulated each SNR for 1000 replicates, using 600 trials per replicate, a blocksize of 15, and 8 stimulation sites as described in table 2.

Table 5. Different state noise levels. ${b_{{\text{diff}}}}$ represents the difference between the most effective and second most effective stimulation site.

$\sigma$	0.5	0.4	0.3	0.2	0.1	0.05	0.04	0.03	0.02	0.01	0.005
$\frac{{{b_{diff}}}}{\sigma }$	0.06	0.075	0.1	0.15	0.3	0.6	0.75	1	1.5	3	6
dB	−24.44	−22.5	−20	−16.48	−10.46	−4.44	−2.5	0	3.5	9.54	15.56

2.11.5. Discriminability of different stimulation effects

Even the best optimization algorithm can still converge away from the global optimum, particularly if the measurement (RT) is noisy. This might be mitigated by an ensemble approach, where optimization is performed repeatedly from different initial conditions, or simply by allowing the random choices inherent in the algorithms to proceed from different seeds. Clinically, this could be realized by redoing the optimization with the same patient on multiple days, as has been suggested in [34]. We simulated this ensemble case, by repeatedly re-simulating the same generator and different sensor model. We tested ensembles of 1, 2, 3, 4, 5, 6 and 7 repetitions over 1000 generator-sensor pairs selected as described above, with 600 trials per replicate, blocksize 15, stimulation problem of size 2, 4, 6 and 8, and noise level as defined in the 'Grid search for model noise processes' section. The selection for computing accuracy was then based on a hard majority vote across repetitions.

Modeling includes estimation of parameters using the COMPASS EM algorithm. The full parameter estimation problem given in (3) has many unknowns, which in practice often leads to convergence on a local maximum. This is particularly problematic when estimating the state noise W, leading to the simplifications in (4) and (5). To enable those simplifications, we performed a grid search to identify reasonable fixed/cross-patient values for W, seeking values that minimize the deviance between actual and predicted data during the expectation step of EM.

Figure 2 shows the deviance for the state-space model of (1) and (2), fitted using different noise variances for ${x_{{\text{base}}}}$ and ${x_{{\text{conflict}}}}$. A region of values corresponding to [4, 19] (${\sigma _1}$ ≡ [0.0014,0.25]) for ${x_{{\text{base}}}}$ and [2, 34] (${\sigma _2}$ ≡ [0.000007629,0.5]) for ${x_{{\text{conflict}}}}$ produced small deviance on average (figure 2(A)). Due to the unavailability of ground truth about optimal scale values, the specific best point within this region cannot be determined. In practice, any W drawn from within the low-deviance zones would likely work well. We verified that this region of low deviance is present in single-participant data (examples in figures 2(B) and (C)). In inspecting those data, we noted that some participants (e.g. S1 shown in figure 2(B)) have a narrower range of acceptable W values. We thus used this narrower range for all further simulations described below. This corresponded to scale [7, 16] (${\sigma _1}$ ≡ [0.0039, 0.0884]) for ${x_{{\text{base}}}}$ and [13, 27] (${\sigma _2}$ ≡ [0.00008,0.0110]) for ${x_{{\text{conflict}}}}$.

Zoom In Zoom Out Reset image size

A good model fit using (4) and (5) acts as an approximation of actual participants receiving stimulation, which means that a good model would generate the same distribution as that of the participants' (S1–S6) data. However, empirically due to the assumptions for (4) and (5), the estimated model requires additional model tuning for exactly matching the distribution of any individual participant. We validated that the generator models, with and without tuning, produced plausible behavior.

Without any tuning (simple fitting to empirical data via COMPASS), the RTs emitted by the generator follow the expected gamma distribution (figures 3(A) and S1). While not a good model for S1 specifically, the generator model output thus follows the expected statistics of actual patient data, i.e. could be used as an example of a new, unobserved patient. To verify this, we fit a gamma distribution to the simulated RT, then calculated the empirical distribution of the RTs expected given that fit. The KS distance between the generator-model RT's pdf and the pdf of a gamma distribution corresponded to p > 0.05 in 100% of the simulation runs across patients. As further verification, we showed that interference trials have higher RTs by the amount that would be expected during MSIT (figure 3(B)). We also verified that, if necessary, a specific patient's data could be matched through the generator model. Hand-tuning parameters A, $\mu$, w, v and $\alpha$ generated simulations that almost perfectly matched with the empirical observations (figures 3(A) and S1; table S1). The calculated KS distance between the empirical and the simulated RT (1000 trials) from the tuned model had a p > 0.05 for 98.9, 96.2, 90.4, 100.0, 91.5, 98.6% of the 1000 simulation runs for S1–S6 respectively.

Zoom In Zoom Out Reset image size

The sensor model tracks the generator well (figures 4(A)–(C) and S2) across the simulations, with a median NRMSE of 0.0935 for the actual/estimated RT, 0.0541 for ${x_{{\text{base}}}}$, and 0.0784 for ${x_{{\text{conflict}}}}$. This suggests that the sensor models can accurately track the unobservable generator state, and that this holds true across many random realizations. However, this tracking is not instantaneous. As seen particularly in figure 4(C), the nature of the sensor model, the random initialization of ${x_0},$ and the chosen noise covariances cause the tracking to be temporarily inaccurate for the first 20 trials.

Zoom In Zoom Out Reset image size

Similarly, the model structure means that the tracked state variable ${x_{{\text{base}}}}$ does not change immediately when stimulation is applied. Since we assume a linear dynamical system with a substantial autoregressive component, it takes 5–10 trials for the effect of a continuous perturbation to settle to a new steady state (figure 4(D)). In general, the number of trials to reach steady state is dependent on the stimulation effect strength, i.e. stronger perturbations will take longer, but also will be more easily detected despite small stochastic variations. Regardless, this result argues that the effects of stimulation can only be evaluated over blocks of trials, not on individual trials.

Having validated the generator and sensor model and identified settling time as a consideration, we evaluated which optimization algorithms could best identify optimal stimulation sites for improving MSIT performance (decreasing task RTs). We measured the accuracy (probability of identifying a known global optimum) of multiple algorithms as a function of (1) blocksize for stimulation, (2) number of stimulation sites ('problem size'), (3) number of trials performed, (4) distance between the stimulation effects of the stimulation sites, (5) SNR of stimulation effect relative to state noise, and (6) number of times an algorithm is repeated during an ensemble strategy.

The accuracy of all algorithms improves with increasing blocksize, but plateaus around a blocksize of 10 (figure 5(A) and tables 6, S2). This is visible across all algorithms. This is consistent with the observation in figure 4 that the stimulation effect is not instantaneous, but takes a few trials to manifest in ${x_{{\text{base}}}}$. Further, at these higher blocksizes, UCB1 is generally the best algorithm. When we compared all algorithms across problem sizes and blocksizes, UCB1 was the most likely to converge to the known optimum (tables 6 and S2). Bayesian algorithms were the best for small blocksizes and single trial simulations, but lost this advantage in part because as the blocksize increases, the estimate of ${x_{{\text{base}}}}$tends to the actual ${x_{{\text{base}}}}$. At the same time, most algorithms consistently outperform BF, emphasizing the overall value of more advanced optimization strategies.

Zoom In Zoom Out Reset image size

Accuracy of all algorithms decreases with increases in the number of stimulation sites. This is a consequence of a decrease in the available samples for a fixed convergence criterion (number of trials) and blocksize. Most algorithms outperformed BF in most cases as the problem size (number of sites) increased, and UCB1 in particular consistently dominated BF (figures 5(C) and S3). However, UCB1 also maintained its advantage over the Bayesian and greedy algorithms (tables 6 and S2). The decrease in accuracy affects the Bayesian algorithms more since they are based on inferring the underlying distribution of the data, whereas greedy algorithms comparatively start performing better.

Accuracy increases with increase in the convergence criterion (number of allowed trials). However, in practice, the number of trials is limited, as a participant can only tolerate a certain number of trials/certain time on task on a given day before fatiguing. In our clinical experience, a 64-trial block takes a bit over 5 min, and it is difficult for a participant to perform more than 9–10 such blocks in close succession. Hence, we emphasized 600 trials as an approximation of this upper limit. In this analysis UCB1 again outperformed other algorithms (tables 7 and S3). It specifically maintained its advantage over BF at all convergence criteria above 500 trials (figure 5(B)), again demonstrating the benefits of a more efficiently exploring optimizer. Importantly, the number of trials collected is the sole measure of an algorithm's efficiency. The tested algorithms do vary in the time needed to compute each trial's outcomes, but all are far faster than the available inter-trial time (table S4).

Accuracy depends on the closeness of local optima to the global optimum. When given two stimulation choices to optimize, BF and UCB1 could identify the optimum with probability above chance, even in very hard problems (figure 5(D)). They only fell to chance in the degenerate case b7, when the two sites were identical. In this specific analysis, there was no advantage for UCB1 or other algorithms over BF; both had almost identical performance as discriminability fell. All algorithms are compared in table S5.

Accuracy of the algorithms is sensitive to noise. In other words, if there is more noise in the data generator, it will mask the stimulation effect and the overall accuracy will decrease. We compared algorithm performance across different levels of SNR, defined based on ratio of the stimulation effect to the state noise process variance. UCB1's advantage was less clear in this case; Bayesian and even greedy algorithms could outperform it in specific cases, although UCB1 was still either first or second across problem sizes (tables 8 and S6).

Ensemble optimization further improved performance (figure 6). Even in the most difficult problem of eight stimulation sites, UCB1 maintained accuracy of nearly 80%, as compared to 47% when simulated once. In simpler problems, e.g. four stimulation sites as was originally tested in [37], UCB1's accuracy exceeded 90%. Performance slightly improved as the ensemble process lengthened from 3 to 5 d, but only asymptotically. UCB1 performed better than almost all algorithms over the ensemble process (table S7). It was also possible to improve performance above 90%, even on difficult problems, by accepting a 'good enough' local minimum (e.g. the 2nd best rather than optimal stimulation site). This improvement in accuracy was observed across algorithms. TS in particular improved compared to its accuracy in finding the single global optimum. Greedy algorithms moved closer to the performance of UCB1, but UCB1 still held a slight advantage over all algorithms. All alternate algorithms performed better than BF. (Tables S8 and S9).

Zoom In Zoom Out Reset image size