Ethical Statement

All the animals have been used in accordance to the Italian and European Laws on animal treatment in Scientific Research (Italian Bioethical Committee, Law Decree on the Treatment of Animals in Research, 27 Jan 1992, No. 116, guideline n. 86/609/ European Community). The National Research Council, where the experiments have been performed, adheres to the International Committee on Laboratory Animal Science (ICLAS) on behalf of the United Nations Educational, Scientific and Cultural Organizations (UNESCO), the Council for International Organizations of Medical Sciences (CIOMS) and the International Union of Biological Sciences (IUBS). The research has been approved by the Ministry of Health and classified as “Biella 1/2014” into the files of the Ethical Committee of the University of Milan.

Electrophysiology

Fifty male albino rats (Sprague-Dawley, Charles River, Calco, LC, Italy, 270–350 g) were chosen out of the set of the animals employed in the research.

The rats underwent preliminary barbiturate anaesthesia (50 mg/kg ip) for the surgical experimental preparation. The trachea was cannulated to gain a connection to the anaesthesia-ventilation device. A 16-gauge butterfly was then positioned in the root of the lateral tail vein. A surgical opening was then done over electrodes the skull vault with the removal of the skin the galea capitis muscles and of their fibrous insertions over the parietal bones. The periosteum was then delicately scraped off from the parietal, frontal and occipital bones. The rats were then placed in a stereotaxic frame and the animal were paralyzed by intravenous Gallamine thriethiodide (20 mg/kg/h) injection and connected to the respiratory device delivering (1stroke/s) an Isoflurane® 2.5% 0.4 to 0.8 l/min in Oxygen 0.15–0.2 l/min gaseous mixture. Curarization was maintained stable throughout the whole experiment by Gallamine refracted injections. During the experiment the anaesthesia level was continuously monitored by the EEG recordings by four leads placed, respectively in a fronto-occipital row (in stereotaxic coordinates referenced to the bregmatic point, where F, P and O refer to Frontal, Parietal F + 1.2 mm AP-2.2 mm LL; P1: −1.2 mm, 2.5 LL; P2: 3 mm AP, 2.7 mm LL; P3: 5.5 mm AP, 2.7 mm LL; The fifth, reference lead, was placed on the occipital bone posterior to the Lambda reference point. O: −2 from Lambda, 2 mm LL). Briefly, five partial bone holes have been trephined at those reference points and five external leads were then placed in the bone cradles and connected to the recording system. A hypertonic salted cream was used to ameliorate the electric interface between the leads and the bone. During the experiment the EEG recordings continuously monitored the anaesthesia level.

We simultaneously recorded spiking and local field potential activities were simultaneously recorded from anesthetized rats, by two microelectrode matrices from three stations in the brain: the thalamic ventro-postero-lateral complex nuclei (VPL) and the primary somatosensory (S1) cortex. We also concurrently recorded from four EEG derivations. The neuronal electrophysiological recordings have been obtained contralaterally to the stimulated paw; the concurrent EEG recordings were obtained ipsilaterally to the stimulated paw.

The neuronal recordings were obtained by two matrices of extracellular tungsten or Pt-Ir electrodes were framed in 3 × 3 arrays of single shanks, inter-tip distance 150–200 mm, tip impedance 0.5–1 MΩ (FHC Inc., ME, USA). The coordinates have been estimated from the Paxinos – Watson Stereotaxic Atlas of the Rat Brain88. In detail, the two matrices were placed at −1.2 mm AP, 2.6 mm LL and −6 mm AP, 2.5 mm LL for the somatosensory cortex and the thalamic nuclei respectively. The thalamic regions were targeted with a postero-anterior slant of 25° of the matrix to avoid spatial interferences with the cortical matrix. This obliged to recalculate the depth by a simple geometric correction of the usual measure. The cortical matrix was inserted 400 μm deep at the superior border of the granular layer and then slowly advanced, by an electronically controlled microsteppers (Narashige, Japan), at 10 μm steps for probing the responses of local neurons to exploring sensory stimuli on the contralateral posterior paw, until a clear response was evident and repetitive on at least six out of the 8 recording microelectrodes of the matrix within a final depth of 600 μm. The thalamic probe was advanced electronically by a second electronically driven microstepper (Transvertex, Stockholm) from a starting depth of 4700 μm down to 5800 μm (fastly driven to avoid tissue damage). The concurrent EEG recordings were obtained ipsilaterally to the stimulated paw.

EEG recordings

