Multiparametric non-linear TENS modulation to integrate intuitive sensory feedback

Objective. Transcutaneous electrical nerve stimulation (TENS) has been recently introduced in neurorehabilitation and neuroprosthetics as a promising, non-invasive sensory feedback restoration alternative to implantable neurostimulation. Yet, the adopted stimulation paradigms are typically based on single-parameter modulations (e.g. pulse amplitude (PA), pulse-width (PW) or pulse frequency (PF)). They elicit artificial sensations characterized by a low intensity resolution (e.g. few perceived levels), low naturalness and intuitiveness, hindering the acceptance of this technology. To address these issues, we designed novel multiparametric stimulation paradigms, featuring the simultaneous modulation of multiple parameters, and implemented them in real-time tests of performance when exploited as artificial sensory inputs. Approach. We initially investigated the contribution of PW and PF variations to the perceived sensation magnitude through discrimination tests. Then, we designed three multiparametric stimulation paradigms comparing them with a standard PW linear modulation in terms of evoked sensation naturalness and intensity. The most performant paradigms were then implemented in real-time in a Virtual Reality—TENS platform to assess their ability to provide intuitive somatosensory feedback in a functional task. Main results. Our study highlighted a strong negative correlation between perceived naturalness and intensity: less intense sensations are usually deemed as more similar to natural touch. In addition, we observed that PF and PW changes have a different weight on the perceived sensation intensity. As a result, we adapted the activation charge rate (ACR) equation, proposed for implantable neurostimulation to predict the perceived intensity while co-modulating the PF and charge per pulse, to TENS (ACRT). ACRT allowed to design different multiparametric TENS paradigms with the same absolute perceived intensity. Although not reported as more natural, the multiparametric paradigm, based on sinusoidal PF modulation, resulted being more intuitive and subconsciously integrated than the standard linear one. This allowed subjects to achieve a faster and more accurate functional performance. Significance. Our findings suggest that TENS-based, multiparametric neurostimulation, despite not consciously perceived naturally, can provide integrated and more intuitive somatosensory information, as functionally proved. This could be exploited to design novel encoding strategies able to improve the performance of non-invasive sensory feedback technologies.


