Hippocampal output profoundly impacts the interpretation of tactile input patterns in SI cortical neurons

Summary Due to continuous state variations in neocortical circuits, individual somatosensory cortex (SI) neurons in vivo display a variety of intracellular responses to the exact same spatiotemporal tactile input pattern. To manipulate the internal cortical state, we here used brief electrical stimulation of the output region of the hippocampus, which preceded the delivery of specific tactile afferent input patterns to digit 2 of the anesthetized rat. We find that hippocampal output had a diversified, remarkably strong impact on the intracellular response types displayed by each neuron in the primary SI to each given tactile input pattern. Qualitatively, this impact was comparable to that previously described for cortical output, which was surprising given the widely assumed specific roles of the hippocampus, such as in cortical memory formation. The findings show that hippocampal output can profoundly impact the state-dependent interpretation of tactile inputs and hence influence perception, potentially with affective and semantic components.


Tactile input patterns drive SI neurons via multiple dynamically gated pathways
Hippocampal output differentially gates the transmission in these pathways Hippocampal output is thereby part of cortical anticipation, equivalent to a ''prior'' Its function may be to impact semantic and emotional aspects in tactile perception

INTRODUCTION
Analysis of the human brain has indicated that the neocortex is functionally a heavily interconnected network. 1 Various analyses in the rodent tell a similar story. [2][3][4] We have previously shown that the processing of tactile information also at the cellular level in the primary somatosensory cortex (SI) is dependent on other parts of the cortex. This was shown both by stroke-like lesions and by electrical manipulation of brain activity in neocortical areas with a remote location to the SI. 5,6 Recent anatomical studies have indicated that the extent of connections made by individual corticocortical axons have previously been greatly underestimated 7,8 (see also https://www. janelia.org/project-team/mouselight; https://mouse.braindatacenter.cn/), hence giving structural support for the observations that information in the neocortex is globally integrated. Given that perception to some extent depends on expectation, which through the internal cortical state may exert impact already before the tactile afferent information enters the SI cortex, these findings suggest that tactile perception could be influenced by a globally integrated action across the neocortex.
The hippocampal formation is closely interconnected with the neocortex. The hippocampus receives sensory information from a large number of different neocortical areas feeding, via the entorhinal cortex, into the subiculum. 9 Indeed, the hippocampus is also engaged in processing of tactile input. 10 Through its output via the subiculum, the hippocampus can monosynaptically activate the anterior thalamus 11 and indirectly thereby also widespread areas of the cortex. Subiculum neurons can also have widespread projections, such as to the hypothalamus, regions of the temporal/entorhinal cortices surrounding the hippocampus, and the ventral striatum, 11,12 which also enables the hippocampus to potentially impact the activity across widespread areas of the neocortex.
Using intracellular recordings from SI neurons, we have previously shown that repeated presentations of a given spatiotemporal tactile input pattern from digit 2 combine with the state of the neocortical network to result in a variety of response states in SI neurons. 13 This finding indicated that traditional analysis methods involving response averaging would risk missing essential aspects of how the neocortical neurons process information, related to the perceptual process. Interestingly, the relative orthogonality of these response states suggested the presence of attractor-type dynamics, 13,14 reminiscent of what has also been proposed to characterize the neuron population activity dynamics in the hippocampus. 15, 16 More recently, we iScience Article showed that manipulation of the internal cortical state through electrical activation of a remote neocortical area can change the preferred network solutions, or the types of response states observed in SI neurons, for every given tactile input pattern. 5 Remarkably, this change occurred despite that the remote cortical perturbation did not induce any overt response in the recorded neurons. That finding suggested that the effects on the response states in individual neurons were indirectly mediated via more subtle changes in the state of the neocortex globally, which in turn would alter the number of ''open'' network pathways supplying the recorded neuron with information related to the tactile event as a function of time (i.e., cortical state). Such gating of input pathways and information distribution across the neocortex is likely part of the central mechanisms for generating perception, and alteration of such pathways could be a root cause of illusions. Here, we address the question if also the hippocampal output, generated by electrical activation of its neural output stage in the subiculum, can impact the responses in SI neurons to given spatiotemporal tactile input patterns. We find that the hippocampal activation dramatically altered the responses of the SI neurons to tactile inputs, suggesting that the state of hippocampal activity could have a profound effect on the cortical interpretation, or perception, of such inputs. Perhaps the most surprising aspect is that we found no major qualitative difference between the present hippocampal impact versus the previously explored remote cortical impact, 5 despite that the hippocampus is widely considered to have very different specific roles in neocortical processing and memory.

RESULTS
We made whole-cell patch-clamp recordings from putative pyramidal neurons (see STAR Methods; also 13 ) between layers III-V in the forepaw region of the SI cortex, while stimulating the distal part of the second digit with eight fixed spatiotemporal tactile activation (TA) patterns delivered in pseudo-random order ( Figure 1A), a protocol we have used in numerous previous publications. 5,6,13,[17][18][19] Here, these TA input patterns were provided in isolation or in combination with a preceding brief stimulation of the subiculum/terminal CA1 region of the hippocampus (HIPTA) (Figures 1A-1C). As the study relied on comparisons between individual raw responses with large variance, we needed a high number of repetitions of each stimulus, and we could therefore include only recordings which maintained a high quality for a duration of at least 25 min (N = 9 neurons).
As we have reported earlier, 5,13 each individual neuron in SI cortex displayed a wide variety of response types when the exact same TA stimulus was repeatedly presented ( Figure 1D). This indicates that response averaging would risk missing essential features of these cortical responses as averaging would tend to overemphasize the impact of the peripheral component of the response versus its internally generated components. Since the latter presumably is essential to the perceptual process in the cortex, which we were interested in here, we instead focused on the characterization of the pool of individual raw responses evoked by each stimulus pattern and each stimulus condition. This process involved the step to cluster the responses evoked by the same stimulus pattern. We used a clustering method based on the temporal profile of the evoked responses ( Figure 1D; see also 5 ), where the temporal response profile reflects the temporal structure of the responses in the different subnetworks forwarding input to the recorded neuron. 5 Clustering relied on analysis of spike-free intracellular responses, which reflect the synaptic input to the cell. In each recording we therefore applied a mild hyperpolarization current (in the order of À10 pA) to prevent the neuron from spiking (the occasional few remaining spikes could be blanked in software),  iScience Article whereas the quality of the recording was occasionally verified by injecting depolarizing current to evoke spike responses ( Figure 1E).
Stimulation of the subiculum (HIP) in isolation (iso) typically evoked detectable responses in the average ECoG signal (recorded from the parietal cortex, Figure 1A), which were several order of magnitudes longer than the duration of the stimulation (7 pulses at 333 Hz, i.e., 18 ms), whereas the TA stimulation alone evoked no detectable response in the ECoG ( Figure 1F, left column). In the average intracellular signal of the recorded SI neurons, this HIP response was often much less distinct and typically less long lasting than in the ECoG, although it did induce a clear tendency toward oscillation in the intracellular signal (Figure 1F, right column). In contrast, the average response to the TA stimulus was always clear in the intracellular recordings from the SI neurons (though widely different from one neuron to the other, see 13 ). Notably, the intracellular responses evoked by HIP iso and by the TA pattern were clearly different from the HIPTA response, indicating that the effect on the global cortical state induced by the HIP stimulation was nonlinearly combined with the impact of the TA input to produce a new response. Figure S1 presents a frequency distribution analysis of the intracellular signal across the neurons, which indicated that the HIP stimulation pushed the intracellular signal in the low frequency range to a somewhat higher frequency, i.e., in line with the observation of the oscillatory response components induced by the HIP stimulation in Figure 1F.
Interestingly, the combination of a TA input with the preceding HIP stimulation could generate widely different outcomes across different TA input patterns. Figure 2 illustrates examples of the impacts the Notably, the hippocampal output caused the S10 TA input pattern to completely lose an early excitatory component, whereas a later excitatory component was instead greatly amplified. For the S5 TA input pattern, the early excitatory component was in contrast not eliminated, in fact a bit amplified but also slightly slowed down. Also the impacts on the two other example TA input patterns F5 and F10 were highly different. These results indicate that the TA and the HIP inputs combined almost exclusively nonlinearly. Nonlinear combinations could for example be explained by the fact that multiple intracortical pathways that mediate the TA synaptic input to the recorded neuron were differentially nonlinearly controlled by the given HIP input, even at different time points during a response. The differential effect across the TA input patterns would then moreover be mediated by partly different such pathways as the same HIP input could have widely different impact at the same time point, depending on the TA input.
Profound impact of the HIP stimulation on the TA response types We next made a detailed analysis of the response types evoked in the SI neurons. Different intracellular response types evoked by the same TA input pattern are a striking phenomenon that we have previously reported. 5,13 Response types identified among a set of individual raw responses can be demonstrated with response clustering, and here we explored if these response type clusters formed by the pool of raw responses evoked by a single TA input pattern were impacted by preceding HIP stimulation. Figure 3A illustrates some example response types evoked by the F5 TA input. Figure 3B presents responses to the same TA input pattern but now preceded by the brief HIP stimulation. In this case, the response types evoked by the TA input pattern were profoundly altered compared to Figure 3A. This can for example be seen in that response peaks were eliminated or added at different points in time during the duration of the TA input pattern when it was combined with preceding HIP stimulation.
In agreement with what we have previously reported, 5,13 Figures 3A and 3B show that the differences between the individual intracellular responses evoked by the exact same TA input were very large, to an extent that would imply that response averaging would miss the dominant aspects of the information present in these recordings. Nevertheless, in traditional quantitative terms, all of the 9 neurons recorded from here did display a substantial average EPSP response (>3 mV, peaking no later than 20 ms after the onset of the TA stimulation) for at least one of the TA patterns ( Figures 1A and 2). As we have analyzed in great detail before, 13,19 the different SI neurons displayed different average responses, and response types, to the same given TA inputs (not illustrated).
In order to verify that the differences in response types observed in the above example occurred systematically, we performed a quantitative comparison of the temporal profiles of the responses using a principal-component analysis (PCA)+kNN (k-Nearest Neighbor) classification analysis (see Methods; same approach as that used by Etemadi et al. 5 ). We first compared the clusters evoked by the same TA input pattern but separated based on the condition (i.e., with or without preceding HIP stimulation); hence, in Figure 3 the TA and the HIPTA response types were considered different sets and therefore clustered independently of each other. This analysis showed that the nine clusters of response types evoked by this TA input pattern were systematically distinct, with little confusion (similarity) between the response types (Figure 3C). This was also true for the seven clusters evoked by the F5 HIPTA input condition ( Figure 3D). But the more critical comparison was that between the TA and the HIPTA clusters. Figure 3E shows that all response types, except the TA#5 and the HIPTA#7 response clusters, were well separated from each other.
Whereas the preceding panels only illustrated responses evoked by the F5 example TA input, Figure 3F illustrates that, even when all response clusters for all eight stimulation patterns, across the two input conditions, were compared, the confusion between all of these response clusters was still highly limited. These results hence indicate that the HIP conditioning stimulus profoundly impacted the response types evoked by each given TA input and also confirmed our previous findings that the different response types evoked by each TA pattern is relatively unique compared to the response types evoked by other TA input patterns. 13 Figure 3F also indicates that the HIP condition did not in this respect alter the relationship between the response types evoked by the different TA input patterns, i.e., they were still a relatively unique set for each TA input. This is interesting as it implies that whereas each of the eight TA input patterns typically evoked around 7 response types (as can be read out from Table 1 indicating that the chance level per each pattern was about 15%), under the HIP condition there were still about 7 response types per pattern, but the whole set of response types were also different from those evoked by the corresponding TA input without the HIP condition. This suggests that the HIP stimulation resulted in a whole new set of network solutions to each individual TA input pattern. Black part of the traces represent the time period when there was ongoing TA input (the duration of the actual stimulation is indicated by the gray bars); the preceding gray part of the traces were recorded before the onset of the TA input pattern. (B) Response types evoked in the same cell and using the same tactile input pattern when preceded by the HIP stimulation. The HIP stimulation occurred during the blanked period of the recording. Note that some blanked shock artifacts of the TA input pattern left some remnants, which almost looked like small spikes in the intracellular signals.
(C) Confusion matrix for the nine response types evoked by the F5 tactile input pattern, following a mapping of these response types to principal component space. The presence of a ''non-blue'' diagonal in the matrix indicates that the identified response types, evoked by one and the same stimulation pattern, were highly distinct from each other (as previously analyzed in detail 5,13 ). UC, unclassified responses.

OPEN ACCESS
Whereas the analysis in Figure 3 assumed that the TA and the HIPTA response types were different sets and therefore clustered independently of each other, Figure 4 instead assumes that TA and HIPTA responses were part of the same distribution. Therefore, in this case the clusters were identified without an a priori separation of the responses based on condition. Importantly, the clustering algorithm could still identify clusters which with a very high probability consisted of responses of only one or the other condition. This is illustrated for one example neuron, for two different TA input patterns, F5 and S20, in Figures 4A and 4B. Across all eight TA input patterns, for the same example neuron, this was a consistent finding ( Figure 4C). In our population of recorded neurons (N = 9), almost all tended to a bimodal distribution ( Figure 4D), whereas only two of the neurons had a more Gaussianlike distribution (for example, the neuron represented by brown bars in Figure 4D). In the inset, each neuron is represented individually, and these two outlier neurons are indicated by thicker lines (they were outliers also in the frequency analysis presented in Figure S1). The other 7 neurons had very little confusion between the two conditions, i.e., they had strong bimodal distributions. Overall, these findings underscore that the HIP conditioning stimulation profoundly impacted the response Table 1. Separability of the clusters for each clustering and classification context as measured by the mean F1 score The data in this table are the mean and standard deviation from the entire population of recorded cells. Note that the chance level depends on the number of clusters and varies for each neuron and pattern and cell and that all F1 scores were many multiples above chance level. The HIP suffix refers to TA patterns preceded by HIP stimulation. The ''Pooled'' suffix refers to when response clusters from the TA and HIPTA conditions were compared together. The ''Cell'' context represents classification across all TA and HIPTA clusters for a neuron.

DISCUSSION
Here we showed that brief preceding electrical activation of the subiculum, an input-output hub of the hippocampus, profoundly impacted the responses to tactile input patterns in SI neurons, located a long distance from the subiculum/CA1 terminal region ( Figure 1A). As we have previously described, 5,13 the variety of intracellular response states evoked by each individual tactile input pattern indicates that there is a range of dynamically gated cortical subnetworks supplying the sensory input to each neuron. The cortical network can thereby be described as containing a wide range of network solutions for any given spatiotemporal pattern of sensory input. What we described here is that the specific cortical network solution that applies at a given moment also depends on hippocampal activity. Hence, as further detailed below, hippocampus is part of what defines the potential circuitry mechanisms underlying perception by impacting the predictive or anticipative code in the neocortical circuitry.
Though the relative physical distance between SI and subiculum is large, neurophysiologically the distance is not so large. The subiculum projects monosynaptically to the anterior thalamus. From the thalamus, the subiculum output may reach the SI cortex within just a few synapses. There are also other potential pathways, for example, to hypothalamus and ventral striatum and then to thalamus, which is another potential route by which hippocampal output can impact SI activity indirectly. Even so, the magnitude of the impact iScience Article that we observed here was surprising. The synchronized activation that was achieved by the hippocampal microstimulation likely made the effect more powerful than a normally generated hippocampal output. On the other hand, the observed effects remained for several 100 ms after the termination of the activation of the output. In addition, a previous study showed that there can be a strong correlation between SI neurons in deep cortical layers and hippocampal neurons during oscillatory activity. 20 And, hippocampal neurons can be active at high firing rates, 21,22 not seldom in relative synchrony. 23 This suggests that hippocampal output can continually modify the global cortical state with high potency, including impacting the responses of SI neurons to tactile stimuli.
The hippocampus is by many authors considered to be at the ''top'' of the hierarchy of cortical processing, and it is considered critical for memory consolidation, i.e., to drive synaptic plasticity in the neocortex via replay of previous experiences. 24 In comparison with our previous study of the perturbation of SI tactile processing generated from a remote cortical area, one would thus expect quite different effects from the hippocampal output, for example, qualitatively different impacts on the intracellular signal in the SI neurons. However, in short, the effects we found were instead remarkably comparable to those we obtained from a remote location in the neocortex. 5 In comparison, the hippocampal output seemed potentially more potent in generating long-lasting oscillatory activity (Figures 1 and S1), possibly through a more efficient synchronization of the global thalamic activity, in turn possibly due to the impact of the hippocampus on the anterior thalamus. 11 Otherwise, we observed qualitatively similar nonlinear effects, i.e., that the HIPTA responses were strongly nonlinear combinations of the TA and the HIP responses. Such nonlinear effects are in turn indicative of that the perturbation impacted the gating of transmission of tactile information from wide areas of the neocortex to the SI neurons. 5 But the fact that the impact of hippocampal output here did not qualitatively differ from neocortical output raises the interesting possibility that hippocampal cortex may not exert such a different neocortical impact compared to neocortex itself, in neurophysiological terms, in the time window of 1 h.
Our low-noise intracellular recordings represented the temporal evolution of the combined activity of thousands of neurons being afferent to the recorded neuron, 25 and the specific trajectories observed reflect how those thousands of inputs are differentially gated. As we have previously underscored, the fact that the cluster analysis identified multiple preferred response states does not imply that these response states are fixed entities or the only response states possible. 5,13 The cluster analysis merely indicates that there are some central response features that some responses share, whereas at least one of these features was not shared by the members of the other clusters. The unique, in fact, the orthogonal, nature of these response states ( Figure 3; Table 1) indicates that these features may correspond to multiple independent network solutions, which are dynamically gated as a product of the cortical internal state and how it combines with the external sensory input. The relative orthogonality between these response states, as indicated by the PCA/decoding analysis, indicates that the network in this regard may work according to attractor dynamics, 5,13 conceptually introduced in a recent review. 26

Broader implications of our findings
By focusing on the trial-to-trial response variability, rather than different forms of response averaging, which is otherwise a very common neurophysiological approach, our approach takes the alternative angle of addressing how the hippocampal-driven cortical state would impact the SI neuron responses to given sensory input patterns. Hence, we move beyond the idea of considering the neocortex as essentially driven by external stimuli but instead focus on what degree of response variation can be generated as a consequence of the cortical state. This approach is in agreement with recent perspectives arguing that the paradigm of a sensory-driven cortex is by itself inadequate and that we need to address how changes in activity intrinsic to the brain impact the representation of sensory input, 27 which is in turn considered needed to explain cognition. This approach is well in line with recent neuroanatomical studies indicating that the profusion of global interconnections across the neocortical neuronal networks has been vastly underestimated 7 (https://www.janelia.org/project-team/mouselight; https://mouse.braindatacenter.cn/) and that also the cortico-thalamo-cortical pathways are much more widespread 28 than what is typically accounted for in models of cortical network operation. Such profuse connectivity would seem to inevitably lead to widespread representations of any cortical signals, including those driven by sensory inputs. Indeed, several recent analyses at the cellular level indicate that the cortical representation of tactile information is dependent on globally distributed networks. 6,17,18,29,30 Hence, these earlier findings speak in favor of the cortex being a highly recurrent network, rather than having a predominantly feedforward network architecture. Here we show that this recurrent network seems to include even the hippocampus. iScience Article In a recurrent cortical neural network architecture, it is to be expected that a tactile input pattern, or in principle any other sensory input pattern, would combine with the ever-evolving state of the neocortical network. This state would equal a component of the ''prior,'' ''expectation,'' or ''prediction'' currently residing in the neocortical network. That internal state component will be an important determinant of the resulting activity distribution across all neurons of the cortical circuitry that results when the tactile input pattern arrives-which is a potential definition of a ''percept.'' Our results indicate that the diversity of cortical response states is profoundly impacted by hippocampal output. Whatever information that is being processed in the hippocampus will hence to some extent impact how we perceive a given tactile input, which for example could contribute to semantic and emotional aspects in haptic perception. 31-33

Limitations of the study
As we have extensively analyzed and discussed elsewhere, 5,13 the anesthesia inevitably impacted the recorded responses but not to an extent that would invalidate the main finding that hippocampal output shifted the global cortical state to an extent that it would alter the SI neuron representations, at the integrated level translating to the percept, of given tactile input patterns. Furthermore, the physiological network structure, quantified as the recruitment order of a high number of simultaneously recorded cortical neurons across different sources of network excitation, does not alter with the type of anesthesia used here. 34,35 Also, the type of network attractor dynamics that we observed here is likely closely related to the type of neuron population-level dynamics observed in the hippocampus, 15 which remains stable whether the animal is awake or in deep sleep. 16 Notably, using the same anesthesia as we have used here, highly detailed features of visual cortical processing have been discovered. 36,37

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

Lead contact
For any additional inquiry and information related to materials and resources used in this work Henrik Jö rntell is the lead contact (Henrik.Jorntell@med.lu.se).

Materials availability
This study did not use any new/unique materials and/reagents. d Further information about the analyzed data discussed in this paper will be provided by the lead contact upon demand. Code to perform the analysis will be provided by the lead contact on demand.

EXPERIMENTAL MODEL AND STUDY PARTICIPANT DETAILS
We made recordings from 9 anesthetized male adult Sprague-Dawley rats (weight: 250-380 g, age: 12 G 2 weeks) in acute preparations. Before experiments, animals were maintained by Lund University animal facilities with 12h light/dark condition, 2-3 animals per cage (type 3H), and with free access to food and water.

Institutional permission
Ethical approval for this research was obtained in advance from the local animal ethics committee in Lund/ Malmö (Ethical permit ID: M13193-2017).

Surgical procedures
In order to make acute in vivo recordings, adult Sprague-Dawley rats were initially prepared in the same way as in previous studies. 5,13 Briefly: 1) The animal was sedated by inhaling air mixed with isoflurane gas (3%, for 2 min); 2) To induce general anesthesia, a mixture of ketamine/xylazine (ketamine: 40 mg/kg and xylazine: 4 mg/kg, accordingly) was injected intraperitoneally; 3) an incision in the inguinal area of the hindlimb was made to insert a catheter in the femoral vein for continuous infusion of Ringer acetate and glucose mixed with anesthetic (ketamine and xylazine in a 20:1 ratio, delivered at a rate of 5 mg/kg/ h ketamine). 4) A hemicraniectomy ($2 3 2 mm) was made on the right side of the skull at the area covering the somatosensory cortex (SI) on the right-hand side ($2 3 2 mm), located at (from bregma): Ap: -1.0 -+0.1, ML: 3.0-5.0. iScience Article A second cortical exposure over the visual cortex was made to gain access to the caudal CA1 and the subiculum complex of the hippocampus, located at (from bregma): AP: À6, ML: 3. The stimulation electrode was advanced 5.8 mm at an angle of 8-10 to reach the subiculum complex. For the entire experiment, the EcoG was recorded via a surface electrode placed in the second cortical exposed area. To prevent tissue dehydration, the second exposed cortical area was covered with liquid paraffin oil. To stabilize the tissue and to prevent dehydration, the first cortical exposure was covered in agarose.

