Bayesian optimization of peripheral intraneural stimulation protocols to evoke distal limb movements

Experiments were conducted on three adult female Lewis rats (LEW-ORlj, Charles River Laboratories, France) (weight ∼ 220 g). The animals were group-housed on a 12 h light-dark cycle with food and water ad libitum. The three rats were chronically implanted with a SELINE intraneural electrode [25] (figure 1(a)) in the left sciatic nerve, placed ∼0.5 cm proximally to the bifurcation into the tibial and common peroneal nerves. Wires were connected to a 3D-printed pedestal that was sutured to the fascia of the back-muscles (multifidus muscle and gluteus superficialis), through a 2 × 2 cm piece of surgical mesh (Mersilene mesh, Ethicon Inc., NJ). To acquire intramuscular electromyographic (EMG) activity, bipolar electrodes (Teflon-coated stainless steel wires) were chronically implanted in four muscles of the left hindlimb: tibialis anterior (TA) (responsible for the flexion of the ankle), soleus (SOL), gastrocnemius medialis (GM), and plantaris (PL) (jointly responsible for the extension of the ankle). Wires were routed subcutaneously to a percutaneous connector (Omnetics Connector Corp., MN) that was fixated to the animal's skull with dental cement. All surgeries were performed under aseptic conditions and general anesthesia with isoflurane in oxygen enriched air (1%–2%). Details on the rats involved in the study are specified in supplementary table 1 (available online at stacks.iop.org/JNE/18/066046/mmedia).

Zoom In Zoom Out Reset image size

Data acquired during an experiment of single pulse stimulation of the sciatic nerve (figure 1(b)) in the three rats were used to test the algorithm offline. The experiment was conducted 4–5 weeks after implantation. Rats were slightly sedated (Domitor, 0.0025–0.005 ml kg⁻¹) and positioned at the edge of a table with their legs hanging down. Electrical stimulation was delivered by the IZ2H stimulator (Tucker David Technologies, USA). Biphasic cathodic-first current pulses (1 Hz) with a pulse-width of 40–80 us were applied from each working active site of the SELINE electrode at increasing amplitude values, from subthreshold to saturation of the compound muscle action potential (CMAP) (increments of 30–50 uA, three repetitions per value). Intramuscular EMG signals were amplified (1000×, Differential AC amplifier, AM systems, WA), bandpass filtered (10–2000 Hz) and acquired at ∼25 kHz using the RZ2 Processor (Tucker David Technologies, USA). Hindlimb kinematics was recorded using a high resolution camera.

Experiments were carried out on five adult female macaques (Macaca fascicularis) (4–9 years old, weight between 3.1 and 3.9 kg) that were initially involved in other studies: acute electrophysiology was performed in three animals during a terminal procedure before euthanasia, and chronic electrophysiology was performed in two others. The animals were group-housed in an enriched indoor room and had access to water and food ad libitum. All the monkeys were implanted with a custom-made transverse intrafascicular multichannel electrode (TIME) [24] in the median nerve and two of them were additionally implanted in the radial nerve. In particular, the TIME used in this study (the Mk-TIME) (figure 1(c)) was adapted from a previous design [36] to fit the dimensions of the monkey arm nerves. Implantation of the Mk-TIME followed previously described methods [37]. We inserted the median electrode ∼2 cm proximally to the elbow and the radial electrode ∼2 cm proximally to the epicondyle along the humeral bone. To record intramuscular EMG activity, each animal received the implantation of bipolar electrodes in up to six flexor muscles of the hand: flexor carpi radialis (FCR) and palmaris longus (PL) (wrist flexors), flexor digitorum profundus (FDP) and flexor digitorum superficialis (FDS) (finger flexors), thenar (THE) (thumb flexor), first dorsal interosseous (FDI) (index abductor). The two monkeys with the Mk-TIME in the radial nerve were also implanted in four extensor muscles of the hand: extensor carpi radialis (ECR) (wrist extensor), extensor digitorum communis (EDC) (finger extensor), extensor carpi ulnaris (ECU) (responsible for wrist ulnar deviation), abductor pollicis longus (APL) (thumb extensor). Chronic EMG recordings were performed using pairs of Teflon-coated stainless steel wires, implanted according to the procedure described in [3]. Acute EMG recordings were performed with subcutaneous needles (Ambu® Neuroline). Intrinsic hand muscles were always implanted acutely. All surgeries were performed under aseptic conditions and general anesthesia induced with midazolam (0.1 mg kg⁻¹), methadone (0.2 mg kg⁻¹), and ketamine (10 mg kg⁻¹) and maintained under continuous intravenous infusion of propofol (5 ml kg⁻¹ h⁻¹) and fentanyl (0.2–1.7 ml kg⁻¹ h⁻¹). All the details about the monkeys involved in the study and the implants received are specified in supplementary table 2.

