1. Computational and Systems Biology
  2. Neuroscience
Download icon

Rate and timing of cortical responses driven by separate sensory channels

  1. Hannes P Saal  Is a corresponding author
  2. Michael A Harvey
  3. Sliman J Bensmaia  Is a corresponding author
  1. University of Chicago, United States
  2. Vanderbilt University Medical Center, United States
Research Article
  • Cited 10
  • Views 1,020
  • Annotations
Cite as: eLife 2015;4:e10450 doi: 10.7554/eLife.10450

Abstract

The sense of touch comprises multiple sensory channels that each conveys characteristic signals during interactions with objects. These neural signals must then be integrated in such a way that behaviorally relevant information about the objects is preserved. To understand the process of integration, we implement a simple computational model that describes how the responses of neurons in somatosensory cortex—recorded from awake, behaving monkeys—are shaped by the peripheral input, reconstructed using simulations of neuronal populations that reproduce natural spiking responses in the nerve with millisecond precision. First, we find that the strength of cortical responses is driven by one population of nerve fibers (rapidly adapting) whereas the timing of cortical responses is shaped by the other (Pacinian). Second, we show that input from these sensory channels is integrated in an optimal fashion that exploits the disparate response behaviors of different fiber types.

https://doi.org/10.7554/eLife.10450.001

eLife digest

Our sense of touch depends upon receptors in our skin that send signals to the brain about the objects with which we interact. Different types of touch receptors respond in different ways when we grasp and manipulate objects; for example, by altering the strength of their response or its timing.

Saal et al. have now investigated how neurons in a part of the brain called the somatosensory cortex, which processes touch signals from the hand, respond to the information from the different receptor types. First, recordings were made of the electrical activity of the touch receptors and the neurons in the brain of monkeys. Using this data, Saal et al. built computer models that allow the response of neurons in the brain to be predicted from the responses of the touch receptors.

The models showed that signals from different types of touch receptors shape the response of neurons in the brain in different ways. One receptor type controls how strong a neuron’s response will be, while another one controls the precise timing of the response. Further investigation revealed that this way of combining the signals from the different receptors preserves as much information as possible about objects and thereby helps the brain to process information acquired by touch quickly and efficiently.

Future experiments will examine how touch is represented in two structures that lie between the receptors and the somatosensory cortex: one in the brainstem, the other in a brain region called the thalamus.

https://doi.org/10.7554/eLife.10450.002

Introduction

Perception reflects the seamless integration of signals from a variety of sensory receptors that respond preferentially to different aspects of the environment. The classic example of such integration is color vision, where input from different cone receptors is integrated to extract a specific feature of objects, approximately corresponding to their absorption spectrum (Gegenfurtner, 2003). Signals are also combined across different sensory modalities, and this integration process is often optimized to extract stimulus information (Ernst and Banks, 2002). Given that different sensory channels encode different stimulus features and exhibit different response properties, the central nervous system is faced with the challenge of how to integrate such disparate input signals.

The sense of touch is mediated by at least three main classes of mechanoreceptive afferents, each of which responds to different aspects of skin deformation (Johnson, 2001): slowly-adapting type I (SA1) afferents are most sensitive at low frequencies (<10 Hz), rapidly adapting (RA) afferents at intermediate frequencies (in the so-called flutter range, from 10 to 50 Hz), and Pacinian (PC) afferents at high-frequencies (peaking in sensitivity at around 250 Hz). However, the frequency sensitivities of the different afferent classes overlap considerably (Figure 1A), and natural stimuli typically engage more than one class (Saal and Bensmaia, 2014). Not surprisingly, then, information about most tactile events is carried by signals from several afferent classes (Saal and Bensmaia, 2014; Johansson and Flanagan, 2009; Abraira and Ginty, 2013), and this input is integrated across submodalities even at the first stage of cortical processing, namely primary somatosensory cortex (S1) (Pei et al., 2009; Carter et al., 2014). Importantly, S1 neurons convey information about contacted objects not only in the strength of their responses but also in their timing (Harvey et al., 2013; Zuo et al., 2015). For example, the amplitude of skin oscillations—such as those elicited during the exploration of textured surfaces—is encoded in the strength of the cortical response whereas their frequency is encoded in its timing. To achieve such multiplexing requires that afferent signals be integrated in a principled way.

Figure 1 with 2 supplements see all
Characterizing the transformation between peripheral and cortical responses. 

(A) Frequency–amplitude pairs that elicit responses in rapidly adapting (RA) and Pacinian (PC) afferents. Light blue and orange lines indicate absolute response thresholds for RA and PC afferents (as determined by the peripheral model), respectively, while the dark blue and orange lines denote tuning thresholds (i.e., eliciting one spike on each stimulus cycle). The black line corresponds to the highest amplitude the stimulator could deliver at a given frequency and likely exceeds the maximum amplitude that one is liable to encounter during every day tactile experience. While RA and PC afferents are maximally sensitive at different frequencies, most frequency–amplitude pairs elicit responses from both afferent classes (indicated by the shaded region). (B) Illustration of the model that describes the transformation from peripheral to cortical responses. A broadband noise stimulus activates both RA and PC afferents. RA and PC population responses are each convolved with a temporal filter, then summed and rectified to culminate in the responses of individual S1 neurons. (C) Two examples of vibratory stimuli. (D) Simulated RA (blue) and PC (orange) population firing rates along with spike trains of a subset of neurons in the population. See Figure 1—figure supplement 1 for details on how the peripheral population models were validated. (E) Recorded spikes (black ticks) over five stimulus repeats and model predictions (green traces) for two cortical neurons, whose RA and PC filters (plotting filter magnitude over time) are shown on the left. Cortical neurons differ in their response properties, such as burstiness and temporal precision, and the model captures these broad differences with different RA and PC filters. See Figure 1—figure supplement 2 for assessment of the general predictive power of our model.

https://doi.org/10.7554/eLife.10450.003

Here, we measure the responses of S1 neurons to a wide variety of simple and complex vibrotactile stimuli and reconstruct the responses of afferent populations to these same stimuli. We then implement a simple model that allows us to examine how responses of the different afferent populations combine to drive responses in cortex. First, we find that most S1 neurons integrate information from multiple afferent classes. Second, signals from the different classes drive cortical responses in very different ways: RA afferent input exerts an excitatory influence on S1 responses and is the primary determinant of their firing rates; in contrast, PC input exerts a balanced excitatory and inhibitory influence which sharpens the timing of cortical responses without having a major impact on the rates. Finally, we show that the process of integration maximizes information transmission across the range of stimulation conditions that might be encountered during everyday interactions with objects.

Results

We recorded the responses of 18 peripheral tactile afferents (14 RA and 4 PC) and 118 cortical neurons in areas 3b (49) and 1 (69) of S1 to a variety of vibrotactile stimuli, ranging from simple sinusoids to complex bandpass noise with frequencies between 50 and 800 Hz and amplitudes ranging from sub-micrometer to more than 500 µm. Stimulation conditions spanned those that might be experienced when we scan a textured surface (Bensmaïa and Hollins, 2003; Bensmaïa and Hollins, 2005; Weber et al., 2013; Manfredi et al., 2014). As these stimuli drove RA and PC to the exclusion of SA1 afferents (Muniak et al., 2007), we focused the analysis on these two afferent populations. The goal was to quantitatively predict S1 responses based on afferent responses. Because the stimuli used in the peripheral experiments were analogous but not identical to those used in the cortical experiments, we used simple but precise models of afferent responses (derived from recorded responses, Figure 1—figure supplement 1) to reconstruct RA and PC population responses to the stimuli used in the cortical experiments (Figure 1C,D). To examine how input from these different populations is integrated in cortex, we describe the transformation from the periphery to cortex as a linear filter over the RA and PC population responses whose output is subsequently rectified (based on the linear-nonlinear-Poisson framework (LNP), Figure 1B). This approach allows us to examine the filters estimated for each cortical neuron to assess the contribution (if any) of each afferent population to that neuron’s response and to characterize the dynamics of this integration. To the extent that this simple model captures different response properties of cortical neurons, the respective filters for these neurons will be different (Figure 1E).