REAGENT or RESOURCE SOURCE IDENTIFIER
The anesthetics used were chosen because ketamine has previously been reported to not dramatically alter the neuronal recruitment order in spontaneous activity fluctuations and in stimulation-evoked responses as compared to the awake animal. 34 As we have previously discussed extensively, 13 anesthesia was required in order to achieve identical tactile stimulation patterns (where the electro tactile interface was the key, but would not be accepted in the awake animal) over a sufficiently long period of time. It also served to minimize brain activity noise caused by uncontrollable movements and internal thought processes unrelated to the stimuli. The level of anesthesia was assessed both by regularly verifying the absence of withdrawal reflexes to noxious pinch of the hind paw and by continuously monitoring the irregular presence of sleep spindles mixed with epochs of more desynchronized activity, a characteristic of sleep. 38 The animal was sacrificed by an overdose of pentobarbital at the end of the experiment. For histological evaluations, the animals were perfused with 4% paraformaldehyde (PFA) and the brain was removed for further processing (see below).

Data collection: In vivo neuronal recordings
In order to analyze synaptic inputs, we performed whole cell (intracellular) recordings in the current clamp mode. Patch pipettes were made from borosilicate glass capillaries pulled to 6-10 MU of impedance, using a Sutter Instruments P-97 horizontal puller. The pipettes were back-filled with an electrolyte solution of (in mM): potassium gluconate (135), HEPES (10), KCl (6.0), Mg-ATP (2), EGTA (10). The solution was titrated to pH of 7.35-7.40 using 1 M KOH. Recording signals were amplified using the EPC-800 patchclamp amplifier (HEKA Elektronik) in the current-clamp mode (bandwidth from DC up to 100 kHz). Throughout recording sessions, time continuous data was acquired and digitized at 100 kHz using the CED 1401 mk2 hardware controlled from the Spike2 software (Cambridge Electronic Devices, CED, Cambridge, UK).
The patch electrode was inserted into the SI area 39 and advanced in a stepwise manner with applied positive pressure to prevent blockage of the electrode tip. The pipette was advanced in a slow forward movement (0.3 mm/sec) until a sudden change in the electrode resistance was detected. When a dramatic increase in the magnitude of the isolated spike, evoked or spontaneous, occurred, the positive pressure was switched to a brief negative pressure. Also, to create a more optimal condition for the establishment of a GigaOhm seal, a mild hyperpolarizing current (in the order of À10 pA) was applied to the electrode. Then, to gain access to the intracellular space, episodical negative pressures were applied to the electrode tip. The data recording started after obtaining low noise access to the intracellular environment with a stable signal. The recordings were characterized by stable membrane potential of < -55 mV in down states without any bias current, and with a peak-to-peak value between up and down states of >10 mV, with spike amplitudes of >25 mV in the beginning and end of stimulation protocols. During the stimulation protocols, the cells were prevented from firing with a mild hyperpolarizing current. All neurons recorded were putative pyramidal neurons rather than interneurons based on that they exhibited infrequent bursts of two or three spikes but had an absence of longer bursts or sustained periods of high firing. 35 All recordings were made in the neocortical layers II/III-V (350-1100 mm of depth).