In the case of acute implants, electrophysiological experiments were conducted just after the surgical intervention, keeping the animals under anesthesia with intravenous perfusion of propofol (5 ml kg⁻¹h⁻¹) and fentanyl (0.2–1.7 ml kg⁻¹h⁻¹) and continuously monitoring the vital signs. As for the two chronically implanted animals, Mk-Me was anesthetized during the experimental session analyzed, while Mk-Ol_c was awake, sitting on a chair and slightly sedated (ketamine, 3.5 mg kg⁻¹, NaCl 1:1, 0.1 ml/15 min). During all the sessions the animal's forearm was restrained on an armrest while the hand was either relaxed or held open by elastic bands tied to the fingers. Electrical stimulation was delivered by the IZ2H stimulator (Tucker David Technologies, USA) as bursts of charge-balanced cathodic-first biphasic pulses. Stimulation parameters (channel, burst-length, pulse-width, amplitude and frequency) were set through custom-made MATLAB (MathWorks, Natick MA) or C++ (Visual Studio®, USA) routines. Intramuscular EMG activity was amplified (1000×, PZ5, Tucker David Technologies, USA) and acquired at ∼12 kHz using the RZ2 Processor (Tucker David Technologies, USA). Hand kinematics was recorded using up to three high resolution cameras. Grip force was measured using a custom sensor [38].

Data acquired during an experiment of burst stimulation of the median and radial nerves (figure 1(d)) in n = 4 monkeys (supplementary table 2) were employed to test the algorithm offline. During the experiment, all the working active sites of the median and radial Mk-TIMEs were scanned with 500 ms monopolar stimulation bursts at increasing current intensities. The stimulation frequency was fixed at 50 Hz and the pulse-width was set at either 40 or 80 us depending on the channel impedance. As for the selection of intensities, different procedures were undertaken in the four monkeys. For Mk-Pa, we tested all the amplitude values from movement threshold to movement saturation (increments of 10 uA, three repetitions per value). For the other three animals, the amplitude values tested were manually selected within the functional range in a step that was neither uniform nor consistent, resulting in not completely exhaustive datasets. Each intensity value was applied for 1–6 repetitions in this case. In Mk-Me, two of the tested channels were removed from the dataset. Indeed, they generated the pronation of the forearm, and this type of movement could not be characterized because the EMG activity of the muscle responsible for it (pronator teres) was not recorded (supplementary table 2).

The algorithm was tested online on Mk-Le, acutely implanted with the Mk-TIME in the median nerve. Before the actual test, a training phase was performed to gather preliminary information about the implant just performed. During this phase, we applied a few bursts of stimulation from all the Mk-TIME active sites to determine which ones elicited a visible motor response and which amplitude values corresponded to movement threshold and saturation. Based on this quick inspection, we restricted the algorithm search space to functional channel-amplitude pairs: non-functioning channels were discarded, and for the remaining active sites, the amplitude was restricted to the values between the minimum and maximum determined experimentally (steps of 10 uA). The burst-length, frequency and pulse-width were fixed at 500 ms, 50 Hz and 40 us, respectively. Then the actual algorithm test began. The first burst of each run was delivered from the first channel of the map (figure 1(d)) at minimum amplitude. For the following iterations, the algorithm selected the active site and amplitude value to be tested, and a stimulation burst with these parameters was applied. The recorded EMG signals were processed and the extracted measures were added to the set of previous observations that the algorithm used to choose the subsequent action. Due to the limited time available for the experiment, we terminated each run after 60 iterations.

