NeuroMOSAICS: A collection of neurostimulation datasets - Multi-scale Open-Source Across Interfaces Conditions & Species.

1Introduction¶

Neurostimulation can be applied on a large panel of nervous system disorders and injuries by modulating motor, sensory and cognitive functions through electrical stimulation. This technology offers, amongst others, a promising method to restore motor functions in patients living with paralysis due to neurotrauma, with strong potential for enhanced efficacy through improved spatiotemporal precision Ting et al., 2021James et al., 2018. Neuromodulation systems are becoming increasingly complex Schüpbach et al., 2017Anderson et al., 2018, and this trend is expected to continue and scale further in the future De Ridder et al., 2021. Advances in manufacturing, such as the miniaturization of electrodes enabling higher spatial selectivity, the development of biocompatible designs for targeting areas closer to the nervous system, and the integration of wireless systems providing greater freedom of movement, are key drivers of this evolution Won et al., 2020Jun et al., 2017Tsai et al., 2017. Whether in large channel-count or smaller devices such as current clinical deep brain stimulation (DBS) leads, the effectiveness of neuromodulation depends on system identification and on fine-tuning numerous control parameters to meet individual needs, problems that advanced optimization algorithms have the potential to solve more efficiently than traditional approaches.

1.1A.I. In neurostimulation optimization¶

A wide range of pre-clinical studies, including on epilepsy Pineau et al., 2009 and Parkinson’s disease Gao et al., 2023 are actively exploring how to enhance neuromodulation results using diverse machine learning (ML) methods such as decision trees Gebodh & Bikson, 2021De Andres et al., 2021, random forests Gebodh & Bikson, 2021Shukla et al., 2014, neural networks (NN) Gebodh & Bikson, 2021Shukla et al., 2012 support vector machines (SVM) Gebodh & Bikson, 2021De Andres et al., 2021Medaglia et al., 2020 and Reinforcement Learning (RL) Pineau et al., 2009Gao et al., 2023Yu et al., 2021. However, those ML methods rely on pre-training and require a huge number of iterations to converge to a solution, making them currently impractical to use in an online context. Online neurostimulation optimization means in vivo trial-and-error tuning of a neural interface, where decisions involve selecting from a large number of parameter options, and the constantly evolving environment induces both noise and non-stationarity in the response.

To date, real-time neurostimulation optimization in clinical studies is performed by adaptive stimulation Arlotti et al., 2018Rosa et al., 2015Swann et al., 2018. This closed-loop method traditionally relies on turning on and off DBS and varying the amplitude based on neurophysiological feedback, but this approach is not tailored for a complex multivariate problem such as selecting what specific stimulation pattern to activate, among an immense number of options for varying contact, frequency, current amplitude and pulse width. To offer a solution adapted to optimization of multi parameters, some studies focus on a class of algorithms that does not require training and that relies on parametric models of stimulation properties. Bayesian Optimization based on Gaussian Process (GPBO) can offer an effective solution Bonizzato et al., 2023Choinière et al., 2024Laferriere et al., 2020Sarikhani et al., 2022Aiello et al., 2023. GPBO tests single-input parameter combinations (“queries”), leveraging responses to build an evidence-based approximation (“surrogate” function) that models how stimulus choices affect the desired output. This approach allows to efficiently explore high dimensional spaces in a small number of iterations.

1.2Optimizing stimulation for movement disorders¶

Restoring motor control after neurotrauma or in disease is a very complex neuromodulation challenge since the movement is controlled by distributed brain and spinal networks with variable dynamic states and outputs. Promising studies already show how advanced optimization methods can enhance mobility Bonizzato et al., 2023Grahn et al., 2014, Losanno et al., 2021, but the field of motor neuroprosthetics is still young. Novel exciting demonstrations of motor control restoration feature complex neuroprosthetic systems, with a high degree of control complexity Rowald et al., 2022Wagner et al., 2018. Research to develop optimization agents to fine-tune such systems once they are implanted largely benefits from offline modeling and simulation Capogrosso et al., 2013Rowald et al., 2022.