To test the model’s performance, we compared predicted and recorded cortical responses at different timescales. We performed this analysis on a separate data set (consisting of sinusoidal stimuli), to estimate how well the model would generalize to stimuli that differed considerably from the noise traces used to fit the model. We found that the model captures the cortical responses well on both coarse and fine timescales as assessed by calculating the correlation coefficients between the predicted and recorded responses. Indeed, across neurons, the model captured 67% ( ± 2% standard error of the mean) of the predictable variance (Ramirez et al., 2014) for 1 s bins, and 35% ( ± 2%) for 5 ms bins. In other words, the model was able to predict neuronal responses with considerable accuracy given their inherent variability.

RA and PC input is integrated differently

We found that most cortical neurons in areas 3b and 1 receive input from both RA and PC afferents, as evidenced by the significant contribution of both to the overall prediction (Figure 2A). Furthermore, the linear filters obtained for each afferent population, while similar across cortical neurons, were consistently different from each other (Figure 2B): First, while the magnitude of the PC filters was correlated with that of the RA filters (r = 0.41), PC filters were consistently of lower magnitude (Figure 2C, t(70) = 6.4, p < 0.01, paired t-test), with the RA contribution on average about three times stronger than the PC one. Second, RA filters were much more extended in time (with a mean above 30 ms) than were PC filters (with a mean below 20 ms) (Figure 2D, t(188) = 9.6, p < 0.01, two-sample t-test). Third, while the net contribution of the RA input to the cortical response was positive, PC input was on average balanced; that is, excitatory and suppressive components of the PC filters were of equal magnitude. In fact, for about half of cortical neurons, the PC contribution was net suppressive (Figure 2E); the difference in the net contribution of RA and PC input was significant (t(188) = 16.1, p < 0.01, two-sample t-test). PC filters were also much more diverse and spanned the range of this metric (from -1 to +1), in contrast to RA filters, which, as mentioned above, tended to be net excitatory (>0).

Figure 2 with 1 supplement see all
Rapidly adapting (RA) and Pacinian (PC) input is integrated differently.

(A) Fraction of neurons in areas 3b and 1 that receive input from both RA and PC populations (grey bar) or from only the RA population (blue) as assessed by whether input from the respective class is weighted more highly than might be expected by chance. Most S1 neurons receive input from both RA and PC afferents. (B) Mean RA and PC filters (solid lines) and the 5th and 95th percentile of filter values across all S1 neurons indicated by the light-shaded regions. RA and PC filters are systematically different. (C) Ratios of RA and PC filter magnitudes for neurons integrating information from both afferent classes. RA input tends to be weighted more strongly. (D) Width of RA and PC filters. RA input tends to be integrated over longer timescales than does PC input. Insets show example filters that are close to the population average (over 100 ms). (E) Excitation indices for all RA and PC filters, where 1 denotes purely excitatory input, 0 balanced input, and -1 purely suppressive input. RA filters almost always lead to a net positive excitation of the cortical neurons, while the PC filters are more diverse, on average balanced, and often suppress cortical response. Insets show example filters that are close to the population average (over 100 ms). See Figure 2—figure supplement 1 for validation of the analysis.

https://doi.org/10.7554/eLife.10450.006

RA input drives cortical firing rates

The different filter shapes for RA and PC input reflect profound differences in the way the two afferent populations shape cortical responses. Indeed, the firing rates of individual S1 neurons are almost entirely determined by RA input as evidenced by the fact that a model based on RA input alone can predict S1 firing rates—averaged over 1 s long periods—almost as well as RA and PC responses combined; this was true for responses to both noise and sinusoids (Figure 3). PC input alone, on the other hand, cannot account for the observed cortical firing rates. Indeed, these would be much more frequency-dependent if driven by this population of nerve fibers.

Rapidly adapting (RA) input determines cortical firing rates.

Reconstruction accuracy of the full model compared to one that includes only RA or Pacinian (PC) input as measured by correlation between actual and predicted firing rates (measured over the duration of the trial, namely 1 s). Error bars denote the standard error of the mean across neurons. Firing rates elicited by both noise and sinusoidal stimuli are well predicted by the full model. While the performance of RA-only filters is almost as good as the full model, reconstruction accuracy drops dramatically if only PC input is included.

https://doi.org/10.7554/eLife.10450.008

PC input determines precise spike timing

We have previously shown that the responses of mechanoreceptive afferents are highly precise and repeatable, and their timing conveys information about the frequency composition of skin vibrations (Mackevicius et al., 2012) and about the identity of textures scanned across the skin (Weber et al., 2013) at resolutions on the order of milliseconds. Similarly, the temporal patterning of responses in somatosensory cortex is precise and informative down to a millisecond timescale (Harvey et al., 2013). RA and PC input are thus integrated in such a way to culminate in temporally patterned cortical responses. One might expect that the shorter integration window of the PC filters combined with their strong suppressive components would result in this sensory channel driving more precise spike timing than does its RA counterpart. Indeed, when removing the PC contribution from the model, predictions of spike timing in the cortex became less temporally precise (Figure 4A). We quantified this effect in two ways. First, we filtered the recorded and predicted time-varying firing rates to noise stimuli at different temporal resolutions and examined how much of the variance in the responses could be accounted for based on RA input or PC input alone (Figure 4B). We found that only PC input can account for high frequency fluctuations in the cortical responses of many neurons, while RA input is sufficient at coarser temporal resolutions for almost all cortical neurons (Figure 4C). In fact, the more precise the responses of a cortical neuron were, the more of its response variance was explained by PC input (Figure 4—figure supplement 1). Second, we tested how well the model could account for the well-documented patterning in cortical responses to sinusoidal vibrations: over a wide range of frequencies, cortical neurons produce entrained responses to sinusoidal stimulation of the skin (Harvey et al., 2013; Salinas et al., 2000; Mountcastle et al., 1969). That is, they produce one spike or burst of spikes within a very restricted phase of each stimulus cycle. We examined, then, whether we could predict the degree of entrainment of cortical neurons to sinusoidal stimuli from the RA and PC input. We found that the strength of phase-locking to sinusoids ranging in frequency from 50 to 300 Hz was well captured by the model (yielding a correlation between predicted and measured phase-locking of 0.78). Importantly, while RA input was dominant at the low frequencies, only the PC input could account for the observed phase-locking at 200 Hz and above (Figure 4D).

Figure 4 with 1 supplement see all
Precise spike timing of S1 neurons is driven by Pacinian (PC) input.

(A) Responses to noise stimuli from two example neurons with particularly precise responses. Spiking probabilities derived from the full model (green traces) are much more temporally precise than are those derived from a rapidly adapting (RA)-only model (blue traces), reflecting the importance of PC input in shaping the response at fine timescales. Five repetitions of the recorded cortical spike trains are shown below the model traces in black. (B) Fraction of the explained variance (R2) by either the RA (blue) or PC (orange) model relative to the full model across different temporal resolutions for the two neurons shown in A. PC input is important to explain responses on fine timescales. (C) Fraction of explained variance (R2 fraction) relative to the full model by RA and PC input for all neurons at two different temporal resolutions. While RA input can explain most of the observed variance at coarse temporal resolutions (left panel), PC input is needed to explain the timing of cortical responses at fine temporal resolutions (right panel). See Figure 4—figure supplement 1 for further analysis. (D) Vector strengths predicted by RA (blue) and PC (orange) models relative to their measured counterparts at different frequencies across all neurons. RA input accounts for entrainment to sinusoidal stimuli at low frequencies, while PC input is needed at higher frequencies.

https://doi.org/10.7554/eLife.10450.009

RA and PC input is integrated optimally to process natural tactile scenes

Having established how input from different sensory channels is integrated in somatosensory cortex, we sought to understand why the two channels drive cortical neurons in different ways. We reasoned that the process of integration should reflect differences in the response properties of the input channels. Specifically, PC responses to skin vibrations and scanned textures have been consistently shown to be more temporally precise than their RA counterparts (Weber et al., 2013; Mackevicius et al., 2012). Furthermore, the two afferent classes differ in their frequency sensitivity profile: PC afferents are considerably more sensitive than RA afferents at frequencies above about 150 Hz (see Figure 1A). Finally, we surmised that the process of integration would also reflect the statistics of the stimuli experienced during interactions with objects, as is the case in other sensory modalities (van Hateren, 1992; Field, 1987; Smith and Lewicki, 2006).