BO [32] was identified as a perfectly suitable strategy for our purpose, as it is an ideal tool for optimizing black-box objective functions under query limits. Briefly, BO models the unknown objective function as a random function and assigns it a prior distribution. At each iteration, a new response of the function is observed, and a posterior distribution is built by conditioning the prior on all the data observed so far. An acquisition function is finally applied to the posterior distribution to determine the next query point. The acquisition function trades off the exploration of high-uncertainty regions and the exploitation of high-reward regions to reach the global optimum of the objective with the minimum number of iterations. In practice, it is a matter of exploiting promising regions so as not to waste time on useless observations, but at the same time exploring the sample space sufficiently so as not to get stuck in local optima. To apply this technique, we had to define a scalar objective function, a prior for our measurable quantities, and an acquisition function. We used objective functions built from the aggregation of multiple EMG measures, multi-output Gaussian processes (GPs) [39] as priors, and customized upper confidence bound (UCB) like acquisition functions [32]. Details are specified in the following sections. The functioning of our GP based BO (GP-BO) algorithm is summarized in figure 2.

Zoom In Zoom Out Reset image size

Unlike previous studies [34, 35], here we tested BO to solve a multi-attribute neurostimulation optimization problem. In particular, we characterized and optimized the motor responses evoked by intraneural PNS using features extracted from the EMG signals of the $M$ implanted muscles $\left( {\left[ {{\mu _1},\, \ldots ,\,{\mu _M}} \right]} \right)$. For single pulse stimulation (rat experiments), we computed the peak-to-peak of the CMAP, identified as the segment of the EMG signal from 2 to 14 ms after a stimulation pulse. Values obtained for each muscle were then normalized to the maximum found during the experiment. For burst stimulation (monkey experiments), we calculated the average amplitude of the normalized EMG envelope. Specifically, raw EMG data were band-pass filtered between 50 and 500 Hz (Butterworth, 3rd order). Stimulation artifacts were removed using a 3rd order Savitzky–Golay filter with a smoothing window of 2.5 ms. The linear envelope was then computed by rectification and low-pass filtering at 6 Hz (Butterworth, 3rd order). The signals of each muscle were normalized with respect to (a) the maximal activity obtained during the whole experimental session in the case of the offline tests, (b) the maximal activity obtained during the training phase for the online test. Finally, we calculated the average amplitude of the resulting normalized envelope.

To effectively evoke a specific motor function with neurostimulation, we have to search for parameters that strike an appropriate balance between the response of muscles that promote the movement and muscles that impede it. Indeed, because functionally distinct motoneurons are positioned in close proximity within neural structures, a change of stimulating channel or an increase in pulse amplitude may result in many muscles being recruited together in an unselective manner. We solved this multi-objective optimization problem using the scalarization approach, i.e. we aggregated the EMG-based output variables into a single scalar to optimize in our search. We defined two different objective (or reward) functions ${\text{obj}}\left( {\boldsymbol{{x}}} \right)$ for the offline and online algorithm tests. During the offline tests, the search objective was to maximize the activity of the muscles considered functional for the target movement while minimizing the activity of the undesired muscles. In this way, we aimed to selectively evoke the target motor function and at the same time maximize its strength. We linearized this objective by expressing it as:

$$\begin{equation*}\mathop {{\text{argmax}}}\limits_{{\boldsymbol{{x}}} \in D} {\text{ obj}}\left( {\boldsymbol{{x}}} \right) = \mathop {{\text{argmax}}}\limits_{{\boldsymbol{{x}}} \in D} \sum\limits_{i = 1}^M {w_i}{\mu _i}\left( {\boldsymbol{{x}}} \right)\end{equation*}$$

where ${\boldsymbol{{x}}} = \left[ {{x_1},{x_2}} \right]$ (1 = channel, 2 = amplitude), $D$ is the decision set, ${\mu _i}$ is the EMG response of the $i{\text{th}}$ muscle $(0 \le {\mu _i} \le 1)$, ${w_i} = 1$ if the $i{\text{th}}$ muscle is functional for the target movement and ${w_i} = - 1$ if the $i{\text{th}}$ muscle is nonfunctional for the movement. We chose to use $\left| {{w_i}} \right| = 1\,\forall \,i$ (all the muscles equally contribute to the multi-attribute optimization) to make the procedure as generalizable as possible. Based on their physiological function, we identified which of the implanted muscles were functional for each target movement (supplementary table 3). With this definition of objective, the ideal EMG pattern for a given target movement consists in ${\mu _i} = 1$ for all functional muscles, and ${\mu _i} = 0$ for all nonfunctional muscles. For the online test, we defined the reward function as the Euclidean distance between the evoked and a predefined target EMG pattern:

