Significance When we grasp an object, thousands of tactile nerve fibers become activated and inform us about its physical properties (e.g., shape, size, and texture). Although the properties of individual fibers have been described, our understanding of how object information is encoded in populations of fibers remains primitive. To fill this gap, we have developed a simulation of tactile fibers that incorporates much of what is known about skin mechanics and tactile nerve fibers. We show that simulated fibers match biological ones across a wide range of conditions sampled from the literature. We then show how this simulation can reveal previously unknown ways in which populations of nerve fibers cooperate to convey sensory information and discuss the implications for bionic hands.

Abstract When we grasp and manipulate an object, populations of tactile nerve fibers become activated and convey information about the shape, size, and texture of the object and its motion across the skin. The response properties of tactile fibers have been extensively characterized in single-unit recordings, yielding important insights into how individual fibers encode tactile information. A recurring finding in this extensive body of work is that stimulus information is distributed over many fibers. However, our understanding of population-level representations remains primitive. To fill this gap, we have developed a model to simulate the responses of all tactile fibers innervating the glabrous skin of the hand to any spatiotemporal stimulus applied to the skin. The model first reconstructs the stresses experienced by mechanoreceptors when the skin is deformed and then simulates the spiking response that would be produced in the nerve fiber innervating that receptor. By simulating skin deformations across the palmar surface of the hand and tiling it with receptors at their known densities, we reconstruct the responses of entire populations of nerve fibers. We show that the simulated responses closely match their measured counterparts, down to the precise timing of the evoked spikes, across a wide variety of experimental conditions sampled from the literature. We then conduct three virtual experiments to illustrate how the simulation can provide powerful insights into population coding in touch. Finally, we discuss how the model provides a means to establish naturalistic artificial touch in bionic hands.

The human hand is endowed with thousands of mechanoreceptors of different types distributed across the skin, each innervated by one or more large myelinated nerve fibers (1). These fibers convey detailed information about contact events and provide us with an exquisite sensitivity to the form and surface properties of grasped objects (2, 3). During object manipulation and tactile exploration, the glabrous skin of the hand undergoes complex spatiotemporal mechanical deformations, which in turn, drive very precise spiking responses in individual afferents. Coarse object features, such as edges and corners, are reflected in spatial patterns of activation in slowly adapting type I (SA1) and rapidly adapting (RA) fibers, which are densely packed in the fingertip (3, 4). At the same time, interactions with objects and surfaces elicit high-frequency, low-amplitude surface waves that propagate across the skin of the finger and palm and excite vibration-sensitive Pacinian (PC) afferents all over the hand (5⇓⇓–8).

Recording the activity of tactile nerve fibers in monkeys or humans is technically difficult, is slow, and generally yields responses from a single unit at a time (9, 10). Although such recordings have yielded powerful insights into the neural basis of touch, they provide a limited window into the information that the hand conveys to the brain, which is distributed over thousands of responding fibers.

To fill this gap, we have developed a model with which we can simulate the responses of all mechanoreceptive afferents that innervate the palmar surface of the hand to arbitrary spatiotemporal patterns of skin stimulation, taking into account skin biomechanics and receptor biophysics. Model parameters are derived from spiking data obtained from monkeys and validated by reproducing both the strength of and temporal patterning in the responses of afferents to a wide range of stimuli, measured independently in monkeys and humans by several research groups, including our own. With the model, we can simulate in real time the responses of hundreds of afferents to stimuli of arbitrary complexity. We anticipate that the model will be an important tool in somatosensory research to characterize the peripheral representation of tactile stimuli. The model will also be useful in providing somatosensory feedback through interfaces with the peripheral nerve for use in neuroprosthetic devices by converting the output of touch sensors on the prosthesis into biomimetic afferent responses, which can then be implemented through electrical stimulation (11⇓–13).