We focused on natural textures, which are ideally suited to address this question. When running our finger across a textured surface, high-frequency oscillations are elicited in the skin (BensmaÏa and Hollins, 2003; 2005; Manfredi et al., 2014), which in turn elicit highly temporally patterned and repeatable responses in both RA and PC afferents (Weber et al., 2013). In cortex, neurons respond to skin oscillations with precisely timed spikes that encode the skin deflections with millisecond precision (Harvey et al., 2013). In light of this, how should RA and PC responses to natural tactile stimuli be integrated to be maximally informative about texture-elicited skin oscillations at a fine temporal resolution? To address this question, we simulated RA and PC responses to different textures scanned across the skin and optimized cortical RA and PC filters for information transmission about skin oscillations.

We found that the filters optimized to convey information about texture closely matched the filters derived from measured S1 responses: RA filters were mostly excitatory and temporally broad; PC filters were narrower and comprised distinct suppressive components (Figure 5A). The similarities between the recorded and optimized filters also held up quantitatively. First, we found that most (71%) optimized filters comprised contribution from both RA and PC fibers, while the rest relied on RA input alone, similar to what we found in our cortical data set. Second, we found that the optimal filters typically weighted the RA input more strongly than the PC input, as did the filters derived from cortical responses (Figure 5B). Third, the optimal RA filter widths were wider than PC filter widths and largely overlapped the range of the actual filter widths (Figure 5C). Finally, we found that the excitation/suppression indices of the optimized filters approximated those derived from our cortical data set in that RA input was always net excitatory, while PC input spanned the range between highly suppressive, balanced, and excitatory (Figure 5D).

Rapidly adapting (RA) and Pacinian (PC) input is integrated optimally.

(A) Three examples RA and PC filters, each optimized for a different natural texture (colored traces) and paired with a filter from a cortical neuron that most closely matches it (gray traces). (B–D) Comparison of filter statistics for actual (black) and optimized (lavender), and distorted (green) filters for the RA/PC magnitude ratio (B), filter width (C), and excitation/suppression index (D). For the optimized filters, bars denote the range covered by the filters optimized for different textures. For the actual filters, bars indicate the range between the 10th and the 90th percentile.

https://doi.org/10.7554/eLife.10450.011

To test whether our fitted and optimized filters might be an artifact of our approach rather than reflect the integration properties of S1 neurons, we also optimized filters on a distorted data set. During natural interaction with surfaces, the amplitude of elicited vibrations decreases as their frequency increases, leading to roughly a power–law relationship (Manfredi et al., 2014; Wiertlewski et al., 2011; Delhaye et al., 2012). When we re-optimized the filters to stimuli with a distorted frequency composition, we found that the resulting filters did not match the ones observed in cortex (green bars in Figure 5B–D): The PC signal was weighted more heavily, was integrated over longer timescales, and was mostly excitatory, in contrast to what we observed in both the cortical filters and the ones optimized on our natural texture data set. We conclude that RA and PC integration in S1 is optimized to encode stimuli that are commonly encountered during interactions with objects.

Discussion

We set out to examine how input from different mechanoreceptive afferent classes is integrated in the responses of individual cortical neurons. To this end, we simulated the responses of populations of RA and PC tactile afferents and examined how these are combined to drive S1 responses. This modeling effort yielded several important conclusions. First, the vast majority of cortical neurons receive (indirect) input from both RA and PC afferents. Second, RA and PC afferents drive S1 neurons in fundamentally different ways: RA input determines the rate of the response whereas the principal contribution of PC input is to sharpen the temporal precision of the response. Third, the integration of afferent input is consistent with a process that maximizes the information conveyed in the cortical responses by exploiting the different response properties of RA and PC afferents.

Input from different sensory channels is integrated in the responses of individual S1 neurons

Our results unequivocally show that RA and PC input are integrated in the responses of individual S1 neurons. That is, the responses of S1 neurons cannot be explained from the responses of RA or PC neurons alone. While it has traditionally been assumed that the different tactile submodalities are segregated at least through primary somatosensory cortex (Sur et al., 1981; Sur et al., 1984; Chen et al., 2001; Friedman et al., 2004) (although cortical PC neurons have generally proven elusive, see below), recent evidence suggests that input from different classes of somatosensory afferents is indeed integrated between periphery and cortex (Saal and Bensmaia, 2014; Pei et al., 2009; Carter et al., 2014; Prescott et al., 2014). In this study, stimuli fell in the range of frequencies that excites RA and PC fibers at the exclusion of SA1 fibers, which only respond to low-frequency stimulation. However, SA1 and RA input have been shown to be integrated in the responses of S1 neurons when these were studied with ramp-and-hold skin indentations (Pei et al., 2009) so the convergence of submodality-specific input onto individual S1 neurons seems to be a general phenomenon. However, it is still an open question whether the integration of input from different mechanoreceptors takes place only at the cortical level, or whether they appear already at the subcortical level in either the cuneate nucleus or the thalamus. Convergence of somatosensory submodalities has been observed in the brainstem of rats (Sakurai et al., 2013) and the same is likely true in primates.

The effect of PC input on S1 responses

Our results reconcile some seemingly contradictory findings regarding the contribution of PC afferents to cortical responses. First, vibrotactile coding in general and the role of PC input specifically have been probed almost exclusively with high-frequency sinusoidal stimuli. These stimuli have been shown to evoke only a short transient response in S1, which then quickly decays back to baseline, a phenomenon that can be accounted for by the balance of excitation and suppression in the PC drive. That PC-elicited responses are short and weak might explain why PC-driven cortical neurons are thought to be extremely rare in S1 (Lebedev and Nelson, 1996; Mountcastle et al., 1990). Second, the main effect of PC input is to shape the timing of cortical spikes, which could account for the extremely precise timing of S1 responses at high frequencies (Harvey et al., 2013) that has been shown to contribute to tactile perception (Zuo et al., 2015). Third, the balanced drive of PC afferents explains why the firing rates of S1 neurons are not dependent on stimulus frequency even though PC rates are (Harvey et al., 2013). Finally, the PC input has been found to inhibit cortical responses (Tommerdahl et al., 2005; Tommerdahl et al., 2010), consistent with our finding that PC filters are often net suppressive.

Efficient coding and the integration of sensory channels

Sensory cortex is tuned to the statistics of the natural environment, as has been demonstrated most notably in vision (van Hateren, 1992; Field, 1987) and audition (Smith and Lewicki, 2006), by invoking the efficient coding principle (Barlow et al., 1961). Generally, it has been assumed that the peripheral input stems from a single homogeneous population of neurons (e.g. retinal ganglion cells). However, this is often not the case: in touch, information is carried by multiple separate populations of neurons, each with different sensitivities and response properties. To encode information efficiently, then, the disparate response properties of these different channels need to be taken into account. While both RA and PC afferents respond to most high-frequency skin vibrations, such as those elicited during texture scanning, both afferent classes do so in different ways. First, RA afferents increase their firing rate with increasing stimulus amplitudes (Arabzadeh et al., 2014), and do so in a manner that is relatively independent of stimulus frequency, while PC responses saturate at low stimulus amplitudes and are highly frequency-dependent (Muniak et al., 2007). RA firing rates therefore convey a much more informative signal about the amplitude of tactile stimuli than do PC firing rates. That S1 firing rates depend mostly on RA input exploits the fact that these afferents convey unambiguous information about stimulus amplitude. Second, PC afferents exhibit much higher temporal precision in their spiking patterns than do RA afferents (Mackevicius et al., 2012), especially for high-frequency skin vibrations such as those elicited when we run our fingers across a texture (Weber et al., 2013), and this precise spike timing conveys information about tactile events and shapes perception (Mackevicius et al., 2012). That the timing of S1 responses is driven primarily by the PC input exploits the highly informative nature of the PC timing signal.

Separate sensory channels that convey complementary but overlapping information are commonplace in sensory systems and not limited to the sense of touch; the visual (Field and Chichilnisky, 2007), auditory (Cant and Benson, 2003), gustatory (Zhang et al., 2003), olfactory (Uchida et al., 2014), and vestibular (Goldberg, 2000; Sadeghi et al., 2007) systems all involve many types of afferents. Distributing sensory information across channels has distinct advantages such as parsing the behaviorally relevant range (Dominy and Lucas, 2001), keeping energy expenditure low (Gjorgjieva et al., 2014), and optimizing information transmission in the presence of noise (Kastner et al., 2015). To the extent that these parallel input channels represent information in disparate ways (differing in response latency, adaptation properties, or spiking precision, among others), their integration should reflect and exploit such differences, a process that can be described using the type of model introduced here.