$$\begin{equation*}\mathop {{\text{argmin}}}\limits_{{\boldsymbol{{x}}} \in D} {\text{ obj}}\left( {\boldsymbol{{x}}} \right) = \mathop {{\text{argmin}}}\limits_{{\boldsymbol{{x}}} \in D} \parallel {\boldsymbol{{\mu }}}\left( {\boldsymbol{{x}}} \right) - {{\boldsymbol{{\mu }}}_{{\text{target}}}}\parallel .\end{equation*}$$

We used as target the stimulation-evoked EMG activity obtained in another animal (Mk-Me, figure 6(b)), manually selecting trials in which a visually adequate movement was generated. This target implant was particularly successful in terms of selectivity, as it produced several different movements. We chose this objective relying on the reproducibility of muscle recruitment between implants.

Before defining the model prior, we rearranged the 2D input space denoted by channel and amplitude to facilitate the algorithm in input–output mapping. Channels were represented through a distance-based map (figures 1(b) and (d)) to encode the similarity of muscle recruitment based on spatial proximity. For offline testing, we also adjusted the amplitude dimension. Because the amplitude functional range was quite different from one channel to another and the pulse-width was either 40 or 80 us, we normalized the amplitude values tested for each active site between 0 and 1. Channel impedance variability could otherwise have been inadvertently encoded in our model.

The probabilistic model linking the input space thus defined to the multiple EMG-based output variables was assumed to be a multi-output GP. GPs allow to intuitively exploit the a priori knowledge of input–output mapping through the selection of a mean, a covariance function and their hyperparameters [39]. While GPs were originally employed for single-output scenarios, they then found application in various fields for multi-output (or multi-task) learning problems [40, 41]. Assuming that the outputs are correlated (in our case some muscles are expected to have related responses to stimulation), multi-task learning is based on the idea that by exploiting the information shared between the tasks, we can improve the prediction compared to modelling the tasks individually. The key is then to develop admissible multi-output covariances that capture not only the similarity between the inputs but also the interaction among the outputs. In this sense, several solutions have been proposed [40, 41].

Here, we compared three different multi-output models. First, we threated the multiple outputs as independent GPs with shared hyperparameters [42]. We will hereafter refer to this model as IGP. By imposing the same hyperparameters some correlation between the outputs is still implicitly assumed to exist. Then, we implemented two widely used approaches that explicitly model the correlation between the outputs. They both belong to the group of separable models, i.e. the multi-output covariance is decomposed into the sum of products between a kernel over the input space and a kernel for the outputs (inputs and outputs are threated separately). Specifically, the second model we used is a particular case of the linear model of coregionalization (LMC), called semiparametric latent factor model (SLFM) [43]. Briefly, the LMC expresses the $M$ outputs as a linear combination of $Q$ independent latent functions: ${f_m}\left( {\boldsymbol{{x}}} \right) = \,\mathop \sum \limits_{q = 1}^Q {a_{m,q}}{u_q}\left( {\boldsymbol{{x}}} \right)$, Here, we chose $Q = M$. The multi-output covariance in this model can be expressed as ${K_M}\left( {{\boldsymbol{{x}}},{\boldsymbol{{x}}}^{\prime}} \right) = \mathop \sum \limits_{q = 1}^Q {A_q}{k_q}\left( {{\boldsymbol{{x}}},{{\boldsymbol{{x}}}^{\prime}}} \right)$, where ${A_q}$ is the correlation (or coregionalization) matrix and ${k_q}\left( {{\boldsymbol{{x}}},{{\boldsymbol{{x}}}^{\prime}}} \right)$ is the covariance of the $q{\text{th}}$ latent function. The LMC-type SLFM defines the coregionalization matrix as: ${A_q} = {{\boldsymbol{{a}}}_q}{\boldsymbol{{a}}}_q^T$. Finally, we tested a particular case of the intrinsic coregionalization model (ICM), called multi-task GP (MTGP) model [44]. The ICM is a simplified version of the LMC as it employs only a single latent function. The ICM-type MTGP parameterizes ${A_1}$ based on the Cholesky decomposition ${A_1} = L{L^T}$. We compared the IGP, MTGP e SLFM models in our offline tests, while for the online test we solely used the IGP approach.

