Johanna Norrlid
Jonas M. D. Enander
Hannes Mogensen
Henrik Jörntell- Neural Basis of Sensorimotor Control, Department of Experimental Medical Science, Lund University, Lund, Sweden
The brain has a never-ending internal activity, whose spatiotemporal evolution interacts with external inputs to constrain their impact on brain activity and thereby how we perceive them. We used reproducible touch-related spatiotemporal sensory inputs and recorded intracellularly from rat (Sprague-Dawley, male) neocortical neurons to characterize this interaction. The synaptic responses, or the summed input of the networks connected to the neuron, varied greatly to repeated presentations of the same tactile input pattern delivered to the tip of digit 2. Surprisingly, however, these responses tended to sort into a set of specific time-evolving response types, unique for each neuron. Further, using a set of eight such tactile input patterns, we found each neuron to exhibit a set of specific response types for each input provided. Response types were not determined by the global cortical state, but instead likely depended on the time-varying state of the specific subnetworks connected to each neuron. The fact that some types of responses recurred indicates that the cortical network had a non-continuous landscape of solutions for these tactile inputs. Therefore, our data suggest that sensory inputs combine with the internal dynamics of the brain networks, thereby causing them to fall into one of the multiple possible perceptual attractor states. The neuron-specific instantiations of response types we observed suggest that the subnetworks connected to each neuron represent different components of those attractor states. Our results indicate that the impact of cortical internal states on external inputs is substantially more richly resolvable than previously shown.
Introduction
Behavioral, mental, and perceptual functions of the neocortex depend on its internal state. An unresolved issue is how internal states affect the responses to a given external input. A state in the brain can be described as the combination of activity in all of its neurons (Spanne and Jorntell, 2015). Since the neurons make synaptic connections with each other, their activities are not independent, which is reflected in reports of constrained “realms” of possible response combinations in populations of neurons (Luczak et al., 2009; Golub et al., 2018). External inputs to the neocortical circuitry, which generate spatiotemporal patterns of activation arising in the multitude of sensors throughout our bodies, further constrain the space of possible neuronal responses (Luczak et al., 2009). An important aspect of perception, and the foundation of illusions, is that a response is determined not only by the quality of the sensory input but also depends on the physiological circuitry structure, which results from the internal state of the cortex at the time the input arises (Arieli et al., 1996; Fiser et al., 2004; Curto et al., 2009; Berkes et al., 2011). However, although the internal state is a high-dimensional latent state (Stringer et al., 2019a,b) more detailed physiological characterization of the state-dependent influence on cortical circuitry responses to naturalistic sensory input is so far lacking. This is not surprising given that a direct demonstration would require a precise estimate of what the experimental subject is thinking and how that thinking is instantiated in the circuitry, at the time of the stimulus delivery. Otherwise, the internally generated constraints become an uncontrolled variable, which will be interpreted as internal system noise in sensory-evoked responses. Here we aimed to deduce information about the character of the interactions between the time-evolving cortical internal state and external inputs, based on an analysis of the detailed nature of the cortical neuron responses evoked by inputs consisting of several alternative fixed spatiotemporal patterns of tactile sensory activation, each delivered at a high number of repetitions.
It is difficult to achieve exactly reproducible sensor activation patterns in a living organism where the sensors are located in compliant or movable tissue (Hayward et al., 2014) and their exact location or tuning relative to external stimuli are moreover subject to uncontrollable efferent control by the brain, for example via muscles. Therefore, we used electrical intracutaneous stimulation to deliver a set of reproducible but also perceptual rich sensory input patterns, mimicking active touch, to the brain (Oddo et al., 2017; Figure 1A). In order to minimize system noise caused by uncontrollable movements and internal thought processes unrelated to the stimuli (see “Materials and Methods” section for further details), and in order to make the rats accept long–term stimulation of the skin electrodes, we used light anesthesia. The intracellularly recorded signal represents the integrated synaptic inputs from tens of thousands of neurons (out of 25,000,000 neurons in the rat neocortex; Bandeira et al., 2009). As these neurons, by being connected to the recorded neuron, are potentially representing activity in different subnetworks, the intracellular signal is a read-out of the instantiation of the current brain state that is specific to the subnetworks connected to that neuron. Here we find that each out of eight pre-defined tactile input patterns generates a wide range of response types that are specific to each neocortical neuron. The findings suggest a multi-component but structured and subnetwork-specific interplay between internal cortical states and sensor input patterns.
Figure 1. Electrical skin stimulation and general properties of the intracellular responses. (A) Locations of the four pairs of intracutaneous needle electrodes (“channels”) inserted in digit 2 to stimulate the tactile primary afferents. The eight spatiotemporal patterns were re-used from a previous publication (Oddo et al., 2017). (B) Left: neuronal recordings were made in an exposed cortical area of 4 by 2 mm, located in the center of the photo. An electrocorticography (ECoG) surface electrode was placed on the brain. Right: a schematic illustrating the location of the exposed cortical area (red, dashed rectangle). (C) Left: a stained layer III pyramidal neuron from one of the recordings. The white dashed rectangle represents the zoomed-in area shown to the right. (D) Intracellular traces illustrating the spontaneous activity and responses evoked by the stimulation patterns (gray horizontal bars). “+” indicates occasional spikes. The arrow indicates a spontaneous activity pattern resembling an evoked response. Colored lines indicate the times of individual stimulation pulses in pattern F(inf) and F5, respectively. Recording from Neuron #7. Stimulus artifacts are truncated for clarity.
Materials and Methods
Experimental Design.
Ethical Approval
All animal experiment procedures in the present study were in accordance with institutional guidelines and were approved in advance by the Local Ethics Committee of Lund, Sweden (permit ID M118-13 and M13193-2017). All experiments were made using acute preparations under general anesthesia.
Surgical Procedures
Adult male Sprague Dawley rats of male sex (N = 10, weight 240–383 g, Taconic) were briefly sedated with isoflurane (3%, 1–2 min), anesthetized with an intraperitoneal injection (40 mg/kg ketamine, 4 mg/kg xylazine) and maintained under anesthesia with a continuous infusion (ketamine and xylazine in a 20:1 ratio, approximately 5 mg/kg ketamine/h) administered through an intravenous (IV) catheter inserted into the right femoral vein. A hemicraniectomy (approximately 4 by 2 mm) exposed the area of the right somatosensory cortex (Figure 1B). An electrocorticography (ECoG) electrode was positioned on the surface of the brain in order to continuously monitor the depth of anesthesia by ensuring the presence of sleep spindles mixed with epochs of more desynchronized activity, a characteristic of sleep (Niedermeyer and da Silva, 2005). The level of anesthesia was additionally characterized by an absence of withdrawal reflexes in response to noxious pinches of the hind paw. The type of anesthesia used here has little disruptive effect on the physiological network structure at short time spans (in the order of 100 s of ms) as judged by the preservation of the order of neuronal recruitment of neocortical neurons in spontaneous brain activity fluctuations (Upstates, recordings obtained using multi-electrode arrays in the rat) and stimulus-evoked responses (Luczak and Bartho, 2012). Anesthesia drags down the overall activity in the neocortical network (Constantinople and Bruno, 2011), though, and in general, can be expected to make those networks function with a lower degree of precision. Nevertheless, for the present study, the method of stimulus delivery (see below) would not be accepted by the awake animal, and meeting the requirement for long–term intracellular recordings was further facilitated by the anesthesia. To create further mechanical stability, and to protect the brain tissue from dehydration, an agarose gel was applied to cover the exposed part of the cortex. After finishing the neuronal recordings the animal was sacrificed with pentobarbital (140 mg/kg IV).
Artificial Touch Inputs
In order to achieve as realistic spatiotemporal patterns of tactile sensor activation as possible, while preserving the aim of a high reproducibility of the patterns, we used an artificial fingertip equipped with four neuromorphic sensors to generate a set of spatiotemporal patterns of skin activation to be used in the electrical interface with the rat skin. The procedures have previously been described in greater detail and the patterns used here were the same as before (Oddo et al., 2017; Genna et al., 2018; Enander et al., 2019). As discussed in this previous work, the artificial fingertip allowed synthesizing spatiotemporal patterns of skin sensor activation at quasi-natural firing rates that follow a natural overall temporal modulation, or “envelope”, that the biological skin sensors are known to display under dynamic indentation (Jenmalm et al., 2003). This aspect is important because the circuitry of the cortex can be expected to have experienced many events with similar envelopes of tactile afferent activity and is therefore likely to have adapted its circuitry structure to effectively process variations of that overall pattern of activity modulation (Berkes et al., 2011; Oddo et al., 2017). By delivering this input electrically to the primary afferents at local skin sites we bypassed the potentially variable step of skin sensor activation that occurs even with highly controlled mechanical skin stimulation. Hence, the decoding capacity of cortical neurons could be studied in relative isolation from such variability of skin sensor activation. Just like natural tactile inputs, the input we provided can be expected to be distributed and processed through the neuronal networks in the cuneate nucleus, thalamus, and neocortical circuitry before it reaches the neurons we recorded from. Hence, the measured responses are bound to reflect at least in part the natural processing mechanisms of the brain. Accordingly, in humans, electrical nerve stimulation with a much lower resolution than in the present set of experiments are known to generate sensory impressions that are in part perceived as unnatural but also generate diversified and meaningful tactile percepts (Tan et al., 2014; Oddo et al., 2016).
The probes by which the sensory activation patterns of the artificial fingertip were generated and the explanation of their names were thoroughly introduced in Oddo et al. (2017), but a brief explanation follows. The bionic fingertip was moved against objects/probes of different roundness. The numbers within the labels of each pattern (Figure 1A) are exponents in an exponential function that approximately correspond to the radii of the curvatures of the different objects/probes. A low number indicates a higher sharpness, a high number indicates a lower sharpness. Therefore, “infinity” means a flat probe. The letter of the label indicates the adaptive tuning of the biomorphic sensors (S, F; for slow and fast, respectively). These spatiotemporal stimulation patterns were delivered via four pairs (or channels) of intracutaneous needle electrodes (Figure 1A) in a pre-defined random order, where each pattern lasted for less than 340 ms and the consecutive deliveries of the patterns were separated by 1.8 s in order to allow a relaxation of the cortical activity between consecutive deliveries of stimulation patterns (Oddo et al., 2017). Each of the eight patterns was delivered 100 times, except for four out of the total 13 neurons for which the recording was lost after 36, 47, 50, and 80 repetitions, respectively. In addition, for each of the four stimulation channels used, we delivered the same number (up to 100) of repetitions of isolated single-pulse stimulations. These isolated single-pulse stimulations were delivered in chunks of five stimulations from the same channel separated by 300 ms from each other. Thus, for each channel, there were 20 such chunks intermixed with the stimulation patterns in random order. Each channel was stimulated at 0.5 mA with a pulse width of 0.14 ms, which is higher than the threshold of about 0.2 mA reported for tactile afferents using this type of electrocutaneous stimulation (Bengtsson et al., 2013), but lower than the threshold for activating nociceptive afferents (Ekerot et al., 1987).
Neural Recordings
We made recordings in the region of the primary somatosensory cortex (S1) of the forepaw, as estimated by the focus of the local field potentials (measured between layers III and V, corresponding to the depths of maximum field potential negativity recorded in each track). The coordinates of this region were −1.0–1.0 mm relative to bregma and 3.0–5.0 mm lateral to the midline (Figure 1B). Individual neurons were recorded with patch-clamp pipettes in the intracellular, whole-cell current-clamp mode. Patch-clamp pipettes were pulled from borosilicate glass capillaries to 6–15 MOhm using a Sutter Instruments (Novato, CA) P-97 horizontal puller. The composition of the electrolyte solution in the patch pipettes was (in mM) potassium-gluconate (135), HEPES (10), KCl (6.0), Mg-ATP (2), EGTA (10). The solution was titrated to 7.35–7.40 using 1 M KOH. During slow advancement of the recording electrode (approximately 0.3 μm per second) made with positive pressure applied, electrode tip resistance and responses evoked by electrical skin stimulation were continuously monitored to identify encounters with neurons. Once encountered, the positive pressure was changed to negative pressure, and a weak hyperpolarizing current was applied with the aim of obtaining a GigaOhm seal on the neuron. Successful access to the intracellular signal of the neuron, following additional negative pressure once the seal was established, was followed by the release of pressure and the start of the data collection. Using a weak hyperpolarizing bias current, neurons were prevented from spiking. All intracellular data were digitized at 100 kHz using CED 1401 mk2 hardware and Spike2 software (Cambridge Electronics Devices, CED, Cambridge, UK). The criteria used for inclusion of an intracellular recording, or the time period of such a recording to be included in the analysis, were a stable membrane potential of < −55 mV in Down states, a spike amplitude of >25 mV before and after the termination of the protocol and a peak-to-peak difference between the Up and Down states of >10 mV. All neurons recorded were putatively located within layer III-V based on the recording depth measured from the cortical surface (Narayanan et al., 2017). For identification of neuron identity, in addition to depth, we used the nature of the firing during spontaneous activity (i.e., if the neuron was fast-spiking, bursting, and what duration and intensity of bursts the neuron displayed). Before entering the intracellular mode all neurons recorded here exhibited infrequent bursts of two or three spikes but had an absence of longer bursts or sustained periods of high firing. Based on this criterion, they were considered to be pyramidal cells rather than interneurons (Luczak et al., 2009). Three out of the 13 successfully recorded neurons were also stained with neurobiotin and histologically recovered. They were thereby confirmed to be pyramidal neurons located in layer III (Figure 1C).
Statistical Analysis
Post-processing—General
The neuronal recording signal was imported from Spike2 to Matlab (2016a, Mathworks), where it was low-pass filtered using a moving average of over 50 μs, i.e., five samples width given a 100 kHz sampling rate. Stimulation artifacts were removed using a combination of adaptive filtering and blanking of artifacts. Occasional remaining intracellular spikes were removed using adaptive filtering, with a recursive fitting algorithm that created a generic spike shape for the neuron (Mogensen et al., 2019), which could be subtracted from the signal at all occurrences of the spike. This allowed for excitatory postsynaptic potential (EPSP)-like events to be detected also when the membrane potential was influenced by spiking activity. Since all evoked responses were analyzed by visual inspection, there was a quality check of the intracellular signal throughout the recording period.
Post-processing—Clustering Algorithm for Separation of Response Types
On visual inspection, the time-voltage curves of the intracellular responses to repetitions of a given stimulation pattern appeared to fall into certain classes, or types, where the differences between responses of different types appeared to be larger than the variability of different responses within any given type (Figure 2A).
Figure 2. Different response types evoked by the same stimulation pattern. (A) Each panel i-iii shows, for a sample cell (Neuron #5), examples of qualitatively different superimposed raw data traces, evoked by the same stimulation pattern (S5). Note that the stimulation pattern (vertical colored lines) outlasted the traces. Stimulus artifacts are reduced for clarity. (B) The six identified response types for all responses evoked by the 100 repetitions of this stimulation pattern in this cell are shown separately. The mean response of each type is indicated by the thick blue curve, the threshold difference by the light blue span, and the raw traces are indicated in gray, alongside the stimulation pattern. (C) Sensitivity analysis of the parameter settings (for “Threshold” and “Overshoot”, respectively, see “Materials and Methods” section) of our response type classification/clustering method illustrated for two different neurons. The impact of the parameter settings on the fraction of responses ending up as ungrouped (proportion ungrouped) and on the number of response types identified vs. the maximal number of response types identifiable within the parameter space (proportion grouped) are illustrated for both sample neurons (Neuron #5 is also illustrated in panels A,B). Note that the parameter settings used for the main analysis are indicated by the cross-hairs. (D) The Kruskal Wallis test result of the separation between the members of each response type from the members of other response types evoked by the same stimulation pattern in the same neuron, plotted as a separate time curve for each stimulation pattern (N = 104; gray traces). The blue trace represents the average across all stimulation patterns. The corresponding plot for the H statistic is shown in the diagram below. (E) Time evolution of the specificity of the responses of each response type. The black curve is the fraction of all response types where specificity could be detected (at p < 0.05, Kruskal-Wallis test). The red curve shows the corresponding fraction for responses with shuffled response type labels (within each stimulation pattern, for each cell). Gray bars show the durations of the eight stimulation patterns.
To quantitatively evaluate this potential clustering, we used an unsupervised clustering method. This clustering method was merely based on the amplitudes of the intracellular membrane potential at each sample time point. The clustering of different responses was based on the accumulated differences/similarities of these amplitude values across the whole series of sample time points that constituted the response. Thereby, it is based on clear neurophysiological metrics and in this respect differs from traditional machine learning algorithms, in which it is typically not clear which metrics underlie the clustering.
The clustering method defined both the number of response types and the number of members/individual responses that belonged to each response type. The clustering method was also capable of reporting zero response types if no response types could be identified according to the criteria of the clustering method, and hence it ran a low risk of clustering noise compared to some machine learning algorithms. The following procedure was used to sort the intracellular response curves into types:
1. For each cell, the 350 ms time-voltage curves from the onset of stimulation for each of the (up to) 100 repetitions of a given stimulation pattern were compared.
2. To remove high-frequency fluctuations the responses were low-pass filtered (with a 1 ms wide moving average) and re-sampled to 1,000 Hz. In order to focus on the temporal shape of the responses, a moving average of 100 ms was subtracted from the signal to remove its offset and the response voltage was then normalized to a 1.0–0.0 range based on the highest peaks and deepest troughs for each 350 ms time-voltage curve. This was made to ensure that the method captured the shape of the response curve over time.
3. For each neuron/stimulation pattern, each time-voltage curve was subjected to pairwise comparisons against all other time-voltage curves evoked by the same stimulation pattern in the same neuron. The evaluation of similarity between each pair of time-voltage curves was done on a per-sample time point basis (i.e., 350 sample points per time-voltage curve), where the difference for each sample point was calculated. If the difference for a sample point was below a threshold value (this “Threshold difference” was set to 0.13 normalized units, see Figure 2C for the sensitivity analysis of how this threshold was chosen), the difference was set to zero for that sample point. If the difference was above the threshold value, the overshoot was calculated as the absolute value of the actual difference above the threshold value for that sample point. If the mean overshoot across the 350 samples fell below a threshold (this “Overshoot threshold” was set to 0.08 normalized units, see Figure 2C) the two time-voltage curves were classified as the same response type. This procedure helped reducing sensitivity to high-amplitude transient membrane potential shifts while emphasizing persistent low amplitude differences between time-voltage curves in the classification/clustering.
4. This procedure was repeated so that all time-voltage curves were compared with all others of the neuron/stimulation pattern (i.e., typically 100 responses). This resulted in the identification of several types of responses. The members of the response type with the largest number of members were removed from the set of 100 responses, and the sorting (steps 1–4) was repeated until there were no remaining responses left to sort.
5. Last, for the set of detected response types, for each response type with less than five members, the members were categorized as belonging to the “ungrouped” response type and that response type was no longer a valid response type.
Sensitivity Analysis for the Selection of Parameter Values for the Clustering Method
The clustering method described above is unsupervised but depends on the choice of parameter values for the Threshold difference and the Overshoot threshold (step 3 above). The Threshold difference parameter is the absolute value of the distance from the voltage of the reference response (in normalized units) on a per-sample point basis (Figure 2B, light blue area). The Overshoot threshold is the mean of that distance across all sample points.
To test the sensitivity of the method to the choice of parameter values, we explored the resulting outcome (number of response types identified, fraction of the responses that fell into the ungrouped category, and the accuracy by which the identified response types could be identified) across parts of the parameter space for two neurons’ sets of responses to different stimulation patterns. This exploration is visualized in Figures 2C, 3B, which suggested that our choices of parameter values were located in the middle of a smooth landscape of outcomes and that our clustering algorithm was not “brittle” in the sense that gradual changes in the two parameters Threshold and Overshoot did not result in dramatically different outcomes.
Figure 3. Specificity of the different response types evoked by the same stimulation patterns in the same neuron. (A) Confusion matrix for the members of each identified response type for the sample stimulation pattern and neuron illustrated in Figures 2A,B. Decoding accuracy for each response type is indicated by the color code indicated in the calibration bar. (B) Sensitivity analysis of the parameter settings of our response type clustering method illustrated for two different neurons, here with respect to the average response type decoding accuracy (see A) but across all stimulation patterns. (C) Confusion matrix for all response types identified for all stimulation patterns for this sample neuron. Same color code as in (A). (D) To illustrate the specificity of the response types recorded across all eight stimulation patterns, this plot indicates the decoding accuracy across all stimulation patterns for each neuron. The figure also displays the chance level for each neuron, and the decoding accuracy when the responses recorded for each neuron were shuffled.
This is further illustrated in Supplementary Figure 1, where we analyze the consistency of the members of the response clusters for adjacent parameter values. The clusters defined by a parameter setting for an element in a parameter grid (Supplementary Figure 1) are compared against the clusters defined by each of the eight neighboring elements (with single-step increments or decrements in one or both parameters) in turn. Each comparison reports the proportion of the responses/members of one cluster in the center element that are cluster members also of a single cluster in the neighboring element, i.e., if the center element has one cluster with members [0, 1] and a neighboring element has two clusters with members [0, 1, 2] and [3, 4, 5] that comparison would yield the “agreement” value 1.0 since all members of the center element are contained within the first cluster of the neighbor element. The reverse comparison would instead have yielded an agreement value of 0.33 (0.66 for the first cluster and 0 for the second cluster, 0.33 on average). Then the comparison is completed for all eight surrounding elements (except for edges and corners of the grid) and a mean “cluster agreement” value is calculated for the center element. Supplementary Figure 1 also reports the max “cluster agreement”, which is the mean of the best cluster agreement values for all eight comparisons.
Statistical Evaluation of the Identified Response Types
We first tested whether the response types of each stimulation pattern, as identified by the clustering method above, were different from each other. For each sample point (350 sample points per response, see step 2 above), we tested if the distributions of the amplitudes in each cluster/response type were different from each other using the Kruskal-Wallis test. For each set of response types for each stimulation pattern, this procedure was repeated for each sample point up to 1,200 ms after the onset of the stimulation pattern, i.e., for a much longer time than the duration of the stimulation pattern, which yielded a time curve of the p-value that could be visualized individually for each of the 104 stimulation patterns (eight patterns for 13 neurons) in the whole data set (Figure 2D). By defining a threshold p-value of < 0.05, we could also display the results of this analysis as a time curve of the fraction of the response types fulfilling this criterion (Figure 2E). In the latter analysis, we also tested whether a shuffling of the response type labels between all the responses evoked by a specific stimulation pattern in a specific neuron resulted in a collapse of this fraction down to the 5% level, which is the expected value for randomly distributed responses (Figure 2E). This shuffled analysis also functioned as a control for the multiple comparisons problem.
Evaluation of the Specificity of the Identified Response Types for Each Stimulation Pattern in the Same Neuron
We also evaluated the separability of the whole membrane voltage-time curve vectors for each of the identified response types, for each stimulation pattern in each neuron separately, as well as their separability vs. the “ungrouped” responses, using a combination of principal component analysis (PCA) and k-nearest neighbors (kNN)-classification (Figure 3A).
We started by using the mean membrane voltage-time curve vector of each response type (see step 2 above) to calculate the Principal Components (PCs). The number of PCs used was the number required to explain at least 95% of the total variance of the mean signals. Finally, we used the principal component coefficients to transform each individually recorded membrane voltage-time curve vector from the time domain to the principal component domain, reducing the dimensionality of each response from M = 350 (sample time points) to N = [1–6] (PCs).
The position of each response in PC space was then used for verification of the accuracy of the classification/type identity of each response. In order to determine the separability of the detected response types and their possible confusion with the “ungrouped” responses, we used a kNN classification procedure. Half of the responses were randomly selected as the training set. For each response belonging to the test set, we identified the four closest responses in the training set by calculating the Euclidean distance in PC space. The response was classified as belonging to the same response type as the relative majority of the four neighbors, where the classification either matched the response type assigned to it by the clustering method or not. We performed 40 iterations of the classification procedure, each with a different training set, and averaged the fraction of correctly classified responses in each iteration to get the mean correct classification value for each response type. The results of the ensuing kNN decoding are visualized in confusion matrices, which were also used to extract the mean decoding accuracy and the F1 score.
Evaluation of the Specificity of the Response Types Across All Stimulation Patterns in the Same Neuron
To evaluate if the identified response types were specific to the stimulation pattern, we again used PCA and kNN-classification. Here, the PCs were calculated based on the mean membrane voltage- time curve vector of each response type across all eight stimulation patterns. Because the total number of responses considered in this case was in the order of 800 rather than 100 (as above), the dimensionality of the responses was reduced to N = [8–40] PCs, rather than [1–6] PCs. Also, the number of neighbors evaluated for the kNN classification were in this case nine (N = 9).
Evaluation of the Specificity of the Response Types Across the Same Stimulation Pattern in Different Neurons
We also used PCA and kNN decoding to evaluate if the response types were specific to the neuron, for each specific stimulation pattern (see for example Figure 4B), using the same type of approach as above. Here, the PCs were calculated based on the mean membrane voltage-time curve vector of each response type, for a given stimulation pattern, across all thirteen neurons.
Figure 4. Different response types evoked by the same stimulation pattern in different neurons. (A) Averages of intracellular responses of neuron #5 (N = 100 responses) and neuron #11 (N = 80 responses) to stimulation pattern S20 (top). Vertical colored lines indicate the individual stimulation pulses of that stimulation pattern. Similar display for the average responses to two other input patterns (middle, bottom). The responses of Neuron #5 to pattern S5 are also illustrated in Figures 2A,B, 3A. (B) Confusion matrix for all response types identified for all neurons for stimulation pattern F∞. Neurons (N#) recorded in the same experiments are indicated by * and **, respectively. T# indicates the response type.
Brain State Segmentation
During each neuronal recording, a parallel ECoG signal was recorded at a sample rate of 1 kHz from the surface electrode placed on the surface of the cortex nearby the recording site (Figure 1B). To segment the recording into epochs occurring during synchronized vs. desynchronized ECoG states, as previously shown (Enander et al., 2019), the spectral density of the ECoG was calculated with a segment length of 1,000 ms, an overlap of 125 ms, and a constant (mean) detrending. The spectral density of Delta, Theta, and Alpha bands (0–12 Hz) was summed for each segment and the compound value was used for the remainder of the analysis. A desynchronized segment of ECoG was assumed to occur when the compound spectral density dropped below the compound spectral density median for at least two segments in sequence. Therefore, at the onset of every stimulus presentation, and thereby for each response type, there was an identified ECoG state. The fraction of stimulus presentations that occurred within the desynchronized state relative to the total number of presentations was used for statistical comparisons (Figure 5). As the probability of a brain state is expected to influence the observed ratios of occurrences of events within each brain state the statistical method used to make the comparisons was the paired t-test.
Figure 5. No preferred ECoG state for specific response types. (A) Two sample responses belonging to the same response type evoked during synchronized and desynchronized ECoG states. Top two traces are time-compressed ECoG recordings. The pink time zone is expanded below to illustrate the corresponding intracellular responses of a recorded neuron, superimposed. Note the rising Upstate starting before the onset of the stimulation (gray arrow) for the gray trace (synchronized ECoG). (B) Two sample responses from the same neuron and stimulation pattern as in (A), but from response Type 1. Both traces occurred during synchronized ECoG states. The red trace had a rising Upstate that started (red arrow) at almost the exact same time as in the gray trace in (A). (C,D) The fraction of responses occurring in the desynchronized ECoG state for the responses evoked by each stimulation pattern is shown in red, for each neuron, and for each member of a response type (blue). The fraction for each response type is hence paired with the corresponding fraction for the stimulation pattern. The bar chart in (B) illustrates the results of this investigation only for response types with more than nine members (nine responses). In (D), the same display but for all response types identified (including those with less than 10 members, N = 494).
Post-processing—Responses Evoked by Isolated Single-Pulse Stimulation Pulses
As described above, isolated single-pulse stimulations were delivered for each stimulation channel separately in sequences of five repetitions separated from each other by 0.3 s, whereas different such chunks of five pulse stimulations (N = 20 per channel) were randomly intermixed with the full stimulation patterns. We analyzed the responses to such isolated single-pulse stimulations with respect to the onset latency time, the time-to-peak, and the response amplitude, using a tailor-made point-and-click user interface (Figure 6, Table 2).
Figure 6. Random variability of the isolated single-pulse responses. (A) Six superimposed raw intracellular traces evoked by isolated single-pulse stimulations to channel #3, recording from Neuron #7. Pairs of arrows indicate the response onset latency time, the time-to-peak and the peak amplitude for three sample responses. (B) 3D plot of the measured parameters for all responses evoked by isolated single-pulse stimulation of channel #3 in the same neuron. (C) Similar display for all four channels used. (D) The matrix presents the outcome of a statistical investigation of whether the rank order of the single-pulse stimulations within the five pulse chunks (each pulse separated by an interval of 0.3 s) had a statistically significant impact on the response amplitude, by comparing the 20 responses of any one group (the first pulse, for example) to those of each other group (the third pulse, for example) in turn (Wilcoxon rank sum test for pairwise comparisons, p-values are indicated according to the color scale). (E) For the investigation of which an example is shown in (D), this histogram summarizes the result (the distribution of p values) across all stimulation channels and all neurons (binned).
Table 1. Number of response types and the performance measures for the response type separation for each stimulation pattern separately, based on principal component analysis (PCA) and k-nearest neighbors (kNN) analysis.
Table 2. Quantification of single-pulse responses in the population of neurons.
Post-processing—Responses Evoked by Individual Stimulation Pulses Within Patterns
We also performed an analysis of the responses to the individual pulses within the different stimulation patterns (Figure 7). We constrained the analysis to the following set of individual stimulation pulses: Each of the eight spatiotemporal stimulation patterns consisted of 5–33 individual stimulation pulses and a total sum of 152 pulses altogether in the eight patterns used (Figure 1A). The time between subsequent stimulation pulses within the stimulation patterns varied in the span of 1–123 ms. This means that in some cases, the intervals between subsequent stimulation pulses were too short to identify which of the pulses generated the recorded response. Since the scope of this part of the analysis was to investigate the response to specific in-pattern stimulation pulses, only stimulation pulses that were temporally segregated from previous and subsequent pulses by at least 10 ms were included (as the average response latency time was 11 ms, Table 2). Based on this selection criterion, 52 out of the total 152 pulses were included in the analysis of the responses evoked by the individual stimulation pulses within stimulation patterns.
Figure 7. Response metrics differences between neurons for isolated within-pattern stimulation pulses. (A,B) Examples of the variation of the responses to single-pulse stimulation of the same channel depending on the pulse position within the stimulation pattern and the neuron (shown for Neuron#5 and Neuron#11). Five superimposed raw data traces are shown for each neuron for the responses to stimulation of ch#3 at various time points (as indicated by pu# and arrows below traces) within stimulation pattern S20. (C) Average EPSP amplitude, indicated as multiples of the average EPSP amplitude to the isolated single-pulse stimulation of the corresponding channel in the same neuron (Figure 6), shown for each neuron and each analyzed within-pattern stimulation pulse (N = 52), corresponding to a total of 56,762 responses. For each analyzed stimulation pulse, the stimulus pattern and the sequential position of the pulse within that pattern (pu#) for each respective stimulation channel (ch#) are indicated at the bottom (organized according to their order of occurrence within the stimulation patterns). White entries indicate that the response fraction did not surpass the spontaneous level of EPSP events by more than two standard deviations according to the automatic EPSP detection method and the response was hence discarded from further analysis. (D) Similar display as in (A) , but for the relative response latency time, i.e., the latency expressed as multiples of the average response latency time to isolated single-pulse stimulation of the corresponding channel. (E,F) Pairwise comparisons of differences in EPSP amplitudes between all neurons. For two sample stimulation pulse positions (S10, ch#3, pu#6; and F5, ch#3, pu#4, respectively) we analyzed whether the changes in relative EPSP amplitude (compared to the neuron’s isolated single-pulse response) differed between neurons. The p-values of the Wilcoxon rank sum test for each pairwise comparison of the raw EPSP-like responses are reported as a color code in the matrix (N = 56 comparisons).
Responses to the individual stimulation pulses were analyzed both manually and automatically. The manual part of the analysis consisted of a visually guided definition of the onset latency, amplitude height, and latency-to-peak using the tailor-made point-and-click user interface (as in Figures 6A–C). The automatic part consisted of detection of EPSP-like events using tailored template matching routines—its sole purpose was to identify if EPSPs evoked by a particular stimulation pulse were so infrequent that they were at risk of not surpassing the spontaneous occurrence of similar EPSP events (in the recording times in between the presentation of the stimulation patterns), in which case they were to be excluded from further analysis. EPSP templates consisted of a series of 5–20 time-voltage thresholds with individually variable voltage variance and were defined manually for each neuron based on a large sample of EPSP-like events (≫100) occurring after in-pattern stimulation pulses. They were visually confirmed to not omit mid to large EPSP-like events (>3 mV) that occurred spontaneously at randomly sampled time points throughout the duration of the recordings. The response fraction of a neuron to each repetition of a stimulation pulse was defined as the number of repetitions evoking an EPSP-like event, as judged by the automated EPSP identification in the time range 4–18 ms after the pulse onset, divided by the total number of repetitions of that pulse. The baseline activity of that same EPSP-like event was determined by counting its spontaneous activity in time bins of 14 ms width (i.e., the same width as the response time window) for 12 consecutive bins preceding the onset of the stimulation pattern. As each stimulation pattern, as a rule, was repeated 100 times, we typically obtained a total of 1,200 such bins. The response fraction of the spontaneous activity was obtained by taking the average activity across these 1,200 bins. The threshold activity for the EPSP template, i.e., the response fraction that an in-pattern stimulation pulse needed to exceed in order to qualify as an evoked rather than spontaneous response, was defined as the mean plus two standard deviations of the response fraction of the spontaneous activity. If the response fraction was below the threshold activity, or if there were less than five manually detected EPSPs, the response of that in-pattern stimulation pulse was considered not significant and was discarded from further analysis.
Statistical Analysis Summary
For statistical evaluation of the identified response types, we used the Kruskal-Wallis test (Figures 2D,E). For pairwise tests of response fractions occurring under desynchronized brain states (Figure 5), we used paired t-test. For pairwise comparisons of EPSP-like responses (Figures 6, 7), we used the Wilcoxon rank-sum test for pairwise comparisons as the data was not following a normal distribution.
Results
Large Variations in Neocortical Internal States and Responses
While providing specific spatiotemporal tactile afferent activation patterns to the skin (Figure 1A), we made intracellular, whole-cell patch-clamp in vivo recordings from single neocortical pyramidal cells [putative layer III-V pyramidal neurons in the somatosensory cortex (S1), Figure 1B], three of which were morphologically verified to be layer III pyramidal cells (Figure 1C). The time-varying states of the subnetworks connected to the recorded neuron generated a rich background of spontaneous Up and Down states, mixed with episodes of intermediate states, against which the responses evoked by the sensor input patterns could often stand out as distinctly different (Figure 1D). Responses evoked by the same stimulation pattern appeared to be impacted by the preceding state, as reflected in the spontaneous activity (Figure 1D, top two traces). In addition, in some cases, the spontaneous activity could even resemble the responses evoked by the stimulation (Figure 1D, indicated by an arrow in the bottom trace), in agreement with previous observations (Berkes et al., 2011).
Distinct Types of Responses on Repeated Application of Identical Tactile Input Patterns
Given a consistent input pattern (Figure 1A), the response to any element in the stimulation pattern might still be modulated by cortical state changes, internally generated and/or impacted by prior elements of the pattern. Indeed, repeated delivery of one specific stimulation pattern generated multiple different responses that on visual inspection tended to sort into a few different categories, or types (Figure 2A). In order to explore this separability, we applied a clustering method that sorted responses on basis of their membrane potential at each sample time point. In the illustrated case, this method indicated that the majority of the responses evoked by this specific stimulation pattern were separable into six recurring response types (Figure 2B). For the same neuron, the other seven stimulation patterns had three-six response types each. Across all neurons recorded, the responses of the intracellular membrane potential evoked by each stimulation pattern were on average divisible into 3.8 + /−1.3 response types (Table 1).
For the clustering method we used, the set thresholds of allowed variability of the membrane potential were naturally critical factors for the outcome. A parameter sensitivity analysis showed how different parameter settings affected the outcome in two different sample neurons (Figure 2C). The aim for the parameter settings was to have as few ungrouped responses and as many groups as possible. However, across the population of neurons, there was no single setting that would maximize these two factors across the whole population of neurons (Figure 2C). As shown in Figure 2C, the parameter settings that were chosen were located in the center region of a relatively smooth landscape of possible outcomes. This was reassuring, as it showed that our clustering algorithm was not unstable in the sense that gradual changes in the two parameters Threshold and Overshoot did not result in dramatically different outcomes (see also Supplementary Figure 1).
We next investigated if the clustered response types were different from each other using the Kruskal-Wallis test. We compared the responses that belonged to each response type to the responses evoked by the same stimulation pattern in the same neuron but belonging to the other response types. Whereas the Kruskal-Wallis test results for individual response types could fluctuate across different sample time points (Figure 2D), it was clear that the probability of H0 was substantially lowered during the time window that the stimulation pattern was on 0–350 (ms, approximately). When we instead plotted the fraction of the response types that were statistically different (at p < 0.05) from the other responses evoked by the same respective stimulation pattern, this separability was true for more than 60% of the response types across every time point for the duration of the stimulation patterns (Figure 2E). When the stimulation patterns ceased, however, the separability dropped relatively rapidly (within 100–150 ms) to chance levels (Figure 2E). Furthermore, when the response type labels were shuffled among the responses evoked by the same stimulation pattern in the same neuron, the separability collapsed to chance level (0.05; red trace in Figure 2E) for all time points.
As another independent verification method, we next used Principal Component Analysis (PCA) to find out to what extent the shapes of the full duration traces of the response types were distinctly different from each other. This method provided a measure of the accuracy, or the distinctness of separation of the individual responses, as summarized in the confusion matrix for the sample stimulation pattern in this neuron (Figure 3A). The mean accuracy across the different response types, indicated by the diagonal, was on average 56.7% for this stimulation pattern (including ungrouped responses; the chance level was = 1/7 = 14.3%). Across the population of recorded cells, the accuracy of the separation of the responses into the identified response types was generally above 60% for each of the eight stimulation patterns (Table 1), verifying that the identified response types were composed of a set of responses that to a relatively high degree were orthogonal to the responses belonging to other response types. A sensitivity analysis of the parameter choice for the clustering method showed that the chosen parameter values did not provide the highest possible decoding accuracy (Figure 3B), but was instead a balance between this metric, the number of ungrouped responses, and the number of response types identified (Figure 2C).
Whereas the majority of the evoked responses could be classified as belonging to one of the response types (Table 1), “ungrouped” responses (Figure 3A) were by definition a much broader and heterogeneous class and consequently had a much greater risk of confusion (i.e., a low value in the diagonal and a higher prevalence of well-above-zero values outside the diagonal, Figure 3A). To evaluate the “ungrouped” responses we used the F1 score. A high F1 score value indicates a small risk of confusion with any of the specific classes of responses. In Figure 3A, the F1 score for ungrouped responses was 0.32. This value does suggest some degree of confusion, even though the confusion matrix indicates that there wasn’t any specific response type that the ungrouped responses were confused with. Across the dataset, the F1 score was 0.48 (Table 1), indicating that the ungrouped responses were overall well separated from the responses classified as belonging to the defined response types, and therefore potentially represented additional response types that our method could not reliably identify due to the limited number of repetitions (100 or less) available for each input pattern.
We also used PCA to compare the response types evoked by all eight stimulation patterns (Figure 3C). As each of the eight stimulation patterns on average evoked about four responses types (Table 1), a large number of comparisons were made in this investigation (and hence the analysis results for pattern S5 in Figure 3A, for example, cannot be expected to be preserved, as there were so many more responses to be compared with, which substantially increases the probability of confusion/misclassification for each response). In this case, the accuracy of the response types (ungrouped responses were excluded here) was on average 57.6%. This is a remarkably high specificity, considering that the chance level for this comparison was only 2.9% (i.e., 1/34, for the comparison involving 34 response types). Across all cells recorded, this cross-stimulation pattern response type accuracy was 38.2 +/−21.0%. Hence, in each individual neuron, a majority of the evoked responses for a single stimulation pattern could be divided into distinct types that were separable with high accuracy from responses not belonging to that type, regardless of whether they were evoked by the same or by other stimulation patterns. In order to confirm that the reported decoding accuracy was not due to chance, we shuffled the stimulation pattern labels between the responses recorded in the same neurons. In this case, the decoding accuracy dropped substantially (Figure 3D). The fact that it did not drop all the way to chance level was likely due to the circumstance that many of the evoked responses, just like the spontaneous activity (Figure 1D), did share some common features.
Uniqueness of the Response Types Between Neurons
We have previously found that the intracellular responses evoked by a given stimulation pattern are different between neurons (Oddo et al., 2017), which on average was the case also in the present set of recordings (Figure 4A). Here, the issue was instead if the various response types evoked by the same stimulation pattern were distinct from each other across the recorded neurons, which would indicate differences in the subnetworks connected to them. Figure 4B illustrates a confusion matrix of a sample stimulation pattern, for which the accuracy of response types separation was 42.7% (chance = 1/48 = 2.1%) across the 13 neurons (some neurons were recorded in the same experiment, as indicated). Across all eight stimulation patterns, this accuracy was 39.5 +/−3.5%. Moreover, in the sample illustration (Figure 4B), we found that 45 out of 48, or 93.8%, of the response types were separable, from the responses of all other types (i.e., decoding accuracy higher than chance). Across all stimulation patterns, 44.1 +/−2.9 response types, or 89.6 +/−3.0% of each response type identified, were similarly separable from all other response types evoked by the same stimulation pattern in different cells.
Response Types Were Not Associated With Specific ECoG States
Given that anesthesia tends to increase the duration and depth of synchronized cortical activity states (Constantinople and Bruno, 2011), during which responses could possibly become more stereotyped, one suspicion could be that some response types occurred only during synchronized cortical activity whereas other response types occurred only in the desynchronized state. However, comparing the prevailing ECoG state across individual responses of the same response type demonstrated that this was not a rule (Figure 5A). A response evoked in the middle of a desynchronized ECoG state could be classified to the same response type (Type 0) as another response evoked in the middle of a synchronized ECoG state, even in cases where the response was evoked in the middle of a starting Upstate in the neuron (Figure 5A; note that Up-Down states tend to increase with synchronized states but are not equivalent to such states; Poulet and Petersen, 2008). To further show that a rising Upstate alone was not a strong determinant of the response type, a sample trace of another response type (Type 1) is illustrated in Figure 5B, alongside another trace also belonging to Type 1. These two examples illustrate that the presence or absence of an Upstate initiation was not in itself a decisive predictor of which response type was going to be evoked.
We next systematically compared the relative probability of the desynchronized ECoG state for each stimulation pattern (a total of 8 times 13 neurons = 104 data points) with the corresponding probability for each of the 494 response types identified (Figures 5C,D). The issue addressed was if some of the response types were more likely to be associated with a desynchronized or a synchronized ECoG state. In Figure 5C, we only analyzed the distribution for response types with more than nine individual responses/members. In Figure 5D, we also included response types with fewer than 10 members, which resulted in a separate small peak at 0.0 and a few very high values for the type-separated responses. These outlier bars could be explained by chance omission or addition of one response in the desynchronized state, which could greatly impact the ratio for response types with fewer than 10 members as the overall probability of being in the desynchronized ECoG state was around 20%. Nevertheless, under both conditions, the response types were not distributed significantly differently from the overall responses evoked by the same stimulation pattern according to the paired t-test (p = 0.76 and p = 0.68, respectively). Hence, the ECoG indicator of the global brain state was not predictive of the response type, suggesting that the response types were not an effect merely of this global brain state but rather the result of an interplay between the input and the dynamics of subnetwork-specific internal state evolution as the stimulus presentation unfolded.
Neuron-Specific Responses to the Individual Pulses That Composed the Stimulation Patterns
We also performed a separate analysis of the responses evoked by the individual stimulation pulses, of which each stimulation pattern was composed. A striking feature of the responses to isolated single-pulse stimulations, delivered in between the stimulation patterns (analyzed in Figures 1–5), was their highly variable nature (Figures 6A–C, Table 2), in contrast to the known consistency of the primary sensory afferent activation of this type of electrocutaneous stimulation (Bengtsson et al., 2013). A statistical analysis moreover suggested that there was no order-dependence of the response amplitudes evoked by the isolated single-pulse stimulations when delivered in chunks of five pulses separated by 0.3 ms (Figures 6D,E). In contrast, for the same single-pulse stimulations, when they were part of one of the stimulation patterns, the responses displayed a relative specificity that depended on the position of the pulse within the stimulation pattern (Figure 7), and which moreover differed between neurons (Figures 7A,B). Systematic analysis indicated that each neuron constituted a unique pattern within the “mosaic” of response parameter variations across the whole set of “within-pattern” stimulation pulses (Figures 7C,D). This was quantified using a pairwise statistical comparison across all neurons, illustrated in Figures 7E,F for two sample within-pattern stimulation pulses. Across all 52 stimulation pulse positions tested (see “Materials and Methods” section), statistically significant differences between the neurons (at p < 0.05, Wilcoxon rank-sum test for pairwise comparisons; the analysis only included responses that were statistically significant, i.e., non-white entries in Figures 7C,D) were found in 65% of the comparisons of the amplitudes, 78% of the comparisons of the time-to-peaks and 86% of the comparisons of the response latency times. Thus, for the majority of the within-pattern stimulation pulses, the changes in response amplitude, time-to-peak, and response latency relative to the isolated single-pulse stimulation were statistically different between the neurons, corroborating the results from Figure 4B that each neuron responded with partly unique time-voltage curves across the duration of the stimulation patterns. Altogether, the results in Figure 7 show that the outcome of the response type identification analysis cannot be explained by a generic sequence of paired-pulse depression or facilitation phenomena, since these would affect all neurons in the same way.
Discussion
Our results show that the responses evoked by tactile sensory input patterns fall into a limited subset of preferred response states that are specific to each input pattern and each cortical neuron. The relative orthogonality of these recurring response types (Figure 3A, Table 1), the identification of which was made possible by the high resolution of the recording and precise repeatability of the stimulation method, indicates that they resulted from specific response states of the cortical networks, rather than being noise or specific patterns of spontaneous activity unrelated to the stimuli. The finding that the response types were relatively unique even when all the response types evoked across all the eight stimulation patterns were compared (Figure 3C) is another strong piece of evidence against that possibility. Moreover, the observation that the response types are relatively unique for each neuron (Figure 4B) indicates the existence of multiple, parallel, neuron-specific subnetworks, which could reflect specific aspects of the current global internal brain state (Figure 8). This is conceivable since the global internal states subsampled by each neuron, in turn, are likely to encompass large parts of the neocortex (Spanne and Jorntell, 2015; Enander et al., 2019; Stringer et al., 2019a,b; Wahlbom et al., 2019), as underscored by the estimate that any neuron in the neocortex connects to any other neocortical neuron with synaptic linkages within no more than five neurons on average (Arbib et al., 1998). The existence and regulation of internal brain states can be expected to be essential for forming the percepts (or illusions), generated by given sensor activation patterns (Geldard and Sherrick, 1972; Robles-De-La-Torre and Hayward, 2001). Our results indicate that the neocortical network physiology has a non-uniform, non-continuous landscape of solutions, which may cause the response types observed to arise as the result of an adaptable time-evolving match between various subnetworks’ internal states and the actual sensory input pattern received.
Figure 8. Summary of findings. Recorded neocortical neurons receive a high number of synaptic inputs, each of which could potentially be controlled by a different subnetwork. Each subnetwork may be gated by other, relatively independent, events in the global neocortical circuitry. When the same tactile afferent stimulation pattern is provided at different points in time, these subnetworks may have different initial states, which will impact which subnetworks are gated in or out, which in turn will affect the appearance of the intracellular response (two traces shown in the top right corner). Circles are neurons, triangles are synapses, lines are axons and variable connections are drawn by dashed lines.
Methodological Considerations
Our response types should not be regarded as fixed entities, for example the specific types detected could grow or shrink with altered parameter settings of the clustering method (see Figure 2C) and there was a degree of response variability within each identified response type (Figure 2B). However, as each response type could subsequently be demonstrated to have a degree of orthogonality to the responses belonging to other types (Figure 3A, Table 1), our data show that the dynamics of the neocortical network created a multitude of alternative, non-continuous response states to any given sensory input pattern. It is highly unlikely that our identified average number of four response types per stimulation pattern represents anything near an actual upper limit—if we had had the possibility to repeat each stimulation pattern millions of times, rather than one hundred, a higher number of response types would most likely have been possible to identify. Nevertheless, as shown in our sensitivity analysis, where we tested hundreds of different parameter settings for our clustering algorithm, almost every clustering setting with more than one identified response type indicated that the identified clusters were objectively separable using PCA + kNN (Figure 3B). Hence, our analysis indicates the existence of discontinuity, or local minima, in the response states of the cortical network.
Alternative approaches to clustering, like machine learning networks or k-means, were not explored here, as these approaches commonly require a target aim of how many clusters should be sought for, or can suffer from instability and high sensitivity to the exact choice of hyperparameters. Our clustering method was instead designed to identify how many clusters could be observed based on clearly defined metrics of the neurophysiological recording data.
How many possible response types exist and what specific membrane voltage curves they would create are more than likely to depend on the condition. Hence, it would not be surprising if they are instantiated differently in the awake behaving animal, where they may be varying even further with the “state of mind”, and differently still when there is an accurate prediction of an expected sensory input being made. However, apart from being methodologically nearly impossible to explore in the awake animal (due to the requirement of achieving identical stimuli over a long time), a correct analysis of such data in this respect would have required a precise estimate of what the experimental subject is thinking at the time of the stimulus delivery, which is today not theoretically possible.
To what extent could our results be a “product” of the anesthetics? Anesthesia increases the probability of the thalamocortical system to enter episodic, coordinated oscillatory modes (Amzica and Steriade, 1998; Constantinople and Bruno, 2011), or synchronized states, even though they are certainly not uncommon in the awake state (Poulet and Petersen, 2008; Bennett et al., 2013; Petersen and Crochet, 2013; synchronized states are widely observed also in awake humans; Sachdev et al., 2015). However, in between such episodes, there are periods of desynchronized activity, or states, which do not differ between the awake and the anesthetized animal, though anesthesia could be expected to lead to an overall activity reduction within desynchronized states (Constantinople and Bruno, 2011). Importantly, the response types were not specifically associated with any of these two opposite global states (Figure 5), which is a strong argument against the possibility that the response types are a peculiarity created by the increased time spent in a synchronized state that is a consequence of the anesthesia. Hence, whereas our recordings do not reflect natural integrated thought processes, they can still provide important information about underlying governing principles of information processing, principles which would be much harder to disentangle in awake recordings. However, there are data from awake behaving animals that are compatible with the presence of multi-structure cortical states, for example, that the activity of a population of neocortical neurons can follow different trajectories of interdependencies depending on the motor task performed (Churchland et al., 2012; Gallego et al., 2017; Russo et al., 2020). Distributed multi-stable states are also compatible with the dynamic routing of brain activity in humans (Finger et al., 2019).
Relation to Previous Literature
Previous analyses of state-conditioned intracellular signal variations evoked by somatosensory inputs have largely focused on the binary issue of Up and Down states using single-shot or brief inputs. The seemingly paradoxical finding that inputs provided in an excited, Up, state result in much lower responses than the same input in a Downstate (Petersen et al., 2003) was in a detailed conductance-level analysis indicated to be partly due to the shunting effect caused by the Upstate being associated with a high level of excitatory synaptic input but also in parallel a high level of inhibitory synaptic input (Hasenstaub et al., 2007). Moreover, it was observed that stimuli during a Downstate could often evoke disproportionately large responses, by inducing transitions to an Upstate. More recent studies suggest such effects to be inducible because the neocortex is a system with criticality phenomena, where small changes at critical states can cause avalanche effects of increased excitability that spread widely through the cortical network (Wright and Wessel, 2017; Johnson et al., 2019). Although our data indicate a much more fine-grained subdivision of cortical states than the binary dichotomy between Up and Down states, such avalanche effects could be the underlying reason for the response types we observed here. If each subnetwork for example exhibits partly independent criticality, where each individual synaptic input to a neuron could in theory represent the activity state of partly independent subnetworks, this would be compatible with the fact that different neurons exhibited different response types to the same inputs. Moreover, the data suggest that different individual neurons are connected to unique combinations of such subnetworks (Figure 8).
Implications for Our View on the Neocortical Mode of Operation
Although our data naturally do not allow identification of the full perceptual processes, they explore the physiological brain mechanisms that support such processes, which likely correspond to mechanistic underpinnings of expectations or predictions (Loeb and Fishel, 2014). The findings suggest that the multidimensional latent state defined by the large populations of neurons across the neocortex (Stringer et al., 2019b) works according to attractor-like dynamics (Ringach, 2009) with multiple attractors for each given input. The interaction between the input pattern and the cortical state at the moment of stimulus onset would thus cause the cortical network to fall into one specific out of many possible attractors, or series of such attractors. Since different neurons would be coupled to different subnetworks of this global network, the same attractor would potentially impact different neurons in different ways, as our data indicate (Figure 4B). The specificity of the response types observed in individual neurons would thus be local subnetwork-instantiations of the time-evolving input-updated brain-wide state estimations of world and body, hence fundamental mechanisms for forming rich perception, and illusions. We expect these principles to reflect a general computational strategy used by the neocortex across all sensory systems.
Data Availability Statement
The original contributions presented in the study are publicly available. This data can be found here: https://doi.org/10.6084/m9.figshare.14681340.v1. Further requests can be directed to the corresponding author.
Ethics Statement
The animal study was reviewed and approved by Local Ethics Committee of Lund, Sweden.
Author Contributions
JN and HJ planned the study and conducted the experiments. JN, JMDE, and HJ designed the analysis. JN, JMDE, and HM conducted the analysis. JN and HJ wrote the article. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the EU Grant FET 829186 ph-coding (Predictive Haptic COding Devices In Next Generation interfaces) and the Swedish Research Council (project grant no. K2014-63X-14780-12-3).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Jerry Loeb (USC Los Angeles) and Matthias Kohler (TU Munich) for pre-reviewing our manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2021.677568/full#supplementary-material.
SUPPLEMENTARY FIGURE 1 | Stability of identified clusters across the full range of parameter settings. (A) Mean agreement between surrounding elements. The position of each center element was defined by the values of its “Threshold” and “Overshoot” parameters. The color code indicates the average number/proportion of response members of the center element that fell into the same clusters as in each of the eight surrounding elements, across all clusters identified. Data for Neuron#5, pattern S5 (same as in Figures 2, 3). (B) Maximal agreement between surrounding elements. Same analysis as in (B), but instead the highest level of agreement between two clusters is illustrated. The color code here instead indicates the average max agreement between the central element and its eight surrounding elements.
References
Amzica, F., and Steriade, M. (1998). Electrophysiological correlates of sleep delta waves. Electroencephalogr. Clin. Neurophysiol. 107, 69–83. doi: 10.1016/s0013-4694(98)00051-0
Arbib, M. A., Numbers, R. L., Érdi, P., Szentágothai, J., Szentagothai, J., and Szentagothai, A. (1998). Neural Organization: Structure, Function and Dynamics. Cambridge, MA: MIT press.
Arieli, A., Sterkin, A., Grinvald, A., and Aertsen, A. (1996). Dynamics of ongoing activity: explanation of the large variability in evoked cortical responses. Science 273, 1868–1871. doi: 10.1126/science.273.5283.1868
Bandeira, F., Lent, R., and Herculano-Houzel, S. (2009). Changing numbers of neuronal and non-neuronal cells underlie postnatal brain growth in the rat. Proc. Natl. Acad. Sci. U S A 106, 14108–14113. doi: 10.1073/pnas.0804650106
Bengtsson, F., Brasselet, R., Johansson, R. S., Arleo, A., and Jorntell, H. (2013). Integration of sensory quanta in cuneate nucleus neurons in vivo. PLoS One 8:e56630. doi: 10.1371/journal.pone.0056630
Bennett, C., Arroyo, S., and Hestrin, S. (2013). Subthreshold mechanisms underlying state-dependent modulation of visual responses. Neuron 80, 350–357. doi: 10.1016/j.neuron.2013.08.007
Berkes, P., Orban, G., Lengyel, M., and Fiser, J. (2011). Spontaneous cortical activity reveals hallmarks of an optimal internal model of the environment. Science 331, 83–87. doi: 10.1126/science.1195870
Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu, S. I., et al. (2012). Neural population dynamics during reaching. Nature 487, 51–56. doi: 10.1038/nature11129
Constantinople, C. M., and Bruno, R. M. (2011). Effects and mechanisms of wakefulness on local cortical networks. Neuron 69, 1061–1068. doi: 10.1016/j.neuron.2011.02.040
Curto, C., Sakata, S., Marguet, S., Itskov, V., and Harris, K. D. (2009). A simple model of cortical dynamics explains variability and state dependence of sensory responses in urethane-anesthetized auditory cortex. J. Neurosci. 29, 10600–10612. doi: 10.1523/JNEUROSCI.2053-09.2009
Ekerot, C. F., Gustavsson, P., Oscarsson, O., and Schouenborg, J. (1987). Climbing fibres projecting to cat cerebellar anterior lobe activated by cutaneous A and C fibres. J. Physiol. 386, 529–538. doi: 10.1113/jphysiol.1987.sp016549
Enander, J.M.D., Spanne, A., Mazzoni, A., Bengtsson, F., Oddo, C. M., and Jorntell, H. (2019). Ubiquitous neocortical decoding of tactile input patterns. Front. Cell. Neurosci. 13:140. doi: 10.3389/fncel.2019.00140
Finger, H., Gast, R., Gerloff, C., Engel, A. K., and Konig, P. (2019). Probing neural networks for dynamic switches of communication pathways. PLoS Comput. Biol. 15:e1007551. doi: 10.1371/journal.pcbi.1007551
Fiser, J., Chiu, C., and Weliky, M. (2004). Small modulation of ongoing cortical dynamics by sensory input during natural vision. Nature 431, 573–578. doi: 10.1038/nature02907
Gallego, J. A., Perich, M. G., Miller, L. E., and Solla, S. A. (2017). Neural manifolds for the control of movement. Neuron 94, 978–984. doi: 10.1016/j.neuron.2017.05.025
Geldard, F. A., and Sherrick, C. E. (1972). The cutaneous “rabbit”: a perceptual illusion. Science 178, 178–179. doi: 10.1126/science.178.4057.178
Genna, C., Oddo, C. M., Mazzoni, A., Wahlbom, A., Micera, S., and Jorntell, H. (2018). Bilateral tactile input patterns decoded at comparable levels but different time scales in neocortical neurons. J. Neurosci. 38, 3669–3679. doi: 10.1523/JNEUROSCI.2891-17.2018
Golub, M. D., Sadtler, P. T., Oby, E. R., Quick, K. M., Ryu, S. I., Tyler-Kabara, E. C., et al. (2018). Learning by neural reassociation. Nat. Neurosci. 21, 607–616. doi: 10.1038/s41593-018-0095-3
Hasenstaub, A., Sachdev, R. N., and McCormick, D. A. (2007). State changes rapidly modulate cortical neuronal responsiveness. J. Neurosci. 27, 9607–9622. doi: 10.1523/JNEUROSCI.2184-07.2007
Hayward, V., Terekhov, A. V., Wong, S. C., Geborek, P., Bengtsson, F., and Jorntell, H. (2014). Spatio-temporal skin strain distributions evoke low variability spike responses in cuneate neurons. J. R. Soc. Interface 11:20131015. doi: 10.1098/rsif.2013.1015
Jenmalm, P., Birznieks, I., Goodwin, A. W., and Johansson, R. S. (2003). Influence of object shape on responses of human tactile afferents under conditions characteristic of manipulation. Eur. J. Neurosci. 18, 164–176. doi: 10.1046/j.1460-9568.2003.02721.x
Johnson, J. K., Wright, N. C., Xia, J., and Wessel, R. (2019). Single-cell membrane potential fluctuations evince network scale-freeness and quasicriticality. J. Neurosci. 39, 4738–4759. doi: 10.1523/JNEUROSCI.3163-18.2019
Loeb, G. E., and Fishel, J. A. (2014). Bayesian action&perception: representing the world in the brain. Front. Neurosci. 8:341. doi: 10.3389/fnins.2014.00341
Luczak, A., and Bartho, P. (2012). Consistent sequential activity across diverse forms of UP states under ketamine anesthesia. Eur. J. Neurosci. 36, 2830–2838. doi: 10.1111/j.1460-9568.2012.08201.x
Luczak, A., Bartho, P., and Harris, K. D. (2009). Spontaneous events outline the realm of possible sensory responses in neocortical populations. Neuron 62, 413–425. doi: 10.1016/j.neuron.2009.03.014
Mogensen, H., Norrlid, J., Enander, J.M.D., Wahlbom, A., and Jorntell, H. (2019). Absence of repetitive correlation patterns between pairs of adjacent neocortical neurons in vivo. Front. Neural Circuits 13:48. doi: 10.3389/fncir.2019.00048
Narayanan, R. T., Udvary, D., and Oberlaender, M. (2017). Cell type-specific structural organization of the six layers in rat barrel cortex. Front. Neuroanat. 11:91. doi: 10.3389/fnana.2017.00091
Niedermeyer, E., and da Silva, F. H. L. (2005). Electroencephalography: Basic Principles, Clinical Applications and Related Fields. Philadelphia, PA: Lippincott Williams & Wilkins.
Oddo, C. M., Mazzoni, A., Spanne, A., Enander, J. M., Mogensen, H., Bengtsson, F., et al. (2017). Artificial spatiotemporal touch inputs reveal complementary decoding in neocortical neurons. Sci. Rep. 8:45898. doi: 10.1038/srep45898
Oddo, C. M., Raspopovic, S., Artoni, F., Mazzoni, A., Spigler, G., Petrini, F., et al. (2016). Intraneural stimulation elicits discrimination of textural features by artificial fingertip in intact and amputee humans. eLife 5:e09148. doi: 10.7554/eLife.09148
Petersen, C. C., and Crochet, S. (2013). Synaptic computation and sensory processing in neocortical layer 2/3. Neuron 78, 28–48. doi: 10.1016/j.neuron.2013.03.020
Petersen, C. C., Hahn, T. T., Mehta, M., Grinvald, A., and Sakmann, B. (2003). Interaction of sensory responses with spontaneous depolarization in layer 2/3 barrel cortex. Proc. Natl. Acad. Sci. U S A 100, 13638–13643. doi: 10.1073/pnas.2235811100
Poulet, J. F., and Petersen, C. C. (2008). Internal brain state regulates membrane potential synchrony in barrel cortex of behaving mice. Nature 454, 881–885. doi: 10.1038/nature07150
Ringach, D. L. (2009). Spontaneous and driven cortical activity: implications for computation. Curr. Opin. Neurobiol. 19, 439–444. doi: 10.1016/j.conb.2009.07.005
Robles-De-La-Torre, G., and Hayward, V. (2001). Force can overcome object geometry in the perception of shape through active touch. Nature 412, 445–448. doi: 10.1038/35086588
Russo, A. A., Khajeh, R., Bittner, S. R., Perkins, S. M., Cunningham, J. P., Abbott, L. F., et al. (2020). Neural trajectories in the supplementary motor area and motor cortex exhibit distinct geometries, compatible with different classes of computation. Neuron 107, 745.e6–758.e6. doi: 10.1016/j.neuron.2020.05.020
Sachdev, R. N., Gaspard, N., Gerrard, J. L., Hirsch, L. J., Spencer, D. D., and Zaveri, H. P. (2015). Delta rhythm in wakefulness: evidence from intracranial recordings in human beings. J. Neurophysiol. 114, 1248–1254. doi: 10.1152/jn.00249.2015
Spanne, A., and Jorntell, H. (2015). Questioning the role of sparse coding in the brain. Trends Neurosci. 38, 417–427. doi: 10.1016/j.tins.2015.05.005
Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D. (2019a). High-dimensional geometry of population responses in visual cortex. Nature 571, 361–365. doi: 10.1038/s41586-019-1346-5
Stringer, C., Pachitariu, M., Steinmetz, N., Reddy, C. B., Carandini, M., and Harris, K. D. (2019b). Spontaneous behaviors drive multidimensional, brainwide activity. Science 364:255. doi: 10.1126/science.aav7893
Tan, D. W., Schiefer, M. A., Keith, M. W., Anderson, J. R., Tyler, J., and Tyler, D. J. (2014). A neural interface provides long-term stable natural touch perception. Sci. Transl. Med. 6:257ra138. doi: 10.1126/scitranslmed.3008669
Wahlbom, A., Enander, J. M. D., Bengtsson, F., and Jorntell, H. (2019). Focal neocortical lesions impair distant neuronal information processing. J. Physiol. 597, 4357–4371. doi: 10.1113/JP277717
Keywords: tactile, in vivo intracellular, neocortex, cortical state, synaptic input
Citation: Norrlid J, Enander JMD, Mogensen H and Jörntell H (2021) Multi-structure Cortical States Deduced From Intracellular Representations of Fixed Tactile Input Patterns. Front. Cell. Neurosci. 15:677568. doi: 10.3389/fncel.2021.677568
Received: 07 March 2021; Accepted: 18 May 2021;
Published: 14 June 2021.
Edited by:
Tycho Hoogland, Erasmus Medical Center, NetherlandsReviewed by:
William W. Lytton, SUNY Downstate Medical Center, United StatesDaniele Linaro, Politecnico di Milano, Italy
Copyright © 2021 Norrlid, Enander, Mogensen and Jörntell. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Henrik Jörntell, aGVucmlrLmpvcm50ZWxsQG1lZC5sdS5zZQ==
† These authors have contributed equally to this work and share first authorship