Materials and methods

Stimuli

The vibrotactile stimuli have been described previously, for both the peripheral (Muniak et al., 2007) and the cortical (Harvey et al., 2013) experiments. In short, tactile stimuli were delivered to the distal pads of the digits using a stainless steel probe diameter driven by a shaker motor. The shaker motor was calibrated before each experimental run such that stimuli were highly accurate and repeatable. We delivered three types of stimuli that varied in spectral complexity, including simple sinusoids, diharmonic stimuli (only used for the peripheral model fitting), and bandpass noise. The ranges of frequencies and amplitudes spanned the range experienced in everyday tactile experience. While the exact frequencies and amplitudes differed somewhat between the peripheral and cortical data sets, the two stimulus sets were largely overlapping. We were thus able to develop peripheral spiking models that simulated the peripheral responses to the stimuli used in the cortical data set (see below). Individual stimuli in this data set were 1 s long and were repeated five times (frozen noise), yielding a total of about 22 minutes (16 for noise and 6 for the sinusoids) of data from each neuron.

Electrophysiology

Cortical experiments

Experimental procedures have been previously described (Harvey et al., 2013) and are only summarized here. All procedures were in accordance with the rules and regulations of the University of Chicago Animal Care and Use Committee. Extracellular recordings were made in the postcentral gyri in four hemispheres of two awake, behaving Rhesus macaques (Macaca mulatta). Recordings were obtained from neurons in areas 3b and 1 that met the following criteria: (1) action potentials were well isolated from the background noise, (2) the RF of the neuron was on the glabrous skin, and (3) the neuron was clearly driven by light cutaneous stimulation. In total, recordings were obtained from 49 area 3b and 69 area 1 neurons. Some of these data have been previously published (Harvey et al., 2013).

Peripheral experiments

Experimental procedures have been previously described (Muniak et al., 2007). All experimental protocols complied with the guidelines of the Johns Hopkins University Animal Care and Use Committee and the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Single unit recordings were made from the ulnar and median nerves of two anesthetized Rhesus macaques using standard methods (Talbot et al., 1968). Recordings were obtained from 14 RA and 4 PC afferents. This data set has been previously used to develop peripheral spiking models (Kim et al., 2010; Dong et al., 2013).

Peripheral model

We fit integrate-and-fire (IF) models to the measured responses of individual afferents. Our IF model is similar to an earlier one (Kim et al., 2010), with two main differences. First, we did not employ complex temporal filters with a large number of weights, but instead restricted the model input to six dimensions (positive and negative rectified stimulus positions, velocities and accelerations) plus a delay parameter, as this parameterization has been shown to be sufficient to reproduce afferent responses (Dong et al., 2013). Second, we included an additional parameter, namely a smoothing factor that controls the width of a Gaussian window that is convolved with the stimulus trace initially. Using this parameter allows us to model the fact that the PC afferents’ sensitivity decreases for frequencies higher than ~300 Hz (Muniak et al., 2007; Mountcastle et al., 1972).

For model fitting, we used the van Rossum spike distance (van Rossum, 2001) between the recorded and model-predicted responses as a cost function and then optimized the model parameters using the patternsearch function in Matlab (The Mathworks, Inc., Natick, MA) using different starting positions. We optimized the model on the noise stimuli first and then alternated between the sinusoidal and diharmonic stimulus sets while decreasing the cost function time constant until the fits did not improve further. In total, we fit single afferent models to 14 RA and 4 PC afferents. RA and PC population responses were generated by averaging over the responses of all RA and PC individual models, respectively. The models have been extensively validated, as detailed in previous studies (Kim et al., 2010; Dong et al., 2013), and captured both the strength and timing of peripheral afferents (Figure 1—figure supplement 1).

One potential concern is that the reconstructed RA and PC population activity might not accurately reflect the peripheral input on which a given cortical neuron's response relies. While individual peripheral responses are very stereotyped and therefore simple to model, individual afferents still differ in terms of sensitivity, and these differences might be reflected in the cortical responses. To address this issue, we generated different RA and PC populations, using different peripheral models, and also varied the size of the peripheral population used to drive the cortical responses. We found that the recovered filters were robust to these changes.

Modeling the transformation between periphery and cortex

We modeled the transformation between the peripheral RA and PC population responses and single unit S1 responses using the LNP framework (Schwartz et al., 2006). LNP models have long been fruitfully used for modeling the responses of cortical neurons across several synapses, often even including the sensory receptors. The model describes the process of integration as follows: First, the peripheral RA and PC population responses are each convolved with a linear filter. Second, the resulting responses are summed and rectified to yield time-varying spiking probabilities. Third, spike trains are generated according to the spiking probability in each time bin. The full model (see also Figure 1B) can thus be expressed as follows:

(1) p(st)=r(kRA×sRA,t+kPC×sPC,t)

where p(st) denotes the probability of a cortical spike at time t, kRA and kPC are vectors containing the RA and PC filters, sRA,t and sPC,t are vectors containing the binned RA and PC population responses within a time window 100 ms prior to time t, and r() denotes the rectifying nonlinearity.

To estimate the linear filters kRA and kPC, we used reverse correlation of the cortical responses onto the RA and PC z-scored population responses covering a time window of 100 ms preceding the cortical response with a bin width of 1 ms. For this procedure, we exclusively used the data from the bandpass noise stimuli as this method can only recover unbiased filters if the RA and PC population responses are uncorrelated with each other over time. Even though our original vibrotactile stimuli exhibited little autocorrelation, the evoked RA and PC population responses were weakly correlated both across time and with each other, which could affect the filter estimation. We addressed this problem using a two-pronged approach. First, we verified that across our entire data set, temporal and population correlations were low (see below). Second, to minimize the impact of any remaining correlations, we calculated the autocorrelation matrix of the RA and PC population input and used its regularized inverse to correct the obtained filters (ridge regression):

(2) kRA,PC=[STRA,PC,T×SRA,PC,T+a*I]-1×STRA,PC,T×sC,

where kRA,PC,T is a vector containing both the RA and PC filters, SRA,PC,T is a matrix where each row corresponds to a time t and contains the peripheral population responses sRA and sPC for the 100 ms preceding t, a is the regularization parameter, I is the identity matrix, sC is a binary vector containing the responses of a cortical neuron (0 for no spike, 1 for spike) for all times t, and []-1 denotes the pseudoinverse.

We varied the regularization parameter a and also tried other correction methods, such as decomposing the autocorrelation matrix using singular-value-decomposition and only preserving strong components (Sripati et al., 2006); we obtained essentially the same results, and thus conclude that the approach is robust and the recovered linear filters accurate (see further analysis below).

The nonlinearity was a piecewise linear function, chosen over a smoother sigmoidal function, because the firing rates of S1 neurons increase steeply as soon as their respective response threshold is crossed (Harvey et al., 2013):

(3) r(x)=b1, if x<0, r(x)=b1+b2*x, if x>0 and b1+b2*x<1, and r(x)=1, otherwise.

To estimate the parameters b1 and b2, we used an optimization algorithm (lsqcurvefit in Matlab) to fit the time-varying spiking probabilities p(st) such that they matched the observed spiking probabilities to the extent possible.

Analysis of linear filters

To measure the magnitude of a filter, we summed its absolute value over time. To determine the temporal width of a filter, we first determined a noise threshold, set to three times the standard deviation of the filter values 80–100 ms before the cortical response; we then counted the total time that each filter’s absolute value was above that threshold. Finally, to measure excitation/suppression indices, we divided the sum of the positive and negative components of a filter by its overall magnitude. When analyzing the impact of removing input from one fiber class, we refitted the filters and nonlinearity using only input from the remaining class.

Testing the robustness and accuracy of filter estimation