To effectively train and test complex ML models or simulate online experimental neurostimulation sessions, a diverse range of neurostimulation datasets is essential Wang et al., 2023Bonizzato et al., 2023. Such datasets must account for variability across subjects, species, and neural interface targets. For example, rodents provide the advantage of generating large datasets due to relatively restrained costs, while NHPs offer a closer approximation to human brain complexity, making both species invaluable for research. The characteristics of any single dataset can vary widely due to environmental factors and experimental design, leading to significant differences not only between rodent and NHP data but also across neural targets such as the spinal cord and cortex; even within subjects of the same species.

Although synthetic data can serve as an initial benchmark for algorithm performance, it falls short in addressing the complexities inherent to in vivo data. Many neurophysiological mechanisms, such as the biological variability of motor responses, remain poorly modeled. This variability causes, among other things, inconsistencies in neural output to identical stimulation parameters and is compounded by noise introduced by experimental setups. These challenges make in vivo validation important, not only for assessing algorithm functionality but also for uncovering patterns of responses to neural stimuli and deepening our understanding of neural circuits response dynamics and non-stationarities. However, expanding the collection of neurostimulation datasets faces significant barriers, including high costs, technical complexity, and ethical considerations.

The recent surge in optimization algorithms within this field opens new possibilities for neurostimulation applications. To further push the boundaries of AI in neurostimulation, it is essential to broaden the dataset pool, thereby supporting a wider range of innovative use cases.

Our goal is to develop NeuroMOSAICS (Neurostimulation Multi-scale Open-Source Across Interface, Conditions & Species), a collection of open-access datasets comprising a wide range of different neurostimulation scenarios: in rats and non-human primates (NHP), in sedated and awake animals, in intact and injured subjects, in small and large input spaces optimization, and in both cortical and spinal neurostimulation interventions. NeuroMOSAICS aims to support the development of neurostimulation optimization algorithms. For each dataset, we highlight common challenges that any optimization algorithm must address when applied in this context. Our goal is not only to provide a series of datasets, but also to shed light on these recurring challenges for anyone interested in advancing this field and accelerating the development of related technologies. This collection is an attempt to regroup and valorize all neurostimulation characterization datasets opened by the community. Beyond enabling the training and evaluation of data-demanding algorithms, this collaborative initiative also aims to make valuable knowledge more accessible, thereby encouraging the ML community to design solutions tailored to neurostimulation problems.

This living white paper intends to evolve in time with new incoming open-source datasets. New datasets will be added as we target additional unmet needs in neurostimulation optimization. Moreover, we aim to incorporate complex datasets that establish benchmarks and open avenues for novel algorithmic solutions. Datasets may also be introduced in response to community demand or as invited contributions. We welcome opportunities to review contributions from groups exploring high-dimensional optimization, to collaboratively grow this evolving body of work.

2Materials and methods¶

In this section we describe how each featured dataset has been collected. We wrote this section with a focus on clarity and a level of readership accessible to those without a background in neuroscience.

2.1Dataset 1 – Cortically evoked muscle responses: NHP forelimb and rat hindlimb¶

This dataset contains limb EMG responses to cortical stimulation in 4 male NHPs (2 macaques, 2 capuchins) and 6 female rats. The dataset, available on OSF, was collected to develop a Gaussian Process Bayesian Optimization (GPBO) framework, an algorithm designed to optimize neurostimulation parameters for maximizing MEPs. Cross-species data can challenge the robustness and adaptability of optimization algorithms, testing their ability to generalize across scales and from simpler to more complex cortical organizations, that approximate better the human brain.

The detailed protocol for surgery and experiment is available on M. Bonizzato et al. Bonizzato et al., 2023Choinière et al., 2024. The rats were implanted with 32-microelectrode arrays in the hindlimb sensorimotor cortex and the NHPs were implanted with 96-electrode arrays implanted in the hand area of the primary motor cortex (M1). For rats, the EMGs were implanted in leg muscles. For NHPs, the EMGs were implanted in hand and finger muscles. When the stimulation protocol was running, the two capuchin NHP were under sedation and the two macaques were awake. In the next sections of this preprint, we will refer to Dataset 1 as the NHP results; however, the rat outputs follow the same underlying concept.

