Small-World Networks (SWNs) represent a fundamental model for the comprehension of many complex man-made and biological networks. In the central nervous system, SWN models have been shown to fit well both anatomical and functional maps at the macroscopic level. However, the functional microscopic level, where the nodes of a network are represented by single neurons, is still poorly understood. At this level, although recent evidences suggest that functional connection graphs exhibit small-world organization, it is not known whether and how these maps, potentially distributed in multiple brain regions, change across different conditions, such as spontaneous and stimulus-evoked activities. We addressed these questions by analyzing the data from simultaneous multi-array extracellular recordings in three brain regions of rats, diversely involved in somatosensory information processing: the ventropostero-lateral thalamic nuclei, the primary somatosensory cortex and the centro-median thalamic nuclei. From both spike and Local Field Potential (LFP) recordings, we estimated the functional connection graphs by using the Normalized Compression Similarity for spikes and the Phase Synchrony for LFPs. Then, by using graph-theoretical statistics, we characterized the functional topology both during spontaneous activity and sensory stimulation. Our main results show that: (i) spikes and LFPs show SWN organization during spontaneous activity; (ii) after stimulation onset, while substantial functional graph reconfigurations occur both in spike and LFPs, small-worldness is nonetheless preserved; (iii) the stimulus triggers a significant increase of inter-area LFP connections without modifying the topology of intra-area functional connections. Finally, investigating computationally the functional substrate that supports the observed phenomena, we found that (iv) the fundamental concept of cell assemblies, transient groups of activating neurons, can be described by small-world networks. Our results suggest that activity of neurons from multiple areas of the rat somatosensory system contributes to the integration of local computations arisen in distributed functional cell assemblies according to the principles of SWNs.

Cell assemblies (sequences of neuronal activations), seem to represent a functional unit of information processing. However, it remains unclear how groups of neurons may organize their activity during information processing, working as a sole functional unit. One prominent principle in complex network theory is covered by small-world networks, in which each node is easily reachable by each other and organized in highly dense clusters. Small-world networks have been already observed on large scales in human and primate brain areas while their presence at the neuronal level remains unclear. The aim of this work was to investigate the possibility that functional, related neural populations, encompassing multiple brain regions, could be organized in small-world networks. We investigated the coherent neuronal activity among multiple rat brain regions involved in somatosensory information processing. We found that the recorded neuronal populations represented small-world networks and that these topologies were maintained during stimulations. Furthermore, by using simulations to explore the hidden substrates supporting the observed topological features, we inferred that small-world networks represent a plausible topology for cell assemblies. This work suggests that the coherent activity of neurons from multiple brain areas promotes the integration of local computations, the functional principle of small-world networks.

Funding: The research has been partially funded by PON project 01-01297(Virtualab) funds. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Copyright: © 2013 Zippo et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Our results confirm the presence of SWNs in crucial stations of the rat somatosensory system, during spontaneous activity and provide evidence that, during stimuli, small-world organization principles are preserved in spite of massive topological reconfigurations. Given that the hidden functional substrate which supports the observed networks remains elusive, we computationally investigated the hypothesis that the Hebbian structure of cell assembly is organized in small-world networks. By means of these computer simulations, we showed that the observed conditions may be the expression of a large-scale functional substrate of cell assemblies represented by small-world topologies, where node membership to each group is variable. This last result was further supported by verifying the computational models either with random or with lattice networks that did not produce consistencies with experimental data.

In the present work we address these points by means of spike and LFP simultaneous multi-array recordings from the ventropostero-lateral thalamic nucleus (VPL), the centro-median thalamic nucleus (CM) and, the primary somatosensory cortex (S1) [24] – [26] . We analyzed VPL, CM and S1 because these regions are actively involved in the tactile information processing in somatosensory system of mammalian brains. Because spikes and LFPs represent a dualistic picture of neurons and local neuronal population states, we decided to test if small-world topologies emerged from these two entities. In order to estimate the strength of functional connections in spiking activity, we introduced a novel function that allows us to detect long-range temporal dependencies which occur in thalamocortical interactions (Normalized Compression Similarity, NCS) that are otherwise unrecognizable by correlation analyses (see Materials and Methods ), the LFP functional connections were instead estimated by a standard measure of phase synchrony.

In the last 15 years, thanks to a substantial advancement in the complex network theory, Small-World Networks (SWNs) emerged as a paradigmatic model and provided the analytical tools to explain a large set of complex networks in the most diverse scientific areas [6] , [7] . Furthermore, it has been suggested that small-world networks represent optimized structures accomplishing adequate balance between information transfer efficiency and reliability [7] – [12] . Following this trend, large networks of brain areas, studied by imaging techniques, have been shown to functionally organize as SWNs [13] , [14] . Although recent works provided evidence of SWNs in small local neuronal populations, the functional connection graphs emerging at the microscopic level of single neuron are still poorly understood [15] – [23] . Specifically it is not clear whether functional groups of neurons from multiple brain regions display a SWN organization.