Introduction
Somatosensory loss is an invalidating pathological condition arising from multiple causes (e.g. peripheral neuropathies, stroke, limb amputation), which affects a huge portion of population (e.g. 93 ′ 000 yearly amputations only in the European countries) [1][2][3]. Sensory loss impairs the interaction with the external environment, impacting the simplest daily life activities, such as grasping a glass of water or walking [4][5][6][7][8][9]. It has been proven that the lack of autonomy and phantom limb pain resulting from this condition [3,4,[10][11][12][13][14][15][16], can be significantly reduced through the adoption of neuroprosthetic devices providing somatosensory feedback [14,17,18]. Indeed, tactile feedback restoration has proven to enhance manual dexterity [9,11,[19][20][21][22], improve mobility, speed and postural stability while walking [4,12,[23][24][25], and reduce cognitive load, positively impacting on the device acceptance and embodiment [2,5,[26][27][28][29]. Importantly, the stimulations used to provide tactile feedback must reproduce as faithfully as possible the natural neural codes [30][31][32][33][34][35][36][37]. Indeed, the optimal artificial sensory feedback should be somatotopic and homologous meaning that the elicited sensation should match the stimulus quality and the stimulus location (e.g. at the heel-strike, we should perceive a pressure sensation directly under the heel) [2,20,38,39]. Somatotopic and homologous, touch-like sensations have been successfully evoked using invasive electrical stimulation of peripheral nerves [11,19,[40][41][42][43][44]. In order to achieve an homologous and informative artificial feedback on neuroprosthetic devices, the development of sensory encoding algorithms is required [45][46][47]. In standard linear encodings, a single stimulation parameter (e.g. pulse width (PW), pulse frequency (PF) or pulse amplitude (PA)) is modulated proportionally and linearly to the sensor output [40]. In particular, it has been shown that the linear modulation of the injected charge allows to encode a higher number of intensity levels, while PF modulation would seem to play a relevant role in defining the quality of the sensation [40,[48][49][50]. Graczyk et al [51] introduced a quantitynamely activation charge rate (ACR)-able to group all the stimulation parameters and directly related to the perceived sensation magnitude. This single stimulation quantity could predict the perceived intensity combining the PF and charge per pulse. The ACR indicates the charge per second above threshold that is provided by the stimulation [51]. Promising results in terms of evoked sensation quality were achieved through the implementation of sophisticated multiparameters stimulation paradigms. Tan et al showed that the sinusoidal modulation of PW was able to elicit a sensation similar to natural touch, probably by favouring a more asynchronous fiber activation [11]. Similarly, new biomimetic approaches have been introduced, based on the hypothesis that mimicking the natural encoding strategies [47,52,53] could bring to significant gains in sensation naturalness (i.e. defined as an electrically-evoked sensation closer to a touch sensation felt with the same skin location of the intact limb). Valle et al and George et al, showed improved manual dexterity and prosthesis embodiment in upper limb amputees [19,41,54], using biomimetic PF modulation integrated into a multiparametric (based on co-modulation of both charge and PF) stimulation scheme. Pivotal advantage of implants is that they elicit focused somatotopic sensations, as reported by patients, and observed in their brain signals [55]. Yet, these implantable devices suffer from a range of issues, starting from necessary surgery, inevitable biological changes at the interface [56], and their long-term stability [2].
Non-invasive approaches (e.g. sensory remapping through electro-tactile stimulation) have proven to be less intuitive than neural interfaces, due to their inability to evoke somatotopic sensations [57]. To address this, transcutaneous electrical nerve stimulation (TENS) has been introduced as intermediate solution between remapping and neural implants for non-invasive, somatotopic tactile feedback restoration. Indeed, it has been successfully tested on upper and lower limb amputees as low-cost, low-risk, temporally stable alternative to invasive interfaces [6,10,[57][58][59][60][61]. However, mostly single-parameter stimulations have been implemented evoking low natural percepts and low informativeness. For instance, object stiffness and size recognition were performed by modulating either only the PA, as in Vargas et al [59,61] or only PW as in D'Anna et al and Shin et al [6,62]. Few more complex, neuromorphic modulations have been implemented in TENS, without leading to significant improvements in sensation quality [47]. Moreover, the elicited sensations are characterized by a low intensity sensitivity (e.g. few encoded intensity levels) [6,62,63]. Neural stimulation codes, resembling more closely the cumulative asynchronous fibers activation during a sensory process (thus more natural neural activation), demonstrated higher naturalness and functionality in implants, therefore implementing these for TENS holds promise to increase the user experience and resulting functionality. While implementing these novel paradigms, the charge injected to subjects is very variable, and if a paradigm is perceived as more intense, this could eventually bias how the subject judges its naturalness resulting in misleading results. Therefore, it is important to inquiry the intensity influence. This leads to three open questions in the field, which could affect the usage and the acceptance of devices based on TENS: The aim of the study was to evaluate the possible benefits of non-linear multiparameter modulation (PW and PF), inspired from invasive approaches [11,19,41,54], and recent modelling [64], on the naturalness, intuitiveness and functional usability, of the TENS-elicited sensations. To this aim, it was firstly necessary to ensure reliable and robust comparison across paradigms, by avoiding design biases such as differences in the perceived intensity. We therefore evaluated the contribution of each modulated stimulation parameter to the perceived sensation intensity and formulated the TENS ACR T equation to design paradigms with the same absolute perceived intensity while co-modulating multiple parameters. Although the paradigms were absolutely perceived as equally intense, we unveiled a strong, intrinsic negative correlation between perceived naturalness and intensity when evaluating the paradigms in a forced comparison. Furthermore, despite not explicitly perceived as more natural, the sinusoidal PF paradigm (Multi-param1), proved to encode tactile information more intuitively and subconsciously integrated. Then, it was implemented in a VR (virtual reality)-TENS platform to assess its objective performance and intuitive integration in providing real-time somatosensory information in a functional task.

Study design
The study design is summarized in figure S1. Novel TENS paradigms were applied targeting the tibial nerve, and in particular the plantar branch that innervates the most distal part of the foot sole ( figure 1(A)). While implementing these novel paradigms, the charge injected to subjects is variable, and if a paradigm is perceived as more intense, this could eventually bias how the subject judges its naturalness resulting in misleading results. Therefore, the purpose was to design equally intense paradigms and to do so, we firstly investigated the relationship between the ACR, proposed by Graczyk et al [51], and perceived sensation intensity on 56 healthy subjects (Intensity Discrimination experiment) (figure 2(A)). Then, we designed three non-linear multiparameter and a standard single parameter stimulation paradigm using an adapted ACR function to ensure same intensity while co-modulating different parameters. These novel TENS paradigms were evaluated in terms of sensation naturalness (e.g. the sensation elicited by the pressure exerted by a small, smooth object on the limb), intensity and intuitiveness (e.g. ability to encode specific tactile information, as a skin indentation [41]) in 14 subjects (Naturalness and Intuitiveness rating task). The paradigms were compared by assigning absolute naturalness ratings to each individual paradigm and by performing a direct forced comparison between encoding strategies in a two-alternative forced choice (2AFC) task ( figure 3(A)). The intuitiveness was tested by mapping the sensation to the representation of a specific tactile stimuli. Finally, we implemented a unique electro-visual platform enabling to disentangle the contributions of different sensory cues and designed experimental tasks that require the accurate use of sensory feedback (Closed-loop VR-TENS functional task, figure 1(C)). About 12 subjects were asked to carry out a functional task under different visual (visual, blurred visual and no-visual feedback) and tactile (no-tactile, tactile complex paradigm, tactile standard paradigm) feedback conditions (figure 4(A)).

Sensation calibration procedure
A sensation characterization procedure was performed at the beginning of each experiment to obtain the subject-specific charge values [65]. Electrical nerve stimulation (TENS) was delivered using Rehamove 3 (Hasomed GmBh, Germany), which is a CE-approved portable and programmable four channels surface stimulator with transcutaneous adhesive electrodes (Circular Electrode Pads 25 mm, TensCare, UK). Stimulation trains of biphasic, balanced, cathodic-first square pulses have been used in this study ( figure 1(A)). The aim of the sensation characterization was to elicit a somatotopic sensation in the distal plantar area by targeting the tibial nerve. The exact electrode placement is subject-specific, and it was determined with an empirical procedure, by arranging the electrode positions around the medial malleolus. The followed protocol was to place the electrodes behind and below the medial malleolus as shown in figure 1(A), since in this area the tibial nerve is more superficial and therefore closer to the electrodes. Stimulation was delivered to the subjects and electrodes were adjusted according to their feedback on somatotopic and in loco sensations. The electrode arrangement was selected as the one eliciting a non-painful somatotopic sensation while eliciting an in-loco (under the electrodes) sensation perceived with a lower intensity than the somatotopic one. Once the electrodes were optimally placed, they were secured with medical tape for the whole duration of the experiment, so to be consistent across different paradigms. The somatotopic charge threshold also presents a high inter-subject variability. Two different levels, corresponding to the minimum charge able to elicit a somatotopic sensation (somatotopic perceptual threshold) and the maximum charge corresponding to a strong sensation (higher threshold) were The ACR T profile was designed to encode a 2 s skin indentation having the same intensity for all the paradigms, while PW and PF are modulated in a different way for each paradigm. Standard is non-multiparametric and characterized by constant PF, Multiparametric 1 has a sinusoidal PF modulation, Multiparametric 2 has a biomimetic PF modulation and Multiparametric 3 has a sinusoidal PW modulation. (C) Hardware and setup of the functional task carried out in a VR-TENS environment. The foot motion was tracked and allowed the subject to perform actions in VR, according to which, the TENS stimulation was modulated. To have a stimulation intensity increasing proportionally to the rising of the water, each water level is linked to a different ACR T (for Multiparametric 1) or PW (for Standard). The ACR T and PW values are equally spaced and, in the Multiparametric 1 paradigm both PF and PW are sinusoidally modulated over time.
obtained. These two values were associated, respectively, to the level 2 and level 8 in an intensity scale ranging from 0 (no sensation) to 10 (painfully strong sensation). To do so, PA and PW ramps at a constant PF of 80 Hz (as in [25]) were performed. Firstly, three ramps with fixed PW at 200 µs and increasing PA (starting at 3 mA and steps of 1 mA) were performed, saving the PA value corresponding to the subject's somatotopic sensation perceptual threshold. Then, the average value between the ones obtained was computed and used as constant PA in the PW ramps. Three PW ramps were delivered starting from a value of 100 µs and increasing with steps of 20 µs: the subjects were asked to report a somatotopic sensation intensity of 2 out of 10 corresponding to a certain PW value and an intensity of 8 out of 10 where 10 was considered the maximum threshold. The results of the three ramps were averaged and the minimum and maximum charge values were computed for the individual subject, as the product between the threshold PA and the PW values at level 2 and level 8 respectively.

Intensity discrimination experiment
To determine whether the ACR, as defined by Graczyk et al for implantable neural interfaces [51], can be used as a reliable reference quantity for sensation intensity encoding also using TENS, an intensity discrimination experiment was carried out on 56 healthy subjects. The ACR Weber Fraction (WF) was computed to analyse the effect of ACR variations on perceived sensation intensity (figure 2(A)). The ACR is defined as [51]: Then, eight equally-spaced ACR values, spanning the range between 25% and 175% of ACR ref were computed. To perform the intensity discrimination, couples of stimuli composed by ACR ref and one of the other ACR values (or ACR ref itself), were presented to the subjects with the following timing: 2 s first stimulus-2 s pause-2 s second stimulus. Following a double-blinded 2AFC paradigm, the subjects were asked to indicate which of the two stimuli was perceived as stronger. Each stimuli couple was sent eight times and both the couples order and the order of the stimuli within one couple were randomized. To examine the impact of PF and PW on sensation intensity the experiment was repeated three times: the stimuli ACR values were always the same, but the ACR variations were once obtained by changing the stimulus PW only, once the PF only, and once by varying both PW and PF of the same amount.
Once the subjects' answers were recorded, psychometric curves were fit to the experimental data. The psychometric curve represents the probability of judging the test stimulus as stronger than the reference. Just noticeable difference (JND) values were obtained from the psychometric curves as the difference between 75% and 50% probability. To better compare discriminability across different stimulation conditions, we computed the WF, which is the JND divided by the reference stimulus, i.e. the ACR value corresponding to a 50% probability [25,51]. Smaller WFs indicate a higher sensitivity, i.e. a smaller difference in input (PW or PF) is needed to perceive a difference in intensity.

Naturalness and intuitiveness rating task
We defined four stimulation paradigms characterized by different parameters modulations to compare them in terms of perceived sensation naturalness. The paradigms were designed to have the same ACR profile, and thus the same perceived intensity. The ACR profile was defined to encode a skin indentation lasting 2 s: 0.5 s pressure increase (loading), 1 s steady pressure (sustained pressure), 0.5 s pressure decrease (un-loading) (figure 3(A)), reminding the trapezoid press-and-hold pattern proposed by Shin et al [62].
After the analysis of the intensity discrimination results the ACR T equation was modified as: and this new formula was adopted in all the experiments of this study. The ACR T value in the flat region of the ACR T profile was the one computed in correspondence of the PW value PW mean so that the perceived intensity would be rated around 5 out of 10.

Sensory encoding strategies design
The paradigms were designed to modulate different stimulation parameters (figure 1(B)): (1) Standard: constant PF at 80 Hz (constant PF) and the proportional modulation of PW, used as non-multiparametric reference stimulation. (2) Multiparametric 1: sinusoidally-oscillating PF (sinusoidal PF). This type of modulation was chosen since it may induce an higher asynchronicity (thus more natural neural activation) in the nervous afferents response [11]. Furthermore, according to previous studies, PF plays an important role in determining the quality of sensation [48,50].
(3) Multiparametric 2: This biomimetic approach was developed based on the hypothesis that reproducing natural fiber activation patterns can increase the sensation naturalness. Indeed, it modulated PF according to the biomimetic firing rate profile (biomimetic PF), extracted using a realistic in-silico model of foot cutaneous afferents. The model, called FootSim [64] emulates the response of all the cutaneous afferents innervating the foot sole to any type of mechanical skin deformation. It allows to populate the virtual foot sole with some or all the tactile afferents classes, to choose the features (location, duration, contact force profile) of the mechanical stimulus applied on the foot sole and then to export the related temporal and spatial afferent fibers activation. The aggregate population firing rate response was then used to inform the design of the multiparametric paradigm. It was generated by populating the virtual foot sole with SA1, SA2, FA1 and FA2 mechanoreceptors and by applying a virtual stimulus (3 mm deep indentation of a 150 mm diameter round pin in correspondence of the medial metatarsal) (figure S2(A)). The obtained cumulative firing rate profile (figure S2(B)) was used to generate the PF profile. However, since the biomimetic approach was defined for invasive stimulations (delivered directly through electrode implants in close contact to the nervous system), some modifications were made in order to consider the tissue between the electrode and the target nerves. Firstly, it was imposed that the minimum frequency would not be lower than 40 Hz (so that the subject would not perceive the single electrical pulses and would feel the stimulation as being more continuous, basing on what observed by Ng et al [30]) and that the higher frequency would not overcome 400 Hz (limit defined by the stimulator manufacturer).
In all the paradigms, the PA was kept constant at the somatotopic perceptual threshold value found in the characterization phase.
For the Standard, multiparametric 1 and multiparametric 2 paradigms, after generating the ideal indentation ACR T profile and the PF profile, the PW profile to be used for the stimulation was computed by inverting the ACR T equation: Moreover, a constraint, according to which the computed PW values would not be higher than the subject's PW max was added, otherwise the stimulation could have been over the maximum threshold. Because of this, the resulting ACR T profile (delivered stimulation profile), computed using the PW profile modified with these constraints, would slightly differ from the ideal one. To ensure that this would not result in a significant difference between the paradigms in terms of intensity, it was checked that the difference between the real and ideal ACR T profile was not higher than the WF (found in the intensity discrimination task). All the paradigms satisfied this condition. The same workflow was followed for the implementation of multiparametric 3, starting from the ideal ACR T profile and the PW profile, computing the PF profile and checking that the variations between real and ideal ACR T profile would be lower than the WF.

Sensation characterization evoked by different stimulation paradigms
The experiment was carried out on 14 healthy subjects. After the characterization phase the subjects were firstly asked to report the sensations evoked by the four paradigms. The subjects were asked to fill a form indicating the perceived foot sole area (projective field), the intensity of both the in-loco (under the electrodes) and the somatotopic sensation by rating it between 0 (no sensation) and 10 (strong sensation. Then, they were asked to rate the naturalness of the somatotopic sensation between 0 (unnatural, e.g. paraesthetic sensations) and 10 (fully natural touch sensation), similar to the definition already proposed in [11,41,51,[66][67][68][69]. The naturalness was defined as 'how close the artificial sensation is to a natural, non-painful sensation elicited by the pressure exerted by a smooth object on the limb (e.g. a small pedal pressed with the foot or the pressure of that area of the foot on the floor)' . Since TENS can elicit different types of sensations, each paradigm was also qualitatively described by choosing sensation descriptors (light touch, pressure, vibration, tingling, slow pulsation, stinging, electricity, pain) [65]. Among these descriptors, it is desirable to avoid uncomfortable and unpleasant sensations such as stinging, electricity and pain. TENS can elicit other types of pleasant and informative sensations like vibration, tingling, slow pulsation, which are however not homologous to 'pressure exerted by a smooth object on the limb. ′ Finally, the subjects selected the most fitting dynamical representation of the sensation among some short videos representing different interactions between the foot sole and an external object (sustained constant pressure, indentation ramp-and-hold, vibration, skin stretch, other) to evaluate whether the encoded sensory information was correctly and intuitively conveyed by the artificial tactile stimulation (intuitiveness) (figure 3(A)).

2AFC comparison
Then, each paradigm was directly compared to the others in a 2AFC task. All the possible paradigm couples were formed and delivered to the subjects in a double-blinded fashion. Each two-paradigms couple was repeated eight times and the order of the couples and the paradigms order within one couple was randomized. For each couple the subjects were asked to indicate which paradigm had felt as more natural and which had felt as more intense. For each couple one point was added or subtracted whether the paradigm was judged as more or less intense/natural, respectively (figure 3(A)).

Closed-loop VR-TENS functional task
Since the somatosensory feedback provided by the stimulation should be exploited by end users to properly modulate the exerted force in functional daily life activities, we purposely-designed a customized VR-TENS environment featuring three visual feedback modalities and tactile feedback. The aim was to analyse the intuitive, subconscious and objective contribution of different encoding strategies. To this aim, it was necessary a meticulously designed functional paradigm, able to finely disentangle the influence of different sensory cues. We developed and tested a novel multisensory platform, thanks to the combined use of VR and TENS, allowing to disentangle different sensory cues (e.g. vision and touch). This enabled a high spatial and temporal precision of the visual and TENS artificial touch cues. At the same time, it created a realistic task mimicking real-life activities, such as a subject pressing a brake pedal, in highly controlled experimental conditions. The VR hardware used was an HTC Vive system (HTC, China). The coding platform was Unity 3D which allows VR programming in C# language.
While performing the task the subject was wearing the VR headset and sitting on a chair with their foot above the floor level and supported by a box ( figure 1(C)). The movements of the subject's right foot were tracked by a motion VIVE tracker placed on the foot and the interaction with the game happened by pressing a button on the controller that the subject was holding in their hand. Somatosensory feedback was provided by TENS eliciting a somatotopic plantar distal sensation. The stimulation intensity was modulated according to the subject's actions. The communication between the VR and the stimulator was mediated by a dynamic link library which allows to write on a Pipe Server (Boost C++ Libraries). The sinusoidal co-modulation schemes of PW and PF were pre-programmed in the stimulator firmware. At the beginning of a trial, as soon as the subject pressed the virtual button, a message was sent to the PipeServer (and therefore the stimulator), resulting in the initiation of the stimulation at the desired amplitude (the lowest). Then, the stimulation would continue following the multiparameter predefined scheme until the VR would send an update. When a new VR event (e.g. a new water level) was generated, a message string was sent to the stimulator through the PipeServer allowing to select the required amplitude for the preprogrammed sinusoidal modulation ( figure 1(C)). The message contained the instructions to change the offset values around which PF and PW would oscillate. Updates from VR to stimulation could occur at a frequency higher than 60 Hz, allowing to adapt the stimulation to the subject's actions in real-time.
The VR scenario consisted of an open space where the subject was sitting on a wooden bridge with his right foot close to a virtual button. The subject was facing a water column with four highlighted levels and the aim of each trial was that of pushing the virtual button to make the water in the water column rise to match a target level indicated at the beginning of each trial. IR cameras tracked the subject's foot motion to determine if and how much they were pressing the button. The water level and the intensity of the stimulation increased proportionally to the pressure exerted on the button: for example, if the button was not pressed at all, the stimulation was absent and the water was at the base of the column, if the button was completely pressed the water reached the top of the column and the stimulation was delivered at maximum intensity ( figure 4(B)). The end of each trial was determined by the subject pressing a button on the VR controller when they thought they had reached the target water level.

Real-time TENS intensity encoding
A multiparametric and standard stimulation paradigms were implemented in the VR-TENS platform to evaluate the differences in terms of information richness and intuitiveness.
Both paradigms were designed to encode four distinct intensity levels corresponding to the four water levels in the column ( figure 1(C)).
For the multiparametric paradigm each water level corresponded to a different ACR T value. Those were obtained by computing the subject's ACR T1 (correspondent to PW min of the calibration phase) and ACR T4 (correspondent to PW max of the calibration phase) and by dividing the interval between these 2 values into 3 equally spaced parts. Since the multiparametric paradigm featuring a sinusoidal PF modulation (Multiparametric 1) had achieved the best intuitiveness score in the naturalness rating experiment, that was the one implemented in the VR-TENS environment. At each level the PA was kept constant at the perceptual threshold while both PW and PF were sinusoidally modulated. Thus, the PF was modulated according to the formula: The chosen offsets were 60, 80, 100 and 140 Hz and they were chosen so that the difference between the log 2 (offset) of one level and the log 2 (offset) of the next one was always the same. After having defined the PF profiles the PW profiles were computed by inverting the ACR T equation.
For the Standard paradigm, instead, the sensation intensity was modulated by varying the PW between four equally spaced values in the interval between PW min and PW max , while keeping the PF constant at 80 Hz.

Visual feedback conditions
Different types of visual and tactile feedback (no stimulation, stimulation with Standard or Multiparametric 1 paradigm) and different combinations of the two, were provided to the subjects while performing the functional task ( figure 4(A)).
Three different visual feedback conditions were tested: • Visual feedback: water clearly visible.
• No visual feedback: invisible water.
• Blurred visual feedback: water was visible, but the screen was blurred (as in [70]).

Functional task protocol
The experimental protocol was composed by the following steps: (1) Sensation calibration procedure (as described in 2.2). (2) Functional task in the VR-TENS environment.
In the first round of experiments 3 blocks of 160 trials each were performed by the subjects. Each block was characterized by (I) no stimulation (II) Multiparametric 1 stimulation, (III) Standard stimulation.
The blocks were presented in random order and in each of them 50% of the trials was characterized by clearly visible water and the others by invisible water.
In the second round of experiments, the vision was always blurred and 3 blocks (no stimulation, Multiparametric 1 and Standard) of 80 trials each were performed. In both rounds of experiments, between one block and another the subjects were asked to fill-in a questionnaire regarding the sense of embodiment and task intuitiveness with the different feedback conditions. The questionnaires were adapted from [71].

Evaluation metrics
The performances were evaluated with objective and subjective measures.
Error (Err): absolute difference between the target and the selected water level ( figure 4(A)). It indicates how precise the subject was in reaching the target level. It was computed as a percentage of the total water column height.
Execution time (ET): time needed to complete one trial, defined as the time passing between the moment in which the subject starts to push the virtual button and the moment in which they press the button on the controller ( figure 4(A)).
Consistency across force trajectories. In order to be able to examine how similar the water trajectories of different trials were, the standard deviation between the water trajectories was computed. The consistency was obtained as: 1 − avg (std of trajectories) .
Non-Smoothness of trajectories: the jerk was computed as the standard deviation of the velocity vector of the trajectories, since a high velocity indicates high variations in the slope of the trajectory and, thus, more oscillations.
Embodiment, intuitiveness and self-confidence scores: computed from the questionnaires. The scores were adapted from the work of Gonzalez-Franco et al [71].
Level discrimination accuracy: subjects' accuracy in level discrimination was computed by considering the percentage of successful trials for each target level. For each final water level, the 'Actual level' was obtained as the closest target level (the four levels coloured in dark blue on the water column). A trial was considered as successful if the actual level matched the level which had been set as target for that trial. The percentages of successful trials for each target level were represented in confusion matrices and the overall accuracy was computed.
Visuo-tactile integration: sensory integration was assessed for both paradigms by comparing the performance measures obtained when having only blurred visual feedback (no stimulation) and only tactile feedback (water completely invisible) with the ones yielded when both tactile and blurred visual feedback were provided.

Participants
All the experiments were performed on healthy subjects of both sexes, between 20 and 30 years old. Some of them were already accustomed to peripheral nerve transcutaneous electrical stimulation. The experiments were approved by the ETH Zurich's ethics commission (EK 2019-N-97) and all the subjects signed an informed consent form. The experiments were performed in accordance with the proposal approved by the ETH Zurich's ethics commission and in accordance with the declaration of Helsinki.

Statistical analysis
Data processing and statistical analysis were performed using Python 3.7 with the following libraries: Pandas, NumPy, Matplotlib, SciPy. Resulting bar plots show mean values and the standard error of the mean. Asterisks on plots indicate the statistical significance level characterized by p < 0.05 ( * ) and p < 0.01 ( * * ). The normality of the distributions has been checked using a Shapiro test and the homogeneity of variance has been checked with a Levene test. If distribution normality and uniform variance were not rejected, one-way repeated measures ANOVA test was used, otherwise a Kruskal-Wallis test was used. Post-hoc correction for multiple groups comparison was performed with Tukey's HSD test.

Results
In this work, we firstly evaluated the contribution of different stimulation parameters (PW and PF) to the ACR and therefore the perceived sensation magnitude. Psychometric curves were fit from the 2AFC task. Then, JND and WF for each modulation strategies were computed. Secondly, four paradigms (three multiparametric and one non-multiparametric) were characterized in terms of naturalness and intensity of elicited sensation, both independently and by direct comparison in a 2AFC task. They were also compared in terms of sensation intuitiveness, evaluating their ability to convey specific tactile information. Then, we implemented the multiparametric stimulation encoding yielding the best intuitiveness results in a VR-TENS environment to test its performance in a functional task. The usefulness and informativeness of the feedback provided by a multiparametric paradigm (Multiparametric 1) was tested against tactile feedback provided by a single-parameter stimulation (Standard) under different visual feedback conditions: visual feedback, no visual feedback, and blurred vision.

PF and PW do not have the same impact on perceived intensity
The psychometric curves are shown in figure 2(B). JND was evaluated as the ACR range between 0.5 and 0.75. The direct comparison between the three different conditions is shown in figure 2(A). The WF when changing only PW (WF = 0.077) was statistically significantly lower than PW-PF co-modulation (WF = 0.200, p < 0.05) and only PF modulation (WF = 0.430, p < 0.01). The modified version of ACR function for non-invasive TENS stimulation is shown in figure 2(C). The ACR T function has a logarithmic relationship with PF. Figure 1(B) shows an example of each different paradigm (standard and multiparametrics), with the modulation of PW, PF and the resulting ACR T function representing a skin indentation. The paradigms were compared in terms of naturalness, perceived intensity and intuitiveness (e.g. type of perceived stimulus). Figure 3 Although it was not explicitly judged as more natural, the intuitiveness scores indicated that the paradigm Multiparametric 1 achieved a better performance than the others in conveying the skin indentation information ( figure 3(C)). It yielded a correct indentation recognition 43% of the times (against a 21% correct recognition with Standard, 7% with Multiparametric 2 and 29% with Multiparametric 3, chance level: 20%). The whole range of sensations that we elicited with TENS for each paradigm is shown in figure 3(D).

Multiparametric paradigms do not increase perceived naturalness
When moving to direct paradigms comparison in the 2AFC task ( figure 3(E)), the score distributions analysis did not show any statistical difference between Standard and the two multiparametric paradigms Multiparametric 1 (p = 0.9434) and Multiparametric 2 (p = 0.2619) in terms of naturalness. Multiparametric 3 had, instead, significantly lower naturalness scores than the other three paradigms (p < 0.01) and significantly higher intensity scores (p < 0.01), despite being calibrated with the same intensity. Considering the intensity scores, Multiparametric 2 also had significantly higher scores than Standard (p < 0.05). Analysing the 2AFC naturalness were compared in terms of naturalness and intensity. The paradigms were designed with the same ACR T profile so that they would have the same intensity and they would encode a skin indentation. The paradigms were firstly characterized one by one by filling in questionnaires and then directly compared following a 2AFC paradigm. (B) Questionnaire results: intensity and naturalness ratings distributions. (C) Intuitiveness choices: the paradigms were designed to encode a skin indentation, so a better intuitiveness performance was shown by the paradigm having the highest indentation occurrence (Multiparametric 1). (D) Occurrence of sensation descriptors associated to each paradigm. (E) Two-AFC paradigm comparison results: naturalness and intensity scores distributions and correlation between perceived intensity and naturalness (p < 0.05 ( * ), p < 0.01 ( * * )). N = 14 subjects. scores with the intensity scores, a strong negative correlation between perceived intensity and naturalness emerged: R 2 =0.994 with p = 0.00266.

A multiparametric, sinusoidal PF modulation provides more intuitive somatosensory feedback
To test whether multiparametric stimulation paradigms can be more intuitive than standard ones, as suggested by the outcome of the study on naturalness, the Multiparametric 1 and the Standard stimulation paradigm were implemented in real-time in a VR-TENS environment to provide somatosensory feedback to subjects performing a functional ecological task.
Whereas the two paradigms did not have a statistically different impact on the performance when visual feedback was presented, providing a multiparametric tactile feedback resulted in a significantly better performance when the water was not visible (figure 4(C)). Indeed, when the water was not visible the Multiparametric 1 feedback yielded a significantly lower (p < 0.01) Error than the Standard (average of 11.51 ± 0.44% for Multiparametric 1 and 13.94 ± 0.59% for Standard) while requiring approximately the same time to complete a trial (average Execution Time of 1.49 ± 0.042 s for Multiparametric 1 and 1.47 ± 0.039 s for Standard). The Multiparametric 1 stimulation also resulted in a higher consistency among trajectories ( figure  S4(A)), with average value of 0.496 ± 0.01 against 0.409 ± 0.008 of Standard, whose value was similar to the one achieved when no stimulation was provided (0.416 ± 0.007). The Standard feedback, however, led to a significantly lower jerk of trajectories (0.332 ± 0.009 against 0.378 ± 0.009 of Multiparametric 1 and 0.403 ± 0.011 in case of no stimulation).
When subjects' vision was blurred, the Multiparametric 1 stimulation yielded a significantly more precise and faster performance than the Standard paradigm, showing an average error of 7.98 ± 0.39% (while it was 9.79 ± 0.47% for Standard and 11.17 ± 0.58% with no stimulation) and an average Execution Time of 0.87 ± 0.031 s (against 1.2 ± 0.045 s for Standard and 0.99 ± 0.038 s without stimulation). Moreover, the Multiparametric 1 paradigm led to a slightly higher overall accuracy in level discrimination (0.786 while it was 0.729 for the Standard), and the correlation between the used stimulation paradigm and the accuracy was confirmed with a Fisher's exact test (p = 0.03065) ( figure 5(A)). From the trajectories analysis, instead, no significant differences between Multiparametric 1 and Standard performance were highlighted in the case of blurred vision ( figure S4(A)).

Multiparametric paradigms allow for visuo-tactile integration
A visuo-tactile integration analysis was performed for both Multiparametric 1 and Standard feedback, to verify whether, in presence of both visual and tactile feedback, the sensory information was integrated by the subject and exploited to achieve a better bimodal performance. If the subject is correctly integrating the sensory information, the performance should improve when both tactile and blurred visual feedback are provided (bimodal condition) with respect to the unimodal conditions (only blurred vision and only touch). The Multiparametric 1 paradigm showed a significantly lower (p < 0.01) Error (8.23 ± 0.45%) (figure 5(B)) and higher smoothness (0.33 ± 0.009%) of trajectories ( figure S5(B)) when combined with visual feedback with respect to both the unimodal conditions, while the Execution Time in the bimodal condition (0.88 ± 0.036 s) was not significantly different than the unimodal visual one (0.99 ± 0.044 s) ( figure 5(B)). The consistency of trajectories did not show statistical differences among the three conditions ( figure S5(B)).
The Standard paradigm was, instead, not integrated by the subjects. This is firstly shown by a significantly higher Execution Time in bimodal condition (1.23 ± 0.052 s) with respect to the unimodal visual condition (0.99 ± 0.044 s) and secondly, by the fact that the bimodal condition did not lead to significant improvements (p > 0.05) in the performance precision (9.67 ± 0.44%) when compared to the visual feedback condition (11.23 ± 0.66%) ( figure 5(B)). Considering the trajectories, the bimodal condition displayed a significantly higher consistency (0.54 ± 0.009) than both unimodal conditions but also higher jerk (0.32 ± 0.011) with respect to the unimodal touch condition (0.27 ± 0.009) ( figure S5(B)).

Discussion
Since the last decade, electrical stimulation of peripheral nerves is the state-of-the-art technology employed for somatosensory feedback restoration in patients affected by sensory loss [2,44]. Invasive technologies (e.g. neural implants) have been widely investigated. Many studies focused on developing multiparameter stimulation patterns with the final goal of improving the naturalness of the stimulation, i.e. how much the artificial feedback faithfully reproduces the aimed tactile percept (e.g. touching an object with the intact limb, walking on the floor with the intact limb) [11,12,19,40,42,72]. Conversely, there is still room for advancements for TENS in terms of sensation quality and richness [6,10,57,58]. In fact, it has been shown that simple stimulation paradigms, such as the ones proposed by Shin et al [62], D'Anna et al [6] and Vargas et al [59,61], elicited sensations characterized by low intensity resolution, naturalness and intuitiveness.
To address this issue, we implemented multiparametric non-linear stimulation paradigms, featuring the simultaneous modulation of multiple parameters. To be able to reliably tune the sensation intensity while modulating multiple stimulation parameters, we evaluated the relationship between ACR and sensation strength resolution. The ACR was introduced by Graczyk et al [51], to group all the neural stimulation parameters in a single reference quantity of sensation magnitude. When introducing the ACR, Graczyk et al [51] adopted a linear relationship between PF and sensation intensity, resulting from a direct nerve stimulation in two upper-limb amputees. Other studies reported instead increase of perceived intensity as a logarithmic function of PF [11,49,73,74]. In our results, the co-modulation of PF and PW yielded statistically different JNDs and WFs. This result shows that the discriminability of two electrical stimuli is lower when varying the PF with respect to PW (e.g. the intensity variation achieved with a small PW percentage change can be obtained only with a much higher PF percentage change). The possible underlying reason may be that, while using TENS, the capacitive components of the skin and other tissue layers generate a frequencydependent non-constant impedance [6,63,[75][76][77], causing the non-linear relationship between PF and perceived intensity. Thus, we formulated a new ACR metric for TENS stimulation (ACR T ), respecting the logarithmic relationship PF-perceived intensity. Despite not preserving the ACR physical meaning (e.g. total charge delivered per second), this reformulation allowed to map ACR T variations with TENS perceived intensity variations resulting from PF and PW modulation. Therefore, the adaptation of ACR formulation was a necessary step to obtain equallyintense paradigms while modulating different parameters simultaneously (PF or PW). The reformulation ensured that PF variations would lead to lower ACR variations. The ACR T formulation as sensation intensity reference was tested in discrimination tests and functional assessments.
Since in invasive applications the naturalness of evoked sensations had been significantly improved by the simultaneous modulation of multiple stimulation parameters [19,41], we designed three multiparametric stimulation paradigms. Multiparametric 1 and Multiparametric 3 showed a sinusoidal modulation of PF and PW respectively, based on the hypothesis that it could lead to a more asynchronous fiber activation [11]. Multiparametric 2 was developed following a biomimetic approach [41,44,54]. The paradigms were compared with a standard PW linear modulation in terms of naturalness, through absolute questionnaires and a 2AFC task. To ensure that the naturalness judgement would not be biased by a different sensation intensity, the four paradigms were designed to have the same ACR T profile, and consequently, the same perceived intensity. As expected, all paradigms were perceived with the same intensity when evaluated independently from the others. Furthermore, the naturalness ratings did not show any statistical difference between paradigms, indicating that, differently from what was observed for invasive stimulations [19,41], a multiparametric parameter modulation in TENS does not improve the sensation naturalness or quality. This is probably due to the electrode-skin interface and the underlying tissue layers that affect the current propagation and the nerve activation [75,76,78]. In particular, the biomimetic paradigm may have failed since it was based on simulations run with a model tailored for invasive stimulation. Thus, a more faithful TENS biomimetic approach could be implemented with a model taking into account also the electrical properties of the tissue layers interposed between electrodes and nerve fibers [78]. Such approach could be characterized by a PF profile based on the fibers firing rate patterns and a charge profile generated according to the desired fiber recruitment. Furthermore, a low sensation naturalness may be caused by the impossibility of selectively targeting specific nerve fibers, as also suggested by Ortiz-Catalan et al [67]. Interestingly, although they were independently judged as equally intense, the paradigms showed difference in perceived intensity and naturalness when directly compared between each other. In fact, the trial results displayed a strong negative correlation between naturalness and intensity scores, meaning that the most intense paradigm was also considered as the least natural. This suggests that higher stimulation intensities will probably generate more electrical or unpleasant sensations, as also observed in invasive stimulation by Ortiz-Catalan et al [67] and in electro-cutaneous stimulation [45]. However, it is important to consider that intensity modulation is pivotal for performing a functional task requiring sensory feedback, such as the fine pressure of a pedal. Therefore, the tradeoff between absolute naturalness homologous to the input (pressure on an object) and intensity range of elicited sensations (and consequent intensity resolution) must be considered. For instance, sensations like vibration or tingling, perceived as qualitatively different, i.e. not homologous to pressure, could still be incorporated in the sensory motor scheme and used during a functional task, especially when homologous sensations cannot be elicited. Thus, sensation naturalness is a multidimensional concept (e.g. dependent on multiple factors). Despite not being judged as more natural, the multiparametric paradigm characterized by sinusoidal PF modulation (Multiparametric 1) allowed to correctly recognize the encoded tactile information more often than the other paradigms. This result is important considering that, regardless its conscious judgement, artificial feedback can be considered as natural, or intuitive, if it is promptly processed and exploited to effectively interact with the surrounding environment. This leads to a decrease of the cognitive effort required by the usage of neuroprosthetic devices [4,11,27,47], in particular when visual feedback is not available or unreliable, playing a pivotal role in enhancing the device acceptance.
The ease of processing the feedback information provided by the stimulation was evaluated through a functional assessment in a VR-TENS environment. In fact, while naturalness and intuitiveness rating task presented more subjective and descriptive results, the functional assessment allowed us to quantify those results through objective measures (e.g. accuracy, time). The higher intuitiveness of Multiparametric 1 was confirmed by the significant improvement in the functional performance with respect to the Standard encoding paradigm. In agreement with the results of Rock and Victor, we observed that when normal visual feedback was provided, vision was prevalent over touch [79]. However, when visual feedback was absent or impaired, the increase of the weight given to tactile inputs [80] hence highlighted the impact of the chosen paradigm on the task outcome: Multiparametric 1 resulted in a faster, more confident and more accurate performance showed by a statistically lower Error and Execution Time, higher trajectory consistency and higher level discrimination accuracy. Furthermore, the Multiparametric 1 paradigm also achieved a visuo-tactile integration (bimodal versus unimodal conditions) [70]. Instead, the addition of the Standard stimulation to blurred visual feedback significantly slowed down the performance, signalling an increased required mental effort and unsuccessful sensory integration of this paradigm with vision.
Therefore, although we assessed that it is difficult to achieve highly natural evoked sensation with TENS even with more sophisticated multiparametric stimulation (in striking contrast to those tested with implantable technology), a stimulation based on sinusoidal PF, resulted to be more intuitive and successfully integrated with visual feedback in a functional task. This indicates that a mixed stimulation encompassing non-linear multiparametric modulation is able to convey real-time, detailed information about the feedback dynamics, resulting in a better functional outcome in closed-loop with respect to the standard paradigms. Possible reason for it could be due to the fact that non-linear multiparametric sinusoidal modulation resembles more closely the cumulative asynchronous fibers activation in the nervous afferents response [11]. The main limitation of non-linear multiparametric paradigms relies in more complex programming; indeed, the simultaneous sinusoidal modulation of two parameters requires highly programmable stimulators with fine modulation of all stimulation parameters. The VR-TENS platform allowed to disentangle different sensory cues (e.g. vision and touch) [81] and to evaluate the different contribution of each input in a modular environment. At the same time, it allowed to create realistic tasks mimicking real-life activities, such as a neuropathic patient suffering by sensory loss pressing a brake pedal. This hence suggests that TENS could be a valid low-cost, low-risk alternative to implanted neuroprosthetic devices in providing intuitive sensory restoration in functional activities. The implementation of multiparametric stimulation paradigms could enhance the feedback intuitiveness and usefulness particularly when other sensory inputs are not available.
A limitation of this work is the fact that our results were obtained on healthy subjects. In fact, given the significant performance improvements brought by multiparametric TENS feedback to healthy subjects, it would be interesting to evaluate whether its benefits could be even more outstanding in sensory loss patients. Moreover, the intuitiveness of the multiparametric feedback could be tested also in dual tasks, in which the subject has to focus on performing multiple activities at the same time.

Conclusion
This work has been focusing on the design, implementation and testing of innovative, sophisticated neurostimulation paradigms for somatosensory feedback restoration in neuroprosthetic devices based on TENS. These paradigms aimed to improve the naturalness, richness and informativeness of the TENSevoked tactile feedback. The relationship between stimulation parameters and perceived magnitude has been studied. The ACR, previously defined for invasive stimulation, has been adapted accounting for a lower contribution of PF and used to normalize the intensity perceived among different paradigms. Despite not consciously perceived as more natural, multiparametric stimulation paradigms were better integrated by the users, showing higher intuitiveness and visuo-tactile integration during functional tasks in VR. This suggests that TENS could be successfully integrated in non-invasive neuroprosthetic devices to assist sensory loss patients in motor activities, improving their autonomy and life quality. Moreover, the implementation of multiparametric stimulation paradigms could enhance the feedback intuitiveness, lowering the required cognitive effort and leading to higher device acceptance and embodiment.

Data availability statement
The data cannot be made publicly available upon publication because they contain sensitive personal information. The data that support the findings of this study are available upon reasonable request from the authors. the figures; N K developed the FootSim model; G A and M R developed the algorithms and the software; G V supervised the analyses, discussed the results, reviewed the manuscript; S R designed the study, supervised the analyses, discussed the results and reviewed the manuscript.