Sensitive FBG-based e-skin

The experimented e-skin (Extended Data Fig. 5) had size and shape similar to those of the human forearm. To mimic biomechanical compliance, a soft silicone polymer (Dragon Skin 10 medium, Smooth-On) was used to realize a flexible patch to embed an optical fibre with 21 FBG sensors. It was fabricated by a three-step process¹¹: (1) a first silicone pouring in a custom mould to obtain a bottom layer with a groove; (2) the embedding of the fibre along the grooved path; (3) a second silicone pouring to cover the fibre itself and provide robustness to the patch. The e-skin was then attached to a custom three-dimensionally printed plastic support to guarantee stability. The optically sensitive elements were FBGs (length, 10 mm) distributed along the designed path running throughout the e-skin resemble the sensitivity features of the human forearm, that is, with an increasing density from the proximal to distal parts. These gratings are interference patterns inscribed within the core of a single optical fibre. When broadband spectrum light, sent by means of an optoelectronic interrogator, passes through the fibre cable and hits an FBG, it is partly back-reflected to the input point, as a spectral peak centred at a characteristic wavelength, λ_B:

where n_eff is the effective refractive index of the optical fibre and Λ is the pitch of the specific grating⁸⁰. The λ_B value of the e-skin FBGs were, hence, defined by fabrication and ranged between 1,520 nm and 1,580 nm with a step of 3 nm. Wavelength shifts occur when a deformation is applied to the sensor¹¹, thereby allowing the encoding of tactile stimuli. Considering this feature and to reduce the effects of potential prestrains that could have impacted the residual deformability of the sensors in response to external tactile stimuli, the FBGs were arranged on straight fibre segments.

Indentation data collection

The training and validation of the proposed model required the collection of a dataset of indentations applied over the e-skin surface. For this purpose, an automatized indentation protocol was executed, with 1,846 target points to touch, computed as the centroids of the triangles of a random mesh of the skin surface. A bimanual robotic platform, consisting of two anthropomorphic arms (Racer-5-0.80, Comau), was programmed to cooperate during the indentation protocol (Extended Data Fig. 5): one robot performed the tactile tasks, whereas the other held and oriented the e-skin. The indenting robot reached each target site perpendicularly, touching the skin with a spring-like hemispherical probe attached to the end-effector. Force feedback measured by means of a load cell (Nano 43, ATI Industrial Automation), mounted at the base of the indenter, was provided to control the robot to release the contact when exceeding a threshold (2.5 N) and then fly to a new target point. The control code consisted of dedicated routines (LabVIEW, National Instruments) running on both an industrial controller (IC-3173, National Instruments), to send commands to the robots, and a PC, for managing and checking the experiments. During the indentation sessions, force and FBG signals, the latter read by means of an optical interrogator (FBG Scan 904, FBGS Technologies), were collected at a 100-Hz rate, and saved into text files together with the Cartesian coordinates of the contact sites.

Assessment of receptive fields of e-skin FBG sensors

The receptive fields of the FBG sensors populating the skin were represented as two-dimensional contour maps depicting the spatial distribution of sensitivity to contact force. Hotspot regions in dark blue indicate the maximum sensitivity (Fig. 1b and Extended Data Fig. 1). To obtain the FBG receptive fields, the wavelength variation signals (Δλ) from each FBG sensor were averaged over 500 ms captured on the highest force plateau, that is, the second force level, and divided by its average value. The obtained results were then normalized with respect to the highest FBG sensitivity value, to have sensitivity maps ranging from 0 to 1. This procedure was performed using all the available indentation samples.

First-order neurons: type II PA model