Neurons in the brain form highly composite networks sustained by a complex and variously distributed thread of synapses. Although the detailed anatomical connections are still under investigation [1] it has been shown that, during in-vivo recordings, active neurons form functional assemblies that do not necessarily depend on the underlying anatomical connectivity [2] . Therefore, while anatomical connectivity is stable for relatively long times, the functional connections are highly dynamic and depend on the particular tasks triggered by internal and external events. Critically, the organization principles that control the functional connections among single neurons and small neuronal populations are still poorly understood [3] . The early hypotheses [4] , that distributed groups of neurons can operate as a functional unit through coordinated activities, is now supported in observations of transient cell assemblies over different cortical regions [5] . It remains, however, unclear how these dispersed neurons may gather in a functional unit for information processing.

From our results we conclude that neural activity may be structured in overlapping functional groups, each of them organized as a small-world network. Finally, we computed the admissible values of small-world networks per number of nodes that ranged within and , keeping consistency with the and values observed in experiments.

At each run we randomly extracted 100 units (comparable with the number of units that we could simultaneously record during the experiments) and we computed the network statistics. We observed that small-world statistics underwent a bimodal behavior ( Table 5 ), increasing for , peaking when and decreasing with up to . At peak we found comparable statistics values with those measured experimentally both on spiking ( , ranksum test) and LFP activities ( , ranksum test). To support these findings, we evaluated two opposite null hypotheses replacing small-world networks with ring lattices and random networks. For the former we used the same algorithm of Watts and Strogatz setting p (the probability of edge rewiring) equal to 0 while for the latter we set . In both cases, we did not find consistency with the experimental observations.

For these purposes, we implemented a large-scale simulation of neurons that embedded small-world networks. Then we estimated small-worldness by systematically varying . We simulated a neural network of neurons, a biological relevant neuronal ensembles exhibiting a magnitude about five fold the estimated size for a rat cortical hypercolumn [38] . The size of each functional group was chosen by sampling from a discrete uniform distribution . This distribution and its parameters are consistent with the findings of [39] . For further details about the implementation see the relative Materials and Methods section.

We thus investigated the hypothesis that neuronal networks are arranged in functional groups, recalling the concept of cell assembly, where the membership to each assembly is variable [4] and assembly neurons can be spatially distributed. To test this hypothesis we developed an artificial neuronal network where nodes were arranged in dynamical groups and emitted spikes following small-world criteria. We then sampled the activity of small ensembles of nodes and compared the network statistics from the synthetic spiking activity with those obtained by recordings. Once verified such hypothesis, we estimated how many functional groups can coexist within a neural network while keeping consistency with experimental observations.

In the previous sections we provided evidence, both at a pre- and post-synaptic level, that cortical and subcortical neurons can be functionally organized as small-world networks. Our results are consistent with recent findings from large scale cortical recordings [10] , [36] , [37] and suggest the existence of sparse functional networks composed of neuron assemblies that encompass multiple cortical and subcortical areas. It is of utmost importance to consider that each functional network can recruit its units not necessarily in close proximity [5] , [27] , it remains although unclear what kind of functional substrate supports the observed networks and specifically, how large-scale networks of neurons are organized such that small-world configurations can be observed sampling from a small subset of nodes.

Moreover, by analyzing the modularity index of the community structures, we found that graphs in pre-stimulus conditions had larger Q indices ( , ) than graphs in post-stimulus conditions ( , ; , ranksum test). In conclusion, the incoming stimulus forced delegated networks to merge the active communities indicating that more neural groups expressed coordinated oscillations.

By analyzing the distribution of node betweenness centrality, we found that graphs in pre-stimulus conditions ( , ) had smaller betweenness centrality ( , ranksum test) than graphs in post-stimulus conditions ( , ) indicating that stimuli cause nodes to be more recruited. Furthermore, we found that the betweenness centrality was not equally distributed over the recorded regions (VLP, S1, CM) both in pre- and in post-stimulus conditions: when performing an ANOVA-1-way test on the distributions of centrality in the three regions, we found that they had not the same mean ( in pre-stimulus, in post-stimulus). Furthermore, the Figure 5C shows that the centrality in CM decreased ( , ranksum test), while it increased for VPL populations ( , ranksum test) and decreased for S1 neurons ( , ranksum test). This suggested that neural populations in VPL received greater network load in stimulus processing.

(A) Synchrony of LFP phases in each recording site during spontaneous, spike responsive and spike non-responsive configurations. Responsive stimuli increase the global and decrease the local synchronies. (B) Average node degree distributions of functional graphs (LFPs) extracted by pre- and post-stimulus conditions. (C) Average betweenness centrality balance over the three recorded regions (VPL, S1, CM) in both conditions.