The EEG electrodes were placed along the stereotactic coordinates (in front-back order, Bregma as relative zero AP, zero ML reference point) as follows: FrontalAnterior (FA) +3 mm MedioLateral (ML) −2 mm, ParietalAnterior (PA) −0.15 mm, ML −2.8 mm, MiddleParietal (MP) −3 mm, ML. −3 mm, PosteriorParietal (PP) −6 mm, ML. −2.8 mm, reference in Occipital bone −9 mm, ML. 2 mm. (121). For the analyses, we selected preferentially the EEG data from the second derivation placed over the sensory cortex mirroring the contralateral somatosensory primary cortex where the neuronal recordings were obtained.

Electrophysiological extracellular recordings

For the electrophysiological recordings, two holes were drilled on the skull of 3 mm2 for the cortical and the thalamic matrix accesses. The holes were drilled centered respectively on the cortical access centered around a reference point at −1.5 mm AP and 2.5 mm ML, and the electrode was driven around 450 to 800 micrometer deep by an electronically controlled microsteppers (Narashige, Japan). The thalamic access hole was centered at −6 mm AP and −2.5 mm ML. The thalamic matrix was inserted with a slant at 25° and driven at least at 5500 micrometers in depth and then advanced electronically by a second electronically driven microstepper (Transvertex, Stockholm) until full responses were observed to peripheral test stimuli. Fast thalamic and cortical responses to light tactile stimuli with a brush-test on the sciatic innervation field (the plantar aspect of the left hindlimb) were the anatomo-functional acceptance criteria for acquisition. All the experimental blocks were organized with periods of ongoing activity recordings lasting around 20 min. and not more to preserve at most the data homoscedasticity, in the additional stable conditions of gaseous anesthesia. After a cycle of spontaneous and stimulated activities was completed, we repeated twice the original recordings. Then we advanced in depth the electrodes, 20 mm and 100 mm respectively for the cortical probe and the thalamic matrix ensemble, to reach an adjacent region, then checking again with the test stimulus the responsiveness of the newly recorded regions. In positive cases we repeated the recording cycle as above. We recorded from five to six stations in progressive steps for each animal.

For signal amplification and data recordings a 32 channel Cheetah Data Acquisition Hardware was used (Neuralynx, MT, USA, sampling frequency 32 kHz). Electrophysiological signals were digitized and recorded with bandpass at 6 kHz and 300 Hz for spikes and at 475 Hz–1 Hz for the EEG. The data stored were analysed off-line both by Matlab and by locally developed software. A histological confirmation of the placement of the electrodes was then obtained on brain coronal sections stained with cresyl-violet.

Preliminary data analyses

For signal amplification and data recordings a 40 channel Cheetah Data Acquisition Hardware was used (Neuralynx, MT, USA, sampling frequency 32 kHz). Electrophysiological signals were digitized and recorded with bandpasses at 6 kHz/300 Hz for spikes, 180 Hz/1 Hz for Local Field Potentials. The data stored were analyzed off-line both using Matlab and by locally developed software. The neural firing rates had a mean of 31.4 Hz with standard deviation of 26.8 Hz.

After the recordings the LFPs were downsampled to 0.5 KHz. We used for filtering the same techniques described in. After filtering and downsampling, the spike contamination of LFP signals was null avoiding further spike removal techniques. The spikes were extracted and sorted by using the Wave_clus MATLAB toolbox. Sorted cells with average rates below 4 Hz and above 100 Hz were excluded from the analysis. Furthermore, neurons resulted from sorting which had improbable inter-spike-interval distributions were discarded as well. Recorded neurons were uniformly distributed over the recording matrices and every electrode show distinct neural activity otherwise the matrix was repositioned. At the end of this process, we collected a total of 391 neurons (56 ± 17 in each experiment) out from the set of the acquired signals.

The timestamps of spike occurrences were represented by binary sequences where 1’s labeled a spike. We considered time bins of 1 ms thus avoiding occurrence of multiple spikes within the same bin.

Finally, we split each sequence into fixed-length (100 ms) overlapping windows (Fig. 1), thus obtaining an ordered set of equal length windows.

In order to discriminate groups of recording by the firing rates, we verified the intraclass consistency of the recorded spiking activity within each experimental condition. To this purpose we performed a one-way analysis of variance on the four sets of experiments. No significant difference within each class was observed.

Functional connections by spike-train similarities

Interactions between neurons can generate very complex, time-delayed and asymmetric patterns especially in the thalamocortical circuitry. In this work, we reproposed a framework successfully applied in a similar context wherein Markov stochastic models model spike trains89.