The 21 FBG signals were the input of the PA models. Each sensor signal was resampled at a frequency of 1 kHz to establish a regular firing dynamics consistent with the discretization of the Izhikevich neuronal model⁷⁴. The FBG signal derivatives were computed to obtain the dynamic components of the indentations, that is, the loading and unloading phases of the varying force levels applied to the e-skin. The absolute values of both raw and differentiated signals were calculated to feed the Izhikevich neurons as input current I, to mimic the adaptation dynamics and firing response of type II mechanoreceptors, that is, Ruffini endings (SAII) and Pacini corpuscles (FAII)³⁸. To simulate the physiology of type II PAs in terms of the ratio between slowly and fast-adapting units, that is, two SAIIs to one FAII³⁶, their diverse sensitivity and the different-sized overlapping receptive fields, the raw absolute FBG signals and their derivatives were multiplied by four and two distinct gains. Hence, six mechanoreceptor models were considered for each FBG sensor (Fig. 2), resulting in 126 input currents for the first-order neuron models. Extended Data Table 1 shows the gains used for the SAII and FAII models. The SAII gains were selected through some pilot tests and empirical analyses, where we assessed whether the increase in gain was generating an increase in the firing rate within the physiological range documented in the literature⁴³. Specifically, the average ISIs (that is, the inverse of the average firing rates) corresponding to each force range were computed, and their logarithmic values were fitted to the logarithm of the stimulus intensity percentage (force values normalized by the maximum of 4 N; Fig. 2h). The results were then compared with background neurophysiological findings⁴³, particularly to the ISI characterization of SAII units. The FAII gains were selected following the same empirical approach, that is, by evaluating the spiking activity at load transients. It is worth noticing that a general procedure for establishing the gains of the input currents of the Izhikevich model for artificial mechanoreceptor simulation has not yet been established in the literature. Thus, several studies, such as ref. ⁵⁸, reported this empirical method of gain selection according to the specific application.

The Izhikevich model consists of a system of differential equations that can be solved via the Euler method⁷⁴. Equation (2) describes the dynamics of the membrane potential v of the neuron; equation (3) represents the recovery variable u, responsible for the repolarization of the neural membrane; equation (4) describes the membrane potential restoration and the recovery variable update when v reaches the 30-mV firing threshold:

where A, B, and C are standard variables of the Izhikevich model; t represents time; a is the timescale of u; b is the sensitivity of u to the membrane potential; c is the membrane resting potential value; and d modulates the dynamics of the after-spike reset of the recovery variable u. For the implemented mechanoreceptors, we choose the parameters Gain_n (Extended Data Table 1) that reproduced a regular firing behaviour.

Encoding of contact intensity by SNN PAs

The activity of type II mechanoreceptors was analysed to investigate relationships with the intensity of the applied load. For this purpose, the perpendicular force profile was extracted for all the 1,846 indentations. For each indentation, SAII and FAII units were considered separately. First, their corresponding firing rates were calculated on moving time windows of 100 ms, with an overlap of 50 ms. For each window, the firing rate was extracted as follows:

where the subscript X indicates the type of mechanoreceptor (SAII or FAII), w_t is the duration of the time window, N_spikesX is the number of spikes of the X units in w_t and N_{neuronXactive} is the number of active neurons of each type in the selected window. The resulting firing rates were compared with the corresponding force profile (Fig. 2a–f). Then, the SAII firing rates were grouped for force levels, ranging from 0 to 4 N, with a step of 0.5 N. The corresponding distributions were analysed and the relevant statistics computed (median values and first and third quartiles; Fig. 2g and Supplementary Table 1). For each load-intensity level, the mean SAII firing rate value was extracted to further derive the relationship with the applied loads. For that, the linear Pearson correlation coefficient ρ was calculated. Then, the first-order polynomial that fitted the logarithm of the SAII average firing rates to the logarithm of the percentage of the applied load, normalized by 4 N, was extrapolated. The R² coefficient was computed to assess the goodness of the linear fit (Fig. 2h).

Second-order neurons: CN model

The second layer of the SNN was composed of 1,036 CNs, modelled with the mathematical implementation based on the exponential integrate-and-fire approach²⁹. This model reproduces the complete dynamics of the membrane potential of CNs, as described in equation (5), along with a detailed modelling of the activity of low-threshold voltage-gated calcium channels and calcium-activated potassium channels:

In equation (5), C_m is the capacitance of the neural membrane; V_m is the CN membrane potential; I_L is the leak current of the neuron (equation (6)); I_spike is the spike current that recreates the action potential onset and the fast neuron depolarization (equation (7)); I_ion is the ionic current resulting from the summation of the currents of voltage-gated calcium channels (I_Ca) and calcium-activated potassium channels (I_K; equation (8)); I_ext is the external current that can be injected into the neuron (in this study, it is equal to 0); and I_syn is the synaptic current (equation (9)), where each synapse (i) is activated by a PA.

where w_exc,i is the excitatory synaptic weight; t* is the time at which a spike occurs and w_inh is the inhibitory synaptic weight. I_syn can be, hence, expressed as the sum of the excitatory and inhibitory synaptic currents of all the synapses of the individual PAs. The definitions of the other variables and their respective values are reported in Supplementary Table 6.