Subsequently, we compared pre- and post-stimulus functional graphs by computing the networks statistics in 50 and 100 ms windows before and after the onset of tactile stimuli. A time of 50 ms covers most of fast neural oscillations while a time of 100 ms includes slower ones. We found that both conditions (pre- and post-stimulus) exhibited small-world properties and the statistics were not different ( Table 4 , in all comparisons, ranksum test). In addition, the window size per se did not change the small-worldness values. Interestingly, the tactile stimulus onset triggers substantial increases in the number of inter-area connections but does not have a significant effect on the intra-area connections ( Figure 5A ). We concluded that small-world statistics are preserved in LFPs during both spontaneous and evoked activities. Thereafter, we analyzed the extracted graphs in pre- and post-stimulus conditions by their node degree distributions that were not significantly different ( Figure 5B , , ranksum test).

We asked whether these observations on a single recording could generalize to our full dataset. In order to test for this possibility we first divided our LFP recordings into two classes: responsive and non-responsive. LFP recordings were considered responsive when the average evoked firing rate, measured from the same recording channels (see Materials and Methods ), was larger than the mean plus 5 times the standard deviation of the basal firing rate. We found that the number of intra-area (local) functional connections was preserved across the three conditions of spontaneous, non-responsive and responsive LFP activity ( for CM, VPL and S1, ranksum test). However, during responsive LFP recordings, the average number of inter-area (or global) functional connections was substantially larger for all possible inter-area combinations ( for CM-VPL, CM-S1, VPL-S1, ranksum test).

Functional connections are disposed on the LFP recording sites. Green nodes represent CM channels, red nodes represent VPL channels and blue nodes represent S1 channels. (A) Before a tactile stimulation, LFPs are tightly coupled among the LFP channel of the same brain areas. (B) After an effective tactile stimulation, LFPs broke their inter-site synchronies and established cross-site phase couplings.

For evoked activity recordings, stimulus occurrence triggered significant topological reconfigurations in functional graphs, as shown in Figure 4 (i.e. a representative case). The picture reports a functional connection graph estimated before ( Figure 4A ) and after ( Figure 4B ) the stimulus onset. During spontaneous activity, CM (green), VPL (blue) and S1 (red) showed tight intra-area and poor inter-area synchronizations. After the stimulus onset inter-area synchronizations emerged while intra-area synchronizations were maintained roughly constant.

For spontaneous activity recordings, LFP functional connection graphs exhibited time-invariant small-world properties ( Table 3 ). Even in this analysis, increasingly larger windows makes smaller values thus the threshold was decreased when windows became larger. Time windows larger than 1 second showed no small-world network configurations.

In order to estimate the functional connections between pairs of channels we used the phase synchrony measure (see Materials and Methods ) [34] . Phase synchrony is more suitable than NCS for continuous signals and it has been widely applied to EEG and LFP analyses [35] . It is normalized in the range , where unity occurs when phase coupling is exact and zero indicates that it is null.

Our second aim was to estimate the topology of functional connection graphs obtained from the LFP activity recorded simultaneously to the spiking activity. We estimated the LFP functional connection maps by following the same sequence of analyses used for spiking activity. First, we estimated the functional connections between all possible electrode pairs in order to generate the adjacency matrix, then we binarized this matrix by using a variable threshold (see Materials and Methods ) and finally we computed the network statistics.

By analyzing the the modularity index (Q) of the community structures, we found that graphs in pre-stimulus conditions had larger Q indices ( , ) than graphs in post-stimulus conditions ( , ; , ranksum test). Hence the incoming stimulus information forced the functional networks to reduce the modularity indicating that more neurons were involved in the stimulus representation.

Subsequently, we analyzed the extracted graphs in pre- and post-stimulus conditions by their node degree distributions which were not significantly different ( Figure 3B , , ranksum test). Moreover, by analyzing the distribution of node betweenness centrality, we found that graphs in pre-stimulus conditions (mean, , standard deviation, ) had smaller betweenness centrality ( , ranksum test) than graphs in post-stimulus conditions ( , ) indicating that stimuli induced greater centralization in nodes. Moreover, we found that the betweenness centrality was not equally distributed over the three regions (VPL, S1, CM) both in pre- and in post-stimulus conditions. By means of an ANOVA-1-way test over the distributions of centrality in the three regions, we found that means were significantly different ( in pre-stimulus, in post-stimulus). Furthermore, Figure 3C shows that the betweenness centrality of neurons from CM increased ( , ranksum test), from VPL it decreased ( , ranksum test) and from S1 it increased ( , ranksum test). Since the betweenness centrality generally refers to the network node load [33] , the last result suggested that neurons from S1 received more load in the processing of stimulus information.

(A) Correlation between the standard deviation of difference matrices versus the spike responsiveness in each stimulus session. Correlations were computed by least-square regression (red lines, R = 0.637). (B) Average node degree distributions of functional graphs (spikes) extracted by pre- and post-stimulus conditions. (C) Average betweenness centrality balance over the three recorded regions (VPL, S1, CM) in both conditions.

