Autonomous optimization of neuroprosthetic stimulation parameters that drive the motor cortex and spinal cord outputs in rats and monkeys

Abstract

Neural stimulation can alleviate paralysis and sensory deficits. Novel high-density neural interfaces can enable refined and multipronged neurostimulation interventions. To achieve this, it is essential to develop algorithmic frameworks capable of handling optimization in large parameter spaces. Here, we leveraged an algorithmic class, Gaussian-process (GP)-based Bayesian optimization (BO), to solve this problem. We show that GP-BO efficiently explores the neurostimulation space, outperforming other search strategies after testing only a fraction of the possible combinations. Through a series of real-time multi-dimensional neurostimulation experiments, we demonstrate optimization across diverse biological targets (brain, spinal cord), animal models (rats, non-human primates), in healthy subjects, and in neuroprosthetic intervention after injury, for both immediate and continual learning over multiple sessions. GP-BO can embed and improve "prior" expert/clinical knowledge to dramatically enhance its performance. These results advocate for broader establishment of learning agents as structural elements of neuroprosthetic design, enabling personalization and maximization of therapeutic effectiveness.

Keywords: Bayesian optimization; artificial intelligence; black-box optimization; brain-computer interface; machine learning; neural interfaces; neuromodulation; neurostimulation; neurotechnology; precision medicine.

Graphical abstract

Figure 1

A versatile learning agent for neurostimulation optimization (A) Neurostimulation of motor regions such as M1 causes muscle responses. These controllable electrical signals delivered to the brain can feature complex spatiotemporal patterns. (B) A GP-BO-based algorithm iteratively searches for the most effective input x, capable of eliciting a desirable output y. In doing so, it balances “exploration” (testing unknown inputs) and “exploitation” (pursuing performant solutions). The trade-off between the two is determined by the hyper-parameter k. (C) Intracortical microstimulation (ICMS) evokes selective movements, each featuring a characteristic pattern of EMG activity. Left, stimulus causing a flexion movement and EMG pattern associated. Center, stimulus causing an extension movement and EMG pattern associated. Right, pictures of both evoked movements. (D) An objective function can be cast over the movement, represented as a scalar value. This number encodes one representative feature of EMG activation or kinematic output (f: X → Y, where X is a multi-dimensional space of stimulation parameters and Y is a movement feature, such as peak EMG or a step height), or an arbitrary combination of several biomarkers. The example shows a linear combination using weights w_i. BO has the objective of maximizing (or minimizing) this value. (E) Typical operation of the proposed algorithm consists in a sequence of queries following a pattern that gradually converges over “hotspots.” These are locations in the input space where the objective function is maximal. Two weighted maps of estimated value and uncertainty take part in the acquisition process. (F) GP-BO is a flexible class of learning agents, which features the possibility to inject prior knowledge (e.g., known effective range of values) to further accelerate the search, as well as extracting posterior understanding of the neurostimulation problem (i.e., after testing several subjects, the most effective range of value is explicitly accessible to the user and can be used as prior). See also Figures S1–S3.

Figure 2

Online cortical search for the most performant electrodes to evoke hindlimb EMG responses in rats (A) In sedated rats, we searched for the most performant electrode evoking leg movements using a 32-channel electrode array implanted in the leg sensorimotor cortex using algorithmic hyper-parameters previously optimized offline (see Figure S1). (B and C) The algorithm efficiently determined an estimated performance map across active sites of the implant for two antagonist muscles. Top: average response (reference value map), an empirical reference built by pooling all trials acquired during the entire recording session (i.e., data acquired during all repetitions of the GP-BO, extensive [random] and greedy searches). The value metric was the size of the motor evoked potential (MEP), reported as % of maximal EMG response (average response from the best electrode). Middle: Average estimated performance map (exploration) and comparison with the extensive (random) search. GP-BO search intensified over optimal solutions. The GP-BO algorithm developed an earlier knowledge of the position of performant hotspots. This advantage was maintained throughout the execution of the search. Bottom: Average number of stimuli delivered from each electrode (exploitation) and comparison with the extensive (random) search. The algorithm intensified stimulation of performant hotspots, while the extensive (random) search equally sampled the whole search space. This means that by the end of the 32 queries, in addition to a better knowledge of the hotspots, the algorithm used approximately half of its queries to deliver an effective intervention. (D and E) Average performance for exploration (D) and exploitation (E) over all n = 10 experiments collected in four rats. Data displayed as mean ± SEM. During a series of 32 queries, the algorithm ranked above the extensive (random) search both in exploration and exploitation performance. Because of the random query procedure, the extensive search exploitation performance consists of a flat line. The ripple along this line is due to this performance being estimated using a finite number of experiments (32.5 ± 6.3 runs). See also Figures S4 and S5.