2.2Dataset 2 – Forelimb muscle responses to spinal stimulation in rats¶

This dataset contains EMG responses to spinal cord stimulation in the forelimbs of 4 female Long Evans rats. EMG electrodes were implanted bilaterally in the Biceps, Triceps, Flexor Carpi Radialis, Deltoid, and Dorso muscles. Spinal cord stimulation was delivered via a 64-electrode array placed epidurally.

• Experiment Description

A 64-electrode array was placed epidurally on the cervical spinal cord (C3–C6) in 9 distinct placements across 4 sedated rats. In each animal, up to three array positions were tested within this spinal segment. EMG activity was recorded from muscles in both forelimbs in response to single-pulse stimulation delivered through the array.

• Surgical Procedure

The experiments were terminal. Rats were anesthetized with urethane (1 g/kg) prior to surgery. After 30 minutes, they were transitioned to isoflurane anesthesia (∼2%). Isoflurane was withdrawn prior to stimulation, and the rat maintained under the initial dose of urethane. Vital signs, temperature, and hydration were continuously monitored.

EMG Implants EMG microwires were implanted in the Biceps, Triceps, Flexor Carpi Radialis, Deltoid, and Dorso muscles on both sides. Correct placement was confirmed by electrical stimulation through the wires and observation of evoked movements. Wires were inserted transcutaneously.

Spinal Cord Implants The dura mater was exposed from C3 to C6. Initially, a single electrode was used to determine muscle response thresholds. A 64-electrode array was then placed epidurally at up to three different locations.

• Data Collection

For each placement, 10 repetitions of single-pulse stimulation were delivered to each electrode in random order using Tucker-Davis Technologies (TDT) equipment. Stimulation amplitude was set to the threshold value identified during single-electrode testing.

2.3Processing¶

The raw signals in both datasets have been rectified.

3Results¶

In this section, we frequently refer to the heatmaps in each Figure. Each heatmap corresponds to a specific muscle from a subject and displays the average quantification of motor evoked potentials (MEPs) across different electrodes. Those electrode arrays allow us to map the topographic organization of the implanted motor region. Typically, the best-performing electrodes cluster spatially. Each query refers to an individual stimulation.

3.1Dataset 1¶

Dataset 1 was collected during in vivo sessions that involved repeated stimulations of electrode locations in M1. These stimulation sessions were initially conducted on sedated capuchins and on awake rats and macaques. The average MEP responses across spatial locations show a topographic organization of M1 such that nearby cortical locations have similar MEP patterns (Figure Figure 1 A.). While identifying the optimal electrode in this condition may seem straightforward, it presents significant challenges. The heatmaps represent averages of MEP, but individual responses are influenced by excitability, noise factors, and possibly non-stationarity, leading to considerable variability. This variability can obscure important details at the single-trial level. In Figure Figure 1 B., the upper-left panel displays the average response across multiple stimulations delivered through a specific electrode, providing a global overview of the evoked activity. In contrast, the two other panels offer a more detailed perspective : the distribution of peak amplitudes reveals the consistency of responses across trials, while the stack of individual traces allows visualization of how responses evolve over time. Given the goal of minimizing the number of queries, the task becomes finding the right balance between using enough queries for a robust solution and limiting the search to minimize stimulation.

The clear topographic organization of the responses observed in sedated NHPs become considerably noisier and more selective in awake NHPs (see Figure Figure 2 A.). These NHP implants were older, and the quality of the electrical interface at the EMG wire-muscle contact can degrade over time, potentially contributing to variations in noise levels in the recorded muscle responses across subjects. We can also identify two ranges of electrodes (column 8 and 9 in Figure Figure 2 A.) that consistently exhibit very low signal amplitudes across all muscles. This suggests that these portions of the electrode matrix may not have been in full contact with the tissue. Additionally, the NHPs were awake, which may have further contaminated MEPs with spontaneous muscle activations, which are possibly present even when the condition of no muscle activity pre-stimulation is enforced. Consequently, this dataset provides an opportunity to test optimization approaches on data that differ significantly from the previous dataset.