We used the function Normalized Compression Similarity (NCS), formally defined as: given that x and y are two spike trains, the NCS is equal to:

where the C function represents the compressed sequence length and is the sequence concatenation operator (e.g. 0101 · 101 = 0101101). If NCS(x, y) is close to 1, the sequences x and y are considered similar. If close to 0, the sequences are strongly dissimilar.

We evaluated the NCS function on time windows (100 ms length) of the recorded spiking activity assuming that relative high values of similarity corresponded to actual functional connections (Fig. 1E).

Functional connections by LFP phase synchrony

LFPs are low frequency signals reflecting a wide range of synaptic events. In this work we investigated the synchrony of LFP phases originated in different recording sites during spontaneous and tactile evoked activities. We measured phase synchronies between two recorded LFP sequences (x and y) by the following function

where e is Napier’s constant, H is Hilbert Transform, arg is the argument function and i is the imaginary unit90. The Hilbert transform and the argument were computed with, respectively, the hilbert and the angle Matlab functions. When γ (x, y) is equal to 1 (0), then x and y are perfectly synchronous (asynchronous).

Complex Brain Network

By using the NCS and γ functions, we estimated the functional connections of the recorded neuronal networks. We first split each recorded sequence into 100 ms time windows (Fig. 1) and then we computed the adjacency matrix for all neurons or electrodes. The resulting matrices exhibited values in the unitary interval. The functional connections extracted from extracellular recordings are the counterpart of non-oriented graphs.

Since in previous work we noticed that variable thresholds equal to a higher percentile of the weight distribution and vary between 0.2 and 0.8, did not affect result we selected the 75 percentile.

For the analysis of these graphs, we introduced a set of common statistics from the Complex Network Theory able to detect possible matches between the extracted graphs and prominent topologies like small-world networks (Table 1)91,92. A small-world network is generally obtained by evolving a basic ring lattice graph, where each node is connected to their K neighbours. The chosen neighbourhood involves typically much less nodes than the total node number N. Randomly adding and removing edges from the starting graph achieve the graph evolution. The resulting graph has many, typically small, quasi-complete subgraphs (cliques) where each node is connected to every other node in the clique. Furthermore, small-world networks exhibit short average distances between nodes46,47.

From a functional perspective, small-world networks can express two important information-processing features: information integration and segregation. Functional segregation recruits specialized processing within densely interconnected nodes (cliques). Functional integration combines information processed in distributed nodes or cliques. These network properties can be measured by two statistics: the clustering coefficient (C) and the characteristic path length (L). The former measures how close the neighbours of a node are to being a clique. The latter estimates the average shortest path length in the graph, i.e. how much the nodes are accessible. Both measures, implemented in a Matlab toolbox, were used for our network analyses (clustering_coef_bu.m, charpath.m)93.

In complex network theory, several graph measures take specific meaning only if they are compared to the same graphs subject to randomization or latticization (often called null networks). Both procedures guarantee that the node degree distributions of the original graphs were preserved. We computed, by using the Matlab function randmio_und.m, the randomized version of our graphs and we estimated Cr and Lr. These null network values are required to verify the small-world nature of the graphs.

The functional graphs obtained by our analysis were further characterized to study the information flow. For this aim, we computed a measure of centrality (betweenness) for graph nodes, an estimate of the number of shortest paths from all vertices to all others that pass through that node. Because it can be interpreted as a measure of the load of a node within the network, the distribution of node centrality highlights how the information flow is balanced within graphs.

Ultimately, we analyzed networks that evolved in time dropping and recruiting nodes and connections and networks from different experimental conditions. Such a methodology requires the discussion of potential issues.

First, unconnected nodes were rare but could occur after adjacency matrices were binarized. For this reason, we removed graphs in which less than the 99% of nodes were connected.

Second, network statistics were applied on network with different sizes (for spiking activity) because the recording sessions returned a variable number of active neurons. However, by analyzing the observed variance of network size we concluded that C and L couldn’t be significantly affected by our network size changes. Significant changes appeared for synthetic networks that increased their size by orders of magnitude. However, we discarded graphs that were outliers (beyond 5th and 95th percentile) of the node, edge and density number distributions in order to obtain a better homogeneity49. In the work, we refer to these two conditions as admissibility criteria. Samples of functional graphs (represented by their weighted and binarized adjacency matrices) from different groups are shown in Fig. 7.