Figure 3

Online cortical search for the most performant electrodes to evoke arm muscles EMG responses in NHP (A) We tested the efficacy of the GP-BO algorithm in a sedated capuchin monkey. We searched for the most performant electrode to evoke EMG responses in muscles driving hand/wrist movements using a 96-channel electrode array chronically implanted in M1. The value metric was the size of the motor evoked potential (MEP). (B and C) During a series of 96 queries, the algorithm ranked above the extensive (random) search both in exploration and exploitation performance. Data displayed as mean ± SEM over all n = 4 experiments collected in 1 NHP. (D and E) The algorithm efficiently determined an estimated performance map across active sites of the implant for two antagonist muscles. Left: Average response (reference value map). Top right: Average estimated performance map (exploration) and comparison with extensive (random) search. The algorithm features an earlier knowledge of the position of performant hotspots, and this better knowledge was maintained throughout execution. Bottom right: Number of stimuli delivered to performant hotspots (exploitation) and comparison with extensive (random) search. Algorithmic search intensified over optimal hotspots, while the extensive search evenly queried the entire space. This meant better knowledge of the true optima (exploration) and earlier access to effective intervention (exploitation). (F) We tested the efficacy of the algorithm in two awake macaque monkeys sitting quietly. A 96-channel electrode array was chronically implanted in the hand representation of M1. (G and H) During series of 96 queries targeting n = 4 muscles, the algorithm performed similarly to what was obtained in the sedated capuchin, ranking above the extensive (random) search both in exploration and exploitation performance. Data displayed as mean ± SEM. See also Figures S4–S6.

Figure 4

Injecting prior knowledge immediately increases performance (A) GP-BO strategies can incorporate existing expertise provided by the user, in this case embodied as a prior on the process mean. A “tabula rasa” approach refers to an empty mean prior. In this experiment, a prior was derived from indirect knowledge obtained from a previous study. Cortical motor maps of ankle flexion movements were collected using visual inspection of responses evoked with intracortical microstimulation (ICMS) in 25 rats. Average results were compiled in a 2D map of estimated probability of evoking ankle flexion (normalized 0%–100%). Across all these maps, the rostral (top) part of the of the array was more likely to generate ankle flexion. (B) We performed online experiments in seven rats in which the algorithm had to find the electrode that evoked the greatest response in the ankle flexor muscle. The value metric was the size of the motor evoked potential (MEP). In one rat, an additional experiment was done for the knee flexor muscle (muscle number 8). These maps had high inter-subject variability. (C and D) In an online demonstration, injecting priors immediately boosted exploration and exploitation performance. GP-BO outperformed extensive and greedy searches, which were provided with the same prior information. Data displayed as mean ± SEM over n = 8 experiments. See also Figure S7.

Figure 5

Continual learning over an evolving neuroprosthetic scenario (A) Schematic representation of the continual learning framework. GP-BO features an explicit representation of prior (before optimizing) and posterior (after optimizing) knowledge, which allows designing continual learning frameworks, where posterior information from a previous session is injected as prior information for a new one. (B) From the first to the third week after SCI, in repeated sessions we sought to identify the electrode providing the most effective activation of leg flexor muscles within a cortical array. The value metric was the size of the motor evoked potential (MEP). The panels show the average response (reference value map) for one rat. As expected, in this post-SCI scenario, neuroplasticity and recovery caused a continuous evolution of cortical representation of movements. (C) We tested the capacity of GP-BO to find the best electrode along recovery, using a very limited number of eight queries in each session, for a quick tuning of electrode choice that would maximize the duration of treatment time. On the first session we provided an initial prior consisting of the average responses across the arrays of the 25 rats of Figure 4A collected 1 week post-SCI. In subsequent sessions, the output of the previous session was used as prior. In this evolving scenario and with such a limited number of queries, GP-BO steadily maintained higher performance than extensive (prior-ranked) and greedy searches. (D) In the same rats we ran an additional online experiment in parallel during the third week post-SCI. To simulate the loss of the best electrode in the array, one of the stimulator outputs was removed and we ran another series of eight queries in these conditions. (E) The loss caused an immediate drop of performance, where GP-BO allowed a rapid recovery of baseline performance. (F) Although, due to the very limited number of queries, this experimental scenario was not designed to study exploitation properties, GP-BO displayed a better exploitation than all alternative methods. All data displayed as mean ± SEM over n = 8 experiments.

Figure 6