Discussion Limitations of the Model. Although the model faithfully reproduces key response properties of tactile afferents, it is also subject to a number of limitations. First, we are able to simulate SA1, RA, and PC but not slowly adapting type II (SA2) fibers. Indeed, the models are based on electrophysiological data obtained from rhesus macaques, which are devoid of SA2 afferents (43). Second, the models were derived from responses to stimuli that were indented into the skin and therefore, do not incorporate lateral sliding and the concomitant shear forces. Rather, scanning is mimicked by moving an “image” of the indentation pattern across the skin frame by frame. Note, however, that tangential forces are often highly correlated with normal ones during sliding (44), and therefore, this approximation is sufficient under most circumstances. The model also does not incorporate the effect of fingerprints on texture responses (7) or the onset of slip (45). Tangential sliding and its onset are governed by the friction between the skin and object surface, which is a complex mechanical problem that is not yet fully understood (46). As our understanding of how tangential forces shape afferent responses improves, this aspect can be incorporated into the model. Third, the skin mechanics model treats the skin as a flat surface, when in reality, it is not. The 3D shape of the skin matters during large deformations of the fingertip. For example, pressing the fingerpad on a flat surface causes the skin on the side of the fingertip to bulge out, which in turn, causes receptors located there to respond (47, 48). Such complicated mechanical effects can be replicated using finite element mechanical models (49) but not using the continuum mechanics (CM) model adopted here. To the extent that friction is a critical feature of a stimulus—for example, when sliding a finger across a smooth, sticky surface—or that the finger geometry plays a critical role in the interaction between skin and stimulus—as in the example of high-force loading described above—the accuracy is compromised. Under most circumstances, the model will capture the essential elements of the nerves’ response. Fourth, SA1 and RA afferents have been shown to exhibit complex RFs with multiple hotspots (24, 50), because single afferents innervate multiple mechanoreceptors (43, 51). Afferent branching is not currently implemented in the model; rather, RFs comprise a single hotspot and are isotropic. However, the model can be readily elaborated to accommodate afferent branching, particularly when we better understand the mechanisms by which input from the various mechanoreceptors is integrated to drive the afferent response (52). Applications. The proposed model can be a powerful tool to investigate the sense of touch. Indeed, recording the responses of human or monkey afferents is technically challenging and only yields responses from a single fiber at a time. Even multielectrode arrays only yield responses from a sparse sample of afferents (53) or the aggregate activity of a large number of fibers (54, 55). The model allows us to simulate the responses of entire populations of afferents to arbitrarily complex stimuli. These simulated responses can be used to (i) test peripheral neural codes (as illustrated above with the two simulated experiments), (ii) assess the degree to which the results from psychophysical experiments can be accounted for based on population response (56), and (ii) investigate how sensory representations are transformed as one ascends the somatosensory neuraxis (cf. ref. 57). Also, the model will be valuable in neuroprosthetic applications that aim to restore the sense of touch through peripheral nerve interfaces (58⇓–60). Indeed, the output of sensors on the prosthesis during object contact can be used as input to the simulated afferents, which then provide an accurate representation of how the nerve from an intact hand would respond to object contact. This biomimetic pattern can then be effected in the residual nerve of an amputee by delivering spatially and temporally patterned stimulation pulse trains designed to evoke the desired nerve activation (11, 12, 17). The implementation is computationally efficient and can run in real time (SI Methods). In the context of neuroprosthetics, the model can also be used as a benchmark against which to compare simpler encoding algorithms and determine the circumstances under which their behavior diverges from that of nerves in an intact arm.