As for the input covariances, in all the three scenarios we used a composite kernel formed as the product of two standard Matérn 52 [39], one in the channel space and the other in the amplitude space. A standard Matérn 52 is a stationary distance-based kernel defined as:

$$\begin{align} k\left( {{\boldsymbol{{x}}},{\boldsymbol{{x^{\prime}}}}} \right) &= {\sigma ^{\,2}}\left( {1 + \sqrt 5 r + 5/3{r^{\,2}}} \right){e^{ - \sqrt 5 r}},\,\,\\ r &= \,\sqrt {{{\left( {{\boldsymbol{{x}}} - {\boldsymbol{{x^{\prime}}}}} \right)}^T}{{{\Sigma }}^{ - 1}}\left( {{\boldsymbol{{x}}} - {\boldsymbol{{x^{\prime}}}}} \right)},\\ {{\Sigma }} &= {\text{diag}}\left( {l_1^{\,2}, \ldots l_n^{\,2}} \right) \end{align}$$

with $n$ being the dimension of ${\boldsymbol{{x}}}$ (here $n = 2$). The Matérn class was chosen because it can depict rough behaviors [39]. The product of the two Matérn along the two dimensions implies that two stimuli ${\boldsymbol{{x}}}$ and ${\boldsymbol{{x^{\prime}}}}$ elicit a similar muscle response only if both their channels and amplitude values are close. Regarding the mean function of the multi-output GP, we used a constant mean for the offline tests and a zero mean for the online test.

The hyperparameters of the models tested (such as the constant $c$ for the mean function, the length-scales ${\boldsymbol{{l}}} = \left[ {{l_{C1}},{l_{C2}},{l_{A1}},{l_{A2}}} \right]$ and the signal variance ${\sigma ^{\,2}}$ for the input covariances, and the parameters of the coregionalization matrix for the MTGP and SLFM models) were tuned on each animal training set (see section 2.3.6): we found the optimal hyperparameters to maximize the marginal likelihood of the training data using the conjugate gradient method [45]. The values obtained for the IGP model are shown in supplementary table 4 for reference.

As mentioned in section 2.3.1, BO updates the probabilistic model at each iteration by conditioning the prior distribution on previous observations and thus makes predictions about future responses. In the case of multi-output GP models, if we condition the prior ${\boldsymbol{{f}}}\left( {\boldsymbol{{x}}} \right)$ on a finite number of responses ${{\boldsymbol{{y}}}_{t - 1}} = \,{\left[ {{{\boldsymbol{{y}}}_1}, \ldots ,{{\boldsymbol{{y}}}_{t - 1}}} \right]^T}$ corresponding to actions ${{\boldsymbol{{X}}}_{t - 1}} = \,{\left[ {{{\boldsymbol{{x}}}_1}, \ldots ,{{\boldsymbol{{x}}}_{t - 1}}} \right]^T}$ and we assume that the noise is i.i.d. Gaussian, we get a Gaussian posterior distribution [41, 46]:

$$\begin{equation*}{\boldsymbol{{f}}}\left( {\boldsymbol{{x}}} \right)\,|\,{{\boldsymbol{{y}}}_{t - 1}}\,\sim \,\mathcal{N}\left( {{{\boldsymbol{{m}}}_{t - 1}}\left( {\boldsymbol{{x}}} \right),{{\mathop \sum \nolimits }_{t - 1}}\left( {{\boldsymbol{{x}}},{\boldsymbol{{x}}}} \right)} \right)\end{equation*}$$

where

$$\begin{equation*}{{\boldsymbol{{m}}}_{t - 1}} = { }{K_M}{\left( {{\boldsymbol{{x}}},{{\boldsymbol{{X}}}_{t - 1}}} \right)^T}{\left[ {{K_M}\left( {{{\boldsymbol{{X}}}_{t - 1}},{{\boldsymbol{{X}}}_{t - 1}}} \right) + {{\mathop \sum \nolimits _M}}} \right]^{ - 1}}{{\boldsymbol{{y}}}_{t - 1}},\end{equation*}$$