Recovering linear filters accurately from data that exhibit correlations is notoriously difficult and error prone. In our case, one of the main concerns is whether cross-correlations between RA and PC responses might introduce biases in the recovered filters that would lead us to conclude that both populations contribute input, when in reality only one does. Additionally, the RA and PC responses are autocorrelated, which might bias our estimates of the filter width, or their net excitation or suppression, and thereby erroneously yield different filters for the two populations (see Figure 2—figure supplement 1A for the RA/PC covariance matrix). To test whether our method is prone to such biases, we simulated a population of 243 neurons with a variety of different linear filters that differed in width, magnitude, excitation/suppression, temporal offset, and received RA input, PC input, or both. We then used these filters to generate simulated spiking data from the noise vibrations. We set the firing rates to span the range encountered in the measured cortical responses. We then tested whether we could accurately recover the linear filters from the simulated responses and whether there were any biases in the measures we used to quantify these filters. First, we found that the recovered filters matched the filters used to generate the simulated data well (Figure 2—figure supplement 1B,C). Second, we verified that our procedure for determining whether neurons received RA input, PC input, or both, worked correctly in the vast majority of cases: 90% of the simulated neurons that received input from both RA and PC populations were identified as such. Importantly, only 3% of simulated neurons that received input from only one of the peripheral populations were misclassified as receiving input from both (Figure 2—figure supplement 1D). Finally, we also checked whether there might be biases when estimating filter metrics such as overall magnitude, filter width, and excitation/suppression indices in RA as compared to PC filters that might explain the observed differences. We found that these biases were generally very small and in no case could account for the effect sizes observed in the cortical data set (see comparison in Figure 2—figure supplement 1E). In summary, our procedure recovers most filter shapes accurately and the small biases that remain cannot account for the observed effect sizes.

Timing analysis

To assess the contribution of RA and PC input at different timescales, we bandpass-filtered the measured and predicted cortical responses to four different frequency ranges, corresponding to different temporal resolutions (3, 5, 13, and 66 ms) and then calculated the correlation coefficients for each neuron. To assess spike timing for the sinusoidal data set, we calculated the vector strength (Goldberg and Brown, 1969) (a measure of how precisely a neuron’s response is aligned with a sinusoidal stimulus) of the recorded and predicted cortical responses.

Optimizing filters by maximizing information

We first simulated the responses of populations of RA and PC fibers to 55 natural textures (Manfredi et al., 2014). In brief, textured surfaces were scanned over the fingertip of human subjects at different scanning speeds, while elicited vibrations were recorded using a Laser Doppler vibrometer (Polytec OFV-3001 with OFV 311 sensor head, Polytec, Inc., Irvine, CA). For each texture, we collected 240 (8 subjects × 3 speeds × 10 repetitions) 500-ms long vibration traces. The recorded traces were bandpass-filtered to between 50 and 800 Hz and then served as input to the peripheral RA and PC spiking models. To adjust for the distance-dependent decay of the vibrations and the size of the contact location, we added a global scaling factor for each the RA and PC traces which was set to ensure that the mean simulated RA and PC firing rates matched those recorded using the same set of textures (Weber et al., 2013). Indeed, across textures the recorded and simulated firing rates matched well (r = 0.83 for RA and r = 0.92 for PC responses), indicating that our simulated responses capture actual texture-elicited responses well. Some textures, mainly coarser one, elicit responses not only in RA and PC afferents, but also in SA1 afferents; we excluded textures with an average SA1 contribution >10 Hz from further analysis. Since we were interested in tactile environments that cause both RA and PC afferents to respond, we also excluded textures where the RA or PC average firing rate was below 20 Hz. In total, this left 25 textures, whose elicited vibrations where in the high-frequency range and excited RA and PC but not SA1 afferents.

Next, we optimized linear filters for each texture individually by maximizing the mutual information I(S,V) between the responses of our model (S) and the texture-elicited vibration traces (V):

(4) I(S,V)=SVp(s,v)log[p(s,v)/p(s)p(v)]

The probability of a spike in an S1 neuron, p(s), can be calculated from the LNP model described above. To estimate p(v), that is, the distribution of stimulus values, we divided the stimulus trace into 1 ms long bins and calculated the absolute value of the deflection amplitude at each time point. The resulting values were divided into 50 bins and p(v) was set to the relative frequencies with which each of the 50 values appeared. We introduced a delay of 15 ms to mimic the response delays observed in cortical responses. For optimization, we calculated the gradient of the mutual information I(S,V) analytically and then used a standard constrained optimization method (fmincon in Matlab). We restricted the average firing probability, Σt p(st)/T, of the optimized neurons to values observed in our recorded data set to prevent neurons from responding at unnaturally high rates (Tkacik et al., 2010). We found that filters optimized for different firing rates differed in a scale factor, but were otherwise identical. We also used a constant set of parameters for the static nonlinearity; again, we found that these did not affect the shape of the optimized filters but rather changed their scale.

Finally, we ran our procedure for filter optimization on a distorted data set. Specifically, we numerically computed the second derivative (accelerations) of the original texture traces described above and re-normalized them to the magnitude of the original traces before simulating peripheral responses to these stimuli and then optimizing the filters. The differentiated traces contain the same frequencies as the original traces, but the higher frequencies are weighted more strongly.

References

  1. 1
  2. 2
  3. 3
    Possible principles underlying the transformations of sensory messages
    1. HB Barlow
    (1961)
    In: WA, Rosenblith, editors. Sensory communication. Cambridge, MA: MIT Press.. pp. 217–234.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
    Fine-scale organization of SI (area 3b) in the squirrel monkey revealed with intrinsic optical imaging
    1. LM Chen
    2. RM Friedman
    3. BM Ramsden
    4. RH LaMotte
    5. AW Roe
    (2001)
    Journal of Neurophysiology 86:3011–3029.
  9. 9
    Texture-induced vibrations in the forearm during tactile exploration
    1. B Delhaye
    2. V Hayward
    3. P Lefèvre
    4. J-L Thonnard
    (2012)
    Frontiers in Behavioral Neuroscience, 6, 10.3389/fnbeh.2012.00037.
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
    Modality maps within primate somatosensory cortex
    1. RM Friedman
    2. LM Chen
    3. AW Roe
    (2004)
    Proceedings of the National Academy of Sciences of the United States of America 101:12724–12729.
    https://doi.org/10.1073/pnas.0404884101
  16. 16
  17. 17
  18. 18
    Response of binaural neurons of dog superior olivary complex to dichotic tonal stimuli: some physiological mechanisms of sound localization
    1. JM Goldberg
    2. PB Brown
    (1969)
    Journal of Neurophysiology 32:613–636.
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
    Critical and maximally informative encoding between neural populations in the retina
    1. DB Kastner
    2. SA Baccus
    3. TO Sharpee
    (2015)
    Proceedings of the National Academy of Sciences of the United States of America 112:2533–2538.
    https://doi.org/10.1073/pnas.1418092112
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
    Cortical neuronal mechanisms in flutter-vibration studied in unanesthetized monkeys. neuronal periodicity and frequency discrimination
    1. VB Mountcastle
    2. WH Talbot
    3. H Sakata
    4. J Hyvärinen
    (1969)
    Journal of Neurophysiology 32:452–484.
  29. 29
    Detection thresholds for stimuli in humans and monkeys: comparison with threshold events in mechanoreceptive afferent nerve fibers innervating the monkey hand
    1. VB Mountcastle
    2. RH LaMotte
    3. G Carli
    (1972)
    Journal of Neurophysiology 35:122–136.
  30. 30
    Frequency discrimination in the sense of flutter: psychophysical measurements correlated with postcentral events in behaving monkeys
    1. VB Mountcastle
    2. MA Steinmetz
    3. R Romo
    (1990)
    The Journal of Neuroscience 10:3032–3044.
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
    Periodicity and firing rate as candidate neural codes for the frequency of vibrotactile stimuli
    1. E Salinas
    2. A Hernandez
    3. A Zainos
    4. R Romo
    (2000)
    The Journal of Neuroscience 20:5503–5515.
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
    Modular distribution of neurons with slowly adapting and rapidly adapting responses in area 3b of somatosensory cortex in monkeys
    1. M Sur
    2. JT Wall
    3. JH Kaas
    (1984)
    Journal of Neurophysiology 51:724–744.
  44. 44
    The sense of flutter-vibration: comparison of the human capacity with response patterns of mechanoreceptive afferents from the monkey hand
    1. WH Talbot
    2. I Darian-Smith
    3. HH Kornhuber
    4. VB Mountcastle
    (1968)
    Journal of Neurophysiology 31:301–334.
  45. 45
    Optimal population coding by noisy spiking neurons
    1. G Tkacik
    2. JS Prentice
    3. V Balasubramanian
    4. E Schneidman
    (2010)
    Proceedings of the National Academy of Sciences of the United States of America 107:14419–14424.
    https://doi.org/10.1073/pnas.1004906107
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
    Stability of rapidly adapting afferent entrainment vs responsivity
    1. BL Whitsel
    2. EF Kelly
    3. KA Delemos
    4. M Xu
    5. PM Quibrera
    (2000)
    Somatosensory & Motor Research 17:13–31.
  53. 53
    On the 1/f noise and non-integer harmonic decay of the interaction of a finger sliding on flat and sinusoidal surfaces
    1. M Wiertlewski
    2. C Hudin
    3. V Hayward
    (2011)
    IEEE World Haptics Conference.
  54. 54
  55. 55