High selectivity becomes particularly problematic when the search space is expanded (such as by adding additional parameters like stimulation amplitude and frequency) and even more so in noisy conditions. If it makes it harder to identify the best electrode, the global maximum of the function to optimize, and increases the risk of getting stuck in local maxima during optimization. For example, while searching for the electrode that grants maximal response in the First Dorsal Interosseous muscle (see Figure Figure 2 A.), one could easily get stuck at coordinates (7,7) on the heatmap. If the algorithm hasn’t yet queried the optimal electrode at coordinates (4,1), and instead queries (7,7) at a given query q, it will produce the best result thus far. The algorithm would then continue to exploit the local maximum around (7,7) without knowing that a better option exists unless it explores very specific electrodes (4,1) or (5,1) (see Figure Figure 2 B.). With around 50 electrodes to query in a limited number of stimulations, the algorithm’s exploration strategy becomes critical to avoid local maxima. This issue of local maxima and high selectivity is also observed in awake rats, as seen in the Left Medial Gastrocnemius available in Dataset 1 on OSF. However, the smaller array (32 electrodes, compared to 96 for NHPs) reduces the complexity of the search.

The temporal aspect of the responses can play an important role in the nature of the response. In any scenario from this dataset, the stronger MEP are peaking at the same time as we can observe in the stacked MEP where deeper blue lines are vertically crossing all the trials.

Considering noise (such as muscle activations unrelated to stimulation), some responses may convey misleading information about the overall response pattern. A clear example is provided by the EDC muscle in Figure Figure 2. The response map indicates a high activity level for electrode (0,6). However, upon closer examination of the stacked EMG traces, it becomes evident that this apparent response is driven almost entirely by a single trial, the 18th stimulation, which shows an isolated increase in MEP amplitude. Notably, this sharp deflection is also visible in other channels, such as electrode(7,2), but occurs at a time (~170 ms post-stimulation) that is unrelated to the actual stimulation event. This suggests the presence of a transient artefact, underscoring the importance of examining individual trials alongside average responses to avoid misinterpretation.

Identifying these outliers in real time is challenging, but it is valuable to do so afterward as it raises awareness of their potential impact on optimization processes. The parameter is_valid classifies outliers in two categories: those that can be excluded due to excessive EMG activity prior to stimulation and those that cannot be excluded, where no prior EMG activity was observed, but the response amplitude is unusually high (such as the 18th trial of electrode (0,6) in ECD muscle, Figure Figure 2). Accurate online processing can easily reject the first of the two classes, while the second is unlikely to be rejected in a real use case. The first category is highlighted in red in the stack MEP figures. The mean and standard deviation of responses in Dataset 1 still include these outliers. A data scientist wishing to remove these spurious responses from the optimization “ground truth”, or objective value, must though recalculate the mean responses.

3.2Dataset 2¶

This dataset was collected from four anesthetized rats, resulting in smooth topographic organization of the spinal cord with high variance for the average responses. Similar to capuchin in Dataset 1, the relaxed muscle state in sedated animals, with no intention of motor control or reflex interference, contributes to the smooth responses. However, unlike Dataset 1, this dataset involves spinal cord (spinal cord) stimulation, leading to some key differences in the results. One notable distinction is that the high-response electrodes are more widely distributed spatially, which may be attributed to the different spatial recruitment properties of spinal cord stimulation compared to motor cortex stimulation. Because dorsal spinal cord stimulation is, from a circuit perspective, closer to the muscle target than M1 stimulation, there is less variability in the responses. Indeed, less synapses (and thus less pulses) lead to a higher repeatability of the responses. This very clear repeatability (a straight deep blue line across each trial in stack MEP figures) makes the identification of outliers easier, and they are also listed in a parameter is_valid and highlighted in red in the stack MEP figures. In this case, the mean and the standard deviation were computed without those outliers.

