Neuronal functional connectivity among multiple areas of the rat somatosensory system during spontaneous and evoked activities
Antonio G. Zippo, Riccardo Storchi, Giuliana Gelsomino, Sara Nencini, Gian Carlo Caramenti, Maurizio Valente, Gabriele E. M. Biella
11 Neuronal functional connectivity among multiple areas of therat somatosensory system during spontaneous and evokedactivities
Antonio G. Zippo , Riccardo Storchi , Giuliana Gelsomino , Sara Nencini , Gian Carlo Caramenti ,Maurizio Valente , Gabriele Eliseo M. Biella , ∗ , ∗ [email protected] Abstract
Small-World Networks (SWNs) represent a fundamental model for the comprehension of many complexman-made and biological networks. In the central nervous system, SWN models have been shown to fitwell both anatomical and functional maps at the macroscopic level. However the functional microscopiclevel, where the nodes of a network are composed of single neurons, is still poorly understood. At this level,although recent evidences suggest that functional connectivity maps exhibit small-world organization, itis not known whether and how these maps, potentially distributed in multiple brain regions, changeacross different conditions, such as spontaneous and stimulus-evoked activities.We addressed these questions by simultaneous multi-array extracellular recordings in three brain regionsdiversely involved in somatosensory information processing: the ventropostero-lateral thalamic nuclei(VPL), the primary somatosensory cortex (S1) and the centro-median thalamic nuclei (CM). From bothspike and Local Field Potential (LFP) recordings, we estimated the functional connectivity maps by usingthe Normalized Compression Similarity (spikes) and the Phase Synchrony (LFPs). Then, by using graph-theoretical statistics (clustering coeffient, characteristic path length, small-worldness), we characterizedthe functional map 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 map reconfigurations occur both in spike and LFPs,small-worldness is nonetheless preserved (iii) The stimulus triggers a significant increase of inter-area LFPconnections without modifying the topology of intra-area functional connections; (iv) Through computersimulations of the fundamental concept of cell assemblies, transient groups of activating neurons can bedescribed by small-world networks.Our results suggest that neural activity among neurons from multiple areas of the rat somatosensorysystem allows for the integration of local computations occurred in distributed functional cell assembliesaccording to the principles of SWNs.
Author Summary
Cell assemblies (sequences of neuronal activations), seem to represent a functional unit of informationprocessing. However, it remains unclear how groups of neurons may organize their activity during in-formation processing, working as a sole functional unit. One prominent principle in complex networktheory is covered by small-world networks, in which nodes are easily reachable by each other and areorganized in highly dense clusters. Small-world networks are already found on large-scales in human andprimate brain areas while their presence at the neuronal level remains unclear. The aim of this work wasto investigate the possibility that functional, related neural populations, encompassing multiple brain re-gions, can be organized in a small-world network. We investigated the coherent neuronal activity amongmultiple rat brain regions involved in somatosensory information processing. We found that the recorded a r X i v : . [ q - b i o . N C ] J a n neuronal populations represented small-world networks maintained during stimulations. Furthermore, byusing computer simulations, we inferred that small-world networks represent a plausible topology for cellassemblies.This work suggests that the coherent activity of neurons from multiple brain areas in somatosensorysystem, which can comprise individual functional units, promotes the integration of local computations,the functional principles of small-world networks. Introduction
Neurons in the brain form highly composite networks sustained by a complex and variously distributedthread of synapses. Although the detailed anatomical connections are still under investigation [1] ithas been shown that, during in-vivo recordings, active neurons form functional assemblies that do notnecessarily depend on the underlying anatomical connectivity [2]. Therefore, while anatomical connec-tivity is stable for relatively long times, the functional connectivity is highly dynamic and depends onthe particular tasks triggered by internal and external events. Critically, the organization principles thatcontrol the functional connectivity among single neurons and small neuronal populations are still poorlyunderstood [3]. The early hypotheses [4], that distributed groups of neurons can operate as a functionalunit through coordinated activities, have now been supported in observations of transient cell assembliesover different cortical regions [5]. It remains, however, unclear how these dispersed neurons may gatherin a functional unit for information processing.In the last 15 years, thanks to a substantial advancement in the complex network theory, Small-WorldNetworks (SWNs) emerged as a paradigmatic model and provided the analytical tools to explain a largeset of complex networks in the most diverse scientific areas [6, 7]. Furthermore, it has been suggestedthat small-world networks represent optimized structures accomplishing adequate balance between in-formation transfer efficiency and reliability [7–12]. Following this trend, large networks of brain areas,when studied from a macroscopic perspective by imaging techniques, have been shown to functionallyorganize as SWNs [13, 14]. Contrarily, the functional connectivity emerging at the microscopic level ofsingle neuron networks is still poorly understood, although recent works provided evidence of SWNs insmall local neuronal populations [15–23]. Specifically it is not clear whether functional groups of neuronsfrom multiple brain regions display a SWN organization.In the present work we address these points by means of spike and LFP simultaneous multi-array record-ings 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 re-gions are actively involved in the tactile information processing in somatosensory system of mammalbrains. Because spikes and LFPs represent a dualistic picture of what (respectively) neurons and localneuronal populations are doing, we decided to test the working hypothesis for both signals. To estimatethe strength of functional connectivity in spiking activity we introduced a novel function that allows usto detect long-range temporal dependencies which occur in thalamocortical interactions (NCS) that areotherwise unrecognizable by correlation analyses. Instead, the LFP functional connectivity was estimatedby a standard measure of phase synchrony.Our results confirm the presence of SWNs in crucial stations of the rat somatosensory system, duringspontaneous activity and provide evidence that, during stimuli, small-world organization principles arepreserved in spite of massive topological reconfigurations. In addition, further results obtained by com-puter simulations showed that the observed functional conditions may be the expression of a large-scalefunctional substrate composed by small neural groups with small-world topologies where node member-ship to each group is dynamical.Consistent with other findings on large-scale brain networks [5, 27–29], our results provide evidence for amodel of functional organization where distinct functional neuronal assemblies are sparse and encompassmultiple brain areas.
Results
Functional Connectivity in Spiking Activity
Our first aim was to estimate the topology of functional connectivity graphs obtained from simultaneousspiking activities of neuronal populations in VPL, CM and S1. Because specific tactile stimulation canpotentially exert a powerful effect on functional connectivity, we set out to evaluate functional connectiv-ity graphs both during spontaneous and evoked activities. In order to compute the graph’s connectivityin the two conditions, we performed pairs of 10 minutes recordings, the first with no stimuli, the secondby delivering fast tactile pulses (see methods) at randomized intervals (500 ±
200 ms) on the five digits.We estimated the functional connectivity in neuronal spiking activity by using the Normalized Compres-sion Similarity (
N CS ) function to detect functional couplings between pair of neurons.
N CS is definedin the range [0 , N CS estimation on allpossible pairs of simultaneously recorded neurons, the resulting adjacency matrices were binarized by athreshold in order to obtain the functional connectivity graphs. Subsequently, the small-world statistics C (clustering coefficient), L (characteristic path length) [6], S and ω were estimated on these graphs [18,30].The first measures the tendency of neurons to segregate into separate sub-networks, the second measuresthe 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 originalnetworks but preserving the node degree distributions. Finally, we further characterized the resultinggraphs by analyzing the node degree distribution, the centrality and the community structure. Specifi-cally, we wondered if stimuli provoked changes in the node degree or the node centrality distribution orif they affected the number of communities.First, we measured the small-worldness statistics in our spontaneous activity recordings and we foundthat S and omega stated that such graphs can be considered small-world networks ( S > ω close to0), irrespective of the time window used for detecting the functional connections (50, 250, 500, 1000 ms;Table 1).To note that increasingly larger windows make smaller N CS values because statistical dependenciesbecome more and more sporadic. Consequently, the threshold decreased when windows became larger.Time windows larger than 1 second reported no small-world network configurations. Then we comparedthe average functional connectivity in time windows of 100 ms duration before and after the stimulusonset. 100 ms was a reasonable time constraint in order to consider most of the thalamocortical inter-actions. The functional connectivity graphs underwent significant topological reconfigurations after thetactile stimulus delivery. As an example, in Fig. 2, the neuron 32 switches from a fully disconnectedstate (Fig. 2A, relative to pre-stimulus condition) to highly connected (hub) state (Fig. 2B) after thestimulus onset. Surprisingly, notwithstanding these profound changes, the small-worldness computed onthe pre- and post-stimulus functional connectivity graphs were not significantly different (
P > .
13 inall comparisons, non-parametric Wilcoxon ranksum test, Table 2). Therefore, the average small-worldproperties were not perturbed by tactile stimulations. Furthermore, we performed a correlation analysisbetween the quantitative changes in functional connectivity graphs and the firing rate variations inducedby tactile stimuli. We found that the changes induced by stimuli were significantly greater (
P < . R = 0 . P < . P = 0 .
48, ranksum test). Moreover,by analyzing the distribution of node centrality, we found that graphs in pre-stimulus conditions (mean, µ = 22 .
84, standard deviation, σ = 1 .
63) had smaller betweenness centrality (
P < . µ = 25 . σ = 2 .
01) indicating that stimuli made more cen-tralized nodes. Moreover, by further analyses we found that the betweenness centrality was not equallydistributed over the three regions (VPL,S1,CM) both in pre- and in post-stimulus conditions: performingan ANOVA-1-way test over the distributions of centrality in the three regions, we found that these hadnot the same mean ( P = 0 .
007 in pre-stimulus, P = 0 .
002 in post-stimulus). Furthermore, the Figure3C shows that the centrality of neurons from CM increased ( P = 0 . P = 0 . P = 0 . µ = 5 . σ = 1 .
45) had less communities (
P < .
01, ranksum test) than graphs in post-stimulus conditions( µ = 3 . σ = 1 . Functional Connectivity in LFPs
Our second aim was to estimate the topology of functional connectivity graphs obtained from the LFPactivity recorded simultaneously to the spiking activity. We estimated the LFP functional connectivitymaps by following the same sequence of analyses used for spiking activity. First, we estimated the func-tional connectivity between all possible electrode pairs in order to generate the adjacency matrix, then webinarized this matrix by using a variable threshold (see Methods) and finally we computed the networkstatistics.In order to estimate the functional connectivity between pairs of channels we used the phase synchronymeasure γ (see Methods) [31]. Phase synchrony is more suitable than NCS for continuous signals and ithas been widely applied to EEG and LFP analyses [32]. It is normalized in the range [0 , γ values thusthe threshold were decreased when windows became larger. Time windows larger than 1 second reportedno small-world network configurations. We wondered if the synchrony was crucial for LFP phases orsmall-worldness can emerge even without tight time constraints. By using NCS rather of gamma func-tion to estimate functional connectivity, small-world networks did not appeared (data not shown).For evoked activity recordings, stimulus occurrence triggered significant topological reconfigurations infunctional graphs, as shown, as a representative case, in Fig. 4. The picture reports a functional con-nectivity graph estimated before (Fig. 4A) and after (Fig. 4B) the stimulus onset. During spontaneousactivity, CM (green), VPL (blue) and S1 (red) showed tight intra-area synchronization and poor inter-areasynchronization. After the stimulus onset inter-area synchronization emerged while intra-area synchro-nization was maintained roughly constant.We asked whether these observations on a single recording could generalize to our full dataset. In orderto 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, measuredfrom the same recording channels (see Methods), was larger than the mean plus 5 times the standarddeviation of the basal firing rate. We found that the number of intra-area (local) functional connectionswas 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 averagenumber of inter-area (or global) functional connections was substantially larger for all possible inter-areacombinations (
P < .
002 for CM-VPL, CM-S1, VPL-S1, ranksum test).Subsequently, we compared pre- and post-stimulus functional graphs by computing the networks statisticsin 50 and 100 ms windows before and after the onset of tactile stimuli. We challenged both windows sizeto establish the possible dependency of diverse time constraints during synchronization event. We foundthat both conditions (pre- and post-stimulus) exhibited small-world properties and the statistics werenot different (Table 4,
P > .
38 in all comparisons, ranksum test) and the window size per se did notchange the small-worldness values. Interestingly, the tactile stimulus onset triggers a substantial increasein the number of inter-area connection but does not have a significant effect on the intra-area connections(Fig. 5A). We concluded that small-world statistics are preserved in LFPs during both spontaneous andevoked activity. Thereafter, we analyzed the extracted graphs in pre- and post-stimulus conditions bytheir node degree distributions that were not significantly different (Figure 5B, P = 0 .
83, ranksum test).By analyzing the distribution of node centrality, we found that graphs in pre-stimulus conditions ( µ =21 . σ = 2 .
17) had smaller betweenness centrality (
P < .
03, ranksum test) than graphs in post-stimulusconditions ( µ = 23 . σ = 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 post-stimulus conditions: performing an ANOVA-1-way test on the distributionsof 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 . P = 0 . P = 0 . µ = 2 . σ = 0 .
71) had less communities (
P < . µ = 2 . σ = 0 . Functional Connectivity in Simulated Spiking Activity
In the previous sections we provided evidence, both at a pre- and post-synaptic level, that cortical andsubcortical neurons can be functionally organized as small-world networks. Our results are consistentwith recent findings from large scale cortical recordings [10, 33, 34] and suggest the existence of sparsefunctional networks composed of neuron assemblies that encompass multiple cortical and subcorticalareas; importantly, each functional network would recruit its units not necessarily in close proximity [5,27].However, it remains unclear what kind of functional substrate supports the observed networks, i.e. howlarge-scale network of neurons are organized such that sampling from a small subset of nodes, we canobserve small-world configurations.To address this question, we investigated the hypothesis that neuronal networks are arranged in functionalgroups, recalling the concept of cell assembly, where the membership to each assembly is dynamical[4] and assembly neurons can be spatially distributed. To evaluate this conjecture we developed anartificial neuronal network where nodes are arranged in dynamical groups and emit spikes followingsmall-world criteria. Subsequently, we sampled the activity of small ensembles of nodes and compare thenetwork statistics from the synthetic spiking activity with those obtained by recordings. Once assertedsuch hypothesis, we estimated how many functional groups can coexist within a neural network keepingconsistency with experimental observations.For these purposes, we implemented a large-scale simulation of N neurons that embeds k small-worldnetworks. Then we estimated small-worldness by systematically varying k . We simulated a neural networkof N = 100000 neurons, about five fold the estimated size for a rat cortical column in the somatosensorysystem [35]. The size of each functional group was chosen by sampling from a discrete uniform distribution ∼ U (50 , N/ k = 0 .
01, peakingwhen k = 0 . k up to 0 .
5. At peak we found comparable statistics values withthose measured experimentally both on spiking ( P = 0 .
60, ranksum test) and LFP activities ( P = 0 . S and ω ) of small-world networks per number of nodes that rangedwithin 0 .
01 and 0 . Discussion
In this paper we show that spiking and LFP activities of neurons, in three stations of the somatosen-sory system of rat brains, present clear signs of small-world functional organization with sub-secondinvariance. Furthermore, we show that this distinctive functional organization persists in the presenceof tactile stimuli, independently of the neural response intensity. Finally, results obtained by means ofcomputational models suggest that small-world networks may represent a consistent and formal modelfor cell assemblies.
Small-world in brain networks
Studies in anatomical brain connectivity discovered that small-world network architectures are a dis-tinctive trait in animals [19, 21, 37], yet with different degree of brain development, including humans[10, 13, 38]. The complementary observation that brain pathologies like epilepsy [39] and schizophre-nia [40] show small-world network disruptions provides further support to the potential interpretive andexplanatory strength of small-word topology. It should be noted that the methodological approaches usedto study the functional connectivity at macroscopic and microscopic levels are slightly different. Studieswith fMRI estimate functional connectivity comparing the BOLD signal of each node to a seed voxelwhereas the functional strength between neurons is computed comparing the electrical activity of all pos-sible neuron couples [41]. Setting apart the technical incongruences and focusing on neuronal functionalconnectivity and its underlying topology, a number of beautiful studies found an intrinsic small-worldtopology in single visual areas of primates [15,17,19,20,22,23] and in neuronal cultures [16]. These elegantapproaches gave access to the richness of the local functional connectivity scaled at the neuronal level.However, none of these studies provided a description of an extended connectivity, involving differentneuronal populations nested in close or far brain regions.Yet, according to the work of Bassett et al. [11], it is suggested that the functional organization in small-word networks may be ubiquitous at different spatial and temporal scales and different experimentalconditions. This organization could represent a stable, and even a necessary, scheme for informationprocessing in brain areas, networks and neurons.
Tactile information processing
In this study our aim was to explore the functional connectivity of different brain areas involved insomato-sensory information processing and to envisage a potential broadening to other brain circuits.As it is known, the precise mechanisms regulating the rich and complex neural thread of brain areasinvolved in the construct of tactile processing are still far from being clarified. However, a response tosensory inputs is well recognized for a number of areas, well including the areas we chose for this study.Namely, several studies showed that mechanoreceptor signals reach both VPL and CM thalamus. Then,the former outputs directly to innervate S1, while the latter outputs are more diffused and addressedto the higher sensorimotor cortical layers [25, 26, 42, 43, 51]. We observed that in the resting state,VPL, CM and S1 showed substantial mutual connectivity and, in addition that, under stimuli, theyestablish an even more integrated architecture enabling fast information exchanges. This functionalconnectivity could be ascribed to the necessity to create larger functional units and nodes with highernetwork load. In accompaniment of these changes, however, the topology of the small-world networksnever degenerated. This dualistic behavior, redistribution of connectivity and permanence of topology,appears an interesting interpretive key, potentially extendable to other complex brain networks and usefulto analyze also pathological signs in brain dynamics, like epilepsies or schizophrenia, where the detectionof disrupted network or loss of topological hallmarks may help novel nosological classifications.
Cell assemblies and the synthetic modell
The analyses from computational models suggest that a cell assembly, i.e. a transient functional unitcomposed by neurons potentially distributed in separated brain regions [4], appear to be organized asa small-world network. The null hypothesis that such groups were random networks does not matchexperimental evidence.Our model assumed that neurons may variably join in concurrent assemblies and we empirically estimatedthat the admissible number of assemblies that a neuron could participate in, covered a range between0.01 and 0.1. It could be objected that the number of shared neurons appears to be relatively high.However, it might be noted that neurons in stimulated conditions probably get a balance between thenumber of assemblies and the inherent metabolic cost [44]. As a consequence, neuron sharing could abatethe metabolic cost of network activations, behaving as a parsimonious limiting factor to an unregulatedgrowth of recruited functional units.Our synthetic model tried to answer a double hypothesis: the first that the topology expressed by oursmall graphs could be a nested and natural expression of a homologous larger population. The reversequestion appears to be if the hypothetical extended topology that we designed in our synthetic modelcould 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 hy-pothetical larger synthetic networks.Our results are also strictly related research suggesting that neuronal oscillations enable selective anddynamic control of distributed functional cell assemblies [5]. We could enrich such a scenario, speculatingthat LFP coherent activities may reflect the integration of local computations which occurred in these dis-tributed cell assemblies. Thus, the small-world topology expressed by synchronized and distributed LFPphases 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 wasindeed very low and, as it comes from other studies, low enough to avoid important suppressive actionsof the neural activity [45]. However, it is implicit that conclusive studies with waking animals are neededto definitely address the existence of these functional topologies (at least) in these three brain regions.We proved the effectiveness of the NCS measure to capture the long-range dependencies and although suchmethod is inspired by related approaches [46–49], it has never been applied before on electrophysiologicaldata. A potential drawback of such approach is represented by its computational cost. Indeed, a numericalincrease of recorded neurons (up to thousands or more), would need a cubic increase of the computationaltime required to calculate the NCS adjacency matrices.Ultimately, inferences from computational models are often risky. In fact, our findings mainly endorsebut not prove the hypotheses. Other topologies, rather than random networks, could deliver consistentresults. A set of other debated topologies (hierarchical modular networks, scale-free networks, etc.) couldbe also considered.
Conclusive notes
In an evolutionary perspective, small-world topologies appear preferentially selected among networktopologies under the natural constraints (efficiency and reliability) of brain expanding complexity in themammal phylogenetic tree [12]. It can satisfy the necessity of internal input integration and grants forbest responses to environmental requirements. Moreover, small-world networks seem coherent with therecent advancements in the physiology of neuronal networks. More punctually, from a functional pointof 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 efficiently participateto the information process.
Materials and Methods
Ethical notes
All the animals have been used in accordance to the Italian and European Laws on animal treatment inScientific 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 UnitedNations Educational, Scientific and Cultural Organizations (UNESCO), the Council for InternationalOrganizations 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 classifiedas “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, 300 −
400 g) were chosen in theset of 11 animals employed in the recordings. All the animals were maintained in a 16 / . . . . − . . × × − . − . − − . − . o slantand driven at least to 5500 µ m in depth and then advanced by a second electronically driven microstepper(AB Transvertex, Stockholm) until responses were observed to peripheral test stimuli. The neuronalrecordings were obtained with two types of matrices, a vertical array devoted to the cortical recordings andtwo planar matrices devoted to the thalamic recordings. The vertical array was a multitrode (MultitrodeType 1, Thomas RECORDING GmbH, Giessen, Germany) with 8 gold contacts (125 µ m contact spacing)with an average impedance of 1 . × − µ m, tip impedance 0 . − Tactile Stimulation −
250 ms (see Fig. 6C). The stimuli semi-randomness was adopted to avoid habituation [53].
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 andrecorded with bandpasses at 6 kHz/300 Hz for spikes, 180 Hz/1 Hz for Local Field Potentials. The datastored were analyzed off-line both using Matlab and by locally developed software. The neural firingrates had a mean of 31 . . . Wave clus
MATLAB toolbox [55]. Sortedcells 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 aswell. Recorded neurons were uniformly distributed over the recording matrices and every electrode showdistinct neural activity otherwise the matrix was repositioned. At the end of this process, we collected atotal of 391 neurons (56 ±
17 in each experiment) out from the set of the acquired signals. Distributionof firing rates and inter-spike intervals is shown in the Supplementary Figures 6A-B.The timestamps of spike occurrences were represented by binary sequences where 1’s labeled a spike. Weconsidered time bins of 1 ms thus avoiding occurrence of multiple spikes within the same bin. Finally, wesplit each sequence into fixed-length (from 50 to 1000 ms) overlapping windows (Fig. 1C), thus obtainingan ordered set of equal length windows.
Functional connectivity by spike-train similarities
Interactions between neurons can generate very complex, time-delayed and asymmetric patterns. Thisrepresents a potential problem in experimental configurations where the communications between distantneurons 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 lasttens of milliseconds [53, 56–58]. In a toy example (shown in Figure 1A-B) it is assumed the existence ofthe interaction between neurons A and B (A → B). The neuron A is recorded and its activation triggersthe long chain that, from site A, produces firing activity to neuron B in another brain region. The0direct path between neurons C and B allows, instead, for synchronous spike patterns well detectable bycorrelation analysis (Fig. 1B).In general, in simultaneous recordings from separated brain sites, it is unlikely to find couples of physicallywired neurons although axonal projections connect the sites. It’s definitely more probable that spikingactivity relations may reflect a complex anatomical substrate, where chains of activations exist betweenthem provoking the observed coherent activity. Such a problem can be solved by mathematical tools ableto model arbitrarily long temporal relationships. In this work, we proposed a novel framework whereinspike trains with arbitrarily long temporal dependencies are modeled by Markov stochastic models. Typi-cally, in regular Markov models, each state depends only on the previous state while higher-order Markovmodels suffer from high state-space computational complexity [59]. Here, we used the Variable-orderMarkov Models (VMMs) [60, 61] because they are able to overcome these limitations. Lossless compres-sion algorithms (LCAs) represent one of the most efficient techniques to estimate VMMs [61]. Within theset of LCAs we chose the Prediction by Partial Matching (PPM) algorithm [62, 63] which is consideredthe best match between prediction accuracy and speed [61]. The last step consists to build a similarityfunction for this kind of spike train stochastic models.In the last 15 years, Vitanyi and colleagues have developed a function, the Normalized Compression Dis-tance (
N CD ) [64–67], which estimates the distance between symbolic sequences. This function performsthe estimation directly by the sizes of the compressed sequences. In fact, can be proved that the betteris the VMM estimation the shorter is the compressed sequence size. In this work, we redefined the
N CD function in order to reverse its assigned values pointing to similarity instead of distance. We called suchfunction the Normalized Compression Similarity (
N CS ). Formally, given that x and y are two neuralsequences (e.g. spike trains), the N CS is defined as follows:
N CS ( x, y ) = 1 − N CD ( x, y ) = C ( x · y ) − min { C ( x ) , C ( y ) } max { C ( x ) , C ( y ) } (1)where the C function represents the compressed sequence length and · is the sequence concatenationoperator (e.g. 0101 ·
101 = 0101101). If
N CS ( x, y ) is close to 1, the sequences x and y are consideredsimilar. If close to 0, the sequences are strongly dissimilar.We therefore evaluated the N CS function on time windows (50-1000 ms) of the recorded (binary) spikingactivity (Fig. 1C-D) assuming that relative high values of similarity corresponded to actual functionalconnections (Fig. 1E). To note that significant
N CS values do not imply significant correlations but theopposite is true. Although the
N CS is asymmetric function, it is not relevant for the causal interactionanalysis and consequently for effective connectivity. Because
N CS has never been used as method toestimate the functional connectivity in neurophysiological data (although many LCAs has been used asinformation measure of synaptic efficacy [46] and spike train similarities [47–49]), our aim was to provethat
N CS effectively captures similar asynchronous spike patterns in comparison to correlation coefficient.Furthermore, we demanded that
N CS did not bias its estimations with respect to independent spike trainsgenerated by simulations.To address these requirements we performed two experiments: i) we evaluated the
N CS and the Pearsoncorrelation coefficient over a set of independent binary spike trains and ii) we evaluated the
N CS and thecorrelation coefficient with couples of binary spike trains (1000 ms long) containing an identical, shortand random spike pattern (50 ms) that is fixed and centered in the first train and drifts (from left toright) in the second trains (see Fig. 7B, overall 37 drifts). The latter procedure creates synthetic spiketrains wherein the common pattern is increasingly distant.In the first experiment we found that the NCS method did not show any significance for independentspike trains as well as for the Pearson coefficient ( µ = − . · − , σ = 0 .
005 for Pearson coefficientand µ = 0 . σ = 0 .
004 for
N CS , Fig. 7A). This guarantees that the
N CS function is unbiased toindependent distributed binary spike trains. In the second experiment, the
N CS proved its efficacy todetect long-range interactions where Pearson coefficient fails. In fact, during the whole shifting epochs,1significance held while Pearson coefficient showed significant only when the two spike patterns were veryspatially close (Fig. 7C, drift number 18-23).
Functional connectivity by LFP phase synchrony
LFPs are low frequency signals reflecting a wide range of synaptic events. In this work we investigatedthe synchrony of LFP phases originated in different recording sites during spontaneous and tactile evokedactivities. We measured phase synchronies between two recorded LFP sequences ( x and y ) by the followingfunction γ ( x, y ) = (cid:12)(cid:12)(cid:12)(cid:68) e i ( arg ( H ( x )) − arg ( H ( y ))) (cid:69)(cid:12)(cid:12)(cid:12) (2)where e is Napier’s constant, H is Hilbert Transform, arg is the argument function and i is the imaginaryunit. The Hilbert transform and the argument were computed with, respectively, the hilbert and the angle Matlab functions [31, 32]. When γ ( x, y ) is equal to 1 (0), then x and y are perfectly synchronous(asynchronous). Complex Brain Network
By using the NCS and γ functions, we estimated the functional connectivity of the recorded neuronalnetworks. We first split each recorded sequence into equal-length time windows (Fig. 1) and then wecomputed the adjacency matrix for all neurons or electrodes. The resulting matrices exhibited valuesin the unitary interval. We repeated the analyses with different window sizes from 50 ms to 1 s (fixingthe maximum dependency order to half of the window size). The functional connectivity extracted fromextracellular recordings can be represented by graphs.A variable threshold (typically equal to a higher percentile of the weight distribution and vary between0 . .
8) selected the strongest connections, thus allowing for the construction of the functionalconnectivity graphs. The choice of a threshold is related to the window length used to estimate theconnectivity strength. Larger windows produce less connections and to keep the network sparsity quiteconstant, smaller thresholds must be chosen.For the analysis of these graphs, we introduced a set of common statistics from the Complex NetworkTheory able to detect possible matches between the extracted graphs and prominent topologies likesmall-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 muchless nodes than the total node number N ( N (cid:29) K (cid:29) ln( N ) (cid:29) C ) and the characteristic path length ( L ) [6,69]. The former measures how closethe neighbors of a node are to being a clique. The latter estimates the average shortest path length in thegraph, i.e. how much the nodes are accessible. Both measures, implemented in a Matlab toolbox [69],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 thesame graphs subject to randomization or latticization (often called null networks) [69]. Both proceduresguarantee that the node degree distributions of the original graphs were preserved. We computed, byusing the Matlab function randmio und.m , the randomized version of our graphs and we estimated C r 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 ω = L r L − CC l [30] state that, respectively, if S > ω is close to 0, the graph can be considereda small-world network. The functional graphs obtained by our analysis were further characterized tostudy the information flow. For this aim, we computed a measure of centrality (betweenness) for graphnodes [70], an estimate of the number of shortest paths from all vertices to all others that pass throughthat node. Because it can be interpreted as a measure of the load of a node within the network [69], thedistribution 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 emergefrom graphs by applying a clustering algorithm to nodes [71]. In this work, communities represent theaggregated functional units under investigation and by counting the number of them we can understandhow node graphs are aggregated in each experimental condition. Ultimately, we analyzed networks thatevolved in time dropping and recruiting nodes and connections and networks from different experimentalconditions. Such a methodology requires the discussion of potential issues [41, 72].First, unconnected nodes are rarely but can occur after adjacency matrices were binarized. For this reason,we removed graphs in which less than the 99% of nodes are connected. Second, network statistics wereapplied on network with different sizes (for spiking activity) because recording sessions return a variablenumber of active neurons. However, by analyzing the observed variance of network size we concludedthat C and L cannot be significantly affected by our network size changes. Significant changes appearfor synthetic networks that increased their size by orders of magnitude. However, we discarded graphsthat are outliers (beyond 5 th and 95 th percentile) of the node, edge and density number distributions inorder to obtain a better homogeneity. 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 [73]. ALayer layout helps to understand the information flow in the network and the specific roles performedby each node. Graphs in the insets are visualized by using the function
GraphPlot of Mathematica thatimplements the spring-electrical algorithm [74] 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 physiologicalevents. Guided by results from recordings, we hypothesized that a cell assembly can be functionallyorganized as a small-world network and we investigated such hypothesis by computer simulations. Asoriginally thought, a cell assembly can be constituted 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 experimentalfacts? (ii) assuming the previous as true, how many assemblies (small-world networks) can exist over aset of neurons keeping consistency with observations?Typical simulation frameworks require a choice of neuron models (Integrate and Fire, Izhikevich, Hodgkinand Huxley, etc.) and of a defined network layout (nodes and connections). These choices can be verycrucial and become even more important if the aim of the study is the functional organization of theunits. For this reason, we proceeded following an unconventional approach assuming that cell assembliesare effectively functionally organized as small-world networks on large-scale networks and, sampling theactivity of a random subset of nodes, we wondered if such activity can elicit small-world organization aswell.3So we first assumed that a set of small-world networks exists over a set of available nodes. If we considereach small-world network as a cell assembly [4], the small-world networks of our model can share theirnodes.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 belongto more than one small-world topologies (Fig. 8A-D).5. According to the typical amount of neocortical neurons in a microcolumn [36], the sizes of small-world networks are randomly sampled by a uniform distribution ∼ U (50 , N/ N is thenumber of nodes.From a formal point of view, the network structure can be interpreted by a Multigraph G = ( V, E ), theset V represents the nodes and the multiset E represents the unordered list of edges inherited by eachsmall-world network.To establish the node functioning we chose a primary criterion claiming that brain processing takes placeby functional segregation and integration [68]. 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 areclosely related. Computations performed in (specialized) peripheral nodes, can be subsequently integratedin central nodes by virtue of the small average shortest path length of small-world networks. Indeed, itcan be easily shown, by using simple computer simulations, that nodes in small-world networks had awide spectrum of centrality values. This had a heavy-tailed distribution in contrast to random networksdistributed 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 thecentrality values estimated by using the centrality. The centrality of nodes can be estimated by severalmeasures (e.g. betweenness, pagerank, node degree, etc.).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 hypoth-esis).2. At most two randomly chosen small-world networks are allowed to run independently in a timewindow. A single neuron can work either alternately or concurrently within different networks ofpre and post-synaptic nodes. This assumption claims that neurons are not exclusive and that theycan partake simultaneously to many (up to 2) tasks.3. Neurons express their spikes within 10 ms time windows. Neurons fire following a centrality criterionbased on the node betweenness centrality. A neuron with high betweenness fires its spikes after alower betweenness neuron. This dynamic picture respects the information integration-segregationparadigm (Fig. 8E-H) [68]. The choice of the window length is arbitrarily and is proportional tothe firing rate.4. Simulation timesteps are set to 1 millisecond. Spike propagation times are uniformly distributedwithin the range [1 , Input : n swn , n nodes , timestep; Output : the binary n nodes -by- timesteps matrix trainsp ← . k ← for i ← to n swn do n ← subsample( n nodes ) ; G ← random watts strogatz( n, k, p ) ; wins ← computeWindows( G ) ; foreach time windows in wins doif . < rand() then centrality ← BetweennessCentrality( G ) ; foreach node A in centrality do trains [ A, wins + centrality ( A )] ← foreach output node B of A do trains [ A, wins + centrality ( A ) + centrality ( B )] ← endendendendend where the function subsample() selects a subset of nodes out of the set of available ones. The function random watts strogatz() returns a random graphs built following the Watts-Strogatz algorithm, thefunction computeWindows() computes the length of the execution windows for the current graph andreturns the list of such windows. The function rand() returns a uniformly generated random numberbetween 0 and 1. At last, the function BetweennessCentrality() computes and returns the centralityof each node.The model has been developed using the Python environment [75]. For the generation of small-wordnetworks we used the networkx package (available at http://networkx.lanl.gov/). Fig. 8I shows a rasterplot obtained by sampling the simulated activity. The darkest blue represents the no-spike event andthe other spike colors are associated with the specific small-word topologies that generated them. Thealgorithm describes a core loop where each small-world network is first randomly created by the libraryroutine watts strogatz graph (probability of rewiring equal to 0 . Acknowledgments
We wish to thank M.E. Giardini, L. Bonini and C. Toschi for helpful suggestions. Finally, we thankKatherine Davis for the revision of the english.5
References
1. Sporns O, Tononi G, K¨otter R, (2005). The Human Connectome: A Structural Description ofHuman Brain. PLoS Computational Biology 1(4): e42.2. Feldt S, Bonifazi P, Cossart R, (2011). Dissecting functional connectivity of neuronal microcircuits:experimental and theoretical insights. Trends in Neuroscience 34(5), 225-236.3. Bassett DS, Greenfield DL, Meyer-Lindenberg A, Weinberger DR, Moore SW, et al. (2010) Ef-ficient Physical Embedding of Topologically Complex Information Processing Networks in Brainsand Computer Circuits. PLoS Comput Biol 6(4): e1000748.4. Hebb D O, (1949). The organization of Behavior. John Wiley & Son Inc. New York.5. Canolty RT, Ganguly K, Kennerley SW, Cadieu C F, Koepsell K, (2010). Oscillatory phase couplingcoordinates anatomically dispersed functional cell assemblies. Proc. Natl. Acad. Sci. USA 107, 17356-17361.6. Watts DJ, and Strogatz SH, (1998). Collective dynamics of ’small-world’ networks. Nature 393,440-42.7. Albert R, Barabasi A, (2002). Statistical Mechanics of Complex Networks. Reviews of ModernPhysics 74,47-97.8. Barrat A, and Weigt M, (2000). On the properties of small-world network models. The Eur. Phys.J., B 13, 547-560.9. Barahona M, Pecora LM, (2002). Synchronization in small-world systems. Physical Reviews Letter89(5), 054101.10. Bassett D, Bullmore E, (2006). Small-World Brain Networks, The Neuroscientist, 12(6):512-523.11. Bassett DS, Meyer-Lindenberg A, Achard S, Duke T, and Bullmore E, (2006). Adaptive re-configuration of fractal small-world human brain functional networks. Proc. Natl Acad. Sci. USA130, 19518-19526.12. V´ertes PE, Alexander-Bloch AF, Gogtay N, Giedd J, Rapoport J, Bullmore ET, (2012) Simplemodels of human brain functional networks. Proc. Natl. Acad. Sci. (USA), 109(15): 5868-73.13. Bullmore E, Sporns O, (2009). Complex brain networks: graph theoretical analysis of structuraland functional systems. Nature Reviews Neuroscience 10, 186-198.14. Horwitz B, (2003). The elusive concept of brain connectivity. Neuroimage 19, 466-470.15. Bettencourt LM, Stephens GJ, Ham MI, Gross GW (2007) Functional structure of cortical neuronalnetworks grown in vitro. Phys Rev E Stat Nonlin Soft Matter Phys 75:021915.16. Downes JH, Hammond MW, Xydas D, Spencer MC, Becerra VM, Warwick K, Whalley BJ, NasutoSJ (2012) Emergence of a small-world functional network in cultured neurons. PLoS Comput Biol8:e1002522.17. Humphries MD, Gurney KN (2006). A means to an end: validating models by fitting experimentaldata, Neurocomputing , 20 (10), 1892-1896.18. Humphries MD, Gurney KN (2008). Network ‘small-world-ness’: a quantitative method for deter-mining canonical network equivalences, PLoS One, 3:e0002051.619. Yu S, Huan D, Singer W, and Nikolic D (2008). A small World of Neuronal Synchrony. CerebralCortex 18, 2891-2901.20. Gerhard F, Pipa G, Lima B, Neuenschwander S, and Gerstener W, (2011) Extraction of networktopology from multi-electrode recordings: is there a small-world effect? Frontiers in ComputationalNeuroscience 5, 4.21. Varshney LR, Chen BL, Paniagua E, Hal DH, Chklovskii DB, (2011). Structural Properties of theCaenorhabditis elegans Neuronal Network. Plos Computational Biology 7, e1001066.22. Takahashi N, Sasaki T, Matsumoto W, Matsuki N, Ikegaya Y (2010). Circuit topology for syn-chronizing neurons in spontaneously active networks. Proc Natl Acad Sci 107(22), 10244-10249.23. Pajevic S, Plenz D (2009) Efficient network reconstruction from dynamical cascades identifiessmall-world topology of neuronal avalanches. PLoS Comput Biol 5:e1000271.24. Van Der Werf YD, Witter MP, Groenewegen HJ, (2002). The intralaminar and midline nucleiof the thalamus. Anatomical and functional evidence for participation in processes of arousal andawareness. Brain Research Reviews, 39:107-140.25. Sherman SM, Guillery RW, (2001). The role of the thalamus in the flow of information to thecortex. Phil. Trans. R.Soc. Lond. B, 357:1695-1708.26. Grossberg S, Versace M, (2008). Spikes, synchrony, and attentive learning by laminar thalamocor-tical circuits. Brain Research, 1218:278-312.27. Fries P, (2005). A Mechanism for Cognitive Dynamics: Neuronal Communication Through Neu-ronal Coherence. Trends Cogn Sci. 9,474:480.28. Siegel M, Donner T, Engel A, (2012). Spectral fingerprints of large-scale neuronal interactions.Nature Reviews Neuroscience, 13:121-134.29. Womelsdorf T, Schoffelen J M, Oostenveld R, Singer W, Desimone R, (2007). Modulation of Neu-ronal Interactions Through Neuronal Synchronization. Science 316, 1609-1612.30. Telesford QK, Joyce KE, Hayasaka S, Burdette JH, Laurienti PJ. (2011) The ubiquity of small-world networks. Brain Connectivity. 2011, 1(5): 367-375.31. Montemurro MA, Rasch M J, Murayama Y, Logothetis NK, Panzeri S, (2008). Phase-of-firingcoding of natural visual stimuli in primary visual cortex, Current Biology, 18: 375-380.32. Pereda E, Quiroga RQ, Bhattacharya J, (2005). Nonlinear multivariate analysis of neurophysiolog-ical signals. Progress in Neurobiology 77, 1-37.33. Zhang N, Rane P, Huang W, Liang Z, Kennedy D, Frazier JA, King J (2010) Mapping resting-statebrain networks in conscious animals. J Neurosci Methods 189:186-196.34. Pawela CP, Biswal BB, Cho YR, Kao DS, Li R, Jones SR, Schulte ML, Matloub HS, Hudetz AG,Hyde JS (2008) Resting-state functional connectivity of the rat brain. Magn Reson Med 59:1021-1029.35. Meyer HS et al., (2010). Number and Laminar Distribution on Neurons in a ThalamocorticalProjection Column of Rat Vibrissal Cortex. Cereb Cortex 20(10), 2277-2286.36. Mountcastle VB, (1997). The columnar organization of the neocortex. Brain 120, 701-722.737. Perin R, Berger TK, Markram H (2011). A synaptic organizing principle for cortical neuronalgroups. Proc Natl Acad Sci USA 108, 5419-24.38. Honey CJ, Ktter R, Breakspear M, Sporns O, (2007). Network structure of cerebral cortex shapesfunctional connectivity on multiple time scales. Proc Natl Acad Sci USA 104, 10240-10245.39. Liao W, Zhang Z, Pan Z, Mantini D, Ding J, (2010). Altered Functional Connectivity and Small-World in Mesial Temporal Lobe Epilepsy. PLoS ONE 5, e8525.40. Micheloyannis S, Pachou E, Stam CJ, Breakspear M, Bitsios P, (2006). Small-world networks anddisturbed functional connectivity in schizophrenia. Schizophr Research 87, 60-66.41. Zalesky A, Fornito A, Bullmore E (2012) On the use of correlation as a measure of networkconnectivity. Neuroimage 60:2096-2106.42. Swadlow HA, Bezdudnaya T, Gusev AG, (2005). Spike timing and synaptic dynamics at the awakethalamocortical synapse. Prog Brain Res 149:91-105.43. Sadikot AF, Rymar VV, (2008). The primate centromedian-parafascicular complex: anatomicalorganization with a note on neuromodulation. Brain Res Bull 78:122-30.44. Fonseca-Azevedo K, Herculano-Houzel S, (2012). Metabolic constraint imposes tradeoff betweenbody size and number of brain neurons in human evolution. PNAS 109(45):18571-18576.45. Wang K, van Meer MP, van der Marel K, van der Toorn A, Xu L, Liu Y, Viergever MA, Jiang T, Di-jkhuizen RM, (2011). Temporal scaling properties and spatial synchronization of spontaneous bloodoxygenation level-dependent (BOLD) signal fluctuations in rat sensorimotor network at differentlevels of isoflurane anesthesia. NMR Biomed 24: 61-67.46. London M, Schreibman A, Hausser M, Larkum ME, Segev I, (2002). The information efficacy of asynapse. Nature Neuroscience 5, 332-340.47. Szczepa´nski J, Amig´o JM, Wajnryb E, Sanchez-Vives MV, (2004). Characterizing spike trains withLempel-Ziv complexity. Neurocomputing 58-60, 79-84.48. Amig´o JM, Szczepa´nski J, Wajnryb E, Sanchez-Vives M, (2004). Estimating the Entropy Rate ofSpike Trains via Lempel-Ziv Complexity. Neural Computation 16, 717-736.49. Zippo AG, (2011). Neuronal Ensamble Modeling and Analysis with Variable Order Markov Models.Ledizioni, Milano, 2011.50. Kohn D, Wickson S, White W, Benson G J, (1997). Anesthesia and Analgesia in LaboratoryAnimals. London: Academic Press.51. Sherman SM, Guillery RW, (2001). Exploring the Thalamus. London: Academic Press.52. Paxinos G, Watson C, (2006). The Rat Brain in Stereotaxic Coordinates, London. Academic Press.53. Storchi R, Zippo AG, Caramenti G, Valente M, Biella GE, (2012). Predicting Spike Occurrence andNeuronal Responsiveness from LFPs in Primary Somatosensory Cortex. PLoS ONE, 7(5): e35850.54. Zanos TP, Mineault PJ, Pack CC, (2011). Removal of spurious correlations between spikes andlocal field potentials, J. Neurophysiol.105, 474-86.55. Quiroga RQ, Nadasdy Z, and Ben-Shaul Y, (2004). Unsupervised spike sorting with wavelets andsuperparamagnetic clustering. Neural Computation 16, 1661-1687.856. Panzeri S, Petersen RS, Schultz SR, Lebedev M, Diamond ME, (2001). The Role of Spiking Timingin the Coding of Stimulus Location in Rat Somatosensory Cortex. Neuron 29, 769-777.57. Latham PE, Lengyel M, (2008). Phase Coding: Spikes Get a Boost from Local Fields. CurrentBiology 18, 349-351.58. Aguilar J, Morales-Botello ML, Foffani G, (2008). Tactile responses of hindpaw, forepaw andwhisker neurons in the thalamic ventrobasal complex of anaesthetized rats. Eur J Neurosci 27(2),378-387.59. Zaki MJ, Carothers CD, Szymanski BK, (2010). VOGUE: A Variable Order Hidden Markov Modelwith Duration Based on Frequent Sequence Mining. ACM Transactions on Knowledge Discovery inData, 4(1):Article 5.60. B¨uhlmann P, Wyner AJ, (1999). Variable Length Markov Chains. The Annals of Statistics 27(2),480-513.61. Begleiter R, El-Yaniv R, Yona G, (2004). On Prediction Using Variable Order Markov Models.Journal of Artificial Intelligence Research, 22,385-421.62. Cleary J, Witten I (1984). Data Compression Using Adaptive Coding and Partial String Matching.IEEE Transactions on Communications, 32, 396-402.63. Teahan W, (1995). Probability estimation for PPM. Proceedings of the New Zealand ComputerScience Research Students’ Conference, University of Waikato, Hamilton, New Zealand.64. Li M, Vitanyi PM, (1997). An Introduction to Kolmogorov Complexity and its Applications.Springer-Verlag, New York, 2nd Edition.65. Bennett CH, Gacs P, Ming Li, Vitanyi PM, Zurek WH, (1998). Information Distance. IEEE Trans-actions on Information Theory 44(4), 1407-1423.66. Chen Z, Francia B, Li M, McKinnon B, Seker A, (2004). Shared information and program plagiarismdetection, IEEE Transactions on Information Theory, 50(7), 1545-1551.67. Cilibrasi R, Vitanyi P, (2005). Clustering by Compression. IEEE Transactions on InformationTheory 51, 1523-1545.68. Tononi G, Sporns O, Edelman GM, (1994). A measure for brain complexity: Relating functionalsegregation and integration in the nervous system. Proc. Natl. Acad. Sci. USA 91, 5033-5037.69. Rubinov M, Sporns O, (2010). Complex network measures of brain connectivity: Uses and inter-pretations. Neuroimage 52, 1059-1069.70. Kintali S, (2008). Betweenness Centrality: Algorithms and Lower Bounds. arXiv:0809.1906v2.71. Newman ME, (2006). Finding community structure in networks using the eigenvectors of matrices.Phys Rev E 74(3), 036104-19.72. Maslov s, Sneppen K, (2002). Specificity and stability in topology of protein networks. Science 296,910-913.73. Sugiyama K, Tagawa S, Toda M, (1981). Methods for visual understanding of hierarchical systemstructures, IEEE Transactions on Systems, Man, and Cybernetics SMC-11 (2), 109-125.74. Fruchterman TM, Reigold EM, (1991) Graph drawing by force directed placement. Software Prac-tice and Experience, 21:1129-1164.75. Lutz M (2010). Programming Python, Fourth Edition. OReally Media.9
Figure Legends
Figure 1. The proposed framework for the estimation of neuronal functional connectivity. (A) A recording session from thalamic and cortical regions. Arrows indicate the effective influenceamong neurons. The electrode tips record the neurons in dark red. (B) The firing patterns of thecortical neuron B produce common firing patterns both with neurons A and C but with different timedelays. In particular, C → B (red spikes) can be easily inferred by correlation analysis instead of A → B (blue spikes) hardly detectable. (C) Recorded signals are processed in overlapping windowslasting 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 thecompressed spike trains (C(A), C(B)) used by the N CS function. Whether
N CS ( A, B ) is greater thana fixed threshold then we can conclude that A → B . (E) An example of functional graph extracted byrecordings in the i th time window. (F) Typical sites of the rat paw for the tactile stimulation.0 Figure 2. Example of stimulus evoked redistribution of spiking functional connections.
Functional weights are redistributed from the pre- (A) to post-stimulus (B) configurations. Red, blueand green nodes indicate neurons respectively from VPL, S1 and CM. Some neurons, not functionallyconnected 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 withhigh degree) with many functional connections with VPL and CM thalamic neurons in (B) while is notinvolved 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.
Tables nodes edges win θ C C l C r L L r S ω ±
14 1117 ±
385 50 0.5 0 . ± .
08 0 . ± .
08 0 . ± .
08 1 . ± .
11 1 . ± .
11 1 . ± .
00 0 . ± . ±
14 947 ±
317 50 0.6 0 . ± .
09 0 . ± .
08 0 . ± .
08 1 . ± .
10 1 . ± .
09 1 . ± .
01 0 . ± . ±
14 1516 ±
216 250 0.3 0 . ± .
02 0 . ± .
01 0 . ± .
02 1 . ± .
14 1 . ± .
06 1 . ± . − . ± . ±
14 1311 ±
246 250 0.4 0 . ± .
04 0 . ± .
03 0 . ± .
03 1 . ± .
12 1 . ± .
07 1 . ± . − . ± . ±
14 1294 ±
228 500 0.3 0 . ± .
05 0 . ± .
04 0 . ± .
05 1 . ± .
38 1 . ± .
08 1 . ± .
04 0 . ± . ±
14 1176 ±
236 500 0.4 0 . ± .
04 0 . ± .
04 0 . ± .
06 1 . ± .
35 1 . ± .
08 1 . ± . − . ± . ±
14 1198 ±
170 1000 0.2 0 . ± .
01 0 . ± .
03 0 . ± .
06 1 . ± .
06 1 . ± .
05 1 . ± . − . ± . ±
14 876 ±
164 1000 0.3 0 . ± .
05 1 . ± .
04 0 . ± .
09 1 . ± .
08 1 . ± .
06 1 . ± .
06 0 . ± . Table 1.
Network statistics for spontaneous spiking activity. pre-stimulusnodes edges win θ C C l C r L L r S ω ±
17 889 ±
297 100 0.6 0 . ± .
08 0 . ± .
05 0 . ± .
04 1 . ± .
11 1 . ± .
09 1 . ± .
09 0 . ± . ±
17 1161 ±
364 100 0.5 0 . ± .
06 0 . ± .
04 0 . ± .
04 1 . ± .
12 1 . ± .
11 1 . ± .
03 0 . ± . ±
17 1356 ±
439 100 0.4 0 . ± .
04 0 . ± .
04 0 . ± . . ± .
13 1 . ± .
13 1 . ± .
01 0 . ± . ±
17 1409 ±
476 100 0.3 0 . ± .
04 0 . ± .
04 0 . ± .
04 1 . ± .
14 1 . ± .
14 1 . ± .
00 0 . ± . ±
17 897 ±
302 100 0.6 0 . ± .
01 0 . ± .
08 0 . ± .
09 1 . ± .
11 1 . ± .
09 1 . ± .
01 0 . ± . ±
17 1159 ±
391 100 0.5 0 . ± .
10 0 . ± .
09 0 . ± .
09 1 . ± .
12 1 . ± .
12 1 . ± .
00 0 . ± . ±
17 1371 ±
480 100 0.4 0 . ± .
09 0 . ± .
09 0 . ± .
09 1 . ± .
14 1 . ± .
14 1 . ± . − . ± . ±
17 1424 ±
517 100 0.3 0 . ± .
09 0 . ± .
09 0 . ± .
09 1 . ± .
16 1 . ± .
16 1 . ± .
00 0 . ± . Table 2.
Network statistics for evoked spiking activity. nodes edges win θ C C l C r L L r S ω
24 186 ±
60 50 0.5 0 . ± .
05 0 . ± .
05 0 . ± .
12 1 . ± .
242 1 . ± .
14 1 . ± .
00 0 . ± . ±
46 50 0.6 0 . ± .
09 0 . ± .
12 0 . ± .
12 2 . ± .
41 1 . ± .
03 1 . ± . − . ± . ±
37 250 0.4 0 . ± .
05 0 . ± .
04 0 . ± .
07 1 . ± .
25 1 . ± .
11 1 . ± . − . ± . ±
25 250 0.5 0 . ± .
09 0 . ± .
13 0 . ± .
08 1 . ± .
41 2 . ± .
29 1 . ± .
06 0 . ± . ±
32 500 0.3 0 . ± .
03 0 . ± .
03 0 . ± .
06 1 . ± .
01 1 . ± .
06 1 . ± . − . ± . ±
25 500 0.4 0 . ± .
05 0 . ± .
04 0 . ± .
06 1 . ± .
23 1 . ± .
08 1 . ± . − . ± . ±
12 1000 0.2 0 . ± .
07 0 . ± .
07 0 . ± .
05 1 . ± .
39 1 . ± .
07 3 . ± .
02 0 . ± . ±
11 1000 0.3 0 . ± .
06 0 . ± .
09 0 . ± .
05 1 . ± .
22 1 . ± .
13 4 . ± .
06 0 . ± . Table 3.
Network statistics for spontaneous LFP activity. pre-stimulusnodes edges win θ C C l C r L L r S ω
24 199 ±
48 50 0.4 0 . ± .
05 0 . ± .
04 0 . ± .
11 1 . ± .
11 1 . ± .
08 1 . ± .
02 0 . ± . ±
47 50 0.5 0 . ± .
06 0 . ± .
07 0 . ± .
12 1 . ± .
23 1 . ± .
01 1 . ± . − . ± . ±
43 100 0.3 0 . ± .
04 0 . ± .
37 0 . ± .
08 1 . ± .
09 1 . ± .
07 1 . ± .
01 0 . ± . ±
44 100 0.4 0 . ± .
05 0 . ± .
05 0 . ± .
10 1 . ± .
19 1 . ± .
11 1 . ± . − . ± . ±
45 50 0.4 0 . ± .
05 0 . ± .
04 0 . ± .
10 1 . ± .
11 1 . ± .
08 1 . ± .
02 0 . ± . ±
46 50 0.5 0 . ± .
05 0 . ± .
07 0 . ± .
11 1 . ± .
26 1 . ± .
14 1 . ± . − . ± . ±
40 100 0.3 0 . ± .
04 0 . ± .
03 0 . ± .
08 1 . ± .
08 1 . ± .
07 1 . ± .
02 0 . ± . ±
42 100 0.4 0 . ± .
05 0 . ± .
05 0 . ± .
10 1 . ± .
18 1 . ± .
29 1 . ± . − . ± . Table 4.
Network statistics for evoked LFP activity.2
Figure 3. Salient facts of functional graphs in spiking activity. (A) Correlation between thestandard 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 degreedistributions of functional graphs (spikes) extracted by pre- and post-stimulus conditions. (C) Averagebetweenness centrality balance over the three recorded regions (VPL, S1, CM) in both conditions.3
Figure 4. Example of LFP phases couplings. Functional connections are disposed on theLFP recording sites.
Green nodes represent CM channels, red nodes represent VPL channels andblue nodes represent S1 channels. (A) Before a tactile stimulation, LFPs are tightly coupled among theLFP channel of the same brain areas. (B) After an effective tactile stimulation, LFPs broke theirinter-site synchronies and established cross-site phase couplings.
Computational model composed by small-world networksnodes edges ρ C C l C r L L r S ω
100 1503 ±
211 0.01 0 . ± .
06 0 . ± .
04 0 . ± .
05 3 . ± .
34 2 . ± .
31 1 . ± . − . ± . ±
192 0.03 0 . ± .
02 0 . ± .
02 0 . ± .
02 1 . ± .
15 1 . ± .
15 1 . ± .
00 0 . ± . ±
181 0.05 0 . ± .
03 0 . ± .
03 0 . ± .
03 1 . ± .
10 1 . ± .
10 1 . ± .
00 0 . ± . ±
184 0.075 0 . ± .
03 0 . ± .
01 0 . ± .
05 1 . ± .
06 1 . ± .
05 3 . ± .
00 0 . ± . ±
156 0.1 0 . ± .
04 0 . ± .
04 0 . ± .
07 1 . ± .
12 1 . ± .
09 2 . ± .
12 0 . ± . ±
169 0.25 0 . ± .
04 0 . ± .
09 0 . ± .
07 2 . ± .
17 2 . ± .
15 3 . ± . − . ± . ±
205 0.5 0 . ± .
05 0 . ± .
07 0 . ± .
06 2 . ± .
17 2 . ± .
14 3 . ± . − . ± . ±
186 0.01 0 . ± .
11 0 . ± .
05 0 . ± .
10 1 . ± .
18 1 . ± .
13 0 . ± .
00 0 . ± . ±
193 0.03 0 . ± .
11 0 . ± .
07 0 . ± .
08 1 . ± .
23 1 . ± .
17 0 . ± .
08 0 . ± . ±
168 0.05 0 . ± .
09 0 . ± .
37 0 . ± .
08 2 . ± .
30 1 . ± .
21 1 . ± .
04 0 . ± . ±
179 0.075 0 . ± .
08 0 . ± .
05 0 . ± .
08 3 . ± .
64 2 . ± .
52 1 . ± .
04 0 . ± . ±
162 0.1 0 . ± .
09 0 . ± .
07 0 . ± .
08 2 . ± .
39 2 . ± .
31 1 . ± .
06 0 . ± . ±
181 0.25 0 . ± .
02 0 . ± .
37 0 . ± .
21 2 . ± .
33 2 . ± .
21 0 . ± .
01 0 . ± . ±
184 0.5 0 . ± .
07 0 . ± .
09 0 . ± .
13 2 . ± .
48 1 . ± .
19 0 . ± . − . ± . Table 5.
Network statistics for the computational model. ρ indicates the density (small-worldnetworks per node) and the winsize is set to 100 ms in all analyses.4 Figure 5. Salient facts of functional graphs in LFP. (A) Synchrony of LFP phases in eachrecording site during spontaneous, spike responsive and spike non-responsive configurations. Responsivestimuli increase the global and decrease the local synchronies. (B) Average node degree distributions offunctional graphs (LFPs) extracted by pre- and post-stimulus conditions . (C) Average betweennesscentrality balance over the three recorded regions (VPL, S1, CM) in both conditions.5
Figure 6. Basilar neurophysiological data. (A) Distributions of firing rate and (B) of Inter-spikeInterval for the representative neurons. (C) Pattern of tactile stimulations by Arduino microcontroller.6
Figure 7. Efficacy of Normalized Compression Similarity (NCS) to detect long-range spikeinteractions. (A) NCS and Pearson coefficient exerted on 100 couples of independent uniformlydistributed binary sequences (1000 bits). Both functions do not show bias. (B) To test the capacity ofNCS to detect significant interactions within time windows of 1000 ms, we fixed in the center of the firsttime-window a binary random pattern (50 ms long). The same pattern was replicated in a sequence of36 time-windows drifting it from the initial to the window end. (C) NCS and Pearson coefficient wereevaluated on these sequences: NCS is able to detect the interaction along the entire drifting processwhile Pearson coefficient is able to return significance only when the reference pattern is almost alignedin both sequences (drifts number 20-23).7