Decision letter

  1. Andrew J King
    Reviewing Editor; University of Oxford, United Kingdom

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for submitting your work entitled "Rate and timing of cortical responses driven by separate sensory channels" for peer review at eLife. Your submission has been favorably evaluated by Eve Marder (Senior editor), a Reviewing editor (Andrew King), and three reviewers, one of whom, Stefano Panzeri, has agreed to reveal his identity.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission.

Summary:

In this study, data-based models were used to examine how information from different mechanoreceptive afferent classes is integrated in the responses of individual cortical neurons. The authors report that cortical activity reflects a combination of rapidly-adapting and Pacinian corpuscle inputs, with the former determining overall firing rate whereas the latter are represented by the precise timing of the activity. Their model results suggest that information about textures is represented in the primary somatosensory cortex by the statistically optimal integration of inputs from these two peripheral receptor channels.

Essential revisions:

Although the three reviewers agree that this is an interesting and important study and that the analysis is rigorous, they have raised the following issues.

1) Are the authors able to say where the inferred convergence of inputs from these two classes of mechanoreceptors takes place – in the cortex or at a subcortical level?

2) Given that model assessment largely revolves around predictive power and predictive power is impacted by reliability, the measurements of reliability mentioned in the Results (second paragraph) should probably be shown. Perhaps the authors could plot reliability vs accuracy. The authors should clarify how they measured reliability.

3) On a related note, factors unrelated to the stimuli, such as the stochasticity of spike generation, lead to a lack of reliability, which the authors comment on. The reader is told how well the model predicts the total variance (about half), but the model is probably performing better than this since some of the variance is due to the lack of reliability in the responses. The authors should estimate how much of the predictable variance (rather than the total) the model accounts for (e.g., Sahani & Linden 2003 NIPS; Ramirez et al. 2014 Nat Neurosci).

4) The authors report (Figure 2A) that most cortical neurons receive "significant contributions" from RA and PC populations. However, if the RA and PC population responses are correlated, a model-fitting algorithm might generate positive weights from a non-connected PC population even though cortical neurons receive input exclusively from RA afferents. It is particularly hard to exclude the possibility that a small (but statistically significant) input might be an artefact. The authors consider this issue in the Materials and methods. They state that the correlations were "low" – but how low? What is the correlation coefficient between the RA and PC population responses as a function of time resolution? Also, the authors use a whitening procedure to attempt to decorrelate (is equation 2 correct? Is the regularizing term not added?), but these procedures are delicate and it is not clear that the correlation issue has been fully dealt with. Further work is needed to address this. It may be useful to use a simulation approach, i.e. simulate cortical neurons with realistic filters that range from single channel to mixed, analyze them with their methods, and test whether their single/mixed nature is accurately recovered. Furthermore, the results (Figure 2A) are based on a statistical test (no details are provided for how this was done or what the critical P value was). Hence, a neuron that counts as convergent might have a strong input from RA together plus a tiny PC input that passes the statistical test, but is functionally irrelevant. Suppose that you define convergent input as where at least 10% of the total input strength comes from each channel, do you still see the same prominent convergence shown in Figure 2A?

5) The authors claim that cortical (time-averaged) firing rates are due to RA, not PC input. They state that their stimuli "drove RA and PC to the exclusion of SA1". On this basis, they consider only RA and PC responses in their study and ignore SA. However, Figures 12 of the supporting reference (Muniak et al., 2007) do show SA responses to both sinusoidal and noise stimuli. Given these data, it does not seem appropriate to ignore the SA afferents. Do the authors have peripheral recordings of SA1 responses? If so, they could apply their model-fitting approach. If they are correct, the weights from the SA population should turn out to be zero.

6) For the analyses shown in Figures 34 where full and RA/PC only models are compared, how were the RA/PC only results obtained? Were the cortical models retrained on only RA/PC input or were the RA/PC contributions in the full model simply set to zero? Given the concern about correlations, the authors should check whether they get the same results using the retraining method if that is not already what they did.

7) The conclusion that the precise timing of S1 neurons is mainly driven by PC input depends on how good the models for the RA/PC afferents are. The authors show a single neuron example (Figure 4A) where removing PC input degrades timing precision, but in Figure 4C (right – the 3 ms resolution case) it is only for a minority of points that PC explains more variance than RA. In fact, for most of the cortical neurons, RA explains more variance. If correct, this seems to undermine the conclusion that PCs drive cortical spike timing. Figure 1—figure supplement 1 shows that the RA afferent model performs poorly at 3 ms resolution. In this case, how is it possible (in Figure 4C right) that they explain so much variance? (It would be useful to define the R squared measure precisely and exactly what "variance" is meant).

8) The authors report that "filters optimized to convey information about texture closely matched the filters derived from measured S1 responses". This is not currently fully convincing. It seems from the Materials and methods that they find the RA and PC filters that maximize mutual information between the stimulus (measured in a 1 ms bin) and the model's response in that bin. They repeat this for 25 different textures. Then they look through these 25 results and pick out the one that looks most like the filters fitted from spike data (this important part of the process is only mentioned in the legend to Figure 5 – it should be in the Materials and methods). There are several concerns with this approach.

A) In classic studies of this kind (e.g. Olshausen & Field, Nature, 2006; Attick, Network, 1992), the optimization is carried out by averaging across a bunch of natural stimuli, not on individual stimuli. It is possible that cherry-picking from 25 separate optimizations could result in spurious matches. It would be much more convincing if the optimization was carried out across the whole set of 25 textures.

B) Even using the texture-picking method, "close" matches are only shown for 3 examples in Figure 5A. Figure 5B-D show similarities between the distribution of width and other parameters between real and optimized filters, but do not show how close these values are on a neuron by neuron basis (or even the means, medians or shapes of these distributions). If there really is systematic, close matching (which would be impressive), there should be a high correlation coefficient between e.g. actual and optimized filter width, etc. across neurons. Some test or comparison with other models is needed to justify the conclusion that filters for the recorded cells are/aren't optimal. Perhaps the authors could distort the textures or shuffle them somehow to see whether they can create "suboptimal" filters.

C) For the mutual information calculation, is it correct that stimulus and response are measured in the same simultaneous time bin? Should response latency not be taken into account?

9) For fitting the afferent models, how did the authors avoid over-fitting? Was cross-validation used when assessing model accuracy in Figure 1figure 1? Why the unusual approach of sequentially fitting to individual stimulus types – why not fit on all 3 stimulus sets at once?

10) It is stated in the first paragraph of the Results section that the stimuli for the cortical experiments were "analogous but not identical" to the peripheral ones. This is potentially a problem, since statistically fitted, nonlinear models cannot be relied on to extrapolate outside their training set. At present, there are insufficient methodological details to understand how different the stimuli really were. Please detail, in the Materials and methods, exactly what the stimuli were in the two cases. Related to this, it is stated in the Materials and methods that the peripheral recordings were done in anesthetized animals. What about the cortical ones?

11) Materials and methods, subsection “Electrophysiology”. The paper includes only 4 PC afferents. The authors need to justify why this small sample size is sufficient or increase the number. Similarly, it is unclear whether multiple peripheral models can be generated using only 4 recorded afferents (subsection “Peripheral model”).

https://doi.org/10.7554/eLife.10450.012

Author response

Essential revisions:

Although the three reviewers agree that this is an interesting and important study and that the analysis is rigorous, they have raised the following issues. 1) Are the authors able to say where the inferred convergence of inputs from these two classes of mechanoreceptors takes place – in the cortex or at a subcortical level?

The model allows us to determine that the integration of signals from different afferent classes is taking place somewhere between the periphery and cortex, but not whether it is taking place in the cuneate nucleus, the thalamus, or in S1 neurons. Not enough is known about the integration properties of the two subcortical structures to make much of an educated guess either. We have now clarified this by adding the following passage to the Discussion section:

"However, it is still an open question whether the integration of input from different mechanoreceptors takes place only at the cortical level, or whether they appear already at the subcortical level in either the cuneate nucleus or the thalamus. Convergence of somatosensory submodalities has been observed in the brainstem of rats (Sakurai et al., 2013) and the same is likely true in primates."

2) Given that model assessment largely revolves around predictive power and predictive power is impacted by reliability, the measurements of reliability mentioned in the Results (second paragraph) should probably be shown. Perhaps the authors could plot reliability vs accuracy. The authors should clarify how they measured reliability.

We now measure reliability in the form of predictable variance (see our response to point 3 for more details) rather than using our own measure and we have revamped the section dealing with model assessment accordingly. A new plot (Figure 1—figure supplement 2A) has been added to show predictable variance as a function of temporal resolution.

3) On a related note, factors unrelated to the stimuli, such as the stochasticity of spike generation, lead to a lack of reliability, which the authors comment on. The reader is told how well the model predicts the total variance (about half), but the model is probably performing better than this since some of the variance is due to the lack of reliability in the responses. The authors should estimate how much of the predictable variance (rather than the total) the model accounts for (e.g. Sahani & Linden 2003 NIPS; Ramirez et al. 2014 Nat Neurosci).

We have switched to using predictable variance for model assessment, which has the added benefit that it makes model assessment more comparable across different temporal resolutions. A new figure (Figure 1—figure supplement 2) showing predictable variance and the fraction of the predictable variance across different temporal resolutions.

The new section on model assessment now reads as follows:

"We found that our model captures the cortical responses well on both coarse and fine timescales as assessed by calculating the correlation coefficients between the predicted and the recorded responses. […] In other words, the model was able to predict neuronal responses with considerable accuracy given their inherent variability."

4) The authors report (Figure 2A) that most cortical neurons receive "significant contributions" from RA and PC populations. However, if the RA and PC population responses are correlated, a model-fitting algorithm might generate positive weights from a non-connected PC population even though cortical neurons receive input exclusively from RA afferents. It is particularly hard to exclude the possibility that a small (but statistically significant) input might be an artefact. The authors consider this issue in the Materials and methods. They state that the correlations were "low" – but how low? What is the correlation coefficient between the RA and PC population responses as a function of time resolution? Also, the authors use a whitening procedure to attempt to decorrelate (is equation 2 correct? Is the regularizing term not added?), but these procedures are delicate and it is not clear that the correlation issue has been fully dealt with. Further work is needed to address this. It may be useful to use a simulation approach, i.e. simulate cortical neurons with realistic filters that range from single channel to mixed, analyze them with their methods, and test whether their single/mixed nature is accurately recovered. Furthermore, the results (Figure 2A) are based on a statistical test (no details are provided for how this was done or what the critical P value was). Hence, a neuron that counts as convergent might have a strong input from RA together plus a tiny PC input that passes the statistical test, but is functionally irrelevant. Suppose that you define convergent input as where at least 10% of the total input strength comes from each channel, do you still see the same prominent convergence shown in Figure 2A?

We followed the excellent suggestion to test our procedure for estimating the filters using simulated neuronal data. We have added a new section in the Materials and methods entitled "Testing the robustness and accuracy of filter estimation" that describes the new tests in detail, along with a new figure (Figure 2—figure supplement 1). In short, this analysis shows that we are able to recover linear filters accurately and that our metrics are only minimally affected by artifacts.

Our criterion for whether a neuron has input from only RAs, only PCs, or both was based on the whether the filters deviated by more than three standard deviations from the noise distribution (measured 100 ms away from the spike) for more than 5 ms. The reviewers are correct that small but significant deviations from the noise distribution (such as introduced by artifacts driven by the RA/PC cross-correlations) might lead to spurious results; indeed, our toy neuron analysis exhibited a number of conditions where this could happen. We thus added another criterion (the one suggested by the reviewers) which states that a neuron can only be classified as receiving input from both classes when the ratio of their filter strengths is not greater than 1:10. In our toy neuron analysis, this additional criterion led to almost no errors in input classification (see new Figure 2—figure supplement 1D). However, this criterion was already met by the neurons in our data set (cf. old Figure 2C) and so the results are unaffected.

Thank you for spotting the typo in equation (2); the regularization term should indeed be added.

The new section in the Materials and methods is reproduced here:

"Testing the robustness and accuracy of filter estimation

Recovering linear filters accurately from data that exhibit correlations is notoriously difficult and error prone. […] In summary, our procedure recovers most filter shapes accurately and the small biases that remain cannot account for the observed effect sizes."

5) The authors claim that cortical (time-averaged) firing rates are due to RA, not PC input. They state that their stimuli "drove RA and PC to the exclusion of SA1". On this basis, they consider only RA and PC responses in their study and ignore SA. However, Figures 12 of the supporting reference (Muniak et al., 2007) do show SA responses to both sinusoidal and noise stimuli. Given these data, it does not seem appropriate to ignore the SA afferents. Do the authors have peripheral recordings of SA1 responses? If so, they could apply their model-fitting approach. If they are correct, the weights from the SA population should turn out to be zero.

Our peripheral data set spanned more frequencies than the cortical data set. Specifically, the lowest frequency used in the cortical data set (and therefore in simulations and analyses of the current paper) was 50 Hz. The figures that the reviewers are referring to both show SA1 responses below 50 Hz and are therefore not applicable for the current study. A better figure to look at in reference (Muniak et al., 2007) would be Figure 4 which shows some (weak) SA1 activation even for frequencies above 50 Hz given rather high stimulus amplitudes (with thresholds generally around just below 100 microns). A valid concern is that these SA1 responses (which are currently not accounted for in the model) could bias the estimation of the RA and PC filters. We therefore ran our filter estimation again, this time restricting our data set to noise traces with an RMS amplitude of <50 microns; this conservative threshold all but excludes any possible contribution of SA1 afferents. While the resulting filters were slightly noisier (because fewer data were used to estimate them), all metrics such as the fraction of neurons with convergent input, the filter magnitude ratios, the filter widths, and the excitation/suppression indices almost exactly matched the ones calculated from the full data set. We thus conclude that any SA1 contribution does not affect our results. While we would have liked to also derive SA1 filters using our model, the fact they fire so sparsely to these stimuli means that their contribution to the S1 firing rates is minimal.

6) For the analyses shown in Figures 34 where full and RA/PC only models are compared, how were the RA/PC only results obtained? Were the cortical models retrained on only RA/PC input or were the RA/PC contributions in the full model simply set to zero? Given the concern about correlations, the authors should check whether they get the same results using the retraining method if that is not already what they did.

We were already using refitted filters for the RA/PC-only models. Compared to the full models, the filters were essentially the same. The only consistent difference was in the fitted nonlinearity, which generally needs to be steeper when only RA or PC input is used in order to achieve the same firing rates as the full model. We have now clarified this in the Materials and methods section.

7) The conclusion that the precise timing of S1 neurons is mainly driven by PC input depends on how good the models for the RA/PC afferents are. The authors show a single neuron example (Figure 4A) where removing PC input degrades timing precision, but in Figure 4C (right – the 3 ms resolution case) it is only for a minority of points that PC explains more variance than RA. In fact, for most of the cortical neurons, RA explains more variance. If correct, this seems to undermine the conclusion that PCs drive cortical spike timing. Figure 1—figure supplement 1 shows that the RA afferent model performs poorly at 3 ms resolution. In this case, how is it possible (in Figure 4C right) that they explain so much variance? (It would be useful to define the R squared measure precisely and exactly what "variance" is meant).

First, the main goal of our analysis was to show that PC input is needed to explain the precise timing of cortical neurons. As such, the important comparison is not whether PC input explains the cortical responses better than RA input, but that RA input alone is not enough to explain the cortical responses at fine temporal resolutions. Second, not all cortical neurons exhibit precise spike timing (see Harvey et al., PloS Biol, 2013). Indeed, cortical neurons that exhibit more precise spike timing also depend more on PC input, as we now show in the new Figure 4—figure supplement 1.

We have added a sentence in the Results section to clarify this point:

"We found that only PC input can account for high frequency fluctuations in the cortical responses of many neurons, while RA input is sufficient at coarser temporal resolutions for almost all cortical neurons (Figure 4C). In fact, the more precise the responses of a cortical neuron were, the more of its response variance was explained by PC input (Figure 4—figure supplement 1)."