Functional weights are redistributed from the pre- (A) to post-stimulus (B) configurations. Red, blue and green nodes indicate neurons respectively from VPL, S1 and CM. Some neurons, not functionally connected in the pre-stimulus graph, are involved in the stimulus processing [3] – [6] , [10] , [24] , [28] , [32] , [34] , [37] , [42] . Conversely, some neurons, previously employed, are excluded in the functional graph [11] . Furthermore, many neurons change their functional roles. For instance, cortical neuron 32 is a hub node (node with high degree) with many functional connections with VPL and CM thalamic neurons in (B) while is not involved in (A). Again, cortical neuron 23 constitutes a small clique with cortical neurons 22 and 40 in (A) while in (B) becomes a small hub node with diverse thalamic and cortical neurons.

Then we compared the average functional connections in time windows of 100 ms duration before and after the stimulus onset. 100 ms was a reasonable time constraint in order to consider most of the thalamocortical interactions. The functional connection graphs underwent significant topological reconfigurations after the tactile stimulus delivery. As an example, in Figure 2 , the neuron 32 switches from a fully disconnected state ( Figure 2A , relative to pre-stimulus condition) to highly connected (hub) state ( Figure 2B ) after the stimulus onset. Surprisingly, notwithstanding these profound changes, the small-worldness computed on the pre- and post-stimulus functional connection graphs were not significantly different ( in all comparisons, non-parametric Wilcoxon ranksum test, Table 2 ). Therefore, the average small-world properties were not perturbed by tactile stimulations. Furthermore, we performed a correlation analysis between the quantitative changes in functional connection graphs and the firing rate variations induced by tactile stimuli. We found that the changes induced by stimuli were significantly greater ( , ranksum test) than those observed in spontaneous activity. This means that stimuli modulated in an effective manner the functional connections between neurons. Namely, the firing rate of neurons was positively correlated with the evoked connection changes ( , , t-test; Figure 3A ) and whenever the recorded neurons were effectively involved in the response to stimuli, these stimuli induced concurrent functional connection changes proportional to the evoked firing rates.

Most of neural activity is supposed to be combined in sub-second time scales [32] . In fact, as expected, for time windows larger than 1 second, functional graphs did not meet the admissible criteria (see Materials and Methods ) and, anyhow, no small-world network configuration was found. This was probably due to the increasingly weaker functional interactions among neurons with such relaxed time delays. Notably, increasing the window size makes the neural interactions harder to detect and thus values became smaller because statistical dependencies become more and more sporadic. Consequently, the selected threshold was smaller when windows became larger.

First, we measured the small-worldness statistics in our spontaneous activity recordings and we found that S and stated that such graphs can be considered small-world networks ( , close to ), irrespective of the time window used for detecting the functional connections (50, 250, 500, 1000 ms; Table 1 ).

Subsequently, the small-world statistics (clustering coefficient), (characteristic path length) [6] , and were estimated on these graphs [18] , [31] . The first measures the tendency of neurons to segregate into separate sub-networks, the second measures the average path length between nodes, the third and the fourth are two measures of small-worldness. The terms and represent normalizing values estimated by algorithms which randomize the original networks but preserving the node degree distributions. Finally, we further characterized the resulting graphs by analyzing the node degree distribution, the betweenness centrality and the community structure. Specifically, we wondered if stimuli provoked changes in the node degree or the node centrality distribution or if they affected the community structures.

(A) A recording session from thalamic and cortical regions. Arrows indicate the effective influence among neurons. The electrode tips record the neurons in dark red. (B) The firing patterns of the cortical neuron B produce common firing patterns both with neurons A and C but with different time delays. In particular, (red spikes) can be easily inferred by correlation analysis instead of (blue spikes) hardly detectable. (C) Recorded signals are processed in overlapping windows lasting hundreds of milliseconds. (D) Spike trains are modeled by VMMs and compressed by LCAs. The functional connectivity strength between the spike trains A and B is estimated by the length of the compressed spike trains (C(A), C(B)) used by the function. Whether is greater than a fixed threshold then we can conclude that . (E) An example of functional graph extracted by recordings in the time window. (F) Typical sites of the rat paw for the tactile stimulation.

We estimated the functional connections in neuronal spiking activity by using the Normalized Compression Similarity ( ) function to detect functional couplings between pair of neurons. is defined in the range , where indicates no interaction and indicates an exact correspondence between the firing patterns of the two neurons considered. This measure was chosen in place of more conventional ones because of its ability to capture both short-range (synchronous) and long-range (delayed) interactions between neurons. Its efficacy was then assessed by analyses on synthetic spike train (see Materials and Methods , Figure 1C–E ). We thus scanned the recorded activity in search for functional relations by using sliding windows of different lengths (50, 250, 500, 1000 ms). The time length of a sliding window defined the largest possible delay at which an interaction could be detected (see Materials and Methods ). After the estimation on all possible pairs of simultaneously recorded neurons, the resulting adjacency matrices were binarized by a threshold in order to obtain the functional connection graphs. Typically, thresholds were chosen in order to select the strongest connections [2] , [30] and, in this study, we used the percentile of the weight distribution. In few cases, to meet the admissibility criteria we chose lower percentiles (see Materials and Methods ).

