A functional spiking neuronal network for tactile sensing pathway to process edge orientation

To obtain deeper insights into the tactile processing pathway from a population-level point of view, we have modeled three stages of the tactile pathway from the periphery to the cortex in response to indentation and scanned edge stimuli at different orientations. Three stages in the tactile pathway are, (1) the first-order neurons which innervate the cutaneous mechanoreceptors, (2) the cuneate nucleus in the midbrain and (3) the cortical neurons of the somatosensory area. In the proposed network, the first layer mimics the spiking patterns generated by the primary afferents. These afferents have complex skin receptive fields. In the second layer, the role of lateral inhibition on projection neurons in the cuneate nucleus is investigated. The third layer acts as a biomimetic decoder consisting of pyramidal and cortical interneurons that correspond to heterogeneous receptive fields with excitatory and inhibitory sub-regions on the skin. In this way, the activity of pyramidal neurons is tuned to the specific edge orientations. By modifying afferent receptive field size, it is observed that the larger receptive fields convey more information about edge orientation in the first spikes of cortical neurons when edge orientation stimuli move across the patch of skin. In addition, the proposed spiking neural model can detect edge orientation at any location on the simulated mechanoreceptor grid with high accuracy. The results of this research advance our knowledge about tactile information processing and can be employed in prosthetic and bio-robotic applications.

To obtain deeper insights into the tactile processing pathway from a population-level point of view, we have modeled three stages of the tactile pathway from the periphery to the cortex in response to indentation and scanned edge stimuli at different orientations. Three stages in the tactile pathway are, (1) the first-order neurons which innervate the cutaneous mechanoreceptors, (2) the cuneate nucleus in the midbrain and (3) the cortical neurons of the somatosensory area. In the proposed network, the first layer mimics the spiking patterns generated by the primary afferents. These afferents have complex skin receptive fields. In the second layer, the role of lateral inhibition on projection neurons in the cuneate nucleus is investigated. The third layer acts as a biomimetic decoder consisting of pyramidal and cortical interneurons that correspond to heterogeneous receptive fields with excitatory and inhibitory sub-regions on the skin. In this way, the activity of pyramidal neurons is tuned to the specific edge orientations. By modifying afferent receptive field size, it is observed that the larger receptive fields convey more information about edge orientation in the first spikes of cortical neurons when edge orientation stimuli move across the patch of skin. In addition, the proposed spiking neural model can detect edge orientation at any location on the simulated mechanoreceptor grid with high accuracy. The results of this research advance our knowledge about tactile information processing and can be employed in prosthetic and bio-robotic applications.
Our ability to touch and feel is usually considered a simple task in daily life. However, there is a complex process behind our sense of touch, from the sensory receptors within the skin to the brain's neuronal activity. When we interact with an object, the activation of different kinds of mechanoreceptors in the skin send various signals regarding object characteristics such as texture, shape, and size. Single-unit recordings of tactile afferents have been used to discover how tactile information is encoded 1,2 . One of the main findings is that stimulus information is distributed over a population of fibers in the form of spike trains. However, our understanding of populationlevel coding is still basic and only a few models have tried to address this important representation [3][4][5] . Indeed, tactile processing involves neural mechanisms that extract high-level features of a stimulus, such as an edge orientation, by integrating information from many low-level inputs 6,7 . The low-level inputs originate from the cutaneous mechanoreceptors which are located in the skin throughout the body and convert skin deformation into spiking responses. The density of mechanoreceptors varies across the body. For example, the density of cutaneous afferents that innervate the human fingertip is about 240 units/cm 2 , whereas in the palm it is only 58 units/cm 28,9 . In general, the hands and lips for primates have the highest density of mechanoreceptors while the legs and back contain the lowest. One of the main features of these first-order neurons is that they have complex receptive fields. Indeed, the afferent's distal axons which innervate the skin form many transduction sites 4,10 . This innervation pattern is a critical peripheral nerve mechanism that allows individual afferents to encode geometric information about tactile stimuli. To the best of our knowledge, the spatial complexity and heterogeneity of innervation patterns have not been integrated in models, such as in the excellent work by Saal et al. 3 .