8) The authors report that "filters optimized to convey information about texture closely matched the filters derived from measured S1 responses". This is not currently fully convincing. It seems from the Materials and methods that they find the RA and PC filters that maximize mutual information between the stimulus (measured in a 1 ms bin) and the model's response in that bin. They repeat this for 25 different textures. Then they look through these 25 results and pick out the one that looks most like the filters fitted from spike data (this important part of the process is only mentioned in the legend to Figure 5 – it should be in the Materials and methods). There are several concerns with this approach.

A) In classic studies of this kind (e.g. Olshausen & Field, Nature, 2006; Attick, Network, 1992), the optimization is carried out by averaging across a bunch of natural stimuli, not on individual stimuli. It is possible that cherry-picking from 25 separate optimizations could result in spurious matches. It would be much more convincing if the optimization was carried out across the whole set of 25 textures.

Olshausen & Field and similar approaches ask how information can be efficiently encoded in a population of neurons, and they optimize for both information transmission and sparseness. We, however, are concerned with the response properties of individual neurons in this paper, and ask whether we can explain the way single neurons integrate information from different inputs by a simple principle of information maximization. The reviewers are right to wonder whether filters similar to the cortical ones might not be trivially obtained by this method, and we have now added a new analysis that we think addresses this question (see below). We have also clarified the rationale for our approach in the Results section:

"We focused on natural textures, which are ideally suited to address this question. […] To address this question, we simulated RA and PC responses to different textures scanned across the skin and optimized cortical RA and PC filters for information transmission about skin oscillations."

B) Even using the texture-picking method, "close" matches are only shown for 3 examples in Figure 5A. Figure 5B-D show similarities between the distribution of width and other parameters between real and optimized filters, but do not show how close these values are on a neuron by neuron basis (or even the means, medians or shapes of these distributions). If there really is systematic, close matching (which would be impressive), there should be a high correlation coefficient between e.g. actual and optimized filter width, etc. across neurons. Some test or comparison with other models is needed to justify the conclusion that filters for the recorded cells are/aren't optimal. Perhaps the authors could distort the textures or shuffle them somehow to see whether they can create "suboptimal" filters.

We have followed the reviewers' suggestions and now derived optimal filters also for a distorted data set, namely one where the usually observed power law between frequency and amplitude of skin oscillations has been broken. We find that the optimal filters derived for these non-natural scenes do not match the cortical ones as well as those derived from natural scenes. Furthermore, we do not want to make any claims that the optimal filters match the cortical ones on a neuron-by-neuron basis, but simply that similar filter properties are observed in the optimal and the actual filters, when we former ones are optimized on a data set containing natural scenes. We have added the following text in the Results section:

"To test whether our fitted and optimized filters might be an artifact of our approach rather than reflect the integration properties of S1 neurons, we also optimized filters on a distorted data set. […] The PC signal was weighed more heavily, was integrated over longer timescales, and was mostly excitatory, in contrast to what we observed in both the cortical filters and the ones optimized on our natural texture data set."

Additionally, we have added the following notes in the Materials and methods section:

"Finally, we ran our procedure for filter optimization on a distorted data set. […] The resulting traces are thus similar but differ in one important aspect (namely frequency composition) from the natural traces."

C) For the mutual information calculation, is it correct that stimulus and response are measured in the same simultaneous time bin? Should response latency not be taken into account?

Our analysis assumes a fixed response latency of 15 ms, comparable to what is observed in the cortical responses. This is now clarified in the Materials and methods.

9) For fitting the afferent models, how did the authors avoid over-fitting? Was cross-validation used when assessing model accuracy in Figure 1figure 1? Why the unusual approach of sequentially fitting to individual stimulus types – why not fit on all 3 stimulus sets at once?

The afferent models comprise few parameters so overfitting is not a serious concern. We have validated the approach extensively in previous published work, in which we systematically cross-validated the fits. Furthermore, the responses of peripheral afferents are very stereotyped and have been extensively described in the literature, so it is easy to spot erroneous response patterns.

Our approach of sequentially alternating between different data sets was implemented mainly to speed up the convergence of the parameters. Similar approaches, where parameters are fitted on alternating subsets of the data, are used in the machine learning community (for example in stochastic gradient descent algorithms).

10) It is stated in the first paragraph of the Results section that the stimuli for the cortical experiments were "analogous but not identical" to the peripheral ones. This is potentially a problem, since statistically fitted, nonlinear models cannot be relied on to extrapolate outside their training set. At present, there are insufficient methodological details to understand how different the stimuli really were. Please detail, in the Materials and methods, exactly what the stimuli were in the two cases. Related to this, it is stated in the Materials and methods that the peripheral recordings were done in anesthetized animals. What about the cortical ones?

Both peripheral and cortical stimulus parameters have been described in previously published papers that are referenced in the Materials and methods section. In both cases, simple sinusoids as well as bandpass noise traces were presented. While the exact frequencies and amplitudes differed somewhat between the peripheral and cortical data sets, there was broad overlap between the two sets. Thus, the peripheral models were only used to simulate responses to stimuli that were reasonably close to those encountered during their training. The models were further validated by assessing how well they could replicate well-known response properties of peripheral afferents, such as phase-locking and entrainment plateaus, some examples of which are shown in Figure 1—figure supplement 1. As mentioned above, we documented a much more systematic validation of these models in previous papers. The cortical recordings were from awake animals. This is now clarified in the text.

We have added the following sentence to the 'Stimuli' section in the Materials and methods:

"While the exact frequencies and amplitudes differed somewhat between the peripheral and cortical data sets, the two stimulus sets were largely overlapping. […] Individual stimuli in this data set were 1 s long and were repeated 5 times (frozen noise), yielding a total of about 22 minutes (16 for noise and 6 for the sinusoids) of data from each neuron."

11) Materials and methods, subsection “Electrophysiology”. The paper includes only 4 PC afferents. The authors need to justify why this small sample size is sufficient or increase the number. Similarly, it is unclear whether multiple peripheral models can be generated using only 4 recorded afferents (subsection “Peripheral model”).

We believe that 4 PC afferent models are sufficient for the following reason. We have more models available for RA afferents and can therefore test how population size affects our prediction accuracy. We find that using only one or two models leads to decreased performance (accounting for 40-50% of explainable variance at a resolution of 100 ms), but that accuracy levels off when at least 4 models are used (57%) and does not increase further with more afferents added (58% using all 14 RA models).

https://doi.org/10.7554/eLife.10450.013

Article and author information

Author details

  1. Hannes P Saal

    Department of Organismal Biology and Anatomy, University of Chicago, Chicago, United States
    Contribution
    HPS, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    hsaal@uchicago.edu
    Competing interests
    The authors declare that no competing interests exist.
  2. Michael A Harvey

    1. Department of Organismal Biology and Anatomy, University of Chicago, Chicago, United States
    2. Department of Neurological Surgery, Vanderbilt University Medical Center, Nashville, United States
    Contribution
    MAH, Conception and design, Acquisition of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  3. Sliman J Bensmaia

    1. Department of Organismal Biology and Anatomy, University of Chicago, Chicago, United States
    2. Committee on Computational Neuroscience, University of Chicago, Chicago, United States
    Contribution
    SJB, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    sliman.bensmaia@gmail.com
    Competing interests
    The authors declare that no competing interests exist.

Funding

National Science Foundation (IOS 1150209)

  • Sliman J Bensmaia

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Erika Dunn-Weiss for getting this work started and Benoit Delhaye and Justin Lieber for comments on a previous version of this manuscript.

Ethics

Animal experimentation: This study was performed in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. All of the animals were handled according to protocol #72042, which was approved by the institutional animal care and use committee (IACUC) of the University of Chicago and the Johns Hopkins University.

Reviewing Editor

  1. Andrew J King, University of Oxford, United Kingdom

Publication history

  1. Received: July 29, 2015
  2. Accepted: December 8, 2015
  3. Accepted Manuscript published: December 9, 2015 (version 1)
  4. Version of Record published: February 2, 2016 (version 2)

Copyright

© 2015, Saal et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,020
    Page views
  • 325
    Downloads
  • 10
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

  1. Further reading

Further reading

    1. Cell Biology
    2. Computational and Systems Biology
    Harry Nunns, Lea Goentoro
    Research Article
    1. Computational and Systems Biology
    2. Ecology
    Farnoush Farahpour et al.
    Research Article Updated