|Home | About | Journals | Submit | Contact Us | Français|
Slowly adapting type I (SA-I) mechanoreceptors encode the edges and curvature of touched objects by generating neural spikes in response to indentation of the skin. Beneath this general input-output relationship, models are of great utility for understanding the sub-processes, as SA-I transduction sites are inaccessible to whole-cell recording. This work develops and validates a SA-I skin-receptor model that combines a finite element model of skin mechanics, a sigmoidal function of transduction, and a leaky integrate-and-fire model of neural dynamics. The model produced a R2=0.80 goodness of fit between predicted and observed firing rates for 3 mm and 5 mm grating stimuli. In addition, modulation indices of predicted firing rates for 3 mm and 5 mm gratings are 0.46 and 0.59 respectively, compared to the 0.71 and 0.72 found in vivo. An analysis of predicted first spikes indicates their latency may also be enhanced by edges, as edge proximity shortened first spike latencies by 26.2 and 41.8 ms for the 3 mm and 5 mm gratings, respectively. The model described here bridges the gap between those models that transform sustained indentation to firing rates and those that transform vibration to spike times.
Slowly adapting type I (SA-I) mechanoreceptors encode curvature and edges, enabling us to grasp and manipulate everyday objects such as buttons and spoons [1–3]. When an object contacts the skin’s surface, mechanical forces propagate through underlying tissue and toward the locations of Merkel cell neurite complexes (MCNC), the transduction sites of SA-I mechanoreceptors. Each SA-I transforms forces local to its transduction site into a series of neural pulses that carry stimulus information along the afferent to the brain. The pulses are stereotyped all or nothing events, referred to as action potentials or spikes.
The series of spikes generated by the SA-I exhibit both dynamic and static phases . Whereas the dynamic phase takes place over the period of indentation ramp-up, the static phase begins once the indentation reaches its end point. The number of spikes generated per second, or firing rate, is greater over the dynamic phase. In either phase, firing rates are greater in response to the indentation of edges as opposed to flat surfaces. This phenomenon, termed edge enhancement, has been observed in vivo.
While these types of in vivo observations match stimuli at the skin with spikes at the SA-I afferent – the gross input-output relationship – other intermediate transformations are not currently observable. For example, we can neither measure the specific forces local to the MCNC nor observe the transformation of those local forces to the timing of elicited spikes. Therefore, researchers have turned to modeling approaches that range from skin mechanics to neural dynamics.
Skin mechanics models, for example, indicate that SA-Is respond more vigorously to edges because of the higher magnitude of stress and strain produced in tissue beneath indented edges. Most skin mechanics models employ either continuum mechanics [4, 5] or finite elements [6, 7]. They transform static indentation or force into a stress or strain quantity (typically strain energy density) and then use a scaling function to obtain the firing rate. These models’ predictions closely match in vivo SA-I firing rates for like stimuli, generally with goodness of fit values between 0.80 and 0.97 [4–7]. However, recent studies contend that changes in the firing rate alone may not fully explain the transmission of tactile information necessary for motor control and perception .
In particular, the timing of individual spikes may enable our rapid reaction time. For example, our motor system responds to features of touched objects within 100 ms, when most afferents have had time to transmit only a single spike . One measure of the neural response derived from that single spike is first spike latency, i.e., the time between stimulus onset and the first spike. First spike latencies and firing rates may carry complementary information, or may carry the same information within different timescales. However, because skin mechanics models produce only firing rate, they are not applicable for investigating first spike latency.
Neural dynamics models, in contrast, afford the prediction of spike times. Examples include the Hodgkin-Huxley model, two-dimensional reductions of the Hodgkin-Huxley model, and leaky integrate-and-fire models [9, 10]. Leaky integrate-and-fire models specialized for vibration have successfully modeled SA-I neural dynamics [11–15]. They use a resistive-capacitive (RC) circuit to calculate the voltage across the SA-I membrane, i.e., membrane potential, as a function of vibration frequency and magnitude. A spike time is recorded when the membrane potential is driven to threshold. The issue is that the input of these specialized leaky integrate-and-fire models is restricted to vibration magnitude and frequency.
This work addresses the problem of sustained indentation by bridging the gap between skin mechanics models that produce only firing rates and neural dynamics models that respond to only vibration. This work describes the design and analysis of a skin-receptor model to predict SA-I spike times elicited by sustained indentation. The model’s ability to predict both firing rate and first spike latency may help in investigating the origination and transmission of edge information, both for rapid behaviors and over longer timescales.
The skin-receptor model developed herein represents skin mechanics with a finite element model, transduction at the SA-I membrane with a sigmoidal function, and neural dynamics with a leaky integrate-and-fire model. The skin-receptor model was fit to published firing rates using response surface methodology and the resulting goodness of fit was evaluated. Additional analysis examined the percentage of first spike latencies predicted during stimulus ramp-up, the modulation of firing rates due to edges, and differences in predicted first spike latencies due to edges.
The skin-receptor model couples sub-models representing skin mechanics, transduction at the SA-I afferent membrane, and neural dynamics.
Skin mechanics transform indentation at the skin’s surface into distributions of stress and strain that propagate through the skin’s layers. To model skin mechanics, a finite element model (FEM) was used to represent a two-dimensional cross-section of fingertip skin, as developed and validated previously . The FEM represents skin microstructures such as the epidermal, dermal, and hypodermal layers of skin as well the undulating interface between the epidermal and dermal layers (Fig. 1). Layers are modeled with a linear elastic material model (Young’s modulus parameters of 136 kPa for the epidermis, 80 kPa for the dermis, and 34 kPa for the hypodermis; Poisson’s ratio of 0.48) [17, 18]. While some models focus upon the skin’s nonlinear and time-dependent properties [19, 20], linear elastic models are sufficient for displacements of 1 mm or less, which amount to less than 10% of the total skin thickness . In addition to skin layers, the FEM represents collagen fibers, fingerprint lines, the nail, and the bone. Boundary conditions are imposed such that nodes along the nail and bone are fully constrained in x and y directions. Indenters that contact the model’s surface are specified as rigid analytic surfaces with a friction coefficient of 0.3 between indenter and skin surfaces . The mesh utilized four-node, bilinear quadrilateral, hybrid with constant pressure elements (ABAQUS type CPE4H). Generalized plane strain elements have been used for fingerpad tissue modeling . The entire mesh contains about 16000 nodes and elements. ABAQUS Standard, version 6.6 was used to compute the deformation at the skin’s surface and the propagation of stress and strain through the layers.
From the resultant distributions of stress and strain, the quantity of strain energy density (SED) is averaged from six elements in the epidermal layer at the trough of one intermediate ridge (Fig 1, Inset), 0.95 mm below the skin’s surface and centered on the model’s x-axis. Although the relative contribution of different measures of stress, strain, and displacement to SA-I receptor responses is not fully known, SED is used here as the proxy measure. Our previous investigations with this same finite element model have shown a close correlation of SED with SA-I firing rate , as have the investigations of others [4–7, 16]. However, these studies also show that the measures of maximum compressive strain and maximum tensile strain are closely correlated to firing rate.
SED (Uo) (Eqn. 1) for nearly incompressible materials is formulated as the energy produced by distortion, i.e., a change in the unit volume, where τoct is the octahedral shear stress and G is the shear modulus of elasticity .
The octahedral shear stress, τoct, is given by equation 2, where σxx, σyy, and σzz are normal stresses and where τxy, τxz, and τyz are shear stresses.
The shear modulus of elasticity, G, is given by equation 3, where ν is Poisson’s ratio and E is Young’s modulus.
The transduction sub-model transforms SED in the vicinity of the Merkel cell neurite complex to current entering the SA-I afferent through its membrane, a quantity known as the receptor current. While the relationship between SED and receptor current is unknown – because the Merkel cell neurite complex is not accessible to whole-cell recording – sensory cells such as hair cells and pain receptors exhibit stimulus-current curves that are sigmoidal [25, 26]. For this reason, this work employs a sigmoidal function (Eqn. 4) where α, γ, and λ specify the shape of the function, and where Uo is the SED that gives rise to the receptor current I.
The model parameters α, γ, andλ are obtained through model fitting, where the difference between model predictions and in vivo firing rates is minimized (Appendix).
A leaky integrate-and-fire model is employed to transform receptor current to spikes. In essence, whenever the membrane potential, u(t), is driven to the spike initiation threshold, , by the receptor current, I, a spike time is recorded and a 1.0 ms refractory period is entered . Algorithmically (Fig. 2), the neural dynamics sub-model simulates time in 0.1 ms increments. For each time increment, the neural dynamics sub-model is either in a refractory period (Fig. 2 “refractory”) and does nothing, or is not in a refractory period (Fig. 2 “!refractory”) and calculates the membrane potential.
To calculate membrane potential in response to receptor current, the SA-I membrane is modeled as a resistive capacitive (RC) circuit. The RC circuit, written in mathematical form in equation 5, defines the dynamics of membrane potential, where τ is the membrane’s RC time constant (Eqn. 6).
Equation 5 can be integrated to obtain equation 7 where is the last spike time, ur is the value that the membrane potential is set to after a spike (reset membrane potential), and s is a term for integrating from time t back to the last spike time .
By specifying the reset membrane potential as 0 and by including an absolute refractory period of 1 ms , equation 7 can be rewritten as equation 8 and equation 9, where is the last spike time adjusted for the absolute refractory period.
The model parameters τ, C, and are determined through model fitting (Appendix).
Four analyses of the skin-receptor model were undertaken. First, the skin-receptor model’s predictions were validated against published SA-I firing rates. Second, the general timing of predicted first spikes was examined. Third, the modulation of firing rates due to edges was quantified. Finally, differences in predicted first spike latencies due to edges were examined.
Model validity was determined by fitting and evaluating the skin-receptor model against primate SA-I firing rates published by Philips and Johnson . Phillips and Johnson recorded these firing rates as various gratings were indented into the fingerpad at multiple locations. Indentation locations were successively incremented by 0.2 mm along the horizontal axis perpendicular to the length of the finger, spanning the receptive field of a single SA-I. At each location, gratings were indented to a depth of 1 mm with a ramp-up time of 50 ms. Elicited spikes were recorded for a 1 second time period and converted to firing rates.
The subset of this data used here consists of firing rates elicited from a single SA-I by gratings with 3 mm and 5 mm central bars (Fig. 1), indented in 27 and 45 successively lateral indentations, respectively. Of these 72 in vivo firing rates, two thirds of those associated with each grating were selected at random to be part of the training set (18 and 30 data points, respectively) with the remaining one third reserved for the test set (9 and 15 data points, respectively).
The skin-receptor model was evaluated with the goodness of fit between model predictions and observed firing rates, calculated for both in-sample and out-of-sample data. The in-sample test determines how well the model can learn and represent the relationship between dependent (firing rate) and independent (location and grating) variables during model fitting, i.e., how well the model predicts data in the training set. The out-of-sample test determines how applicable the learned relationship between dependent (firing rate) and independent (location and grating) variables is for new independent variables, i.e., how well the model predicts data in the test set. Goodness of fit is calculated by equation 10, where observedi is the in vivo firing rate for a particular location-grating combination, and predictedi is the predicted firing rate for the same stimulus.
In vivo, first spikes are generated during the stimulus ramp-up, i.e. during the first 50 ms . Noted earlier, edge information may be carried by these first spikes. To assess the model’s ability to produce spike times, first spike latencies for all location-grating combinations are examined.
Changes in firing rate due to edge proximity have been previously quantified as the modulation index . The modulation index compares maximum, Rmax, and minimum, Rmin, firing rates, for like gratings (Eqn. 11).
Here, modulation indices are calculated for both predicted and in vivo firing rates elicited by the 3 mm and 5 mm gratings. Thus, four modulation indices are calculated.
Edge proximity may cause a change in first spike latency, just as it is known to change firing rate. While the modulation index of firing rates is generally accepted as a physiologically relevant measure, this measure may not be appropriate for first spike latencies. Furthermore, while this work utilizes in vivo firing rates, the associated in vivo first spike latencies are not available. For these reasons, examining the impact of edges on first spike latencies is restricted to differences in first spike latencies predicted for edges as compared to those predicted for flat surfaces. Here, differences in predicted first spike latencies are calculated for both the 3 mm and 5 mm gratings. Thus, two differences are calculated. Plots of first spike latency versus grating location for both the 3 mm and 5 mm grating are also examined.
Presented are the results from model validation, first spike analysis, firing rate modulation analysis, and first spike edge impact analysis.
Comparing predicted and observed firing rates for like stimuli in the training set yields an in-sample R2 of 0.79 for all stimuli combined, and R2 fits of 0.83 and 0.75 for the 3 mm and 5 mm gratings, respectively. In contrast, the out-of-sample R2 was 0.80 for all stimuli combined and 0.78 and 0.82 for the 3 mm and 5 mm gratings, respectively. Predicted firing rates overlay the observed firing rates for the 3 mm and 5 mm gratings (Fig. 3).
Fifty-nine of the 72 firing rates predicted for the 3 mm and 5 mm gratings are greater than zero. Therefore, 59 first spike latencies are available for analysis. Twelve of the 59 first spike latencies, or 20%, predicted for the 3 mm and 5 mm gratings arrive in under 50 ms from stimulus onset (Fig. 4).
Modulation indices are non-zero for both predicted and in vivo firing rates across the 3 mm and 5 mm gratings (Table 1). Compared to predicted firing rates, modulation indices are greater for in vivo firing rates.
For the 3 mm grating, the minimum and maximum first spike latencies are 44.2 ms and 70.4 ms, respectively. Thus, edge proximity shortened predicted first spike latencies by 26.2 ms for the 3 mm grating. For the 5 mm grating, the minimum and maximum first spike latencies are 43.8 ms and 85.6 ms, respectively. Thus, edge proximity shortened predicted first spike latencies by 41.8 ms for the 5 mm grating. This edge impact can be seen in figure 5, where predicted first spike latencies are shorter when an edge is in close proximity of the receptor.
By coupling a finite element model of skin mechanics, a sigmoidal transduction function, and a leaky integrate-and-fire model of neural dynamics, the skin-receptor model predicts spikes elicited by sustained indentation. By predicting spikes, the skin-receptor model facilitates examining multiple measures of the neural response, such as first spike latencies and firing rates. For example, the results indicate that the model can predict the observed firing rates elicited by 3 mm and 5 mm gratings, and that predicted firing rates are enhanced by edges, similar to what is observed in neural recordings. Additionally, the results suggest that first spike latencies are also enhanced by edges.
For all gratings combined, a R2 fit of 0.80 to firing rate data was achieved. The model’s ability to predict the in vivo firing rates lies within the range previously reported (between 0.80 and 0.97). The disparity between the R2 fit achieved by the skin-receptor model and the higher values reported may be due to differences in the nature of the data to which they were fit. For example, while this work fit to neural responses elicited by grating stimuli, a number of skin mechanics models have been fit to neural responses elicited by isolated bar stimuli [6, 7, 16]. In particular, the Gerling model achieved a R2 fit of 0.96 to neural responses elicited by isolated bar stimuli . However when the Gerling model is fit to the neural data for a 3 mm grating, as used in this study, the resultant R2 fit is 0.80, comparable to that of the skin-receptor model.
This may indicate that it is either a) inherently more difficult to fit the in vivo firing rates for the 3 mm grating opposed to the 3 mm bar, or b) generate suitable SED distributions for grating stimuli (Fig. 6). An indication of the latter can be obtained by examining the SED predicted for the 3 mm grating used here and the 3 mm isolated bar used in previous studies. SED predictions for the 3 mm grating are non-symmetric about the central indentation, while the firing rates elicited by the same 3 mm grating are symmetric. This symmetry mismatch is likely responsible for the lower R2 fit reported here, although further analysis is required to confirm this hypothesis.
The skin-receptor model provides a computational test bed for investigating the complex interaction of skin mechanics, transduction, and neural dynamics in the generation of spikes. By predicting spikes, the model allows one to analyze the impact of each sub-model on multiple measures of the neural response. Here, we examined both firing rates and first spike latencies. The model suggests first spike latencies and firing rates are enhanced by edges, indicating that edge information may be available in both. These predictions are in agreement with previous studies, which demonstrated that firing rates are modulated by edges  and that different indentation angles can impact first spike latencies . However, only 20% of predicted first spike latencies, as opposed to the nearly 100% that would be observed in vivo, are generated during the 50 ms stimulus ramp-up. This may stem from the model not differentiating between the dynamic and static phase of the SA-I response, as first spikes are generally part of the dynamic phase.
This shortcoming may be addressed in the future by extending the skin-receptor model so as to exhibit both high firing rates and rapid first spikes in the dynamic phase. As the dynamic phase is sensitive to indentation velocity and magnitude while the static phase is sensitive to indentation magnitude, two superimposed transduction functions, one a function of SED and one a function of change in SED, may allow the model to exhibit the biphasic nature of the SA-I response . Additional phenomena that the model does not currently reproduce include the irregular inter-spike intervals of the static phase and receptive fields with multiple regions of sensitivity .
To predict SA-I spike times in response to sustained indentation at the skin’s surface, the skin-receptor model combines a finite element model of skin mechanics, a sigmoidal function of transduction, and a neural dynamics model. In the analysis of the model, predicted firing rates matched published firing rates with a 0.80 goodness of fit, 20% of predicted first spike latencies arrived within 50 ms, and both predicted first spike latencies and firing rates were enhanced by the presence of edges. This work provides a computational test bed for investigating how the interaction of skin mechanics, transduction, and neural dynamics impact the generation of spikes.
The project described was supported by grants from the Defense Advanced Research Projects Agency (DARPA) and the National Library of Medicine (Grant Number T15LM009462). The content is solely the responsibility of the authors and does not necessarily represent the official views of DARPA, the National Library of Medicine, or the National Institutes of Health. The authors would like to acknowledge E.A. Lumpkin, Baylor College of Medicine, for her insightful feedback throughout this effort.
Model fitting is achieved through response surface methodology (RSM) . The goal of our RSM strategy is to find the combination of values for the model parameters that minimize the sum of squared deviations (ssd) between observed and predicted firing rates for like stimuli in the training set (Eqn. 12), where observedi is the in vivo firing rate for the location-grating combination i, and predictedi is the predicted firing rate for the same location-grating combination.
The RSM process utilized here can be broken into 5 iterative steps.
In the first step, model parameters τ, C, , α, γ, and λ are coded (Eqn. 13), where xk is the coded value corresponding to the actual value of parameter k, ξk.
For the first RSM iteration (i.e., first sequence of 5 iterative steps) the base values, ξkbase, are the initial estimates, and for subsequent RSM iterations they are the values that minimized the ssd in the previous iteration. Coding increments, Δξk, are chosen through trial and error to impact predicted firing rates while remaining relatively small. Most coding increments were less than 5% of the associated base value, and all were less than 12%. Although RSM calculations are carried out in coded variables, those coded variables are converted back to model parameter values (Eqn. 14) when used in the skin-receptor model.
In the second step, model parameters are systematically varied and the ssd (Eqn. 12) is calculated for each variation. Each calculation of ssd for a set of model parameter values is termed a run. Factorial runs are performed where xτ, xc, x, xα, xγ, and xλ each take one of two levels: 1 or −1. Six model parameters and two levels for each gives 64 combinations and therefore 64 factorial runs. These factorial runs are supplemented by one run where the model parameters all take the coded value 0. This yields 65 total runs, each of which relates a different set of model parameter values to a resulting calculation of ssd.
In the third step, linear regression is used to generate a first order approximation of the true relationship between model parameters and ssd (Eqn. 15). The sum of each coded variable, xk, weighted by the regression coefficient, βk, gives an estimate of the ssd, .
In the fourth step, the regression coefficients of the approximation (Eqn. 15) are analyzed to determine the directions and relative magnitudes in which to change model parameters in order to most significantly decrease the ssd. This vector is the path of steepest decent, Δ (Eqn. 16). To calculate the path of steepest descent, the single greatest absolute regression coefficient, βmax, is determined and the change in the associated coded variable, Δxmax, is defined as 1 if βmax is negative and −1 if βmax is positive. Then, each remaining Δxk is calculated from equation 17.
In the fifth step, the skin-receptor model is run with sets of model parameters from successive steps along the path of steepest descent. These runs are completed until the sum of squared deviations between observed and predicted firing rates no longer decreases.
After step 5, the RSM process repeats. These 5 steps are iterated until little or no decrease in ssd is obtained from further iterations. Typically 2 to 4 iterations are required, unless the true relationship being approximated is atypical.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.