Our first aim was to estimate the topology of functional connection graphs obtained from simultaneous spiking activities of neuronal populations in VPL, CM and S1. Because specific tactile stimulation can potentially exert a powerful effect on functional connections, we set out to evaluate functional connection graphs both during spontaneous and evoked activities. In order to compute the graph's connections in the two conditions, we performed pairs of 10 minutes recordings, the first with no stimuli, the second by delivering fast tactile pulses (see Materials and Methods ) at randomized intervals (500 200 ms) on the five digits.

In an evolutionary perspective, small-world topologies appear to be preferentially selected among network topologies under the natural constraints (efficiency and reliability) of brain expanding complexity in the mammalian phylogenetic tree [12] . They can satisfy the necessity of internal input integration and grant for best responses to environmental requirements. Moreover, small-world networks seem coherent with the recent advancements in the physiology of neuronal networks. More specifically, from a functional point of view, small-world networks appear to provide dynamical features, such as communication-through-coherence [5] , , for fast information integration where even far located nodes participate to the information process in an efficient way.

Ultimately, inferences from computational models are often risky. Indeed, our findings mainly endorse but do not prove the hypotheses. Other topologies may deliver consistent results, rather than random or ring lattice networks and a set of other debated topologies (hierarchical modular networks, scale-free networks, etc.) could be also considered.

We proved the effectiveness of the NCS measure to capture the long-range dependencies and although such method has been implemented for related approaches 50 – 53 , it has never been applied before on electrophysiological data. A potential drawback of such approach is represented by its computational cost. Indeed, a numerical increase of recorded neurons (up to thousands or more), would need a cubic increase of the computational time required to calculate the NCS adjacency matrices.

The use of gaseous anesthetics represents a potential limitation of this study. The level of Isoflurane was indeed very low and low enough to avoid important suppressive actions of the neural activity, as acknowledged from other studies [49] . It remains therefore implicit that conclusive studies with awake animals are needed to definitely address the existence of functional topologies (at least) in these three brain regions. Furthermore, the analyzed functional graphs could be compared only to a limited extent, since their basilar features (number of nodes, edges, etc.) were different.

Our results are also strictly related to research suggesting that neuronal oscillations enable selective and dynamic control of distributed functional cell assemblies [5] . We could add to such a scenario, the speculation that LFP coherent activities may reflect the integration of local computations which occurred in these distributed cell assemblies. Thus, the small-world topology expressed by synchronized and distributed LFP phases would support a hierarchical and efficient functional substrate for incorporating cell assemblies.

Our synthetic model tried to meet a double hypothesis: on the one hand that the topology expressed by our small graphs could be a nested and natural expression of a homologous larger population and on the other hand to answer the reverse question concerning the hypothetical extended topology that we designed in our synthetic model and if it could suitably host a subset of units with the described topological properties of our recorded networks. Namely, we argue that there are mutually fitting topological features between our recorded and the hypothetical larger synthetic networks.

A potential objection may be that the number of shared neurons appears to be relatively high. However, it should be noted that neurons in stimulated conditions probably get a balance between the number of assemblies and the inherent metabolic cost [48] . As a consequence, neuron sharing could cut down the metabolic cost of network activations, in accordance with a parsimonious limiting factor to an unregulated growth of recruited functional units.

Our model assumed that neurons may variably join in concurrent assemblies and we empirically estimated that given a population of neurons, the expected number of cell assemblies ranges within .

The analyses from computational models suggest that a cell assembly, i.e. a transient functional unit composed by neurons potentially distributed in separated brain regions [4] , appears to be organized as a small-world network. The null hypotheses that such groups were random or lattice networks does not match experimental evidences.

Along with these changes, however, the topology of the small-world networks never degenerated. This dualistic behavior, redistribution of connections and permanence of topology, appears an interesting interpretive key, which can potentially be extended to other complex brain networks. These features may also be useful to analyze also pathological signs in brain dynamics, like epilepsies or schizophrenia, where the detection of disrupted network or loss of topological hallmarks may help novel nosological classifications.

In this study our aim was to explore the functional connections of different brain areas involved in somato-sensory information processing and to envisage a potential broadening to other brain circuits. The precise mechanisms regulating the rich and complex neural thread of brain areas involved in the construct of tactile processing are still far from being clarified. However, it is possible to claim that a response to sensory inputs is well recognized by a number of areas, including the areas we chose for this study. Namely, several studies showed that mechanoreceptor signals reach both VPL and CM thalamus. The former innervates directly S1, while the latter outputs are more diffused and addressed to higher sensorimotor cortical layers [25] , [26] , [45] – [47] . We observed that in the resting state, VPL, CM and S1 showed substantial mutual connections and, in addition that, under stimuli, they establish an even more integrated architecture enabling fast information exchanges. This functional connectivity could be ascribed to the necessity to create larger functional units and nodes with higher network load. Furthermore, S1 neurons seem to be preferred in stimulus processing indicating that this area hosts the most of information processing among the regions explored.