Intelligent intracortical neuroprosthesis optimizes spatial features of stimulation to alleviate SCI locomotor deficits (A) Rats with SCI were engaged in treadmill walking. Inner loop: EMG activity was processed online with pattern recognition to determine gait phases and trigger phase-coherent cortical stimulation at the onset of foot lift. Outer loop: A red foot marker was real-time tracked to determine step height at each gait cycle. This information was fed to the fixed-ρ learning algorithm, which searched for the optimal cortical stimulation channel (maximizing foot clearance). The value metric was the step height. (B and C) During a series of 32 queries, the algorithm ranked above extensive (random) search both in exploration and exploitation performance. Data displayed as mean ± SEM over n = 4 rats. (D) Improved exploitation performance means delivering effective intervention while still fine-tuning the search: neurostimulation alleviated foot dragging during active optimization, a clear benefit over extensive search. (E) Example trial. Top: Stepwise estimation of performance. Middle: Foot vertical trajectory. Bottom: EMG traces and stimulation timing. +, one step was not detected by real-time tracking, therefore it was ignored and repeated. Initial wide search phases correspond to low exploitation, meaning low foot clearance and presence of dragging. Steps 10–12 corresponded to the first approach to the most performant area of the cortical implant (top right corner). Subsequently, search intensified in the neighboring region and step height was quite consistently close to the 2-cm range. The algorithm converged (star symbol in figure) on the site considered optimal by average response (reference value map): in the final steps dragging was completely alleviated and the rat displayed foot clearance stably above 2 cm in height. (F) The algorithm efficiently determined an estimated performance map for the cortical array. Left: Average response (reference value map). Top right: Average estimated performance map develops earlier (exploration). Bottom right: Average number of stimuli delivered (exploitation). See also Figure S4.

Figure 7

Intelligent neuroprosthesis optimizes multi-dimensional spinal stimulation parameters and alleviates SCI locomotor deficits (A) Rats with SCI were engaged in treadmill walking. Inner loop: Online gait phases detection triggered phase-specific spinal stimulation. Stimulation trains were delivered at one of seven possible timings in relation to foot lift. Outer loop: The variable-ρ learning algorithm optimized spinal stimulation parameters to maximize step height (the value metric), which was real-time tracked using a camera. (B) A 4D input space was searched, including spinal site and stimulus timing, duration, and frequency. The space featured 672 possible combinations, while the algorithm and the extensive (random) search were given 30 steps to optimize parameters. (C) Caps in amplitude and charge were set as operational constraints to both GP-BO and extensive (random) search. As a result, top frequencies are not allowed full amplitude. Longer durations would also be associated with high charges, further reducing the stimulus amplitude available, due to the charge cap (light blue arrow). (D) Maps from online validation. Each of the two submaps is a 2D average projection, corresponding to the region marked with a square on the complementary submap. The darker regions identified by the algorithm correspond to known spatiotemporal characteristics of spinal stimulation and to the imposed operational constraints. (E) Example trial. Top: Stepwise optimization, progressively intensifying over performance “hotspot.” Bottom: Foot trajectory increased over time with algorithmic learning. (F and G) Exploration and exploitation indicators evolution during a series of 30 queries. Data displayed as mean ± SEM over n = 4 rats. (H) Neurostimulation optimization alleviated foot dragging. Bars, means and individual replicates. (I) Top: Average estimated performance map for spinal stimulation (exploration). An accurate representation of spatiotemporal characteristics of spinal stimulation rapidly emerges during algorithmic searches. Bottom: Average density of stimuli delivered (exploitation). See also Figure S4.

Figure 8

Online optimization in a large, high dimensions, input space (A) In sedated rats (n = 6), we searched for the most performant stimulus pattern to evoke responses in the ankle flexor muscle in a 5D space: stimulus position within a bidimensional (2D) 32-channel electrode array implanted in the leg sensorimotor cortex, stimulus frequency, pulse-width, and duration. In total, the input space featured 10,976 options. We compare the performance of GP-BO to benchmark searches when given the opportunity to search only <1% of this space (100 queries). The value metric was the size of the motor evoked potential (MEP) measured at the end of the stimulation burst. (B and C) Exploration and exploitation indicators evolution for spatial parameters during the series of 100 queries. Data displayed as mean ± SEM over n = 6 rats. (D) Optimal choices for frequency and duration are less subject-dependent than pulse-width and stimulus location. Thus, highly variable parameters cannot be well predicted without optimization. Data displayed as mean ± SEM. (E and G) GP-BO was able to infer reliable estimations of subject-specific properties of 5D neurostimulation after a single run of 100 queries. The inferred posterior distributions after a representative single run reliably approximate the body of information collected throughout the entire session. (E) GP-BO personalized cortical location selection (2D) in rats having diverse hotspots. (F) GP-BO correctly identified optimal pulse-width, which highly varied between subjects. (G) GP-BO identified optimal stimulus duration and frequency. Although the two parameters were less variable across subjects, finer details (broader frequency range in rat #2, shift in optimal duration between rats) were correctly modeled. See also Figure S8.