Results
The primary objective of this study is to simulate the tactile pathway from the glabrous skin of the fingertip to the cortical neurons in the somatosensory cortex area 3b in response to edge orientation stimuli. In this way, we discover how the receptive field size modulates information about edge orientation in the spiking pattern of cortical neurons. In the first layer, the SA-I and RA-I afferents are simulated which innervate their mechanoreceptors in a random sampling pattern. Edge stimuli activate a population of mechanoreceptors that send tactile information to the upper layers. The next layer, the CN, integrates and categorizes the coming information. We examine the role of lateral inhibition in the PN population and demonstrate how insufficient inhibitory currents lead to a failure in suppressing skin hypersensitivity. In the last layer, modeling the excitatory and inhibitory receptive fields of the somatosensory cortical neurons shows the classification of edge orientation independent of stimulus location. Two experiments are done; (i) Indented edge experiment: the stimulus vertically indents the simulated patch of skin. (ii) Scanned edge experiment: the stimulus indents and moves across the skin in a specific direction.
First-order neurons. In the first layer, we simulated two afferent populations, namely, SA-I and RA-I. These afferents densely and randomly innervate the skin mechanoreceptors of the human fingertips and encode tactile information through their spike rates and spatiotemporal spiking patterns (Fig. 1A). A set of indentations are applied as the inputs to each afferent model (Fig. 1B) to produce spiking responses (Fig. 1C). The stimulus is the edge at different orientations (from 5°-80°, in 5° increments) indented on the mechanoreceptor grid of the fingertip. The different weights of mechanoreceptors make it possible for afferents to have spatially complex receptive fields (Fig. 1D). Recently, it has been proposed that random innervation of afferents might be a peripheral neural mechanism for extracting geometric features of the touched objects 4 .
The SA-I dynamic model. Myelinated SA-I afferents innervate the skin through unmyelinated fibers (neurites) to form synaptic-like connections with Merkel cells. Indeed, SA-I afferent's end organ is considered as a Merkel cell-neurite complex 27 (2) neurites that produce rapidly adapting firing to mechanical stimulation [29][30][31] . Here, we propose a biologically plausible and functional model to mimic the Merkel cell and neurite response to mechanical indentation stimuli. Unlike previous studies which suggest that the SA-I model responds only to the hold phase of mechanical stimuli [32][33][34] , the new model is composed of a Merkel-cell component (input current) and a neurite component (a derivative of input current). The effect of the SA-I dynamic model on the performance of cortical neurons using spike count analysis is illustrated in Fig. 2. Considering indentation noise in different trials, it is found that the classification performance has a stable characteristic when a biomimetic decoder 35 (see Methods) is used ( Fig. 2A) and consequently it has more robustness. Afferent responses are assessed using Principal Component Analysis (PCA). In this way, considering the first 3-principal components of the feature space (spike count), simulations show that the K-nearest neighbor classifier (KNN) has higher performance for the SA-I dynamic model in comparison with the SA-I static model (Fig. 2B). Furthermore, the RA-I afferent responses against indentation noise are also reported in Fig. 2B. Noteworthy, the number of RA-I afferents is almost double of the SA-I ones (see Methods) and its performance is slightly superior. This figure also illustrates that the initial spiking responses of the SA-I model make rapid and accurate orientation detection.
Decoding the indented and scanned edge orientation. In the indentation experiment, we investigated how edge orientation could be decoded and extracted from the SA-I and RA-I population responses. We consider 2 or www.nature.com/scientificreports/ 3 principal components obtained by applying PCA to the spike counts calculated from firing patterns. In our simulations, the Pacinian afferents have not been considered because they are known to contribute more to texture and vibrotactile encoding 3,8 . In Fig. 3A, as the orientation changes, the edge stimulus lands on the receptive field of different afferents, therefore, information is encoded in the spatiotemporal pattern of activation. Results show that contact stimuli are encoded in the first emitted spikes, and 2-principal components are as informative as 3-principal components which are extracted from the feature space ( Fig. 3B). At the population level, the performance of RA-I is higher than SA-I given the fact that the number of RA-I afferents is greater than SA-I which improves the spatial encoding. Based on the experimental evidence, each type of afferent extracts particular geometric features, and therefore, the combination of both afferents improves the recognition performance. We also simulated the responses of two afferent populations to the scanned edge stimulus on the simulated patch of skin for 16 different orientations (Fig. 4A). The goal is to classify the edge orientation based on the responses at the population-level considering both spike count (number of spikes) and temporal spiking pattern (spike timing). Indeed, afferents randomly innervate the mechanoreceptors and create highly sensitive areas that are non-uniformly distributed within the receptive field 4 . This arrangement acts as a peripheral neural mechanism that allows individual neurons to extract distinct features of the touched objects. When skin is deformed, each afferent produces unique responses and facilitates edge stimulus discrimination. To calculate the temporal distance between spike trains, Victor-Purpura distance (VPd) has been used 36 . Figure 4B shows the neural response when an edge stimulus is moving across the receptive fields of two sampled tactile afferents. In Fig. 4C, as can be seen, the RA-I responses rapidly show the orientation stimuli with high accuracy compared to the SA-I responses. Figure 4C shows the comparative performance analysis for decoding the orientation of the scanned edge to variable time windows. Both spike count (spatial activation pattern) and spike timing-based features are used for orientation detection. Although the classification accuracy for both methods, spike count and spike timing, eventually saturates, the performance curves of the spike timing method (dashed lines) have a delay with respect to the performance curve of the spike count method (solid lines). This www.nature.com/scientificreports/  www.nature.com/scientificreports/ interestingly shows the spatial activation pattern (which is specified by counting the number of spikes) of the early spikes conveys more information compared to the temporal patterns of the afferent responses.
To assess the behavior of the afferents at a single-unit level, one afferent (SA-I and RA-I) is considered which has its receptive field (Fig. 5A). During the experiment, the receptive fields are scanned by 16 edge orientations and each stimulus is repeated 6 times (96 trials in total). Since RA-I afferents have larger receptive fields in comparison to SA-I, the obtained spiking responses are more diverse (Fig. 5B). Figure 5C shows the cross-correlations between the input current of afferents. As can be seen in Fig. 5C, for the SA-I afferents, in lower angles of orientation, correlation coefficients are more similar. Thus, angle discrimination using correlation coefficients in RA-I is easier than SA-I. Figure 5D illustrates the spike timing analysis by computing the VPd metric which illustrates that the RA-I spike responses are more informative than SA-I. Employing the extracted features from Fig. 5C,D, the classification accuracy to classify the edge orientation is reported in Table 1.
Second-order neurons (cuneate nucleus). The second layer of the proposed functional neuronal model considers the projection neurons (PN) and interneurons (IN) to simulate the cuneate nucleus (CN) structure in tactile sensing. The CN neurons mainly contribute to tactile information processing through lateral inhibition which leads to the firing rate filtering 12 . The INs make inhibitory synaptic connections with the PNs (Fig. 6A). We have hypothesized that the lateral inhibition in the CN model prevents the accumulation of noise and preserves spatial information which in turn facilitates the edge detection process. In this way, the impact of CN lateral inhibition on the performance of orientation detection is explored when the edge indents the simulated mechanoreceptor grid. Indeed, we tune the network for the best performance of edge orientation detection. With maximum lateral inhibition, input currents of some PNs decrease significantly. On the other hand, reducing the amount of inhibitory current from INs to the PNs (Fig. 6B,C), leads to an increase in the firing rate of the cortical neurons. A biomimetic decoder is considered by employing the winner-take-all algorithm (a cortical neuron with a maximum firing rate is the winner neuron) (see Methods). Figure 6D illustrates that diminishing the amount of lateral inhibition leads to a decrease in the recognition accuracy of the indented orientations. Indeed, lateral inhibition preserves spatial information and thus facilitates the edge detection process. Firing responses of all tactile processing stages are plotted in Fig. 7A,B, when two edge orientation stimuli indent the mechanoreceptor grid at two levels of inhibition in the CN.
Cortical neurons in the somatosensory area. In the last layer, the responses of cortical neurons to orientation stimuli are simulated. The obtained results are in agreement with the experimental results reported in the literature 35,37 when an edge indents the skin. For simplicity, in our simulations, cortical receptive fields consist of two sub-regions on the skin (excitatory and inhibitory), and they are located close to each other and are tuned for a specific orientation. These receptive fields are distributed over the skin, therefore, edge orientation can be recognized accurately and rapidly by the cortical neurons even when the position of the edge stimulus on the skin is changing. Indeed, cortical neurons in area 3b are sensitive to specific orientations. Specifically, when an edge indents the skin, the input currents of orientation-tuned cortical neurons are raised based on the spatial coincidence between the receptive field sensitive zones of cortical neurons and the edge orientation on the skin. To better understand the function of the cortical neurons, the excitatory and inhibitory currents are plotted (Fig. 8A). It is evident that edge orientation can quickly be decoded using the winner-take-all algorithm. As shown in Fig. 8A, the current amplitudes of the orientation-tuned cortical neurons are increased after 20ms to distinguish edge stimuli, which is compatible with the conduction delay to the cortex in neurophysiological observations 37 . To investigate how the biomimetic decoder recognizes edge orientations indented at different skin positions, we considered three locations with 1.2 mm spacing along the y-axis on the simulated patch of skin. For each trial, the stimulus was applied randomly to one of these locations. The orientation is decoded based on the "winner-take-all" approach (biomimetic decoder). Furthermore, to compare the performance with a common machine learning algorithm, the PCA for feature reduction and the KNN as the classifier are used. The spike count is used as the input features. To explore the evolution of decoding performance over time, the procedure was repeated for 60 progressively expanding time windows, ranging from 0 to 180ms after stimulus onset (Fig. 8B). As can be seen in Fig. 8B, the biomimetic decoder has better performance and achieves 97% accuracy around 180ms . Confusion matrices for both decoders are shown for the 180ms time window (Fig. 8C). Another interesting simulation is to explore the effect of afferents' receptive field size on the orientation recognition by the cortical neurons. The procedure for forming the receptive fields has been described in the Methods section. Here, we examine how the scanned and indented edge stimuli are encoded by the population of cortical spiking neurons when the standard deviation (SD) of the mechanoreceptor distribution is changed (Fig. 8D). Interestingly, it is observed that the larger receptive fields of afferents convey more information about the scanned edge in the first spikes of population-level of cortical neurons for both spike timing and spike count protocols (Fig. 8E). Conversely, in the indented edge experiment, the larger receptive fields make the orientation detection worse (Fig. 8F). The responses of all layers to an edge indentation stimulus are shown and summarized in Movie S1. Figure 9A illustrates sample receptive fields with excitatory and inhibitory sub-regions on the skin. An excitatory (inhibitory) sub-region indicates that any indentation at that location gives rise to an increase (decrease) in neuronal firing. This figure shows that the intensity of neural firing indicates the edge orientation. In this way, as the degree of spatial coincidence between the neuron sensitive zones and tissue deformations caused by edge indentation increases, the neural firing also increases. For a given neuron, some edge orientations show more spatial coincidence than others and therefore yield stronger responses ( Fig. 9B and Movie S2). Raster plots of cortical neurons are depicted for two edge orientations (10° and 40°); the input stimulus is recognized by the highest spiking rate (Fig. 9C). www.nature.com/scientificreports/ www.nature.com/scientificreports/  One of the key features of the human fingertip is its ability to recognize edge orientation. In this way, it was illustrated that the random innervation of the mechanoreceptors by the primary afferents allows the encoding of orientation information through the spatiotemporal spiking pattern. This structure organizes a peripheral neural mechanism for extraction and then transmission of geometric features of the touched objects. The proposed hierarchical spiking neural network successfully discriminated edge orientation stimuli irrespective of edge location. It was shown that using the first spikes of cortical neurons; the orientation of stimuli (scanned or indented edge) was recognizable. The effect of afferent receptive field size was compared in two different experiments (scanned and indented edge). Orientation detection of the scanned edge stimuli in the first spikes of cortical neurons was improved when the afferents' receptive field size was increased. Nevertheless, for the indented edge experiment, the situation was reversed and increasing the size of the afferents receptive field resulted in the reduction of correct detection. The findings showed that the importance of receptive field size depends on the specific tasks and experiments. Recent studies have shown that the main connections in neuronal pathways are formed during the developmental process [38][39][40] . However, the exact cortical dynamics and function have not been studied yet. Here, www.nature.com/scientificreports/ we investigated edge orientation detection through the cortical neurons as a biomimetic classifier. We showed that the intensity of a neuron's response would signal edge orientation because its firing rate would increase with the degree of spatial coincidence between the neuron's highly sensitive zones (excitatory region of receptive field) and the local skin deformations formed by edge indentation. That is, for a given neuron, some edge orientations exhibit more spatial coincidence than others and thus stronger responses are produced. Also, the role of the inhibitory current which forms the lateral inhibition within the cuneate nucleus was studied. Indeed, the simulation results suggest that when lateral inhibition has increased the process of spike filtering is amplified. This leads to the reduction in "noise" within the system and hence the third-order neurons are activated by a strong and consistent signal. This also increases the spatial resolution of the receptive fields and gives them a more distinct border which improves discrimination between two separate points of simultaneous stimulation. Although other forms of lateral inhibition are also observed, the "feedforward" type of lateral inhibition is likely the most significant 41 . Various aspects of tactile sensitivity have been related to different forms of neuronal inhibitory function. Impaired reactions to tactile stimuli in children with autism spectrum disorder (ASD) are frequently reported symptoms. Indeed, impairments in filtering of or adaptation to tactile inputs have been described in ASD 42 . Under the assumption that the inhibitory mechanism is altered in ASD 43,44 , it can be suggested that dysfunction in lateral inhibition of the second layer of tactile processing or malfunction in the formation of the inhibitory sub-regions of the cortical neurons may also have a role. Understanding the specific mechanisms underlying sensory symptoms in ASD is still under investigation which may allow for more specific therapeutic approaches in the future.
The main limitation of the proposed spiking model is the lack of neural recordings for all network layers. Nevertheless, the model is based on the significant literature and published data for model building and validation. The proposed spiking network for a tactile system can be employed in the design and implementation of sensory neuroprostheses applications [45][46][47][48] . Additionally, the broad significance of this work is that the biomimetic tactile sensing and edge encoding are useful in robotic applications for shape recognition and object grasping and palpation [49][50][51] .

