Recalibration of neuromodulation parameters in neural implants with adaptive Bayesian optimization

Objective. Neuromodulation technology holds promise for treating conditions where physiological mechanisms of neural activity have been affected. To make treatments efficient and devices highly effective, neurostimulation protocols must be personalized. The interface between the targeted nervous tissue and the neurotechnology (i.e. human-machine link or neural interface) usually requires constant re-calibration of neuromodulation parameters, due to many different biological and microscale phenomena happening over-time. This adaptation of the optimal stimulation parameters generally involves an expert-mediated re-calibration, with corresponding economic burden, compromised every-day usability and efficacy of the device, and consequent loss of time and increased discomfort of patients going back to clinics to get the device tuned. We aim to construct an adaptable AI-based system, able to compensate for these changes autonomously. Approach. We exploited Gaussian process-based Bayesian optimization (GPBO) methods to re-adjust the neurostimulation parameters in realistic neuroprosthetic data by integrating temporal information into the process to tackle the issue of time variability. To this aim, we built a predictive model able to tune the neuromodulation parameters in two separate crucial scenarios where re-calibration is needed. In the first one, we built a model able to find the optimal active sites in a multichannel electrode, i.e. able to cover a certain function for a neuroprosthesis, which in this specific case was the evoked-sensation location variability. In the second one, we propose an algorithm able to adapt the injected charge required to obtain a functional neural activation (e.g. perceptual threshold variability). By retrospectively collecting the outcomes from the calibration experiments in a human clinical trial utilizing implantable neuromodulation devices, we were able to quantitatively assess our GPBO-based approach in an offline setting. Main results. Our automatic algorithm can successfully adapt neurostimulation parameters to evoked-sensation location changes and to perceptual threshold changes over-time. These findings propose a quick, automatic way to tackle the inevitable variability of neurostimulation parameters over time. Upon validation in other frameworks it increases the usability of this technology through decreasing the time and the cost of the treatment supporting the potential for future widespread use. This work suggests the exploitation of AI-based methods for developing the next generation of ‘smart’ neuromodulation devices.