Yet, Bassett et al. [11] , suggested that the functional organization in small-word networks may be ubiquitous at different spatial and temporal scales and in different experimental conditions. This organization could represent a stable, and even a necessary, scheme for information processing in brain areas, networks and neurons.

Studies in anatomical brain connectivity discovered that small-world network architectures are a distinctive trait in animals [19] , [21] , [40] , including humans [10] , [13] , [41] , yet with different degree of brain development. The complementary observation that brain pathologies like epilepsy [42] and schizophrenia [43] show small-world network disruptions provides further support to the potential interpretive and explanatory strength of small-word topology. It should be noted that the methodological approaches used to study the functional connectivity at macroscopic and microscopic levels are slightly different. Studies with fMRI estimate functional interactions among brain areas comparing the BOLD signal of each area to a seed region. Scientists usually refer to these effects as functional connectivity [2] , [10] , [13] . In this work, instead, we studied networks built from the extracted functional connections computed comparing the electrical activity of all possible neuron couples [44] . Despite the technical incongruences and focusing on neuronal functional connections and the underlying topology, a number of meaningful studies found an intrinsic small-world topology in single visual areas of primates [15] , [17] , [19] , [20] , [22] , [23] and in neuronal cultures [16] . These elegant approaches gave access to the richness of the local functional connection graphs scaled at the neuronal level. However, none of these studies provided a description of an extended connectivity, involving different neuronal populations gathered in close proximity or located in distant brain regions.

In this paper we show that spiking and LFP activities of neurons, in three stations of the somatosensory system of rat brains, present clear signs of small-world functional organization with sub-second invariance. Furthermore, we show that this distinctive functional organization persists in the presence of tactile stimuli, independently of the neural response intensity. Finally, results obtained by means of computational models suggest that small-world networks may represent a consistent and formal model for cell assemblies.

Materials and Methods

Ethical Notes 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). 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). As such, no protocol-specific approval was required. The approval of the Ministry of Health is classified as “Biella 1, 3/2011” into the files of the Ethical Committee of the University of Milan.

Animal Preparation and Stereotaxis Seven male albino rats (Sprague-Dawley, Charles River, Calco, LC, Italy, g) were chosen among the set of 11 animals employed in the recordings. All the animals were maintained in a hour light-dark cycle with access to food and water ad libitum. The rats underwent preliminary barbiturate anesthesia for the surgical experimental preparation. The jugular vein and the trachea were cannulated to gain direct drug delivery access and a connection to the anesthesia-ventilation device. Before the placement of electrodes, rats were paralyzed by intravenous Gallamine thriethiodide ( mg/kg/h) injection and connected to the respiratory device delivering (1 stroke/s) an Isoflurane ( to l/min) and Oxygen ( l/min) gaseous mixture [54]. Curarization was maintained stable throughout the whole experiment by Gallamine refracted injections ( ml of the original solution/h). The anesthesia levels were maintained into ranges which prevented any corneal or retraction reflex (in absence of curarization) with low intensity noxious mechanical stimuli applied on a posterior paw [54]. We chose three areas for the simultaneous neuronal recordings in the left brain: the thalamic median nuclei (in particular the centromedial (CM)), the thalamic ventro-postero-lateral nuclei (VPL) [25], [47] and the primary somatosensory (S1) cortex. Fast tactile stimuli were delivered to the right posterior paw (see Figure 1F) by a suitable electromechanical device. Two holes were drilled on the skull. A mm bone window for the access of the cortical matrix electrodes and a larger bone window ( mm) allowing for the simultaneous insertion of two parallel electrode matrices directed to the thalamic nuclei. The cortical access was set around a reference at mm AP and mm ML on the left [55], and the electrode matrix was driven around to micrometer deep by an electronically controlled microstepper (Narashige, Japan). The thalamic access was centered at the two focus points of mm AP, and mm ML [55]. The electrodes were inserted with a slant and driven at least to m in depth and then advanced by a second electronically driven microstepper (AB Transvertex, Stockholm) until responses were observed to peripheral test stimuli. The neuronal recordings were obtained with two types of matrices, a vertical array devoted to the cortical recordings and two planar matrices devoted to the thalamic recordings. The vertical array was a multitrode (Multitrode Type 1, Thomas RECORDING GmbH, Giessen, Germany) with gold contacts ( m contact spacing) with an average impedance of M . For planar matrices were frames of tungsten or Pt-Ir electrodes, inter-tip distance m, tip impedance M (FHC Inc., ME, USA). Fast thalamic and cortical responses to light tactile stimuli in the plantar aspect of the right hind limb were used as anatomo-functional acceptance criterion for acquisition.