$$\begin{align*}{\mathop \sum \nolimits _{t-1}}\left( {{\boldsymbol{{x}}},{\boldsymbol{{x}}}} \right) &= { }{K_M}\left( {{\boldsymbol{{x}}},{\boldsymbol{{x}}}} \right) - { }{K_M}{\left( {{\boldsymbol{{x}}},{\boldsymbol{{\,}}}{{\boldsymbol{{X}}}_{t - 1}}} \right)^T}\\[9pt] &\quad\times{\left[ {{K_M}\left( {{{\boldsymbol{{X}}}_{t - 1}},{\boldsymbol{{\,}}}{{\boldsymbol{{X}}}_{t - 1}}} \right) + {{\mathop \sum \nolimits _M}}} \right]^{ - 1}}\\[9pt] &\quad\times{K_M}{\left( {{\boldsymbol{{x}}},{\boldsymbol{{\,}}}{{\boldsymbol{{X}}}_{t - 1}}} \right)^T},\end{align*}$$

where ${\mathop \sum \nolimits _M} = {\mathop \sum \nolimits _s} \otimes {I_{t - 1}}$ is a diagonal noise matrix, whose elements are additional hyperparameters that we optimized based on each animal training set (supplementary table 4). Given this posterior, all the points ${\boldsymbol{{x}}} \in D$ ($D\,$ decision set) are assigned a mean vector and a covariance matrix.

For our offline tests, we used UCB-like acquisition functions [32] and determined the next observation ${{\boldsymbol{{x}}}_t}$ as:

$$\begin{equation*}{{\boldsymbol{{x}}}_t} = \mathop {{\text{argmax }}}\limits_{{\boldsymbol{{x}}} \in D} \left( {m_{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right) + { }k\sigma _{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right)} \right)\end{equation*}$$

being $m_{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right)$ and $\sigma _{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right)$ the mean and standard deviation (std) of the posterior distribution of the objective function ${\text{obj}}\left( {\boldsymbol{{x}}} \right)$ at action $t - 1$. The first term of the equation is the exploitation rule: it tends to regions where the expected objective is maximal based on the current data. The second term of the equation is the exploration rule: the search is guided towards high-uncertainty regions. The two processes must be balanced to avoid wasting time and getting stuck in local optima. This balance is determined by the parameter $k\,(k > 0)$. We tuned $k$ to get a final minimum regret (see section 2.4.1) reasonably low and approximately equal across the different models, while having a good convergence rate. The chosen values are shown in supplementary table 4. For the IGP approach, since ${\text{obj}}\left( {\boldsymbol{{x}}} \right)$ was built as a weighted linear combination of independent GPs modeling the EMG responses (see section 2.3.2), its posterior mean and std were computed as $m_{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right) = \sum\limits_{i = 1}^M {w_i}m_{t - 1}^i\left( {\boldsymbol{{x}}} \right)$, $\sigma _{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right) = \sqrt {\sum\limits_{i = 1}^M {{({w_i}\sigma _{t - 1}^i\left( {\boldsymbol{{x}}} \right))}^{\,2}}}$, being $m_{t - 1}^i\left( {\boldsymbol{{x}}} \right)$ and $\sigma _{t - 1}^i\left( {\boldsymbol{{x}}} \right)$ the posterior mean and std of the $i{\text{th}}$ EMG response. In the case of correlated outputs, i.e. for the MTGP and SLFM models, the objective posterior mean and std were calculated as $m_{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right) = \sum\limits_{i = 1}^M {w_i}m_{t - 1}^i\left( {\boldsymbol{{x}}} \right)$, $\sigma _{t - 1}^{{\text{obj}}}\left( {\boldsymbol{{x}}} \right) = \sqrt {\sum\limits_{i = 1}^M {{({w_i}\sigma _{t - 1}^i\left( {\boldsymbol{{x}}} \right))}^{\,2}} + 2\sum\limits_{i = 1}^M \sum\limits_{j = i + 1}^M {w_i}{w_j}co{v_{i,j}}\left( {\boldsymbol{{x}}} \right)}$.

For the online test, we used a readapted version of the UCB decision rule, because the objective was defined as a nonlinear function of independent GPs (see section 2.3.2). Specifically, the next query ${{\boldsymbol{{x}}}_t}$ was determined as:

$$\begin{align*}{{\boldsymbol{{x}}}_t} = \mathop {{\text{argmax }}}\limits_{{\boldsymbol{{x}}} \in D} &\left( - \sqrt {\sum\limits_{i = 1}^M {{\left( {m_{t - 1}^i\left( {\boldsymbol{{x}}} \right) - \mu _{{\text{target}}}^i} \right)}^{\,2}}}\right.\\[9pt] &\quad\left. + k\sqrt {\sum\limits_{i = 1}^M {{(\sigma _{t - 1}^i\left( {\boldsymbol{{x}}} \right))}^{\,2}}} \right).\end{align*}$$

In this case, we used $k = 3$.

Each dataset employed for offline testing was split into a training set, which was used to optimize the model hyperparameters, and a test set, which could be sampled during the algorithm search. The test set was constructed to contain a maximum of two repetitions for each tested combination of channel and amplitude. If more than two repetitions were present, redundant stimuli were added to the training set. For the online test, the dataset acquired during the preliminary phase was used as the training set. The dimension of the final sets and the number of unique stimuli are specified in supplementary tables 1 and 2.

The different versions of the multi-output GP-BO algorithm were implemented in MATLAB (MathWorks, Natick MA) using custom code, with the aid of the GPML toolbox [45] to create the input covariance and learn the model hyperparameters (using the minimize function). A custom MATLAB graphical user interface (supplementary figure 1) was designed for the online test to drive the communication between the stimulator that sent the electrical pulses, the processor that recorded the EMG signals, and the algorithm code.

In the presence of ground truth data, as in our offline tests, BO performance can be quantified using the concept of regret [47]. The instantaneous regret at action ${{\boldsymbol{{x}}}_t}\,$ is defined as $r\left( {{{\boldsymbol{{x}}}_t}} \right) = {\text{ obj}}\left( {{{\boldsymbol{{x}}}^{\boldsymbol{{*}}}}} \right) - {\text{obj}}\left( {{{\boldsymbol{{x}}}_t}} \right)$, where ${{\boldsymbol{{x}}}^{\boldsymbol{{*}}}} = \mathop {{\text{argmax}}}\limits_{{\boldsymbol{{x}}} \in D} {\text{ obj}}\left( {\boldsymbol{{x}}} \right)$ is the optimal solution. The cumulative regret ${R_T} = \,\mathop \sum \limits_{t = 1}^T r\left( {{{\boldsymbol{{x}}}_t}} \right)$ is the sum of instantaneous regrets accumulated over $T$ queries. Based on these definitions, the performance of the algorithm was expressed through the action-by-action evolution of two derived metrics: (a) minimum regret $M{R_T} = \mathop {\min }\limits_{t \in T} r\left( {{{\boldsymbol{{x}}}_t}} \right)$ and (b) time-average regret $A{R_T} = {R_T}/T$. The minimum regret expresses the accuracy of the algorithm in the sense of its ability to find high-reward stimuli, but without evaluating all the observations made. Conversely, the average time regret considers how many of the visited stimuli are successful, and thus describes the actual balance between exploration and exploitation. We also assessed the number of actions required for convergence (detected when the same stimulus was observed at least five times in a row) and how many unique parameter combinations were tested among them (hereafter referred to as unique actions). We compared these values with the size of the decision set, i.e. the number of available unique channel-amplitude pairs.

For the online test, since the optimal solution was unknown, we evaluated the functioning of the algorithm by assessing the dynamics of the minimum and average objective observed so far.

To quantify the movements evoked by the stimulation parameters selected at algorithm convergence, we analyzed the relative kinematics. Hindlimb joints in rats and hand joints in monkeys were manually selected on 2D video frames corresponding to rest and maximal range of movement, using a custom MATLAB routine. 2D kinematic features were then computed as described in supplementary table 5. The values at rest were subtracted from the values at the maximal range of movement.

Grip force was also measured for few stimuli identified at convergence in monkeys. Voltage values were converted in N and normalized with respect to the maximum voluntary contraction (MVC) of Mk-Ol_c squeezing a spherical object [28].

Data are reported as mean ± std or mean ± standard error of the mean (s.e.m). The choice is specified in figure captions. Significance was evaluated using the non-parametric Kruskal–Wallis test with Bonferroni correction for multiple comparisons.