Whereas the clear repeatability from this dataset has a positive impact on the optimization process, the excitability and noise can be, in some cases, particularly challenging. Indeed, the excitability can increase and decrease drastically in a few queries (like for electrode (2,4) in Left Biceps (see Figure Figure 3 B.)), and some electrodes are losing contact on half of the process (like in electrode (4,3) in Left Biceps), making it harder to identify a best electrode and to take a decision on the definition of what is the best electrode. This phenomenon suggests potential future use of this dataset to take into account the non-stationarity of the responses.

The higher variance in this dataset tends to blur the distinction between the best average solution and its neighboring electrodes. This makes it easy to identify the best location area, but also more challenging to point out the precise optimal electrode. As the variability in the responses increases, it becomes more likely that, at query q, the best electrode yields a lower response than nearby ones. For instance, in the Right Biceps muscle (see Figure Figure 3 A.), the peak amplitude distribution from the best electrode (coordinates (2,5)) is centered at 13.1 mV, but there is an outlier stimulation yielding a response of 11.2 mV. Meanwhile, the electrode at coordinates (2,4) has stimulations with outputs around 13.1 mV. If an algorithm repeatedly queries electrode (2,4) and obtains consecutively good motor responses, while a single query of electrode (2,5) produces an outlier of 11.2 mV, the algorithm might incorrectly infer that electrode (2,4) is superior. Although (2,4) is not the global maximum, it remains a strong candidate because it provides stable, high-quality responses and is situated within the highest region of the heatmap.

Another key insight from this dataset is the spatial organization of the spinal cord for motor control. There is a clear lateralization, with the right side of the spinal cord predominantly evoking stronger responses in muscles on the right side of the body, and the left side evoking responses in left-side muscles. Additionally, there appears to be spatial clustering in spinal cord sites evoking motor potential in specific muscles. The spatial organization is non-trivial given that the dorsal aspect of the spinal cord, which we stimulated, is associated with afferent pathways, and motor output depends on reflex responses Capogrosso et al., 2013. This insight has significant implications for optimization in neuromodulation. The function linking stimulation parameters to motor responses is inherently complex and unknown, posing a “black-box” optimization problem. In order to facilitate the optimization task, incorporating prior knowledge Guay-Hottin et al., 2025 about the spatial organization of the spinal cord can significantly enhance the efficiency of optimization algorithms, such as BO, that is designed to leverage prior information to reduce the complexity of the search space and accelerate convergence to optimal solutions Greenhill et al., 2020. For example, integrating the knowledge that stimulating the left (or right) dorsal spinal cord enhances motor responses on the corresponding side of the body allows the algorithm to narrow its spatial exploration. This targeted approach can substantially decrease the number of iterations required to identify the optimal stimulation parameters, improving both the speed and accuracy of the optimization process.

4Conclusion¶

In this living white paper, we introduced NeuroMOSAICS, an open-access, reproducible and diverse dataset collection covering neurostimulation across different species, neural interfaces and search space dimensions. The two scenarios presented in this first version support optimization research in diverse tasks and objectives. The first dataset provides EMG responses from cortical stimulation in both awake and sedated NHPs and rats, while the second dataset focuses on bilateral EMG responses from spinal cord stimulation in sedated rats.

This collection is an attempt to bring together open neurostimulation characterization datasets shared by the community. These datasets, often difficult to collect due to cost, ethical, and technical constraints, are typically shared only on a limited scale. By making them more widely available, we aim to support researchers in developing and validating complex optimization and data-intensive algorithms tailored to neurostimulation challenges. We aim to foster collaboration between neurostimulation researchers and the ML community, encouraging the creation of algorithms better suited to concrete neurostimulation problems.

6CRediT authorship contribution statement¶

• Lison Kardassevitch: Writing – review & editing, Writing – original draft, Visualization, Project administration, Methodology, Formal analysis, Conceptualization.
• Davide Burchielli: Writing – review & editing.
• Numa Dancause: Writing – review & editing, Supervision.
• Marco Bonizzato: Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition, Conceptualization.

7Funding sources¶

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC; RGPIN-2023-04370) and New Frontiers in Research Fund Exploration (NFRFE-2022-00394).

8Declaration of competing interest¶

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.