Tactile Stimulation Controlled stimulation was delivered through a blunted cactus thorn on each of five sites of the rat right hind limb (Figure 1F). The tip was mounted on the dust cap of a speaker and driven through an Arduino microcontroller board (available at http://www.arduino.cc). At the beginning of each stimulation epoch the tip was lightly placed over the skin. Fast ms pressure pulses were applied following a semi-random sequence. Pulses occurred in couplets. The delay between the first pulses of each couplet was set at ms. Every second pulse of each couplet followed the first by a random delay extracted uniformly in the range ms (see Figure 6C). The stimuli semi-randomness was adopted to avoid habituation [56]. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 6. Basilar neurophysiological data. (A) Distributions of firing rate and (B) of Inter-spike Interval for the representative neurons. (C) Pattern of tactile stimulations by Arduino microcontroller. https://doi.org/10.1371/journal.pcbi.1003104.g006

Neurophysiological Recordings and Preliminary Data Analyses For signal amplification and data recordings a channel Cheetah Data Acquisition Hardware was used (Neuralynx, MT, USA, sampling frequency kHz). Electrophysiological signals were digitized and recorded with bandpasses at kHz/ Hz for spikes, 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 Hz with standard deviation of Hz. After the recordings the LFPs were downsampled to KHz. We used for filtering the same techniques described in [57]. After filtering and downsampling, the spike contamination of LFP signals was null avoiding further spike removal techniques [34]. The spikes were extracted and sorted by using the Wave_clus MATLAB toolbox [58]. Sorted cells with average rates below Hz and above 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 neurons ( in each experiment) out from the set of the acquired signals. Distribution of firing rates and inter-spike intervals is shown in the Figures 6A–B. 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 (from to ms) overlapping windows (Figure 1C), thus obtaining an ordered set of equal length windows.

Functional Connections by Spike-train Similarities Interactions between neurons can generate very complex, time-delayed and asymmetric patterns. This represents a potential problem in experimental configurations where the communications between distant neurons are taken into consideration. Standard techniques like correlation analysis are, in many cases, unable to detect such events. Indeed dependencies between neurons from thalamus and cortex can last tens of milliseconds [56], [59]–[61]. In a toy example (shown in Figure 1A–B) the existence of the interaction between neurons A and B is assumed (A B). The neuron A is recorded and its activation triggers the long chain that, from site A, produces firing activity to neuron B in another brain region. The direct path between neurons C and B allows, instead, for synchronous spike patterns well detectable by correlation analysis (Figure 1B). In general, in simultaneous recordings from separated brain sites, it is unlikely to find couples of physically wired neurons although axonal projections connect the sites. It's definitely more probable that spiking activity relations may reflect a complex anatomical substrate, where chains of activations exist between them provoking the observed coherent activity. Such a problem can be solved by mathematical tools able to model arbitrarily long temporal relationships. In this work, we proposed a novel framework wherein spike trains with arbitrarily long temporal dependencies are modeled by Markov stochastic models. Typically, in regular Markov models, each state depends only on the previous state while higher-order Markov models suffer from high state-space computational complexity [62]. Here, we used the Variable-order Markov Models (VMMs) [63], [64] because they are able to overcome these limitations. Lossless compression algorithms (LCAs) represent one of the most efficient techniques to estimate VMMs [64]. Within the set of LCAs we chose the Prediction by Partial Matching (PPM) algorithm [65], [66] which is considered the best match between prediction accuracy and speed [64]. The last step consists to build a similarity function for this kind of spike train stochastic models. In the last 15 years, Vitanyi and colleagues have developed a function, the Normalized Compression Distance ( ) [67]–[70], which estimates the distance between symbolic sequences. This function performs the estimation directly by the sizes of the compressed sequences. In fact, can be proved that the better is the VMM estimation the shorter is the compressed sequence size. In this work, we redefined the function in order to reverse its assigned values pointing to similarity instead of distance. We called such function the Normalized Compression Similarity ( ). Formally, given that and are two neural sequences (e.g. spike trains), the is defined as follows: (1)where the function represents the compressed sequence length and is the sequence concatenation operator (e.g. ). If is close to , the sequences and are considered similar. If close to , the sequences are strongly dissimilar. We therefore evaluated the function on time windows ( ms) of the recorded (binary) spiking activity (Figure 1C–D) assuming that relative high values of similarity corresponded to actual functional connections (Figure 1E). To note that significant values do not imply significant correlations but the opposite is true. Although the is asymmetric function, it is not relevant for the causal interaction analysis and consequently for effective connectivity. Because has never been used as method to estimate the functional connections in neurophysiological data (although many LCAs has been used as information measure of synaptic efficacy [50] and spike train similarities [51]–[53]), our aim was to prove that effectively captures similar asynchronous spike patterns in comparison to standard correlation analyses. Furthermore, we assessed that did not bias its estimations with respect to independent spike trains generated by simulations. To address these requirements we performed two experiments: i) we evaluated the and the Pearson correlation coefficient over a set of independent binary spike trains and ii) we evaluated the , the correlation coefficient and the time lagged cross-correlation with couples of binary spike trains ( ms long) containing an identical, short and random spike pattern ( ms) that is fixed and centered in the first train and drifts (from left to right) in the second trains (see Figure 7B, overall drifts). The latter procedure creates synthetic spike trains wherein the common pattern is increasingly distant. Random spikes were added into each spike sequence to establish the efficacy of each technique in presence of noise. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 7. Efficacy of Normalized Compression Similarity (NCS) to detect long-range spike interactions. (A) NCS, Pearson coefficient exerted on 100 couples of independent uniformly distributed binary sequences (1000 bits). Both functions do not show bias. (B) To test the capacity of NCS to detect significant interactions within time windows of 1000 ms, we fixed in the center of the first time-window a binary random pattern (100 ms long, first plot above). The same pattern was replicated in a sequence of 47 time-windows drifting it from the initial to the window end. NCS, Pearson coefficient and the time lagged cross-correlation (the maximum value) were evaluated on these sequences: NCS is able to detect the interaction along the entire drifting process while Pearson coefficient is able to return significance only when the reference pattern is almost aligned in both sequences (drifts number 23–27). The time lagged cross-correlation detects weak significances in the data because of the spikes that not belong to the pattern. https://doi.org/10.1371/journal.pcbi.1003104.g007 In the first experiment we found that the NCS method did not show any significance for independent spike trains as well as for the Pearson coefficient ( , for Pearson coefficient and , for , Figure 7A). This guarantees that the function is unbiased to independent distributed binary spike trains. In the second experiment, the proved its efficacy to detect long-range interactions where Pearson coefficient and time lagged cross-correlation fails. In fact, during the whole shifting epochs, significance held while Pearson coefficient showed significant only when the two spike patterns were very spatially close (Figure 7C, drift number 23–27). Instead, the time lagged cross-correlation was able to detect weak significances in the data but the spikes that do not belong to the pattern, made the values definitely lows.

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 ( and ) by the following function (2)where is Napier's constant, is Hilbert Transform, is the argument function and is the imaginary unit. The Hilbert transform and the argument were computed with, respectively, the hilbert and the angle Matlab functions [34], [35]. When is equal to ( ), then and 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 equal-length time windows (Figure 1) and then we computed the adjacency matrix for all neurons or electrodes. The resulting matrices exhibited values in the unitary interval. We repeated the analyses with different window sizes from ms to s (fixing the maximum dependency order to half of the window size). The functional connections extracted from extracellular recordings can be represented by graphs. A variable threshold (typically equal to a higher percentile of the weight distribution and vary between the and ) selected the strongest connections, thus allowing for the construction of the functional connection graphs. The choice of a threshold is related to the window length used to estimate the connection strength. Larger windows produce fewer connections and to keep the network sparsity quite constant, smaller thresholds must be chosen. 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. A small-world network is generally obtained by evolving a basic ring lattice graph, where each node is connected to their neighbors. The chosen neighborhood involves typically much less nodes than the total node number ( ). The graph evolution is achieved by randomly adding and removing edges from the starting graph [6], [7]. 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 nodes. From a functional perspective, small-world networks can express two important information processing features: information integration and segregation [13], [71]. 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 ( ) and the characteristic path length ( ) [6], [30]. The former measures how close the neighbors 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 [30], were used for our network analyses (clustering_coef_bu.m, charpath.m). 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) [30]. 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 and ( by latticization, latmio_und.m). These null network values are required to verify the small-world nature of the graphs. In fact, classical and novel measures of small-worldness such as [18] and [31] state that, respectively, if or is close to , the graph can be considered a small-world network. 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 [72], 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 [30], the distribution of node centrality highlights how the information flow is balanced within graphs. For the same purpose we further studied the community structure of our graphs. Communities emerge from graphs by applying a clustering algorithm to nodes [73]. In this work, communities represent the aggregated functional units under investigation and by analyzing them we can understand how node graphs are aggregated in each experimental condition. 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 [44], [74]. 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 and could not 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 and percentile) of the node, edge and density number distributions in order to obtain a better homogeneity. In the work, we refer to these two conditions as admissibility criteria.

Visualization of Graphs Graphs are visualized by using the function LayeredGraphPlot of Mathematica (Wolfram Research, Inc., Mathematica, Version 8.0, Champaign, IL). This function implements the Sugiyama algorithm [75]. A Layer layout helps to understand the information flow in the network and the specific roles performed by each node. Graphs in the insets are visualized by using the function GraphPlot of Mathematica that implements the spring-electrical algorithm [76] that is typically used to draw small-world networks.