Methods The model of peripheral afferents is structured in two distinct sequential stages to capture the structure of the mechanotransduction process (Fig. 1). In the first stage, we compute the deformations experienced by individual receptors given a stimulus applied to the surface of the skin. In the second stage, we compute based on this receptor deformation the spiking response of the nerve fiber using an IF mechanism. Stimulus. A stimulus is defined as the time-varying positions of a set of pins indented into the skin in the direction orthogonal to the skin surface. Each pin makes contact with the skin over a circular area of adjustable radius and can move up and down independent of the other pins. The location of a pin with respect to where it will touch the hand is fixed over time using a coordinate system centered on the fingertip of the index finger (Fig. S6A). A stimulus can be defined by a single pin (corresponding to a circular probe indented into the skin) or multiple pins, which together form a shape (say an edge, a letter, or a curved surface). Pin indentations are independent of each other and can, therefore, form arbitrary spatiotemporal patterns of indentation. Fig. S6. Skin mechanics model. (A) Schematic of the hand used in the model. The reference coordinate system is centered on the center of the index fingertip. Different hand regions are defined that determine the innervation densities of the different afferent types. (B) Skin deflection profile resulting from the indentation of a circular pin (shown in gray; radius = 2 mm). (C, Upper) Quasistatic and (C, Lower) dynamic components resulting from a 200-Hz vibration exerted with a single pin (radius = 2 mm; blue trace) or a grid of evenly spaced pins (spacing = 0.1 mm; radius = 0.05 mm). Receptor depth was set to 0.5 mm. Skin Mechanics. The skin of the fingers and palm is known to have complex, nonlinear mechanics (61⇓–63), with properties that vary widely over different spatial and temporal scales (64, 65). Modeling the precise skin deformation resulting from a given stimulus requires advanced models (66, 67), measurements of individual skin properties (68), many parameters, and considerable computing power. Because we aimed to simulate the responses of whole populations of afferents in quasireal time, we did not attempt to estimate the exact deformation state of the finger but instead, simplified the problem by evaluating mechanical quantities that have been shown to be highly predictive of afferent responses. Accordingly, we focused on two quantities that are strongly associated with neural responses but over different frequency ranges: (i) the local vertical stress based on a quasistatic elastic model of the skin, mainly affecting receptor responses at low frequencies (<100 Hz) (16) and (ii) dynamic variations in pressure propagated across the skin as surface waves, mainly affecting receptor responses at higher frequencies (5, 6). These two quantities are then combined in the spiking model and differentially weighted depending on afferent type. For instance, SA1 afferents are most sensitive to statically indented spatial patterns, and therefore, the quasistatic component is strongly weighted in SA1 models. Conversely, PC afferents do not respond to static indentations but are extremely sensitive to high-frequency vibrations even far away from the center of their RF; therefore, PC models will feature a strong dynamic component. Quasistatic component. The quasistatic component of the model is inspired by an existing CM model of afferent responses (15, 16). CM offers two important benefits: (i) it provides analytical solutions and is, therefore, fast to compute, and (ii) it provides very accurate predictions of the responses of type I afferents (SA1 and RA) to spatial patterns indented into the skin (16). In the CM model, the skin is assumed to be a homogeneous, isotropic, elastic half-space, and the stimulus (i.e., the spatiotemporal pattern of indentations) is applied at the free border of this half-space. The surface of the object is considered to be frictionless and therefore, produce load only in the direction perpendicular to the skin. In the original version of the model, each stimulus pin was described as punctate (that is, with infinitely small diameter). We modified this by defining each pin as a circular punch. This modification has the benefit that large circular pins indenting the skin (i.e., circular probes as commonly used in tactile experiments) can be modeled using a single pin rather than a large number of them, speeding up computation considerably. The deflection produced by a single pin is calculated as follows (Fig. S6B) (69): u ( r ) = f ( r ) ⋅ p , [1] f = { 1 k ( r < r p ) 2 π k sin − 1 r p r ( r > r p ) , [2] k = 2 r p E 1 − ν 2 , [3]where r is the distance from the center of the pin, u is the vertical skin deflection at distance r, p is the force acting on the pin, r p is the radius of the pin, k is the stiffness of the skin viewed by this pin, E is Young’s modulus, and ν denotes Poisson’s ratio. This expression was later generalized to multiple pins (15), adopting the following matrix form: 𝐱 = [ x 1 x 2 ⋮ x n ] 𝐲 = [ y 1 y 2 ⋮ y n ] 𝐮 = [ u 1 u 2 ⋮ u n ] , [4] R i j = ( x i − x j ) 2 + ( y i − y j ) 2 , [5] u = f ( 𝐑 ) 𝐩 , [6]where n is the total number of pins, x and y are the coordinates of each pin, u is a vector containing the indentation depth of each pin, and p is a vector containing the force exerted by each pin. R is a matrix containing the Euclidian distances between pins ( R i j is the distance between pin i and pin j). The stresses acting on each receptor are then obtained in a two-step procedure. First, the indentation depth profile (i.e., the indentation depth of each pin, which constitutes the input of the model) is converted to an indentation force profile (i.e., the force applied by each pin to the surface of the skin). This conversion is achieved by solving the linear system presented in Eq. 6 for p. Because pins can only push into the skin but not pull on it, all pins exerting a negative force (pulling) are removed (because they would not actually be in contact with the skin). The subsystem is then iteratively resolved until all remaining pins exert a positive force (i.e., pushing into the skin) as done in earlier models (16). Second, the stresses acting on a receptor are obtained independently for each pin, and then, their effects are summed by applying the superposition principle. This calculation is made by using the analytical expression of the stresses resulting from the indentation of a circular pin by a force p in the elastic half-space as obtained by Harding and Sneddon (69). The analytical expressions are reported in SI Methods. Note that these expressions assume a particular pressure distribution, and this assumption is not valid when two pins are very close to each other. Nevertheless, the error is negligible if the pin radius is sufficiently small compared with the afferent depth (Fig. S6C). The afferent depth was set according to the literature (51) and is the same for all afferent models (SA1 depth is 0.3 mm, RA is 0.2 mm, and PC is 2 mm). Most stress components have been shown to be highly correlated with afferent responses (16). Here, we chose the vertical component of the stress tensor (i.e., the stress perpendicular to the skin surface), because this component was previously shown to be a very good predictor of both SA1 and RA responses, and it is consistent with the dynamic component described below. Dynamic component. High-frequency vibrations applied to the fingertip, elicited, for example, during making and breaking contact with objects or scanning the finger across a textured surface, travel across the finger and palm and have been recorded as far away as the wrist (5, 6, 8). In the model, we wished to incorporate such surface waves resulting from pin movement. First, the indentation velocity profile is converted into a force variation profile using Eq. 6 as above but after replacing the stiffness coefficient k (Eq. 2) with a viscous coefficient c (c is arbitrarily set to one because of the lack of consistent measurement; therefore, the dynamic component has no unit). Second, the force variation of each pin is propagated to the receptor as a surface wave, with a decay of 1/r (5) and a propagation speed (group velocity) of 8 m/s (as measured in ref. 5 but constant across frequencies for simplicity). As in the quasistatic part, the stress acting on the receptor is obtained by summing the contribution of each pin. IF Model. The output of the skin mechanics model at all receptor locations is fed into IF models that generate spiking responses for each afferent. We fit the parameters of these models on electrophysiological recordings obtained previously (18). These models were similar to earlier ones (57, 70, 71), with some key differences. First, rather than simply using skin indentation depth (along with derivatives) as IF model input, we used the static and dynamic pressure as calculated from the skin mechanics model. This choice of input allowed us to reproduce response properties driven by skin mechanics, such as edge enhancement and surround suppression (Results), as well as effects caused by the propagation of dynamic pressure waves. Second, the models were fit on a wide variety of vibrotactile frequencies, ranging from 1 to 1,000 Hz, in contrast to earlier models, which were fit on much narrower stimulus sets (i.e., either low or high frequencies). Each IF model works as follows (Fig. 1D). The mechanical input (both quasistatic and dynamic) is passed through a low-pass filter (parameter 1) (Fig. 1D) to accommodate the fact that afferents become nonresponsive above a certain (afferent type-specific) stimulation frequency. Then, the dynamic component is differentiated to obtain three signals (quasistatic, dynamic, and differentiated dynamic). Those three inputs can be interpreted as the three terms of a “lumped” mass–spring–damper model, with the quasistatic term being associated with the elastic stress, the dynamic term being associated with the damping, and the differentiated dynamic term being associated with the mass. Next, the three signals are split into positive and negative contributions and rectified, yielding six time-varying inputs (Fig. 1D), which are each multiplied by a weight (parameters 2–7) and summed. Then, the resulting time-varying trace is passed through a saturating nonlinearity, reflecting the fact that afferents’ responses saturate at high intensities (parameter 8). Next, Gaussian noise is added (parameter 9) to mimic the observed stochasticity in afferent responses, which is particularly evident for perithreshold stimuli (72). The resulting trace constitutes the input to a leaky IF model, with “membrane potential” that decays to zero with a characteristic time constant (parameter 10). When the potential hits one, a spike is triggered, the membrane potential is reset, and a postspike inhibitory kernel is added to the potential to build in refractoriness (parameters 11 and 12). The kernel consists of two parts: a fast component, which decays completely within 4 ms, and a slow component, which contributes maximally 8 ms after a spike is triggered and decays within 36 ms (a similar model is shown in ref. 71). In all fitted models, the first component was weighted more heavily than the second component, leading to initially strong but lingering weak inhibition after each spike. Finally, spikes are shifted in time by a small amount to mimic conduction delays (parameter 13), which are largely determined by where the recording was made (proximal vs. distal). An individual IF model, therefore, possesses 13 parameters: the low-pass cutoff frequency, the six weights, the saturation parameter, the noise term, the membrane leak time constant, two parameters to determine postspike inhibition, and the conduction delay. Not all parameters were required for all afferent models. SA1 models did not include a weight for the derivative of the dynamic pressure, because these afferents do not respond to high frequencies. SA1 models also did not include a saturation parameter. RA and PC models did not include weights for the static stress component, because these afferent types do not respond to constant indentations. Thus, SA1 models used 10 parameters, whereas RA and PC models were fit using 11 parameters each. Fitting Procedure. For model fitting, we used as a cost function the van Rossum spike distance (73), which yielded a measure of difference between the recorded and model-predicted spiking responses at a given temporal resolution (set by a time constant). We then optimized the model parameters (apart from the noise term, parameter 9; see below) using the patternsearch function in Matlab (The Mathworks, Inc.) using different starting positions. We used a time constant of 50 ms for the computation of the cost for SA1 and RA afferents and 20 ms for PC afferents. The parameters were both fit on responses to sinusoidal stimuli and bandpass noise traces (18). PC models were fit on all frequencies used in the stimulus set (ranging from 1 to 1,000 Hz). SA1 and RA models were only fit on frequencies between 1 and 150 Hz. After we had zeroed in on good parameter values, we further optimized the parameters using genetic algorithms (ga in Matlab) until they did not improve further. To fit the noise parameter, we first calculated the average spike distance between all pairwise combinations of the five recorded trials for each stimulus in the training data. We then adjusted the noise parameter, such that five simulated runs of the model on the same inputs resulted in the same average pairwise spike distance. In total, we fit single-afferent models to four SA1, nine RA, and four PC afferents. Model Architecture Validation. Our goal was to create a model that could recreate most afferent response properties described in the literature but that was simple enough to allow for fast computation. Furthermore, we wished to avoid overfitting and ensure that the parameters contributed to prediction accuracy. We, therefore, tested a number of simpler models, and we examined which response properties they could successfully recreate and where they failed. Additionally, these tests provide insights about which specific afferent response properties are driven by which model component. Skin mechanics. First, as shown previously in a number of studies (15, 16, 33, 35), the spatial profile of afferent response strength is not accurately predicted by the profile of the skin indentation. Indeed, as predicted by contact mechanics, the distribution of pressure (and therefore, the distribution of stresses inside the skin) is influenced, for example, by edges and the size of the contact area. A contact mechanics component is, therefore, required to incorporate all spatial response properties, such as edge enhancement, surround suppression, and RF sizes. Second, PC fibers are known to respond to vibration far away from their locations, which requires the wave propagation (dynamic) component of the model. However, incorporating skin mechanics is not necessary to achieve high spiking precision (for stimuli located directly above the receptor location), because previous models were able to reproduce precise spiking patterns evoked by vibrating probes without implementing a skin mechanics model (70, 71). Frequency selectivity of different afferent types. Selective weighting of the quasistatic, dynamic, and dynamic derivative inputs will yield models that differ in their frequency response and could, in principle, account for the SA1, RA, and PC frequency response profiles. To test whether a simple model using these inputs might be sufficient to model afferent response properties or whether a more complex spike initiation model is required, we implemented a linear–nonlinear–Poisson (LNP) model (SI Methods has details on the fitting process). This model assumes that spikes are generated by a nonhomogeneous Poisson process and does not include leakage, a spiking threshold, or refractoriness. We found that the LNP model was indeed able to fit and predict firing rates of the vibration datasets, although less accurately than the IF models [training data: R 2 = 0.66 ± 0.14 (SD) and 0.84 ± 0.14 for sinusoidal and noise vibrations, respectively; validation data: R 2 = 0.73 ± 0.16 for diharmonic vibrations] (Results has details on the IF fits). However, the precise spike timing was consistently much worse as shown using metric–space analysis of spikes trains at different timescales (from 0.1 to 100 ms) (Fig. S7 A and B). In particular, we found that the timing precision of the IF models surpassed that of the LNP models at the timescales that have been shown to be particularly informative for the different afferent classes (10 ms for SA1, 5 ms for RA, and 1 ms for PC) (19, 20, 24). An example of this timing precision is the ability of an afferent to phase-lock to sinusoidal stimuli at different frequencies and mechanical noise with different bandpass frequencies (Fig. 2). Fig. S7. Model validation. (A) Comparison of spike timing precision between the IF and LNP models. Shown are raster plots of spikes evoked in the three afferent classes by five repetitions of one of the diharmonic stimuli (Left, SA1; Center, RA; Right, PC). Both measured (black) and simulated spikes (IF in blue and LNP in green) are shown. (B) Normalized spike distance between measured and simulated spike trains as a function of timescale (1/ q ). IF actual distance is in blue; LNP actual distance is in green. Shaded lines show the average distance across all conditions in the diharmonic vibration dataset for each afferent model. Dark lines show averages across afferents. Red traces (DIFF) show the difference between predictions of the two models: the IF model is better than LNP at all timescales but especially at the timescales that have been shown to be most informative for the different afferent classes (10 ms for SA1, 5 ms for RA, and 1 ms for PC). (C–E) Reduced model fits. Removing or clamping IF model parameters during the fitting process leads to impaired accuracy in reproducing well-known afferents response properties as illustrated in C–E. (C) Missing entrainment plateaus for a PC model at three different frequencies when the model was refit without postspike inhibition (parameters 11 and 12). (D) Absolute thresholds that were roughly one-half an order of magnitude too low when a PC model was refit with the decay variable clamped. (E) Excessive firing rates when a PC model was fit without a saturation term. Individual IF parameters. After verifying that the IF model yields more accurate predictions than a simpler spike generation model does, we next evaluated the contribution of individual IF model parameters. Indeed, some of these parameters (for example, the low-pass cutoff, the saturation, and the membrane leak time constant) affect model responses in similar and overlapping ways. To test whether these parameters were necessary to reproduce precise spiking, we fit models that either lacked one of these parameters or had the parameter “clamped” to a certain value during the fitting process (SI Methods). We found that removing any of these parameters considerably reduced the model’s ability to reproduce key response properties; the most notable problems were missing entrainment plateaus (Fig. S7C shows a model with no postspike inhibition), thresholds off by an order of magnitude (Fig. S7D, clamped membrane leak time constant), and excessively high firing rates (Fig. S7E, no saturation). These effects can be explained by the fact that we are fitting complex, frequency-dependent response properties, and precisely matching those properties requires several parameters. Spike Timing. Metric space analysis. To measure the degree to which simulated responses matched their measured counterparts, we used a spike distance metric with an adjustable temporal resolution (74), which we have previously done (19, 20). In this metric, the distance between two spike trains is obtained by computing the lowest possible cost for transforming one train into the other. Adding or deleting a spike incurs a cost of one, whereas moving a spike costs q per unit time. Therefore, by varying the parameter q, we can assess the degree to which the timing of the spikes matches at different levels of temporal precision. Specifically, when q is set to zero, distance is determined solely by the difference in spike count. At higher values of q, moving spikes becomes more costly than adding or removing them, and therefore, the timing of individual spikes increasingly determines spike distance, which remains low if timing is precise. Therefore, the distance increases with q and converges to the sum of spike count in the two spike trains when q tends to infinity (which is equivalent to removing all spikes from the first train and then inserting the spikes from the second train). In all of our timing analyses, the distances between the measured spike trains and the simulated spike trains were measured at 50 log-spaced timescales ranging from 0.1 ms to 1 s (corresponding to q values ranging from 10,000 to 1). The distances obtained were then normalized by the sum of spike counts in the two spike trains, yielding normalized distances that ranged from zero to one, to make them comparable with each other (Figs. S7B and S8 show examples). Fig. S8. Spike timing precision. Difference in normalized distance between simulated (blue) and jittered (green) spike trains to measured spike trains for three afferent classes as a function of temporal resolution. The pairwise difference between the two distances is shown in red. Each row corresponds to a different jitter SD value (1, 2, 5, and 10 ms). Shaded lines show means for individual afferents, and dark lines show means across afferents. The orange vertical lines correspond to the mean distance differences caused by difference in spike count. The red line being below the orange line indicates that the simulated responses are closer to the measured response than the jittered ones. Timing precision of the model. For each experimental condition (from a total of 138 different diharmonic stimuli) and each afferent, pairs of spike trains were used to obtain a distance at every time scale. Pairs consisted of either measured and simulated spike trains or measured responses from two different repetitions with temporal jitter added to the second repetition. Jitter was sampled randomly from a zero-mean Gaussian distribution with a given SD and added to each spike. We tested SDs ranging from 0.5 to 50 ms. The resulting two distances, measured/modeled (shown in blue in Fig. S8) and measured/measured + jittered (shown in green in Fig. S8), were then subtracted to yield a measure of which comparison spike train is closer to the measured data (shown in red in Fig. S8). A positive value indicates that jittered trains are closer, and a negative value indicates that simulated trains are closer. This difference in spike distance was averaged across all timescales above the jitter SD value and computed for each jitter SD value (shown in Fig. 2E). We used a t test to determine whether this averaged difference (corrected for the spike count error) (shown in orange in Fig. S8) was statistically greater than zero, meaning that the model was worse (one-tailed t test, alpha = 0.05). The lowest jitter value for which no significant effect was observed was defined as the precision of the model. Simulated Experiments. Procedures. Edge indentation. Simulated edges (8-mm length and 1.6-mm width) were indented at the center of the index fingertip to a depth of 1 mm for 400 ms, including a 50-ms on ramp, a 300-ms hold phase, and a 50-ms off ramp; 32 equally spaced orientations were tested, each repeated 10 times. In addition to the intrinsic noise in the afferent spiking model, we added experimental variability in the form of Gaussian noise in the exact stimulus position (isotropic with 0.5 mm SD) to mimic jitter in the stimulus presentation as well as small movements of the fingertip. Texture scanning. Three different texture samples were selected from a previously used dataset (7, 20): one fine texture (wool gabardine), one coarse texture (upholstery), and a dot pattern (2-mm interdot spacing). Texture profiles (obtained by profilometry) were used as inputs to the model (with surfaces approximated as rigid). The contact area was defined as circular (with a radius of 4 mm), and the resolution (pin spacing) was set to 0.1 mm. The skin contact area moved across the texture at a speed of 20 mm/s for 0.8 s. The indentation depth was set to 1 mm at the center of contact and followed a circular profile toward the border of contact. Sixteen equally spaced directions were tested, each repeated five times. Gaussian noise was added to the stimulus position (isotropic with 0.5 mm SD) to mimic experimental noise. Data analysis. We used linear discriminant analysis (LDA) to classify orientation or scanning direction (using fitcdiscr in Matlab) using a leave one out cross-validation procedure. We used a spike distance measure (as described above) to compare spike trains at different timescales. The distances between spike trains of different repetition and different orientation were computed for each afferent. The distances were then summed across afferents (by taking the root of the sum of square distances). The final distance matrix was then converted back to Euclidean coordinates. Those coordinates were used as feature vector for the LDA-based classification. The classification analyses were performed over different response time windows, starting after the first spike elicited across the whole population and ending at different times (logarithmically spaced) (Figs. 6C and 7E, x axis tick marks). Two different temporal resolutions were tested (q = 0 and q = 500 s − 1 ).

SI Methods LNP Model. To test whether an IF mechanism was required to reproduce afferent responses, we implemented LNP models and compared them with the full model. The six components from the IF model (rectified quasistatic, dynamic, and dynamic derivative) served as inputs to these models. The same procedure was followed as that described for the IF model (in the text): models were fitted on two datasets (sines plus noise) and then tested on a third (diharmonic). We used a standard sequential method to fit the models. First, a set of six linear filters was obtained using the spike-triggered average computed over a 25-ms window. Second, the linear estimate obtained by filtering and summing six components was passed through a nonlinearity to best fit the actual firing rate (obtained by convolving the spike trains with a 40-ms Gaussian kernel). We found that a nonlinearity of the form a ⌊ x − b ⌋ c (with ⌊ ⌋ denoting rectification) gave the best results. Third, Poisson spikes were generated and then binned at 5 kHz. We used the same metric space analysis as described in Methods to compare the distance between the measured spike trains and the simulated spike train (IF or LNP) at different timescales ranging from 0.1 ms to 1 s. The results are shown in Fig. S7 A and B. Validation of Individual IF Parameters. To avoid overfitting and ensure robust parameter estimates, we verified that individual IF model parameters—specifically, the low-pass cutoff (parameter 1), the saturation level (parameter 8), the time constant/decay (parameter 10), and the postspike inhibition (parameters 11 and 12)—were necessary to improve model accuracy. We removed or clamped individual parameters and then refit the IF models. For this analysis, we focused on PC afferents, because these fibers have the widest frequency sensitivity. The results are shown in Fig. S7 C–E. Benchmarks. The model is computationally efficient enough that a simplified version can be run in real time. Fast runtime is especially important in neuroprosthetic applications, where tactile feedback from a sensorized prosthetic hand would need to be converted into realistic spike trains in quasireal time to be useful to the user. How many pins (defining the complexity of the stimulus) and how many afferents (defining the complexity of the feedback) the model can support in real time depends on the chosen sampling rate for the afferent models. On an Intel i5 quad core running at 3.50 GHz using Matlab 2015b, the model runs in real time with 30 pins and 30 receptors at a sampling rate of 10 kHz, 100 pins and 100 receptors at a sampling rate of 3 kHz, and 300 pins and 300 receptors at a rate of 300 Hz. The model can thus support hundreds of pins and receptors at relatively high temporal resolution in real time. Stresses at the Receptor Caused by a Circular Probe. This section adapts equations from ref. 69. Definitions. r p is the probe radius. P is the vertical force applied by the probe. E is the Young modulus (set to E = 50 kPa), and ν is Poisson’s ratio (set to ν = 0.4 ). x , y are coordinates on the skin surface, and z is the depth of the receptors. Stresses at the receptors a caused by probe p. • Geometry: x = x a − x p y = y a − y p z = z a r = x 2 + y 2 θ = tan − 1 ( y x ) .

• Definitions: ξ = z r p ρ = r r p R = ( ρ 2 + ξ 2 − 1 ) 2 + 4 ξ 2 α = tan − 1 ( 1 ξ ) ϕ = tan − 1 ( 2 ξ ( ρ 2 + ξ 2 − 1 ) ) . • Intermediate integrals: J 1 0 = 1 R sin ( ϕ 2 ) J 2 0 = 1 + ξ 2 R 3 2 sin ( 3 ϕ 2 − α ) J 0 1 = 1 − R ρ sin ( ϕ 2 ) J 1 1 = 1 + ξ 2 ρ R sin ( α − ϕ 2 ) J 2 1 = ρ R 3 2 sin ( 3 ϕ 2 ) .

• Stresses: σ r = P 2 π r p 2 [ J 1 0 − ξ ( J 2 0 − J 1 1 ρ ) − ( 1 − 2 ν ) J 0 1 ρ ] σ r z = P 2 π r p 2 ξ J 2 1 σ θ = P 2 π r p 2 [ 2 ν J 1 0 − ξ J 1 1 ρ + ( 1 − 2 ν ) J 0 1 ρ ] σ z = P 2 π r p 2 ( J 1 0 + ξ J 2 0 ) , with the following if ρ = 0 : J 0 1 ρ = 1 2 ( 1 + ξ 2 ) J 1 1 ρ = ξ ( 1 + ξ 2 ) 2 . • In Cartesian coordinates: σ x = σ r cos 2 ⁡ θ + σ θ sin 2 θ σ y = σ r sin 2 ⁡ θ + σ θ cos 2 θ σ x y = ( σ r − σ θ ) sin ⁡ θ cos θ σ x z = σ r z cos θ σ y z = σ r z sin ⁡ θ .

Acknowledgments We thank Ezra Zigmund for programing help on an earlier version of the model. This work was supported by National Science Foundation Grant IOS 1150209 and the Kimberly Clark Corporation. B.C.R. was supported, in part, by the Conte Center for Computational Neuropsychiatric Genomics (NIH Grant P50MH94267) and the Chicago Biomedical Consortium (L004).

Footnotes Author contributions: H.P.S. and S.J.B. designed research; H.P.S. and B.P.D. performed research; H.P.S., B.P.D., and B.C.R. contributed new reagents/analytic tools; H.P.S. and B.P.D. analyzed data; and H.P.S., B.P.D., and S.J.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1704856114/-/DCSupplemental.