Design: Electrical tactile and cortical stimulations
Once an intracellular recording was obtained, the experiments consisted of delivering tactile stimulation patterns through an electro tactile interface, combined with hippocampal electrical activation. The electrotactile interface consisted of four bipolar pairs of percutaneous stainless steel needle electrodes (isolated except for the tip (0.2-0.5 mm) inserted into the superficial part of the skin of the second digit of the contralateral forepaw). Each pair represented a single channel for providing tactile stimulation ( Figure 1A). Each electrode pair delivered a 400 mA constant current pulse with a duration of 200 msec, which is about two times the threshold for activating tactile afferents using this type of stimulation, 40

QUANTIFICATION AND STATISTICAL ANALYSIS
Post-processing of recorded neuronal data The raw neuronal membrane voltage recordings were post-processed by first defining a stable baseline voltage by applying a least-square linear detrending to 10 s long segments. The recordings were subsequently low pass filtered using a rolling average of 100 bins equating to a boxcar window with a duration of 1 ms. In recordings where occasional spikes appeared, the onset of these spikes were found using a spike shape template and a recursive fitting algorithm. 44 Stimulation artifacts and any spurious action potentials were removed by linear blanking. Finally, the recordings were down sampled to 1 kHz. Overall, the software used in the below analysis was NumPy and SciPy, standard Python libraries in data science.