Methods
Mechanoreceptor and afferent models. To simulate the mechanoreceptor grid, the physiological information from the literature was considered. The SA-I fibers constitute about 25% of all tactile fibers innervating the hand. They are divided into multiple branches near the skin surface to reach the Merkel discs, which are distributed over different fingerprint ridges. This branching and innervation create a receptive field covering an area of about 10 mm 2 on average for each SA-I afferent 8 . Afferent parameters are reported in Table 2. In response to a www.nature.com/scientificreports/ skin indentation, the activated Merkel discs within the receptive field send their partial information via the SA-I fibers which produce sustained spiking responses that slowly decrease over time (Fig. 1E). The RA-I afferents form around 40% of all tactile fibers innervating the hand. They are split into branches to relay signals from multiple Meissner corpuscles. Typically, the RA-I receptive field size is a bit larger than the SA-I receptive field. They respond to mechanoreceptor gird deformation 52,53 , therefore when the skin does not move, no spike is generated.
For modeling and simulation of the tactile afferents, we considered a grid of mechanoreceptors (80*80) as a patch of skin in which individual mechanoreceptors are located at 0.15 mm intervals 5 . A population of 296 firstorder neurons were modeled according to biological observations and were consisted of 100 SA-I and 196 RA-I with overlapping innervation territories. Each primary afferent innervates 28 mechanoreceptors on average 11 and is randomly weighted between 0.1 and 1. This branching leads to having complex receptive fields with multiple hotspots. To create the receptive field of SA-I afferents, we mapped a 10*10 network (100 SA-I neurons) on the 80*80 mechanoreceptor grid to cover all areas. The center of each section from the 10*10 mapped network is considered as a middle point of an individual receptive field. Then, 28 mechanoreceptors were chosen randomly using a Gaussian distribution in the x-axis and y-axis ( σ d in Table 2) and randomly weighted between 0.1 and 1. For RA-I afferents, we mapped a 14*14 network (196 RA-I neurons) on the 80*80 mechanoreceptor grid. The RA-I receptive field was generated similarly to the SA-I afferents. In total, 100 and 196 matrices with dimensions of 80 by 80 for SA-I and RA-I afferents were created, respectively.
The applied force distribution over the skin patch is modeled using a two-dimensional Gaussian function, and thus the output of the individual mechanoreceptor at location x.y is shown as follows 54 : where F is amplitude, (x 0 , y 0 ) is center of pressure, and f x, y is the perceived force for each mechanoreceptor at the location x, y . σ x = σ y = 0.05 are standard deviation. The output of Eq. (1) is shown in Fig. S2A.
The Izhikevich model was used to reproduce the adaptation dynamics of the mechanoreceptors 55 . Indeed, this neuron model is easily able to generate different spiking patterns by simple parameter adjustment 56 . For the Izhikevich model, the membrane potential, v, and the adaptation variable, u, were updated via the following differential equations which were discretized using Euler's method with the discretizing step = 0.1ms: The values of the parameters a, b, c, and d are reported in Table 3. In Eq. (2), for the SA-I model k SA = 1andk RA = 0 and for the RA-I model k SA = 0andk RA = 1.
For the static SA-I model: For the dynamic SA-I model: ( Table 2. Afferent properties used in the simulations 8 .

Afferents
Receptive field size ( mm 2 ) Density/cm 2 Number of afferents on the skin ( 12mm * 12mm) σ d (mm)  where τ r_SA and τ d_SA are the rise and decay time of the proposed SA-I model, respectively. x is an auxiliary variable that is used in the proposed SA-I model.I in is the input current to the network and, τ RA is the time constant for the RA-I model.
Considering neurophysiological observation, the proposed spiking neuronal network of CN is comprised of two subpopulations of PNs and INs. Individual sub-populations are modeled by a random recurrent connection of spiking neurons with a connection probability of 0.2. Primary afferents (PA) make excitatory connections to both PN and IN spiking neuronal networks. Specifically, if a primary afferent excites a PN, neighboring afferents will excite INs to inhibit the activated PN (Fig. 5A). In this way, the concept of lateral inhibition is modeled. The simulated network model has 296 sensory input channels (PAs) that innervate each CN (PN and IN subpopulation). It should be pointed out that this number of afferents is less than the biological observation, which is around 1000s of PAs per CN 57 . To make the decoding problem more biologically relevant, in the simulations, we added intrinsic noise in the simulated pathways and considered the conduction velocities based on the synaptic delay for each layer.
Relying on prior neurophysiological findings, the local circuit in the somatosensory cortex area 3b was mod- where I AMPA andI GABA are the excitatory and inhibitory synaptic currents received by the post-synaptic neuron, respectively. t * ex and t * in are the spike time received from pre-synaptic to the post-synaptic neuron. The latency of the post-synaptic currents is τ L = 1ms 58 . Synaptic constant times are τ dA = 0.4ms , τ rA = 2ms , τ dG = 0.25ms , τ rG = 5ms 58 . In our simulations, these values are multiplied by 7. W k is the synaptic weight of the connection between a pair of neurons.
Cortical receptive fields. Cortical neurons have specific receptive fields that are feature selective. The cortical receptive fields consist of complex spatial excitatory and inhibitory sub-regions on the skin. To create a receptive field for individual PYs, the synaptic plasticity of the selected neuron is activated to strengthen excitatory and inhibitory connections related to the tuned orientation. After that, synaptic weights are fixed for all experiments. For each neuron in the cortical layer, the cumulative spike counts of all afferents are weighted based on their receptive fields and then are summed to determine the cortical neuron response 35 . We allocated a group of 10 neurons per orientation in the cortical neuronal network and their receptive fields were distributed across the skin. Since we simulated 16 different angles, 160 spiking neurons were used for modeling the cortical area. When the excitatory (inhibitory) sub-region is activated on the skin, PY connections become stronger (c-IN inhibits the PYs).
Biomimetic decoder. The biomimetic decoder mimics the behavior of neuronal circuits in the somatosensory cortex 17 . The decoded orientation corresponds to the highest firing neuronal group. The biomimetic decoder gives a scalar value by counting the spikes of each cortical neurons within the neuronal group. www.nature.com/scientificreports/ Experiment simulation. Simulated edges are indented at the center of the patch of skin for 300ms , including a 50ms on-ramp, a 200ms hold phase, and a 50ms off-ramp. Equation (13) simulates an edge to define specific orientation which is determined by parameter θ . The simulated edge vertically indents the simulated patch of skin (80*80 mechanoreceptor grid): xandy are coordination of mesh grid (80*80) and each axis is between − 1 and 1. w = 0.05 is the edge width and finally, θ is the orientation of the simulated edge. The output of Eq. (13) is shown in Fig. S2B for θ = 45 • . To simulate the dynamics of a scanned edge, the speed of moving edge across the patch of skin is 24mm/s (12 mm in 500 ms).
Noise. To improve the biological plausibility of the spiking neural model, a random fluctuation in the membrane potential is added (as intrinsic noise) to the neuron responses 3 . In this way, background activity is included in the population of neurons for each layer. The random noise is specified by a stochastic differential equation using a Gaussian random variable with mean 0 and standard deviation 1. Furthermore, to mimic the small variations in the presentation of stimulus and movement of skin in different trials, jitter in the depth (1 mm ± 0.5 mm SD) has been included. The numbers are used from 3 .
Spike timing analysis at the population level. To address how much spike timing is informative to decode edge orientation, we used the spike train distance metric as defined by VPd which fully has been explained in 59 .

PCA and classification.
Principal component analysis 60 is used to reduce the dimensionality of a feature space (using sklearn in Python). The feature space of simulated neural responses has high dimensionality and therefore is impossible to visualize. PCA gives the best possible representation of a p-dimensional dataset in z dimensions (z < p) by maximizing variance in z dimensions. In our simulations z = 3. Spiking responses are categorized across different time windows, starting from the first spike. For classification of the orientations, a KNN classifier with k = 5 was used and decoding performance was evaluated using fivefold cross-validation. In this way, the process was repeated five times, and each of the five sets was used as the validation set once.