The current I_Ca is described by equation (10) and the current I_K, by equation (11):

where \({\bar{g}}_{{{\rm{Ca}}}}\) and \({\bar{g}}_{{\rm{K}}}\) are the maximum conductance; E_Ca and E_K are the action potentials of the calcium and potassium channels, respectively; and x_Ca,a, x_Ca,i, \({x}_{{{\rm{K}}}_{{{\rm{Ca}}}}}\) and \({x}_{{{\rm{K}}}_{{{\rm{V}}_{\rm{m}}}}}\) are the activity states of the channels. These values are provided in Supplementary Table 6.

In this model, both ion channels (x_Ca,a and x_Ca,i) are considered as sources of the membrane calcium concentration of the CN ([Ca²⁺]). Thus, we proposed that the activity of the total neuron calcium concentration follows equation (12).

The relevant parameters are reported in Supplementary Table 2.

Functional organization of SNN second-order neurons to model a somatotopic map of the e-skin

To functionally mimic the organization of the cuneate nucleus, we generated an e-skin grid mesh with 29 × 38 edges, that is, 28 × 37 subregions. In this mesh, each subregion represented a CN, and its centroid represented the centre of its receptive field. In total, 1,036 CNs were created with circular overlapping receptive fields of 41.67-mm radius each. The criterion for the definition of the size of the CN receptive field was to include at least two FBG sensors in the corresponding sensitive area of the e-skin and, thus, 12 PAs, consistent with the number of dominant PAs that project into a CN found in neurophysiological studies on mammals and simulations of human tactile perception²⁴. This approach allowed us to determine a CN somatotopic map of the e-skin.

Connectivity of SNN and initialization of synaptic weights

In our model, each of the 1,036 CNs was fully connected to the 126 PAs, with excitatory weights ranging from 0 to 1, and to a single inhibitory synapse (Fig. 1a, red triangle), with an IN that grouped the responses of all the 126 PA projections (Fig. 1a, coloured empty triangles). This connectivity enabled the potentiation or depression of all the synaptic weights during the learning process. In this model, the magnitude of the postsynaptic potential projected into a CN for a given sensory input is dependent on both excitatory (w_exc,i) and inhibitory (w_inh) weights, as described in equation (9)^18,29. The presynaptic excitatory weights were initialized as the inverse of the Euclidean distances between the position of the 21 FBG sensors in the skin and the 1,036 centroids of the CNs. Then, the obtained values were rescaled between 0.2 and 1. The initial excitatory weights (w_exc) of the PAs corresponding to the sensors outside the area of the CN receptive fields were instead set to 0 (Fig. 1a, coloured dots). In this way, only the mechanoreceptors inside the CN receptive field region had initial weights greater than 0 (Fig. 1a, coloured filled triangles). All the inhibitory synaptic weights (w_inh) were instead initialized to 0.125.

Synaptic learning protocol of SNN second-order neurons

Inspired by calcium-dependent synaptic plasticity, we implemented a synaptic learning process²⁹, where the total inhibitory weights (w_inh) and presynaptic excitatory weights (w_exc,i) were updated at each stimulus presentation. According to our model, the excitatory weight potentiation of the presynaptic neurons (PAs) occurred when the total calcium activity (\({A}_{{{\rm{Tot}}}}^{{{{\rm{Ca}}}}^{2+}}\)) of the CN, as in equation (13), was strongly correlated with the local calcium activity (\({A}_{{{\rm{Loc}}}_{i}}^{{{{\rm{Ca}}}}^{2+}}\)) of a synapse (PA with the secondary neuron; equation (14)); otherwise, this synapse was depressed:

where \({A}_{{Tot}}^{{{Ca}}^{2+}}\) is the total calcium activity of the CN; k_act = 1 is an arbitrary constant; \({A}_{\rm{Loc}_{i}}^{{\rm{Ca}}^{2+}}\) is the local calcium activity due to a synapse i; τ_r = 4 ms is the rise time; τ_d = 12.5 ms is the decay time; τ_l = 0 ms is the latency time; τ_l = 21 ms is a constant to calculate the ratio; and t* is the time at which a PA spike occurs. To reach supralinearity in the local calcium activity, we used an approximative approach of subtracting an offset, that is, the 75% of the single-pulse activation peak activity, from the local calcium signal, and this resulting value was used as the local calcium activity (\({A}_{{{\rm{Loc}}}}^{{{{\rm{Ca}}}}^{2+}}\)), as proposed in ref. ²⁹.