Definition of evoked responses
The time window included in the analysis started at 5 ms after the onset of the TA stimulation (i.e., for both TA and HippTA stimulation patterns, corresponding to the conduction time from the skin stimulation to the earliest possible arrival of synaptic responses in the neocortex) and ended 400 ms later (some TA patterns lasted almost 350 ms, and responses could sometimes be detected at least 50 ms after the last stimulation pulse).

Clustering method
We used the same method to cluster the responses as we used in a previous publication (Etemadi 2022). This method identifies clusters of similar responses, with respect to their time-voltage curves, evoked by the same stimulation pattern, in an unsupervised manner. Hence, the algorithm was designed to identify clusters where the 'members' of each group were internally similar, while also being distinct relative to the other responses evoked by the same stimulation pattern. Furthermore, the algorithm had to be able to automatically determine the number of clusters that could exist, since this could not be known a priori. iScience Article First, each response was Z score normalized, i.e., the response mean over the whole time period was subtracted from each response and then the remainder of each response was divided by the standard deviation. Secondly, a distance matrix was calculated for the whole set of responses. The distance matrix was populated by calculating the M number of Principal Components (PC) explaining 95% of the variance for N-1 responses. All responses were then fitted to the PCs using the least square method. The resulting PC coefficients were used to position the responses in the M-dimensional PC space. Within this M-dimensional space, the Euclidean distances between a reference response and the rest of responses were calculated. The distances between the reference response and each of the other responses were then put in a matrix and the procedure was repeated until each response had been the reference. Finally, the distances were normalized to a 0-1 range.
From the Euclidean distance calculations above, a linkage matrix was calculated using the hierarchical clustering approach of Ward's minimum variance method. 45 This method sequentially finds the two closest neighboring responses with respect to their Euclidean distance. Then it finds the pair of the second closest neighboring responses and so on until there are no responses left. However, a response may also be closer to a center point between a pair of responses than it is to other responses, in which case the distance to that center point is assigned to the response. This distance is called the cophenetic distance. Sometimes, the shortest distance that can be found is between two such center points. All these distances are used to create the hierarchical structure of responses, where the distances between all such pairs can be compared.
The next step is to identify the cophenetic distance at which the identified clusters are objectively the most separable. The total cophenetic distance of a hierarchy was divided into 100 steps. Each level of cophenetic distance was used as a threshold, which would separate clusters from each other (i.e. a cluster of responses with a cophenetic distance larger than the threshold was considered a unique cluster), whereby an assignment of a cluster id was made. At this point we could perform data shuffling if needed. Data shuffling was achieved by shuffling the leaf position of each response in the hierarchy (the shuffling was performed using the NumPy shuffle method), therefore disrupting the relationship between cluster and response.
We next calculated the separability of the clusters by using a decoding analysis (described in detail in 'Evaluation of cluster separability using decoding analysis' below). The decoding analysis yielded, for each clustering threshold, a measure of the separability of the clusters as the F1-score. With the F1-score, we could apply Gap statistics to identify the cophenetic distance providing the largest cluster separation relatively to the shuffled data. Meaning that for each of the 100 cophenetic distances considered, we obtained an F1 score for the non-shuffled data and another F1 score for the shuffled data, subtracting the F1-score of the shuffled data from the F1-score of the actual data yielded the difference also known as the Gap Statistic. 46 The calculation of the Gap Statistic was repeated three times per clustering threshold and the mean value was reported. Since the specific value of the chosen cophenetic distance defines the number of clusters to be considered, we could identify continuous stretches of cophenetic distances where the number of clusters, and the clusters themselves, remained unchanged. Over each of these continuous stretches we calculated a mean Gap Statistic. The middle value of cophenetic distance of the stretch with the maximum Gap statistic was used to define the cophenetic distance at which the clustering was objectively the most distinct. The clustering obtained at this value was then used for the remainder of the analysis and displays.
While following the above principles, there were some special situations that could arise. A cluster was assumed to be valid only if it had five or more members. All non-valid clusters were put in the same undefined cluster ('unclassified' responses). In addition, if the situation arose that one of the clusters consisted of only one member this response and cluster was excluded from all subsequent analysis and displays.

Evaluation of cluster separability using decoding analysis
We used the same general decoding algorithm to determine the specificity of the response clusters as in several previous publications, 6,17,18,29 including analysis of intracellular recording data with multiple response states to each tactile input pattern. 13 The responses and their respective label were split into a training-and a test set (80%/20% split) stratified by cluster label frequency. The training set was used to train a PCA model explaining 95% of the variance, and the coefficients obtained by fitting the responses with the least-square-method to the principal components were calculated for both the training set and the test set. Finally, a kNN- The decoding analysis was also repeated with the randomly shuffled cluster labels (as described in Clustering method), reported as the ''shuffled'' context. This was done to establish the actual chance level of the data, as opposed to the theoretical chance decoding level (which equals 1 divided by the number of clusters