Introduction
Neuromodulation is a promising approach for treating neurological diseases through its ability to restore motor, sensory and cognitive functions among others [1][2][3][4][5][6][7][8][9][10][11]. It also represents a huge and highly-impactful market, expected to double in size in the near future [12,13]. Neuromodulation is the alteration of neural activity through delivering electrical stimuli directly to a target nervous system area. These methods encompass the need for a bio-compatible neural interface (e.g. implantable active electrodes), connected to a stimulator, which can reliably and effectively deliver custom stimulation patterns to activate or inhibit target regions of the nervous system, such as the brain [14], the spinal cord [1,3,15,16] or the periphery [8,17].
Growing interest is placed into the tailoring of such protocols of stimulation, both for targeting the correct neural structures but also to make the stimulation more effective. The number of parameters to be tuned can exponentially grow depending on the number of active sites (AS) used, and the possibility of multipolar and modulated electrical stimulation. Moreover, inter-subjects and time-variability are important factors to be considered to obtain a personalized, stable and effective neuromodulation [8,[18][19][20]. Thus, innovative methods, based on artificial intelligence (AI) models, are being considered for helping the calibration of such devices, which need an expert technician able to deal with the multidimensional space of possible parameters and choose the optimal combination in a short time [21][22][23][24]. In particular, once the stimulation has been customized, the need for constant control and re-calibration arises from the limited stability of the neural interfaces [25][26][27]. Indeed, the output evoked by a certain stimulation pattern can change with time, leading to an unexpected, less efficient or even dangerous outcome of that neuromodulation protocol [8,26]. Consequently, an expert technician is needed to readapt the stimulation parameters. This would lead to a higher cost, poor usability and efficacy of the device, a very time-consuming process, and consequent patients' distress.
Based on which applications are being considered, re-calibration issues can be derived from either clinical disease progression, alterations at the electrodeneural tissue interface [28], or neuroplasticity [29]. For what concerns the last two, over time variability can be caused by different mechanisms: (i) electrodes' movement (ii) AS' encapsulation due to inflammatory foreign-body reactions which cause instability and changes in impedances [28] (iii) neuroplasticityinduced changes [29].
In adaptive neuromodulation settings for progressive neurodegenerative disorders (where the disease worsens over time), electrical parameters need to be re-adjusted to obtain an optimal clinical outcome [30]. After identifying neural markers that correlate with clinical symptoms (feedback variable), stimulation is tailored in a closed-loop manner to deliver on-demand stimulation that reduces side effects, while increasing therapeutic efficacy. Significant progress in this direction has been made for essential tremor and Parkinson's disease [21,[31][32][33][34][35][36][37][38][39], where with the pioneering study of Little et al the feasibility of an adaptive Deep Brain Stimulation (aDBS) in Parkinson's patients driven by subthalamic nucleus local field potentials in the beta-band was validated [39]. Subsequent studies have shown the superiority of aDBS over continuous DBS for motor symptom improvement, with better control over levodopa-induced dyskinesia. While such systems are finding increasing validation in the neuroscientific community for movement disorders, application in other diseases, such as cognitive and psychiatric ones, is still under development [33], although recent successful applications have been proposed for treatment-resistant major depressive disorder [40] or obsessive compulsive disorder [41]. Similarly, in drug-resistant epilepsy, responsive neurostimulation was proposed [42], where brain stimulation is driven by detection of specific electrocorticography patterns. alternatively, vagal nerve stimulation is triggered by detection of tachycardia [43].
In sensorimotor neuroprosthetic studies [22,44], where the restoration of targeted motor/sensory patterns is obtained through stimulation of motor/sensory structures of the nervous system, the issue of developing a faster and consistent method to improve the efficacy of neuromodulation becomes multifold. While in aDBS, based on the neural marker of interest, stimulation is turned ON or OFF, or amplitude of it is increased/decreased upon need, (sensory-)motor neuroprosthetics require the tuning of many stimulation parameters such as pulse amplitude (PA), pulse frequency, pulse width (PW) and activation of specific sites of stimulation. Depending on which neural interface is being used, such as a Utah Array in the motor/sensory cortex for brain stimulation [45], or polymide-based intraneural electrodes interfacing with the peripheral nervous system [8,46], up to hundreds of AS might need such tuning, with the space of parameters becoming exponentially larger. Indeed, to obtain a certain motor or sensory output at a specific time, co-activations of different targeted muscles or sensory fibers might be needed, for instance agonist-antagonist muscle activation patterns in flexion-extension movements; or in recovering the sensation of grabbing an object which requires more hand areas to be felt. Previous literature in the neuroprosthetic field has experimented the use of Bayesian optimization algorithms for neuromodulation applications in motor prosthesis in animal models [22,44]. In [22], the search for a targeted stimulation pattern of the motor cortex of a nonhuman primate was optimized, by using a hierarchical approach in a multi-electrode application. The authors have exploited the knowledge of the single channel input to output relationship (in this case, stimulation to electromyography (EMG) pattern) and translated it to a more complex, multi-channel level; obtaining a significant improvement on the observations needed to optimize the stimulation protocol with respect to the space of all parameters. Following this, further validation found the optimal stimulation protocol to elicit motor patterns for walking control (in rats) and hand grasping (in monkeys), which are tuned offline reaching a successful performance and converging by exploring only 20% of the search space of all the possible neurostimulation parameters. An online validation in one monkey was achieved [44]. In this application scenario, multiple EMG measures were used as feedback to the algorithm.
Sensory feedback restoration (the subject of the present study) adds a layer of complexity: unlike aDBS or motor prosthetics, where the output can be directly observed by objective motor effects and quantitative measures such as electromyography or local field potentials, sensory feedback is not only subjective but also more complex to follow. It can only be qualitatively assessed by subjects' verbal report at discrete time points and the nature of the feedback does not allow for its continuous tracking with current technology.
These studies provide important benchmark and proof-of-concept towards the usability and efficacy of such algorithms for prosthetics calibration, where the relationship between input and output needs to be learned without any specific knowledge a-priori. Nevertheless, similar approaches able to tackle time variability and the subsequent need of prosthesis recalibration are lacking. For this reason, we decided to exploit Gaussian process-based Bayesian optimization (GPBO) [47,48], (here after referred to as GPBO) for re-calibration purposes. Indeed, GPBO is considered ideal for multi-dimensional problems, where the relationship between the input stimulation parameters and the physiological output is not known a priori. It is applicable in online scenarios of calibration, where the input-output relationship can be experimentally observed. Here, we considered the data from three trans-femoral amputees who were implanted for 3 months with four transversal intrafascicular multichannel electrodes (TIMEs) [49,50] in the tibial branch of the sciatic nerve with the aim of restoring lost tactile and proprioceptive sensations at the level of the missing foot and leg using neuromodulation [2,51,52]. By using this human data, we tested if the AI-based algorithm was able to automatically adapt to the stimulation changes over time maintaining the device stable and effective. We exploited GPBO to find the optimal stimulation protocol, by re-adjusting the previously obtained stimulation. This is done in a progressive manner throughout the days of implantation and by integrating the past information into the predictive process. We built the frameworks for induced threshold and location variations, both translatable to a variety of application fields of neuromodulation where temporal recalibration due to foreign body reactions (FBRs), micro-scale implant's movement, or to neural plasticity is required.
In the first framework (sensation location optimization), we tried to build a model able to identify four AS covering the three functional groups of the foot (frontal, medial, and heel) and leg with the minimum number of stimulation trains possible for a somatosensory neuroprosthesis. The locations elicited by the same AS can change over time [25-27, 53, 54]. Thus, the device needs to be manually re-tuned for optimal functioning. The information gathered from the previous experimental days was integrated into an automatic process as existing observations, with the underlying hypothesis of exploiting, when possible, the stability of some AS.
In the second framework (perceptual threshold optimization), we propose a method able to generalize and adapt to the needed injected charge to electrically elicit an artificial sensation (i.e. perceptual threshold). Indeed, one of the major problems in neuromodulation devices is the need for readjustment of the injected charge values [26][27][28][55][56][57]. This means that continuous recalibration is needed to keep the prosthesis efficient and functional as its intended purpose. More specifically, we wanted to understand how many stimulation trains were needed, when using the GPBO, to find the new optimal perceptual threshold, exploiting the previously functional ones.
By comparing the time needed by these designed models to automatically recalibrate the stimulation parameters, we propose here a novel algorithmic procedure able to tackle the issue of device recalibration [58][59][60]. Interestingly, this study reports how AI-based Bayesian methods could be effectively used to predict the over-time changes in neuromodulation parameters for sensory feedback restoration. Upon tailoring and validation in other application frameworks, this could enhance medical device efficacy and reliability for treating neurological disorders using targeted and adaptive electrical stimulation.

Somatosensory neuroprosthesis based on neural stimulation
Three trans-femoral amputees were included in the study, being the same participants as in [2,61] implanted with four TIME [46] electrodes in the tibial branch of the sciatic nerve for 90 d. Importantly, this is a retrospective study and thus the optimization was done offline.
Each TIME electrode is composed of 14 AS and two grounds (for a total of 56 AS per patient). Trains of bi-phasic symmetric and square pulses were delivered to the nerve, where the PA and the PW were tuned in order to have the subject feel an electrically evoked sensation [2]. The perceptual threshold (i.e. minimum injected charge to evoke a sensation) was saved, together with all the other properties of the reported sensation (i.e. the perceived sensation location, quality and intensity). This sensation characterization was repeated for all the 56 AS.
To characterize the perceptual threshold, an increasing ramp in PA [mA] or PW [µs] was performed by the expert. When the subject felt a distinct sensation, the ramp was stopped and the current value of the injected charge Q = PA·PW [nC] was saved, referred to in literature as the perceptual threshold [62]. While the stimulation was increasing in charge, subjects were asked to report the moment they could feel a distinct sensation, and where this sensation was located. In this study, the charge needed to elicit a distinct sensation (referred to as perceptual threshold, PT) and the tested AS eliciting a certain location were considered as inputs.
The first characterization was defined as the complete characterization of all the functioning AS (those eliciting a sensation with a charge Q < 120 nC maximum injectable charge in the TIME electrode). Since the sensory characterization of all the functioning AS may be performed on multiple days (as being intrinsically very time-consuming), for each functioning AS we gathered the data deriving from the first time it was functionally tried, which does not necessarily correspond to the first day of mapping. Overall, most of the AS were tried within the first week postimplantation. Subjects have been followed up to 90 d post-implantation. Multiple sensation characterizations have been performed throughout this period, with the purpose of collecting and analyzing the time evolution of the evoked sensations and their related stimulation parameters. Indeed, several changes in the optimal stimulation parameters were observed over time [63]. This was due to the implant's micromigrations, biological reactions (i.e. FBR), or sensory adaptation. Thus, the output of the sensation characterization might change regardless of the same combination of input stimulation parameters. With the use of a custom-made mapping GUI, subjects could press on a virtual foot and leg where they were feeling the sensation [62]. The plantar foot was divided into rectangles, and the leg was divided into four main locations [2]. Based on the resolution of the foot, the foot image was divided into a 250 × 450 pixels. Each rectangle's vertex could thus be defined thanks to a system of (x, y) coordinates, representing the width and the height of the image respectively, with x in (0,250) and y in (0,450).
A schematic representation of the framework is shown in figure 1.

Location and intensity: inputs to the algorithm
With the aim of describing the location in a more continuous and intuitive manner, the centroid (i.e. the middle point of the rectangle) of the location was computed as a (x, y) pair. Although the distinction between the medial and the lateral part of the foot is important when considering functional standing and walking, we considered the neuroprosthesis as functional if able to deliver a localized sensation on the frontal, middle, and heel part of the foot in a distinguishable manner. To restore a natural pattern, it is crucial that the subject can perceive in a distinct sensation these three areas [64].
Previous studies exploring sensation restoration in lower limb prostheses reported a functional division of the foot in three main areas: front, mid and heel, with the proprioceptive sensation eventually added as a fourth property of sensation [2,51,52,64,65].
Therefore, only the y-coordinate, spanning from the front to the heel of the foot, was considered as the output able to encode for the sensation location.
Firstly, we confirmed that the level of the injected charge was not influencing the centroid of the perceived location, but rather only on the area. Indeed, it is known that a higher charge is responsible of a more spreading sensation [56]. For this reason, only the AS was considered as the input parameter for analyzing the input-output relationship for the location. Moreover, each AS targets a bundle of fibers innervating different areas of the foot sole and the leg. Different AS could interface different bundles of fibers, thus eliciting different sensation locations. In the same way, due to the position (intra-or extra-fascicular, closer or further from certain fibers) and the potential FBR surrounding an AS, we expect to obtain different perceptual thresholds regardless of their elicited location. For this reason, each AS was considered independently, so that only the charge was considered as the input parameter for the intensity of the felt sensation.

GPBO algorithm
Given the intrinsic variability of the input to output mapping (stimulation parameters/evoked sensation) over time, the need for a continuous recalibration of the neuromodulation arises. This means that, each time, the correct combination of electrical stimulation parameters needs to be found. Recently, the problem of finding the best set of parameters to optimize a certain function is being tackled with AI algorithms. More specifically, GPBO algorithm [47,48] is the state of the art in this context [22,36,66,67]. Automatic re-calibration of neuromodulation parameters. Nervous tissues (at the central or peripheral level) are interfaced with bio-compatible devices able to deliver a custom, targeted pattern of stimulation to induce an inhibitory or excitatory effect on the surrounding neurons. In the present study, a TIME electrode is used to interface the tibial branch of the sciatic nerve of a transfemoral amputee to restore sensations from the missing limb. Each active site of each electrode has to be properly tuned in order to induce a target sensation in terms of location, intensity and type perceived. The time variability of the input-output mapping is a huge problem in neuromodulation: regardless of using the same input parameters, the output is different over time. In somatosensory neuroprostheses, this happens for the elicited sensation's location (which migrates, although the same AS is being used) and the perceptual threshold: i.e. to induce a sensation, a generally increasing amount of charge is needed over time. While the traditional approach relies on a human expert-mediated re-calibration of the prosthesis, which is time consuming and increases the health system cost, we propose an automatized approach, based on a Bayesian Optimization algorithm, able to automatize the recalibration by integrating past information.
Bayesian Optimization is a machine-learning method focused on maximizing (or, without loss of generality, minimizing) an objective function and solving the problem: where: • f (·) is expensive to evaluate, meaning it takes time (or monetary cost) to evaluate what value f has given a choice of parameters x.
• f (·) does not have any known specific property (i.e. concavity, monotony or others) that eases the optimization process (f is black box).
We want to find the global optimum, avoiding local maxima (or minima).
The GPBO algorithm is characterized by two phases: in the 1st phase, the objective function f (·) is modeled with a surrogate. In the second phase an acquisition function suggests where to evaluate the function next (i.e. 'query' or 'observation' point), with the final purpose of maximizing the objective function in the minimum number of evaluations possible. Indeed, not all the possible observations are equivalently useful and only few are promising. The way the objective function is initially modeled is with what is called a Gaussian process prior: a Gaussian distribution with a certain mean µ(x) and covariance Σ(x), being a probability distribution of all the potential values our function f (x) might have at a point x. The covariance (or kernel) chosen for the task at hand is the squared exponential kernel, or radial basis function (RBF) [48]. The RBF kernel is expressed as follows: with σ s being the input signal variance and λ being the length-scale. The signal variance σ s is a measure of the amplitude range of the function, indicating what values the function might span. The length-scale λ indicates how smooth the function is. A large lengthscale indicates that the function needs a big change on the x-axis to have a meaningful change on the yaxis. Another value defines the fit of the Gaussian process, the noise variance σ n , which indicates how much uncertainty around the observed values needs to be considered. This is also measure of how noisy our function is (i.e. how much the output of the function can vary with the same input). According to previous studies, these three parameters have been tuned in a guided manner by analyzing the experimental data at hand [68]. More specifically, the signal variance was defined as The length-scale was experimentally found by analyzing the input-output relationship. A good initialization derives from analyzing the range of the input X, with and λ ∈ [1,10]. λ can be tuned based on the expected linearity of the input-output mapping, with λ = 1 indicating a strong linearity. Finally, the noise variance was defined as the maximum variation of the function value given the same input (i.e. what different outcomes the same input can have). A common initialization for the signal noise derives from the consideration of the expected signal-to-noise ratio κ, setting the signal noise to with κ ∈ [2, 100] . [68]. These hyperparameters are only used as initialization values to build an initial instance of the Gaussian process prior distribution, and are then optimized by the algorithm with the collection of new data points ( figure S2). Importantly, all the hyperparameters have been initialized according to the implant of subject 2, which had the most comprehensive and realistic sensory neuroprosthetic framework, with the underlying hypothesis of similar physiological and computational processes underlying the nature of the two other implants. This means that the initialization of these hyperparameters was the same also for subject 1 and subject 3.
Each time a new observation is made (i.e. the function has been evaluated at x i ), the posterior Gaussian process distribution is updated (posteriorupdating): the mean will interpolate all the evaluations and the confidence interval (containing f (x) with a probability of 95%) will be zero in the observed points.
The reasons why this algorithm was chosen among others are varied and can be divided into the properties of the problem, the algorithm and the neural interface.

Properties of the problem (mapping of electrical parameters to a certain perceived sensation):
• Objective function (i.e. patient's reported sensation to a given stimulation pattern) is expensive to evaluate in terms of time • Black-box and derivative-free problem (we do not expect any specific property of this function a priori) • Parameter space is multi-dimensional and characterized by complex inter-parameter relationships

Properties of the algorithm:
• Suited for online-implementations • Adaptive-scenario of incremental-learning (information is subsequently integrated) • Successfully selected and validated for closely related applications

Properties of the neural interface:
• Inter-subjects variability of sensations • Implant's placement within the nerve is not consistent among implants • Unknown framework of afferent nerve stimulation • Time-variability due to neurophysiological phenomena

Evoked-sensation location search
The purpose of the Bayesian optimization algorithm is to find the optimal set of parameters x that maximizes or minimizes an objective function f (x) in the minimum number of queries (evaluations) possible. In our case, the objective function to minimize can be represented as the difference between a given target sensation, which we will call S * , and the sensation elicited by the electrical stimulation with parameters x, S(x). The problem, we are faced with, is to find the set of optimal parameters, x * , such that |S (x * ) − S (x)| is minimized. When optimizing the sensation location, we aim at finding the optimal set of parameters able to elicit a certain target location L * . Four target locations have been chosen: front, mid, heel and leg (figure 2(A)), which are the ones established in the literature and in previous embodiments in order to have a functional somatosensory prosthesis [2]. More specifically, with the purpose of finding a single global optimum, the most-frontal, the most-backward, the closest to the medium line and the closest to the leg locations have been considered. The Y-coordinate of the centroid was used, as it spans through the heelfront line of the foot (referring to figure 4(C)). For the leg, the Y-coordinate of the centroid was set as equal to 1.5. Thus, the objective function to minimize is where L * is the target location. The objective function changes based on the target pattern; thus, it is being recomputed each time when the target location changes. The space of parameters Ω, with x ∈ Ω, is represented as x = (AS) , with AS being the active site we are referring to. What the acquisition function will point out is thus the optimal choice of AS to activate in order to reach a certain target pattern.
A key factor for the success of the algorithm is the choice of the hyperparameters and the covariance of the Gaussian process (i.e. the kernel). This determines the generalization properties of the algorithm as it stands as a solid knowledge on the input-output relationship when modeling the surrogate function.
A squared exponential kernel has been chosen for the task. This kernel has good generalization properties and is integrable for most of the functions [48].
For the signal variance σ s we analyzed the values that the objective function could take based on the target value. Since the objective function is the absolute difference between the target location S * and the elicited location, this could take a maximal value of 1.5 (if the target was the leg and the actual location was the front, or vice versa). The minimal value was 0, when the two locations were corresponding, leading to a final (approximated) value of σ s = 1 (figure 2(A)). For the length-scale, a visual inspection of the data allowed us to define an appropriate value, able to catch the scale of variability on the xaxis (i.e. of how many AS should we move) in order to have an independent outcome (a new location), defining a value of ℓ = 10 as the optimal one. Finally, for the signal noise σ 2 n , we have set a κ = 2 given the strong noise σ n we expect in the output function, and the low signal-to-noise ratio κ. Based on the above mentioned hyperparameters value, a representation of five samples deriving from the surrogate function is reported in figure 2(B). A graphical representation of the concept of length scale is found in figure 2(C), where arrows are indicating the change in the input-scale needed to obtain independent values on the output-scale.

Data preparation: AS functional resorting
Different type of inputs can be fed into the algorithm: continuous, discrete, or categorical. In the first two cases, close input values are expected to have close output values (always according to the kernel properties), with input values being further away presenting more differentiated behaviors. When considering categorical values instead, inputs are considered as independent. Close AS belonging to the same electrode are not expected to have similar properties, due to the electrode's insertion and the type of fibers they interface. In the same way, AS belonging to different electrodes do not have similar properties. This leads to the conclusion for which no immediate correlation exists when considering physically close AS, and AS could be considered as categorical and independent inputs. On the other hand, we can expect some AS to present a similar behavior based on the output sensation that they elicit (i.e. they interface similar/the same bundle of fibers). Following this principle, the first characterization day, which is hypothesized to be carried out by an expert, was used in order to give a new sorting to the AS. More specifically, AS have been ordered according to the location of the sensation that they elicited. Thus, close AS are AS which elicited close locations, and the input was considered as discrete. This ordering, deriving from the first characterization, has been maintained throughout the whole period of post-implantation follow-up (figure 2(C)) and was defined as an AS resorting.

Time integration of the information in the GPBO
Previously collected information can provide meaningful insights into the stability of a certain active site. Moreover, when searching for a certain target location, AS which are more often eliciting that location are more probably successful. For this reason, we implemented a GPBO algorithm where progressively collected information regarding the elicited location by a certain AS is maintained and integrated into the process, building the posterior distribution. In this way, each time an AS is suggested by the algorithm, regardless of whether it is successful or not, the location elicited by that AS is saved. Thus, AS which are strongly unstable over time, will be characterized by a high uncertainty and noise. On the other hand, AS which are more stable will be characterized by a low uncertainty. The algorithm has been performed on each characterization day by starting from the previously collected information. A representation of the posterior-updating obtained after the first mapping is found in figure 2(D), where observations (centroid elicited by each sorted-AS) are integrated in the Gaussian process, with a target of 0.5. Samples of the surrogate derived from the posterior can also be visualized. The concept behind the GPBO is represented graphically in figure 2(E).

Reset GPBO
Another GPBO framework has been implemented with the aim of understanding the importance and potential added value of integrating past information. In this case, on each characterization day, three stimulations have been performed randomly and a new process has been built, without considering any previous information. This was called 'reset GPBO' (rGPBO).

AS random and spatial resorting
In order to validate our choice of resorting the AS according to their elicited location in the first mapping, we have implemented two further GPBO frameworks: in the first one, AS have been resorted according to the order in which they were tried in the first mapping, and we refer to this resorting as Random resorting (indeed, no information can be gathered in this case on AS similarity). In the second one, AS were resorted according to their spatial relationship, meaning their placement on the electrode shafts. Although our AS resorting procedure leads to a distortion of AS original physical relationship to re-order AS based on their functional relationship (i.e. similar elicited output location), we hypothesize that this method is justified by the neurophysiological setting underlying this application. Indeed, differently from other neuromodulation settings (where physically close AS lead to similar functionality [22,44]) in this case the organization, curvature and branching of fascicles within the nerve does not guarantee that close AS are interfacing similar fibers (in terms of thresholds and projective field). For instance, an AS could be intrafascicular and interface one bundle of fibers innervating the frontal part of the foot, and the nearest AS could be extrafascicular and activating fibers innervating the heel. Finally, since the four electrodes are implanted at different locations in the nerve, there is no relationship between AS belonging to different electrodes.

Assessment
In order to assess the performance of the algorithm, we compared it to the random search algorithm, randomly selecting AS, and to a reset-GPBO, gathering information on the importance of temporal integration of previously collected data. Moreover, we validated the importance of data preparation by comparing the AS functional resorting with the random and the spatial resorting. All the comparisons were made by running post-hoc Tukey-HSD tests with α = 0.05. We compared the number of iterations needed to find all the targeted locations with a stopping criterion of ∆C y < 0.15, repeating the procedure =30 times.

Stability and consistency across AS
We classified the AS into different categories based on a stability index. For each subsequent characterization session, if an AS was kept in the same functional group (either front, heel, mid or leg) then it got a score of 1, else a score of 0. We computed the mean score throughout weeks for each AS and defined as 'unstable' those AS with a score <0.3, 'stable' those with a score >0.6, and 'variable' all the remaining. Also, for each session, we computed the percentage of stable AS, leading to an overall representation of the stability of the implant in terms of perceived location. Importantly, when considering the stability score as a reference for the expert formalization, this has been computed by using the same criterion as in the algorithm (we computed the percentage of AS having a shift in the y-centroid ∆C y < 0.15, rather than belonging to the same functional group). A first consideration regarding the relationship within AS was made, with the specific aim of understanding whether physically close AS (placed on the same electrode shaft and closer together) are more likely to elicit closer locations. Somatotopy is well known and documented at the peripheral level, and this property increases distally: if considering for instance the peroneal nerve innervating the top of the foot, closer fibers in this nerve will most likely innervate closer parts on the foot. At a proximal level though, where subjects have been implanted (tibial branch of the sciatic nerve), the somatotopic fascicle organization and the clustering of fibers innervating the same skin are less prominent [69]. The placement of the electrode, which is strongly dependent on the surgery, variable within subjects and not known a priori also makes it difficult to anticipate what fibers an AS might interface.
Finally, the TIME electrode insertion, which is not straight but rather curved (thanks to its flexible thin film), increases the chance of the electrode to effectively interface different functional groups at once [46]. In order to prove such hypothesis, and finally reject any inter-dependence within AS due to physical vicinity, for all the 4 electrodes implanted in patient 2 we computed the city block distance defined as between all the pairs of the AS (i, j) of the same electrode within the same day and repeated for all the characterization days. With an electrode shaft as represented in figure 4(C) (with two rows and seven columns), the maximal city block distance was . This result clearly derives from the spatial distribution of the AS on the electrode. We correlated this measure with the ℓ 1 norm of the difference of the centroids (i.e. the locations) elicited by the two AS (i, j) (i.e. ∆C y = C y i − C y j ), with the aim of understanding if physically closer AS can be correlated to physically close locations.

Threshold search
Once an AS has been found as the one eliciting the target location, the perceptual threshold charge needs to be found in an adaptive, fast manner. A standard way to proceed is by making ramps of increasing charge (either modulating the PA or the PW) and by asking the subject when a sensation is felt. Although this method is safe, it requires significant time. Moreover, in the same way implants are unstable in terms of the location they elicit, also the perceptual threshold charge is an unstable property, which undergoes significant changes throughout time. These changes need to be addressed, requiring the need to continuously re-calibrate the neuromodulation. With the aim of addressing this issue, a GPBO algorithm has been fitted to this problem. More specifically, the algorithm had to find in the minimum number of iterations possible the closest value to the real perceptual threshold. Importantly, for the intrinsic definition of perceptual threshold, when Q < Q * , with Q * being the perceptual threshold, no sensation can be felt, thus we cannot gather any information regardless of how big the difference of perceived intensity is with respect to the threshold. This leads to our objective function, which is here reported and represented in figure 3(A), where Q * is the target perceptual threshold. In a real application scenario, we do not know a priori the target perceptual threshold, but we can have an online verbal report of the intensity felt by the subject when stimulating with a certain charge. In this specific case, since no stimulations were performed apart from the ones at the perceptual threshold, we built an algorithm whose cost function was based on the perceptual threshold itself, underlying the true value. The space of parameters Ω, with x ∈ Ω, is represented as x = (Q) , with Q being the perceptual threshold charge, in a range [0,120] nC. What the acquisition function will point out is thus the optimal choice of charge to use in order to reach the perceptual threshold. Differently from the location framework, where the observations suggested by the algorithm (and the whole search space) corresponded to experimental data we were provided with (we always knew the elicited location by a certain AS), in this case the algorithm could explore the whole search space. No knowledge was provided on what outcome a certain charge could give, meaning that the observations suggested by the algorithm were not always necessarily included in our experimental data and needed to be extrapolated. Regarding the signal variance σ s , the resulting distribution from our experimental data was of 256.79 ± 359.95, with a median of 121.0. Thus, the value 120 has been chosen as an optimal value for the signal variance ( figure 3(B)). Concerning the length-scale ℓ, the resulting value for the distribution was of 26.69 ± 12.79, with a median value of 32.0. Thus, an optimal value of 30 has been chosen. More specifically, the parameter λ was set to 1, given the high linearity of the underlying objective function ( figure 3(B)). Finally, the signal noise σ n was set to 1.20, given a scaling factor of k = 10 with respect to the signal variance (assuming thus a signal-to-noise ratio equal to 10). Indeed, we expect the function of having a low noise: it is very unlikely that, in a concrete application scenario, the same input charge leads to two different intensities felt [70]. Samples deriving from the surrogate with these hyperparameters are represented in figure 3(C).

Safety factor
The maximum injectable charge with a TIME electrode is 120 nC, and a linear correlation between the injected charge and the perceived intensity has been extensively proven in the literature [2,55,56,70]. A very high perceived intensity could correlate with an unpleasant, painful sensation; therefore a safety factor was introduced in the algorithm search such that the search space for each AS changes based on the previous threshold value. More specifically, we computed the ratio PT (t) /PT (t − 1) for each AS and for each consecutive characterization day. From the distribution of all possible values, which we will refer to as θ, we considered the 0.05 quantile and the 0.95 quantile as the values for the lower and upper bound of the search space ( figure 3(D)). In this way, the new domain was defined separately for each AS and dynamically adjusted on each new characterization day as follows: With θ being the distribution of all the possible ratios PT (t) /PT (t − 1), and Q θ (x) indicating the x-quantile of the distribution. The obtained search space distribution is represented in figure 3(D).

Classical GPBO with previously found optimal threshold
Differently from the location framework, where the information is consequently integrated, in the threshold search a new Gaussian Process is created on each Characterization Day. The previously found, successful perceptual threshold value is tried firstly, and further stimulations suggested by the algorithms are integrated to build the posterior distribution ( figure 3(E)). If this value is correct, the value for that day was found. Else, the domain is narrowed down, based on whether PT(t − 1) is leading to a too strong sensation or to no sensation, towards lower or higher values respectively based on the safety factor. The two possible domain choices are thus The next stimulation tried by the algorithm is the other extreme of the domain, meaning respectively Q = Q θ (0.05) · PT (t − 1) or Q = Q θ (0.95) · PT (t − 1). Based on the intensity of the sensation that this charge delivers, the domain can be enlarged up to its minimum or maximum (for instance, if Q = Q θ (0.95) · PT (t − 1) leads to no sensation, the domain is enlarged up to 120 nC). The search stops when the algorithm reaches a value which is less than ±5 nC from the target value ( figure 3(F)).

Performance assessment
In order to assess the performance of the algorithm, the mean number of stimulations performed (observations) was assessed throughout all characterization days and for all AS. To show robustness in the number of iterations needed by the algorithm, we have run the simulations N = 30 times, obtaining a distribution of number of iterations for each mapping day.
We aimed at understanding, indeed, how many stimulations we needed to deliver for each AS and on each day in order to find the perceptual threshold. Also, since a random search is not plausible in this framework, this method was not considered: indeed, neither when searching for the optimal charge to elicit a sensation, nor when searching for the optimal AS that activates a percept in each area, the expert will likely proceed in a random manner (for example, by randomly selecting charges in the interval between 0 and 120 nC).

Integrating temporal information
Integrating previously collected information can be critical in order to better find the optimal perceptual threshold value. Indeed, if an AS has a stable behavior throughout time, and thus has a perceptual threshold value which only undergoes slight and small variations in time. Then, the search could be facilitated by this information understanding the need to perform more exploitation around that value, without the need to perform exploration in other ranges.

Classification of the perceptual thresholds time series
A k-Means algorithm on the time series characterized by the evolution of the perceptual thresholds of all the AS taken into consideration has been performed as in [63]. The aim was to explore whether AS could be grouped based on their stability behavior with respect to their perceptual threshold change in time. The optimal number of clusters was assessed by the elbow method, choosing the optimal value of k where the decrease in the inertia started to show a linear trend.

Code availability
To support the methodological background of this paper, we are providing the code online (Github link), where data can be processed according to the methodology proposed. More specifically, we are providing the AS resorting procedure (2.4.1) and the codes for location change adaptation and threshold change adaptation.

Adaptation to location change
Patient 2 has been considered for this analysis both in terms of AS and characterization sessions. Noticeably, patient 2 was characterized by a very stable implant, which maintained a high number of AS until the end of the clinical trial (90 d postimplantation). Throughout all the characterization weeks, the foot area was vastly covered by the functioning AS ( figure 4(A)). An imbalance between the AS covering the most frontal part of the foot and the ones covering the heel part of the foot, which was found experimentally in our data, is in line with the previously found proportion of cutaneous afferents innervating the plantar area of the foot [71] (figure S2). Although it was possible to elicit a sensation throughout all the functional areas of the foot, the location reported by the activation of a single AS was not always stable, with some AS eliciting sensations in different parts of the foot throughout the characterization weeks. Out of 54 AS, only four presented a very unstable behavior. 25 AS were showing a stable behavior, belonging to the same functional group in at least 60% of the consecutive mappings, and 25 AS were showing a medium behavior, belonging to the same functional group in 30%-60% of the cases ( figure 4(B)).
We then assessed whether the physical distance of AS was related to the respective difference in elicited locations (i.e. close AS are eliciting close locations in the foot). The Kendall's tau correlation coefficient, used as a quantitative non-parametric measure to assess the numerical correlation between the two variables, was R = 0.102 p = 5.28e −9 confirming the weak inter-dependence between these measures (figure 4(C)). Similar results were also found when considering each electrode separately. The choice of the AS resorting procedure according to the first characterization day was further validated by assessing how closer AS do not experience more similar levels of instability (τ = 0.07, p = 0.001), while AS eliciting closer location on the foot have more similar levels of instability R = 0.37, p = 2.23 · 10 −36 (figure S1). Finally, when comparing how the AS resorting affects GPBO results, we found that AS resorting according to the first characterization leads to better performance than when no sorting is performed or sorting is performed according to physical distances within the electrode (p < 0.001) (figures 4(D) and (E)).
We computed the stability score of each consecutive mapping based on the previous mapping and the first characterization day (figure 4(F)), leading to results of 0.57 ± 0.07 and 0.40 ± 0.07 respectively. We then assessed how our algorithm would perform. We used a Random Search method, the classical rGPBO and the GPBO. Respectively, the mean number of iterations needed for the three methods is: 10.09 ± 5.23, 12.59 ± 4.48 and 6.30 ± 2.28. Only the performance of the GPBO is statistically better than the other two methods (post-hoc Tukey-HSD, p < 0.05) (figure 4(G)). By analyzing how the distribution of the number of stimulations needed to cover all the functional areas could change based on the Mapping Day ( figure 4(H)), we could see how apart from day 36 and day 57, that differed significantly from the remaining days (post-hoc Tukey-HSD, p < 0.05) and required a mean number of 3 stimulations, in the majority of the characterization days, only 1 or 2 stimulations were needed in average to complete the mapping ( figure 4(H)). The summed Euclidean distance of the AS ordering from the first characterization day (as a measure of how different the actual AS sorting is from the initial sorting) was also analyzed. The AS (i.e. as ordered by the first characterization) chosen by the algorithm as the optimal AS are represented in figure 4(I). As visible in the figure, the ordering of the AS has proved useful, as AS that are close to each other are picked, meaning that a significant level of stability of the implant was reached. An illustrative example of the research made by the algorithm is shown in . This is shown for a front/mid and heel target. An example of bad research is provided.

Instability analysis
By analyzing the experimental data collected of patient 2, who had a high implant's stability over time, we computed the percentage of AS which could be considered as stable for the perceptual threshold change with respect to the previously performed mapping. As visible in figure 5(A), the percentage of stable AS reduces drastically at week 3 postimplantation, with a mean value of 17.28 ± 9.49. We computed the mean change in charge ∆Q = |Q (t + 1) − Q (t)| needed for each AS with respect to the previous perceptual threshold, which led to a result of 8.65 ± 10.59 nC on 44 AS. The high standard deviation is a measure of how much AS can have different behaviors. By computing the average PW used during the experiments, being of 60 µs, we obtained that an average 144.13 ± 176.56 µA change with respect to the previous mapping session was present.

Categories of AS stability
The results of the k-means algorithm on time series of perceptual threshold evolution are represented in figure 5(E). The optimal number of clusters (k = 3) was assessed by the elbow method. The algorithm highlighted the presence of three behaviors: a Stable one, characterized by low thresholds (<40 nC) and a smooth dynamic; a Variable one, characterized by thresholds <80 nC and steeper slopes, and an Unstable one, characterized by high thresholds (>80 nC) and more unstable dynamics.

Safety factor
In order to define the safety factor, and thus the multiplicative factors defining the lower and upper bound of the domain, we considered the distribution of the ratios PT (t) /PT (t − 1). Results are shown in figure 3(C). The values of the 0.05 quantile and the 0.95 quantile were respectively Q θ (0.05) 0 = 0.64 and Q θ (0.95) = 2.40, leading to the following domain of research, which was dynamically updated for each AS and for each day: Importantly, although the algorithm research domain has been narrowed down according to the safety factor, the comparison with the expert is still valid as the domain spanned by the algorithm is much larger than the domain spanned by the formulated expert ( figure S3).

Results of the algorithm
All the three patients were considered for this analysis. We analyzed how the algorithm was performing by computing the distribution of the number of stimulations (included the previous threshold value) needed to find the correct threshold on each day. We obtained values of respectively 2.85 ± 1.93, 2.60 ± 2.02 and 4.45 ± 2.87 iterations needed for each patient ( figure 5(B)). An example of the research proposed by the algorithm is shown in figure 5(C), with all the consequent observations made before finding the optimal value. We can observe how, starting from the previous perceptual threshold and after an initial exploration phase, the algorithm chooses in just two steps a perceptual threshold charge which is very close to the target one. When assessing the number of iterations needed to find the optimal perceptual threshold throughout following mapping days in subject 2, we could assess how this number was stable throughout the mapping (post-hoc Tukey-HSD, p > 0.05), highlighting the reliability of the algorithm in finding the optimal threshold in a rather consistent number of stimulations ( figure 5(D)). In order to assess how the algorithm performed in finding the optimal Perceptual Threshold for different types of AS behaviors (as reported by the k-Means algorithm, figure 5(E)), we computed the number of stimulations needed for these different categories of AS and for patient 2. The number of stimulations needed for each category is respectively 1.81 ± 1.55, 2.78 ± 1.87 and 3.47 ± 2.36 for stable, variable and unstable AS (figure 5(F)). A post-hoc Tukey-HSD test showed that stable AS need consistently a lower number of iterations with respect to medium and unstable ones (p < 0.001). This shows that regardless the stability of the AS, the algorithm is always able to find the threshold with a low mean number of iterations in a consistent manner.

Discussion
The issue of time variability of the effective stimulation parameters leads to poor long-term stability of neural implants [8]. Biological tissue reactions and device failures bring the need for a continuous recalibration of neuromodulation devices. Indeed, this is an important aspect to tackle when considering the wide-spread usability and efficacy of such technologies in clinical practice. Although advanced neurostimulation methods are attracting growing interest in the neuromodulation community, with innovative approaches being proposed and validated both offline and online, limited attention has been put on their stability problem. It is as important for proper device functionality and carries the further complexity of considering how past information might be integrated. To address this, we explored how GPBO, a method which is being increasingly exploited in neuromodulation scenarios, can tackle the issue of time variability and the need for re-calibration. Indeed, it can be used as an automatic way to generalize and adapt to input-output changes over time, comparable to an expert technician, while reducing the time needed for use and its cost, while increasing the usability of these devices.
With the data from three transfemoral amputees implanted with neural interfaces for up to 3 months [2,51,64], we presented an offline validation for a method which can serve as a useful tool for further applications of this technology in the neuromodulation field. For instance, one of the major problems in neuromodulation is the need to continuously adapt the injected charge required to elicit the same type of electrically-induced output [9,26,56,72]. In our specific case, the charge needed to induce a sensation in subject (i.e. the perceptual threshold) is known to change throughout time for different mechanisms such as sensory adaptation, fibrotic encapsulation of the active site, interface of different fibers and movement of the AS with respect to the neural structures it targets [8,63,73,74]. This is a critical issue for the proper functionality of any device whose purpose is to modulate neural activity and, if not tackled, can render the interface non-usable or inefficient.
In our method, we assessed if a GPBO was able to model the variability of this Perceptual Threshold and to suggest optimal stimulation parameters, based on those previously effective, and on further observations suggested by the algorithm itself. The goal was to minimize the difference between the automatically chosen perceptual threshold with respect to the target perceptual threshold (the one effectively eliciting a sensation on that mapping). Our method enabled to rely on a very small number of stimulation iterations to obtain the new threshold. Moreover, the safety factor allowed to restrict the stimulation, so that the algorithm is not allowed to explore the entire search space but only a safe range defined from the real experimental data. The approach proved useful in all the subjects, validating the consistency of the method. It also showed the robustness with respect to different type of AS behaviors: from the ones presenting a very stable behavior (with a lower variability of the perceptual threshold evolution) to the very unstable ones (with a very unpredictable behavior).
Unlike the optimization of the perceived location, which is more complex and based on non-trivial input to output relationships, the optimization of the perceptual threshold charge in a uni-dimensional setting can also be achieved successfully with other automatic algorithms, which are not based on AI methods, such as the bisection method, that relies on the theorem of the intermediate value ( figure S4). Nevertheless, these methods are hardly applicable in clinical settings for a variety of reasons. First, they rely on a blind research strategy that requires the use of unnecessarily high charges, which can lead to significant patients' discomfort and lower implantable battery life. Second, they are not suited for multi-dimensional frameworks where a combination of many parameters needs to be optimized, leading to most likely nonlinear and non-convex optimization problems, which we believe will be the future framework for which we propose here the building blocks. Third, the subject's response to stimulation is intrinsically variable and noisy and affected by, among others, habituation, and psychosomatic effects. While the GPBO algorithm can deal with this issue by properly tuning between exploration and exploitation around the collected observations, non-AI based algorithms fail to catch this variability.
A further improvement would be represented by a more detailed exploration of the hyperparameters of the Gaussian Process such as the length-scale, the signal variance and the noise and the kernel type. Although our choices are based on the experimental data, a further hyperparameter optimization would potentially lead to higher performances. Key parameters of the algorithm, such as the acquisition function, could be improved to other functions that have been proven successful in the literature, as the upper confidence bound [22]. Finally, an important aspect would be to better integrate the information on the past evolution of the perceptual threshold of an AS, in order to tune the trade-off between exploration and exploitation. Indeed, if an AS is recognized as stable, more attention should be paid to the integration of temporal information (since the difference between the two values is expected to be small). On the other hand, if an AS is unstable, then a larger exploration is preferred, given that larger changes are expected. A classification module has been implemented with the purpose of understanding how and to what extent it is possible to classify between the three classes, but more studies into this aspect are needed. While some physiological mechanisms behind the change in Perceptual Threshold have been proposed [25,27,63], still a lot remains to be done regarding the change in neural structures being targeted, which also affects significantly and critically the efficiency and usability of the neuromodulation technology being used.
Somatotopy is a known property of the nervous system, which arises due to the rigorous organization of the neural structures where closer neurons have similar properties in terms of their functionality [75]. This is maintained for the visual system [76], the motor [77] and sensory cortices [78] up to the periphery [75]. Because of variable surgical implantation methods, on the interface being implanted and its shape (if compliant with the neural structure or not), and other biological mechanisms, it is not always possible to define a predictable relationship between the vicinity of the AS and the outcome of the stimulation. Indeed, in our specific case, we could prove not only that close AS were not eliciting spatially adjacent locations, but rather no correlation was present. Also, the same AS could undergo changes throughout time post-implantation, even changing the functional group of the location it elicited. We could also confirm that there is no correlation between physically closer AS and their level of stability, while a higher correlation exists between the level of stability of AS eliciting closer locations (figure S1). This might be explained by the fact that if a FBR is ongoing around one AS, this phenomenon could be propagated along that fiber, which could be interfaced with another AS, not necessarily close to the first one. In this way, instability information is not propagated to physically closer AS but rather functionally closer AS. Importantly, for most of the AS that lose their functionality or experience significant degradation, the mechanisms underlying such changes are not to be addressed to an electrode failure but rather to impedance increases and fibrotic encapsulation related to FBRs.
It is crucial to highlight how when an electrode experiences a mechanical failure, this information needs to be propagated to other sites in the same implant, where instead the physical distance between AS matters. These complex neurophysiological mechanisms underlying instability highlight the importance of finding a way to re-calibrate such systems, so that the functionality of the device is always guaranteed regardless the time evolution and changes undergoing at the level of the implant. Our method, a GPBO with consequent integration of temporal information, proved useful for tackling such issue. Indeed, with a lower number of stimulations than other methods such as a reset GPBO or a random search, the algorithm was able to find the ASs needed to elicit the four functional group in a very low number of iterations.
In order to obtain a finer spatial resolution, in the future it could be possible to not only consider the yaxis (posterior to anterior) of the foot, but also the x-axis (medial to lateral) and the z-axis (proximal to distal). Although being able to properly discriminate the heel, the toe or the mid part of the foot already represents a significant restoration of tactile acuity for lower-limb amputees, a more natural sensation could be delivered by encoding for a higher complexity.
Our work presents several limitations. First, the retrospective application of the algorithm in an offline scenario led to a lack of knowledge for the GBPO. Indeed, as we did not explore all the different ranges of injected charge, we could not hypothesize which would be the perceived intensity of sensation at every injected charge, but rather could only give a quantitative measure of how much that charge was different from the reported one. This measure is, on one hand, more precise than the reported intensity (which is much more discretized into a fewer number of levels of intensity perceived) but on the other hand is not provided during a normal re-calibration session (where, clearly, the target value is unknown), meaning that an online validation of such model is needed. It is important to point out that in a realistic scenario we do not have access to the target perceptual threshold (which is the value we are looking for). In this case we would use a subject reported intensity for each new stimulation (new charge value) at the target. We would merge it with the hypothesis of a positive strong correlation between the charge injected and the elicited intensity [56,70,79]. Moreover, the validation of this algorithm in an offline setting results in a less comprehensive assessment of the algorithm's robustness. In an online scenario, it would be possible to further validate the suitability of the neurostimulation parameters' choices, by for instance delivering 2-3 stimulations with fixed electrical parameters to assess the sensation's consistency. Nevertheless, since we are dealing with a stochastic process, we have dealt with the issue of robustness by running the simulation multiple times (N = 30) to assess convergence towards a similar number of steps needed by the algorithm. The importance of validating this approach in an online setting also arises from a more realistic comparison with the neuromodulation expert, which was instead given only a formal derivation in the present study. Indeed, the formalization was based on a retrospective analysis of the stimulations given by neuromodulation experts that participated in the data collection during the clinical trial (which were hand-annotated and examined here). In an online setting, the performance of the expert could instead be established and immediately compared with the algorithm.
Secondly, the safety control layer allows on one hand to have a more controlled and safe way of exploring the range of parameters, but on the other hand is not applicable for those (rare) cases where the perceptual threshold is outside the explorable range of the algorithm. While the expert would manage it by enlarging the explorable range, the algorithm would fail. A trade-off between the two led us to the conclusion that safe research is more important than the potential outliers, given the paucity of those. Third, we lack detailed enough information regarding subject 1 and subject 3. These two subjects, which had a more unstable implant due to catastrophic electrode failures (some electrode stopped working due to mechanical breakage), were also characterized with too few mapping sessions to efficiently validate the algorithm with their data. A further step would be to verify how the model behaves in more unstable scenarios.
A potential approach to tackle this instability is to re-order the AS each time one is used, and its output is different. It is important to consider that each AS in this circumstance is intrinsically independent, given their physical nature and the lack of any relationship between vicinity and outcome. This would mean to consider them as categorical variables, and proceed with a one-hot encoding in the algorithm (at a cost of a very high number of dimensions), losing any possibility of finding a relationship between similar AS. Another aspect to further explore is the extent and importance of temporal information integration. When the same AS (input) has two different outputs (e.g. sensation migrates), this is encoded as noise, rendering the process more and more uncertain. This can be considered as a help for the algorithm which is, then able to discriminate and differentiate between very noisy inputs and less noisy ones. On the other hand, if all the AS have noise, the posterior distribution, when updated, resembles a uniform one and leads to a random-search type of research. Being able to correctly identify what fibers an AS is targeting and how this might change in time might require more complex models, such as hybrid computational modeling [73,80,81], and still requires significant exploration. Further aspects to take into consideration are the possibility of applying multioutput priors and multi-objective optimization [82,83]; where different features are considered jointly instead of independently. In our specific case, this would mean to consider the changes in sensation quality, location and intensity jointly. This would potentially represent a more comprehensive framework of optimization, but it is important to highlight that these three features of sensation are not correlated (i.e. there is no trivial relationship that explains their inter-dependence), and their respective stability in time is not following a similar pattern; thus leading to a harder implementation of such methodologies. The lack of inter-dependence of the two problems also led to the choice of implementing two different algorithmic designs. Firstly, the two objective functions have different properties (with the charge presenting a monotonic relationship to the felt intensity, while the AS presenting a more complex relationship to the location perceived). Moreover, based on our data, there was no trivial relationship between the stability of the location and the stability of the perceptual threshold. Therefore, the two models proposed are different, which enhances the importance of making educated guesses when designing such automated methods. Nevertheless, the need to develop two separate algorithmic instances also enhances the importance of developing a more complex and comprehensive framework for simultaneous optimization of different input parameters. Still, our proposed algorithm is the first, up to our knowledge, to deal with recalibration in sensory feedback restorations; and provides a benchmark for future applications and further developments.
Another barrier encountered by somatosensory neuroprostheses is the lack of a measurable neurophysiological signals: indeed, we based our model on self-reported sensation's properties as verbally indicated by the subject experiencing them. This clearly increases the noise around each evaluation and the unreliability of the output: for instance, it is empirically expected that subjects will respond in a less attentive manner towards the end of the mapping session (which can be very exhausting). This leads to less precise responses, which will affect the performance of any model. In contrast with other neuroprosthetic applications of GPBO (such as for instance the deployment of neuro-stimulation protocols to evoke limb movements), where quantitative and objective measures such as EMG are used as a clear output of the neuromodulation protocol, the absence of purely objective measure in our case affects the automatization of the re-calibration process. Potential solutions could be the recording of electroneurography (ENG) [84] as an objective measure of nerve responses to intraneural stimulation. ENG could be recorded by the implanted TIME electrodes [85], given their usability as both recording and stimulating devices. What is important to highlight though, is that this would require a higher level of processing to properly associate the neural response to the sensation; due to higher complexity underlying perceptive processes with respect to the interpretation of muscle contractions by EMG, for instance.
It is clear that when it comes to onedimensional optimization, data-driven and modelfree approaches could also be exploited. Nevertheless, these approaches (mostly deterministic) fail to catch the probabilistic nature of neurophysiological data, which is instead characterized by an intrinsic noise and output variability, even in the case of equal inputs. For this reason, we believe more complex algorithmic procedures such as the ones proposed here are more tailored for recalibration issues in neuromodulation scenarios. Future challenges to be addressed, specifically for sensory feedback restoration, are the need to build multi-dimensional optimization scenarios where input co-existence and inter-independence is tackled by a proper objective function. While more complex models, could be used (such as reinforcement or deep-learning based models), these require an extensive amount of data that is not yet realistic in human applications. A further challenge for the proposed design would be to modulate exploration and exploitation weights based on the computed instability of the AS of focus, enabling the algorithm to deal with different stability frameworks.
In conclusion, this work presents the successful use of GPBO algorithms for the adaptation of the neurostimulation threshold changes and location changes in sensory feedback restoration, tackling the issue of implant stability and time variability of the input-output mapping in neuroprosthetic applications. The neuroscientific community has extensively explored the mapping between electrical parameters and outcomes in order to speed up first calibration procedures. However, we believe that this study plays an important role integrating this knowledge with the need to re-calibrate neuromodulation devices in the chronic setting. Importantly, it lays the foundations as a benchmark and validation of the GPBO method from which more complexities such as multi-dimensional optimization can be improved. This will allow future applications to develop automatic, fast, and reliable calibration schemes ready to use in clinical settings, helping during the expertmediated recalibration. We are here providing a proof of concept of how this class of algorithms can tackle recalibration issues by integration of previously collected information, together with datadriven model's hyperparameters selection. Automatic re-calibration is the key to future wide-spread use of neurotechnologies since continuous control is crucial for effective and functional neuromodulation.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// github.com/giovannaiello/recalibration_gpbo. Data that supports the findings and the software routines developed for the analysis are available from the corresponding author (S R)