Neuronal Functional Connection Graphs among Multiple Areas of the Rat Somatosensory System during Spontaneous and Evoked Activities

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.


Introduction
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.
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][8][9][10][11][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][16][17][18][19][20][21][22][23]. Specifically it is not clear whether functional groups of neurons from multiple brain regions display a SWN organization.
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][25][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.
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.
Consistent with other findings regarding brain networks at diverse scales [5,[27][28][29], our results provide evidence for a model of functional organization where distinct functional neuronal assemblies are sparse and encompass multiple brain areas.

Functional Connection Graphs in Spiking Activity
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.
We estimated the functional connections in neuronal spiking activity by using the Normalized Compression Similarity (NCS) function to detect functional couplings between pair of neurons. NCS is defined in the range ½0,1, where 0 indicates no interaction and 1 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 NCS 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 75 th percentile of the weight distribution. In few cases, to meet the admissibility criteria we chose lower percentiles (see Materials and Methods).
Subsequently, the small-world statistics C (clustering coefficient), L (characteristic path length) [6], S and v 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 C r and L r 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.
First, we measured the small-worldness statistics in our spontaneous activity recordings and we found that S and v stated that such graphs can be considered small-world networks (Sw1, v close to 0), irrespective of the time window used for detecting the functional connections (50, 250, 500, 1000 ms; Table 1).
Most of neural activity is supposed to be combined in subsecond 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 NCS values became smaller because statistical dependencies become more and more sporadic. Consequently, the selected threshold was smaller when windows became larger.

Author Summary
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.
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 prestimulus condition) to highly connected (hub) state ( Figure 2B) after the stimulus onset. Surprisingly, notwithstanding these profound changes, the small-worldness computed on the preand post-stimulus functional connection graphs were not significantly different (Pw0:13 in all comparisons, non-parametric Wilcoxon ranksum test, Table 2). Therefore, the average smallworld 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 (Pv0:005, 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 (R~0:67, Pv0:001, 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.
Subsequently, we analyzed the extracted graphs in pre-and post-stimulus conditions by their node degree distributions which were not significantly different ( Figure 3B, P~0:48, ranksum test). Moreover, by analyzing the distribution of node betweenness centrality, we found that graphs in pre-stimulus conditions (mean, m~22:84, standard deviation, s~1:63) had smaller betweenness centrality (Pv0:000, ranksum test) than graphs in post-stimulus conditions (m~25:67, s~2:01) 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 (P~0:007 in pre-stimulus, P~0:002 in post-stimulus). Furthermore, Figure 3C shows that the betweenness centrality of neurons from CM increased (P~0:013, ranksum test), from VPL it decreased (P~0:000, ranksum test) and from S1 it increased (P~0:006, 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.
By analyzing the the modularity index (Q) of the community structures, we found that graphs in pre-stimulus conditions had larger Q indices (m~0:48, s~0:11) than graphs in post-stimulus conditions (m~0:24, s~0:01; P~2:04 : 10 {64 , 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.

Functional Connection Graphs in LFPs
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. In order to estimate the functional connections between pairs of channels we used the phase synchrony measure c (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 ½0,1, where unity occurs when phase coupling is exact and zero indicates that it is null.
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 c values thus the threshold was decreased when windows became larger. Time windows larger than 1 second showed no small-world network configurations.
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.
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 (P~0:63 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 (Pv0:002 for CM-VPL, CM-S1, VPL-S1, ranksum test).
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, Pw0:38 in all comparisons, ranksum test). In addition, the window size per se did not change the smallworldness 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, P~0:83, ranksum test).
By analyzing the distribution of node betweenness centrality, we found that graphs in pre-stimulus conditions (m~21:46, s~2:17) had smaller betweenness centrality (Pv0:03, ranksum test) than graphs in post-stimulus conditions (m~23:79, s~1:63) 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 poststimulus 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 (P~0:032 in pre-stimulus, P~0:039 in post-stimulus). Furthermore, the Figure 5C shows that the centrality in CM decreased (P~0:015, ranksum test), while it increased for VPL populations (P~0:033, ranksum test) and decreased for S1 neurons (P~0:013, ranksum test). This suggested that neural populations in VPL received greater network load in stimulus processing.
Moreover, by analyzing the modularity index of the community structures, we found that graphs in pre-stimulus conditions had larger Q indices (m~0:60, s~0:13) than graphs in post-stimulus conditions (m~0:28, s~0:02; P~4:38 : 10 {67 , ranksum test). In conclusion, the incoming stimulus forced delegated networks to merge the active communities indicating that more neural groups expressed coordinated oscillations.

Functional Connection Graphs in Simulated Spiking Activity
In the previous sections we provided evidence, both at a preand 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 Abbreviations: win indicates the window size; h represents the chosen threshold used in matrix binarization; C is the clustering coefficient of the extracted functional graphs; C l is the clustering coefficient computed on the latticizied version of the extracted graphs; C r is the clustering coefficient computed on the randomized version of the extracted graphs; L is the characteristic path length of the extracted graphs; L is the characteristic path length computed on the randomized version of the extracted graphs; S is a small-worldness index equal to C=C r L=L r ; v is a small-worldness index equal to Functional Connections among Multiple Brain Areas PLOS Computational Biology | www.ploscompbiol.org 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.
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. For these purposes, we implemented a large-scale simulation of N neurons that embedded k small-world networks. Then we estimated small-worldness by systematically varying k. We simulated a neural network of N~100000 neurons, a biological relevant neuronal ensembles exhibiting a magnitude about five fold the estimated size for a rat cortical hypercolumn [38]. The 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][4][5][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. doi:10.1371/journal.pcbi.1003104.g002 size of each functional group was chosen by sampling from a discrete uniform distribution *U(50,N=100). 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.
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 k~0:01, peaking when k~0:1 and decreasing with k up to 0:5. At peak we found comparable statistics values with those measured experimentally both on spiking (P~0:60, ranksum test) and LFP activities (P~0:46, 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 p~1. In both cases, we did not find consistency with the experimental observations.
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 0:01 and 0:1, keeping consistency with the S and v values observed in experiments.

Discussion
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, indepen-dently 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.

Small-world in Brain Networks
Studies in anatomical brain connectivity discovered that smallworld 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 smallworld 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.
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  Abbreviations: win indicates the window size; h represents the chosen threshold used in matrix binarization; C is the clustering coefficient of the extracted functional graphs; C l is the clustering coefficient computed on the latticizied version of the extracted graphs; C r is the clustering coefficient computed on the randomized version of the extracted graphs; L is the characteristic path length of the extracted graphs; L is the characteristic path length computed on the randomized version of the extracted graphs; S is a small-worldness index equal to C=C r L=L r ; v is a small-worldness index equal to L r L { C C l . doi:10.1371/journal.pcbi.1003104.t002 organization could represent a stable, and even a necessary, scheme for information processing in brain areas, networks and neurons.

Tactile Information Processing
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][46][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.
Along with these changes, however, the topology of the smallworld 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.

Cell Assemblies and the Synthetic Model
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.
Our model assumed that neurons may variably join in concurrent assemblies and we empirically estimated that given a population of n neurons, the expected number of cell assemblies ranges within ½n=100,n=10.
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 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.
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.

Limitations and Developments
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. Abbreviations: win indicates the window size; h represents the chosen threshold used in matrix binarization; C is the clustering coefficient of the extracted functional graphs; C l is the clustering coefficient computed on the latticizied version of the extracted graphs; C r is the clustering coefficient computed on the randomized version of the extracted graphs; L is the characteristic path length of the extracted graphs; L is the characteristic path length computed on the randomized version of the extracted graphs; S is a small-worldness index equal to C=C r L=L r ; v is a small-worldness index equal to 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][51][52][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.
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.

Conclusive Notes
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,27,29], for fast information integration where even far located nodes participate to the information process in an efficient way.

Animal Preparation and Stereotaxis
Seven male albino rats (Sprague-Dawley, Charles River, Calco, LC, Italy, 300{400 g) were chosen among the set of 11 animals employed in the recordings. All the animals were maintained in a 16=8 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 (20 mg/kg/h) injection and connected to the respiratory device delivering (1 stroke/s) an Isoflurane (2:5% 0:4 to 0:8 l/min) and Oxygen (0:15{0:2 l/min) gaseous mixture [54]. Curarization was maintained stable throughout the whole experiment by Gallamine refracted injections (0:1 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 3|2 mm bone window for the access of the cortical matrix electrodes and a larger bone window (6|2 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 {1:5 mm AP and {2:5 mm ML on the left [55], and the electrode matrix was driven around 450 to 800 micrometer deep by an electronically controlled microstepper (Narashige, Japan). The thalamic access was centered at the two focus points of {6 mm AP, {0:8 and {2:5 mm ML [55]. The electrodes were inserted with a 25 0 slant and driven at least to 5500mm 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 8 gold contacts (125mm contact spacing) with an average impedance of 1:2 MV. For planar matrices were 3|3 frames of tungsten or Pt-Ir electrodes, inter-tip distance 150{200mm, tip impedance 0:5{1 MV (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 anatomofunctional 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 5 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 500 ms. Every second pulse of each couplet followed the first by a random delay extracted uniformly in the range 150{250 ms (see Figure 6C). The stimuli semi-randomness was adopted to avoid habituation [56].

Neurophysiological Recordings and 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  Abbreviations: win indicates the window size; h represents the chosen threshold used in matrix binarization; C is the clustering coefficient of the extracted functional graphs; C l is the clustering coefficient computed on the latticizied version of the extracted graphs; C r is the clustering coefficient computed on the randomized version of the extracted graphs; L is the characteristic path length of the extracted graphs; L is the characteristic path length computed on the randomized version of the extracted graphs; S is a small-worldness index equal to C=C r L=L r ; v is a small-worldness index equal to L r L { C C l . doi:10.1371/journal.pcbi.1003104.t004 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 [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 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. 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 50 to 1000 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, timedelayed 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][60][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  Abbreviations: r indicates the density (small-world networks per node) and the winsize is set to 100 ms in all analyses; C is the clustering coefficient of the extracted functional graphs; C l is the clustering coefficient computed on the latticizied version of the extracted graphs; C r is the clustering coefficient computed on the randomized version of the extracted graphs; L is the characteristic path length of the extracted graphs; L is the characteristic path length computed on the randomized version of the extracted graphs; S is a small-worldness index equal to C=C r L=L r ; v is a small-worldness index equal to L r L { C C l . doi:10.1371/journal.pcbi.1003104.t005 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 (NCD) [67][68][69][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 NCD function in order to reverse its assigned values pointing to similarity instead of distance. We called such function the Normalized Compression Similarity (NCS). Formally, given that x and y are two neural sequences (e.g. spike trains), the NCS is defined as follows: NCS(x,y)~1{NCD(x,y)~C (x : y){minfC(x),C(y)g maxfC(x),C(y)g ð1Þ 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 therefore evaluated the NCS function on time windows (50{1000 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 NCS values do not imply significant correlations but the opposite is true. Although the NCS is asymmetric function, it is not relevant for the causal interaction analysis and consequently for effective connectivity. Because NCS 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][52][53]), our aim was to prove that NCS effectively captures similar asynchronous spike patterns in comparison to standard correlation analyses. Furthermore, we assessed that NCS 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 NCS and the Pearson correlation coefficient over a set of independent binary spike trains and ii) we evaluated the NCS, the correlation coefficient and the time lagged crosscorrelation with couples of binary spike trains (1000 ms long) containing an identical, short and random spike pattern (100 ms) that is fixed and centered in the first train and drifts (from left to right) in the second trains (see Figure 7B, overall 47 drifts). The latter procedure creates synthetic spike trains wherein the common pattern is increasingly distant. Random spikes were added into  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 (m~{1:6 : 10 {5 , s~0:005 for Pearson coefficient and m~0:013, s~0:004 for NCS, Figure 7A). This guarantees that the NCS function is unbiased to independent distributed binary spike trains. In the second experiment, the NCS 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][24][25][26][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 (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 unit. The Hilbert transform and the argument were computed with, respectively, the hilbert and the angle Matlab functions [34,35]. When c(x,y) is equal to 1 (0), then x and y are perfectly synchronous (asynchronous).

Complex Brain Network
By using the NCS and c 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 50 ms to 1 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 0:2 and 0:8) 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 K neighbors. The chosen neighborhood involves typically much less nodes than the total node number N (N&K&ln(N)&1). 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 (C) and the characteristic path length (L) [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 C r and L r (C l 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

S~C
=C r L=L r [18] and v~L r L { C C l [31] state that, respectively, if Sw1 or v is close to 0, 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 C and L could not be significantly affected by our network size changes. Significant changes appeared for synthetic 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][24][25][26][27]. The time lagged cross-correlation detects weak significances in the data because of the spikes that not belong to the pattern. doi:10.1371/journal.pcbi.1003104.g007 networks that increased their size by orders of magnitude. However, we discarded graphs that were outliers (beyond 5 th and 95 th 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 springelectrical algorithm [76] that is typically used to draw small-world networks.

Computational Model
So far, Hebb's idea of cell assemblies represents a fundamental theory supporting important physiological events. Guided by results from recordings, we hypothesized that a cell assembly can be functionally organized as a small-world network and we investigated such hypothesis by computer simulations. As originally thought, a cell assembly can be composed by groups of neurons (even anatomically dispersed) and the membership to each assembly can be dynamical [4].
Specifically, we brought back such considerations in a computational model to investigate two hypotheses: (i) can a network model based on previous assumptions can be consistent with the observed experimental facts? (ii) assuming the previous as true, how many assemblies (small-world networks) can exist over a set of neurons keeping consistency with observations?
Typical simulation frameworks require a choice of neuron models (Integrate and Fire, Izhikevich, Hodgkin and Huxley, etc.) and of a defined network layout (nodes and connections). These choices can be very crucial and become even more important if the aim of the study is the functional organization of the units. For this reason, we proceeded following an unconventional approach assuming that cell assemblies are effectively functionally organized as small-world networks on large-scale networks and, sampling the activity of a random subset of nodes, we wondered if such activity can elicit small-world organization as well.
So we first assumed that a set of small-world networks exists over a set of available nodes. If we consider each small-world network as a cell assembly [4], the small-world networks of our model can share their nodes.
As a whole these facts define the structural property of the model and can be summarized as follow: 1. Each neuron is represented by a node. 2. Functional connections between neurons are represented by edges. 3. Neurons are functionally organized as small-world networks. 4. Many small-world networks are embedded within the simulated network thus each node may belong to more than one smallworld topologies ( Figure 8A-D). 5. According to the typical amount of neocortical neurons in a microcolumn [39], the sizes of small-world networks are randomly sampled by a uniform distribution *U (50,N=100) where N is the number of nodes.
From a formal point of view, the network structure can be interpreted by a Multigraph G~(V ,E), the set V represents the nodes and the multiset E represents the unordered list of edges inherited by each small-world network.
To establish the node functioning we chose a primary criterion claiming that brain processing takes place by functional segregation and integration [71]. This concept, supported by many experimental evidences, explains how segregated specific brain areas work together to produce globally integrated behaviors, cognitions and percepts. We noted that small-world networks and integration-segregation criterion are closely related. The centrality of nodes can be estimated by several measures and we arbitrarily chose the betweenness centrality. Computations performed in (specialized) peripheral nodes, can be subsequently integrated in central nodes by virtue of the small average shortest path length of small-world networks. Indeed, it can be easily shown, by using simple computer simulations, that nodes in smallworld networks had a wide spectrum of centrality values. This had a heavy-tailed distribution in contrast to random networks distributed as a Gaussian-like distribution, i.e. every node had almost the same centrality.
In order to implements these criteria in our model, nodes fires following a rank-order dictated by the centrality values estimated by using the betweenness centrality.
As a whole, these last facts can be summarized as follow: 1. Each small-world network represents the processing of a specific information (the working hypothesis). 2. At most two randomly chosen small-world networks are allowed to run independently in a time window. A single neuron can work either alternately or concurrently within different networks of pre and post-synaptic nodes. This assumption claims that neurons are not exclusive and that they can partake simultaneously in many (up to 2) tasks. 3. Neurons express their spikes within 10 ms time windows.
Neurons fire following a centrality criterion based on the node betweenness centrality. A neuron with high betweenness fires its spikes after a lower betweenness neuron. This dynamic picture respects the information integration-segregation paradigm ( Figure 8E-H) [71]. The choice of the window length is arbitrarily and is proportional to the firing rate. 4. Simulation timesteps are set to 1 millisecond. Spike propagation times are uniformly distributed within the range ½1,3.
The algorithm governing this network is as follows: Input: n_swn, n nodes, timestep; Output: the binary n_nodes-by-timesteps matrix trains pr0.05; kr5; for ir1 to n_swn do nrsubsample(n_nodes); Grrandom_watts_strogatz(n, k, p); winsrcomputeWindows(G); foreach time windows in wins do if 0:01,rand() then centralityrBetweennessCentrality(G); foreach node A in centrality do trains[A, wins+centrality(A)]r1; foreach output node B of A do trains[A, wins+ centrality(A)+ centrality(B)]r1; end end end end end where the function subsample() selects a subset of nodes out of the set of available ones. The function random_watts_strogatz() returns a random graph built following the Watts-Strogatz algorithm, the function computeWindows() computes the length of the execution windows for the current graph and returns the list of such windows. The function rand() returns a uniformly generated random number between 0 and 1. At last, the function BetweennessCentrality() computes and returns the centrality of each node.
The model has been developed using the Python environment [77]. For the generation of small-word networks we used the networkx package (available at http://networkx.lanl.gov/). The Figure 8I shows a raster plot obtained by sampling the simulated activity. The darkest blue represents the no-spike event and the other spike colors are associated with the specific small-word topologies that generated them. The algorithm describes a core loop where each small-world network is first randomly created by the library routine watts_strogatz_graph (probability of rewiring equal to 0:05). Small-world networks take only a fraction of the total node number and several small-world networks can share subsets of their nodes. In a second stage, the generated small-world network expresses its spiking activity following the betweenness centrality as in point 2 of the dynamical assumption. Low level uniformly distributed noise is added to the spike propagation time.
The spikes of the current network occur randomly in equal size time windows (10 timesteps). The spike activity is saved and the loop restarts (source codes are available at http://code.google. com/p/swn-neuronal-networks/).  (16 nodes) showing three small-world networks (B-D) embedded into one network (A). The firing criterion is illustrated in graphs (E-H): each small-world network gives a specific (betweenness) centrality in its nodes. Nodes with low values of centrality (from periphery) spike first, nodes from center, later. (I) A sample of synthetic spiking activity from 100 nodes. Colors are specific for each small-world network. doi:10.1371/journal.pcbi.1003104.g008