In the model, the individual excitatory weights actualization occurred during the presentation of each stimulus during the learning phase, and it is described as the integral of the correlation between local calcium activity (\({A}_{\rm{Loc}}^{{\rm{Ca}}^{2+}}\)) and total calcium activity (\({A}_{\rm{Tot}}^{{\rm{Ca}}^{2+}}\)), following equation (15):

where \({{{\rm{Avg}}}}_{{A}_{{\mathrm{tot}}}^{{{\mathrm{Ca}}}^{2+}}}\) is the average of the last three values of the total calcium activity; Syn_EQ is the synaptic equilibrium defined as a linear function of the total excitatory synaptic weight with a dual slope having point 0 in 10 (decay = 0.04, if ∑w_exc < Syn_EQ; decay = 0.12, if ∑w_exc > Syn_EQ); and K is a constant gain factor defined by the sigmoid function represented in equation (16), with a gain step on the slope of 0.005.

To avoid instabilities in synaptic learning, it was necessary to average the calcium activity; thus, the learning threshold was obtained by multiplying the average of the total calcium activity \({{\rm{Avg}}}_{{A}_{{\mathrm{tot}}}^{{{\mathrm{Ca}}}^{2+}}}\) by the synaptic equilibrium Syn_EQ.

Although the inhibitory synaptic learning is based on the firing rate of the CN calcium channels, w_inh, initially set to 0.125, is responsible for regulating the total activity of the calcium channels in the CN. Thus, the low activity of calcium channels was counteracted with a decrease in w_inh, and vice versa. For the model implemented in this work, w_inh was used to maintain the firing frequency of the calcium channels at a predefined set point of 20 Hz, and the update of w_inh was given by a dual-slope function zeroing at this set point. The synaptic learning protocol is summarized in a pseudo-code in the Supplementary Text.

The synaptic learning process was evaluated through a fourfold validation, such that for each training fold, 1,385 indentation samples were shuffled and presented to the model, thereby totalling 1,385 training cycles. In this step, only 1 s of data was used, including the beginning of the second indentation step (that is, the highest force level for each indentation). These data were the spikes (binary vector) of the first neuronal layer of mechanoreceptors (PAs). All the synaptic learning processes were separately performed for each of the 1,036 CNs.

Data post-processing and SNN performance evaluation

After the synaptic learning, the spike responses of 1,036 CNs were calculated for the indentations of each validation fold and processed to derive the stimulus intensity and location information.

Decoding of stimulus intensity through spiking activity of SNN second-order neurons

To demonstrate the increase in CN spiking activity when progressively higher loads were applied, the cumulative spike numbers inside the receptive field of the stimulated CN were computed for different loading conditions (Extended Data Fig. 2). To this aim, we evaluated the CN spikes in a 1-s window at the beginning of each of the two plateaus of the indentation force profile. In addition, the spatial activations of the CNs for the different force levels were explored. In this regard, the number of active CNs was calculated for the two load conditions and 1,655 of the applied indentations (Fig. 4d). A Mann–Whitney U-test test with a 0.05 significance level was performed to investigate whether the spatially distributed responses of the CNs could encode information on different intensities of the applied tactile stimuli.

SNN performance and temporal resolution for localization of stimuli on the e-skin

Regarding the prediction of the contact position, the spikes of the CNs were first summed over the indentation duration window. The sum of each neuron spike (Nspk_i) was then used for the weighted location estimation, as described in equation (17).

where WL_x,y is the location estimation in the skin x–y plane (considering the two-dimensional projection of the e-skin surface) weighted by the neural activation of the CNs (weighted position); Nspk_i is the number of spikes of a neuron in a time window; i is the CN ID (from 1 to 1,036); and N_iLoc_x,y is the position of the centroid of a second-order neuron in the x–y plane of the modelled functionally organized cuneate nucleus. We estimated the weighted location on a 1-s sliding window with an overlap of 100 ms. To evaluate the location prediction error, we calculated the Euclidean distance between the real position of the indentation on the e-skin and the estimated weighted location (Fig. 5).

In addition, we compared the localization performance of the SNN with those achieved through the weighted average of the 21 FBG sensor wavelengths (FBG WL_x,y), calculated as follows:

where Δλ_i and (x_FBG, y_FBG)_i are the wavelength variation and two-dimensional coordinates of the ith FBG sensor, respectively. The FBG-based weighted average was computed for each timestamp of the corresponding non-zero force values. Then, the average location and Euclidean distance were calculated to obtain a single estimation and a single error value per indentation of the test sets of the four folds (Extended Data Fig. 3). The performance of both neuronal network and FBG weighted average in localizing the test stimuli were also assessed considering ROIs of different sizes. Starting from the centre of the skin, circular areas with an increasing radius-delimited (from 20 mm to 55 mm) portions of the e-skin surface were used to compute the contact estimation error, considering that all the test indentations fell inside the circular region. Median and IQR error values were calculated for each region size, grouping the results of the four cross-validation folds (Fig. 6c and Supplementary Table 4).

We then evaluated the temporal resolution of our system to localize the applied indentations. Specifically, considering that the sensors of the e-skin tracked the dynamics of the applied loading and unloading transients (Fig. 2a–c), we estimated the minimum amount of time needed by the trained SNN to reach its best localization performance. This approach aimed at investigating the stimulus conduction time from the simulated forearm mechanoreceptors, that is, the FBG sensors, to the modelled CNs, designated to localize contacts onto the e-skin. We computed the model localization error (median and IQR) on the second step of the force profile on time windows of different lengths, ranging from 10 ms to 2 s, the latter being the duration of the complete second level of the force profile. The two extreme conditions were statistically compared with a level of significance (α) of 0.05.

Weber two-point discrimination test and analysis

To apply two-point stimuli onto the e-skin through the two-point discrimination test (Supplementary Video 3), a set of custom probes was designed (Fusion 360, Autodesk) and three-dimensionally printed in polylactic acid (Ultimaker S5). The probes featured one or two hemispherical tips (radius, 1.5 mm) with distances ranging from 20 mm to 60 mm, with a step of 5 mm. Each probe was attached to a handle stick by means of a conical spring fixed to a sliding structure, which provided the tool with a translational perpendicular degree of freedom. This structure presented a blocking mechanism so that the tool could be pushed until bottoming out, allowing to control the applied force. Ten experimenters participated in the trials, which consisted of manually indenting the e-skin by means of each of the above-mentioned probes in seven random locations. The single-probe trials were repeated before each of those with two-tip probes. At each indentation, the experimenter had to gently land on the e-skin with the probe, stabilize the posture, push the tool until the descent was not blocked and keep it for 3 s. The protocol instructions and the timing of the different indentation phases were displayed via a custom graphical user interface developed in LabVIEW 2019 (National Instruments). The same software routine allowed us to store the FBG sensor signals, collected by means of an optical interrogator (FBG Scan 904, FBGS Technologies), the corresponding timestamps and the trial information into text files for further analysis. The recordings of the FBG wavelengths were, at first, segmented to focus on the period when the probe was pressed against the e-skin. A Δλ threshold of 0.003 nm was considered to identify the onset of the contact, that is, when at least one sensor signal exceeded it for 200 ms. This corresponded to the beginning of the landing phase. Then, to identify the onset of indentation, a second threshold, that is, twice the average absolute wavelength of all the sensors during the 1 s after the e-skin was touched, was considered. The start of the stimulus corresponded to the moment when the mean value of all the FBG Δλ values exceeded this threshold for at least 500 ms. The resulting indentation portions were zero padded with 500-ms periods to simulate the transition from non-contact to contact states. Given the variability of indentations performed manually by the experimenters, some trials were discarded (Supplementary Table 5). The obtained FBG signals for the selected indentations were fed into the pretrained neuronal model to get the spiking response of the 1,036 CNs. Then, the neuronal response latency was determined within a 200-ms temporal window starting from the first CN spike. These latencies were organized topographically (Fig. 6d,e); hence, a 3 × 3 Gaussian spatial filter (σ = 0.5) was applied. The lowest minimum in the resulting latency map was identified as the first contact point. Then, the occurrence of a second minimum was checked outside a circular region corresponding to the size of the CN receptive field⁴⁴. If two distinct minima were found following this criterion, it was determined that two contact points were applied to the e-skin. The two-point detection rates were then estimated for all the trials, and a psychometric curve was fitted to them. The resulting piecewise logistic function had the following expression:

Finally, the 0.75 probability, that is, when F(x) was 0.75, was considered to determine the e-skin two-point discrimination threshold (Fig. 6f).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.