Stability of spontaneous, correlated activity in mouse auditory cortex.

Neural systems can be modeled as complex networks in which neural elements are represented as nodes linked to one another through structural or functional connections. The resulting network can be analyzed using mathematical tools from network science and graph theory to quantify the system's topological organization and to better understand its function. Here, we used two-photon calcium imaging to record spontaneous activity from the same set of cells in mouse auditory cortex over the course of several weeks. We reconstruct functional networks in which cells are linked to one another by edges weighted according to the correlation of their fluorescence traces. We show that the networks exhibit modular structure across multiple topological scales and that these multi-scale modules unfold as part of a hierarchy. We also show that, on average, network architecture becomes increasingly dissimilar over time, with similarity decaying monotonically with the distance (in time) between sessions. Finally, we show that a small fraction of cells maintain strongly-correlated activity over multiple days, forming a stable temporal core surrounded by a fluctuating and variable periphery. Our work indicates a framework for studying spontaneous activity measured by two-photon calcium imaging using computational methods and graphical models from network science. The methods are flexible and easily extended to additional datasets, opening the possibility of studying cellular level network organization of neural systems and how that organization is modulated by stimuli or altered in models of disease.

Introduction

Distributed and often redundant coding is a hallmark of neural systems [1], providing robustness to single-neuron variability [2] and supporting complexity in the system’s potential behavioral repertoire [3]. A key challenge in understanding this code lies in determining how the nature and strength of correlations between neurons is related to a stimulus [4]. Recent evidence suggests that so-called noise correlations have marked and diverse functions [5], from impacting information encoding and decoding [6-8], to tuning the amount of information present and thus the nature of ensuing cortical representations [9-12]. Correlations in spike trains have also been noted to contain important information about excitability, latency, and synchronization [13-15]. Even apart from task-evoked activity, spontaneous activity and correlations of that activity can profoundly impact cortical responses to a sensory input, thereby playing a critical role in information processing [16, 17].

To better understand the nature of coherent multi-unit interactions both during intrinsic and stimulus-induced processing, it is necessary to have a language in which to study inter-unit interaction patterns. In related work in other species and other spatial scales, network science has proven its utility as just such a candidate language [18]. The notion of a network in its simplest form is akin to the notion of a graph in the field of mathematics known as graph theory [19]. Specifically, an undirected binary graph is composed of nodes, which represent the units of the system, and edges, which link pairs of nodes according to some physical connection, functional relation, or shared feature [20]. This simplest version of a network can also be expanded to include weights on edges, weights on nodes, dynamics on edges, dynamics on nodes, or multiple types of nodes or edges forming a multilayer or multiplex structure [21, 22]. By either the simple or expanded encoding, network models of neural systems seek to distill the most salient organizational features of the system, allowing investigations to focus on how the network topology constrains or supports the system’s function [23]. Importantly, the network modeling approach is flexible in the sense that its components can be redefined at different spatial scales, and is thus equally applicable to cellular data at the microscale as it is to regional data at the large scale [24].

Recent studies have begun to build and characterize network models of cellular activity as measured by calcium imaging [25-31], and have demonstrated their biological relevance across a neural system’s development [32]. For example, one notable study provided initial evidence that immature cells in the developing brain display spontaneous correlation patterns that are characterized by small-world architecture and that critically regulate neural progenitor proliferation [33]. In a mature system (ferret visual cortex), recent evidence suggests that local connections in early cortical circuits can generate structured long-range network correlations that guide the formation of visually evoked distributed functional networks that display striking network modularity [34]. The architecture of correlations in spontaneous activity can be regulated by synaptotagmin [35], modulated by acetylcholine [27], blocked by glutamatergic antagonists [36], and mediated by a combination of intrinsic and circuit mechanisms [36]. Yet, little is known about the conservation or variation of network architecture in spontaneous correlations across different regions of the brain. Moreover, while the activity can be temporally quite precise in a given instance [36], little is known about how patterns of spontaneous activity change over the course of days and weeks after the critical period of development has passed. Understanding the principles of these dynamics is important for understanding the conserved rules that the architecture must obey, as well as the variability that can be exercised to meet the demands of the ever changing internal or external environments.

Here, we take steps to address some of these gaps in knowledge by measuring correlated spontaneous neuronal activity using two-photon calcium imaging, modeling those correlation patterns as networks, and assessing network architecture and dynamics over the course of several weeks. We focus our measurements specifically on mouse auditory cortex because of its rich organizational characteristics, with distributed representations of tone frequency [37], spatially overlapping locations for the representations of pitch and timbre [38], and the capacity for single neurons within the wider network to encode simultaneous stimuli by switching between activity patterns [39]. We choose mouse as our species of interest largely to prepare for future efforts using two-photon optogenetics [40] to perturb the network architecture, with the goal of probing network response to stimulation and validating recently posited theories of network control [41-43]. We begin by testing the hypothesis that networks reconstructed from fluorescence correlations exhibit hierarchical modular structure, and that network modules fluctuate over the timescales of days or weeks. We also test the hypothesis that some units participate in these temporal fluctuations more than others such that the system is best characterized by the existence of a stable temporal core surrounded by a fluctuating and variable periphery. Each of these hypotheses is motivated by prior observations in non-invasive imaging data acquired from humans [44-46], where evidence points to the importance of hierarchical modularity and temporal core-periphery structure for effective cognitive function [18, 45, 47–49]. Thus, collectively our hypotheses are predicated on the notion that neural systems are constrained to display some degree of preservation in network architecture across species, from human to mouse [43, 50, 51], as well as scale invariance, from the level of large-scale areas to the level of small-scale units [52-54].

Results

We recorded spontaneous activity from four awake, head-fixed mice over the course of 5, 5, 6, and 7 sessions spanning between 2 and 4 weeks. Specifically, we used two-photon microscopy to detect changes in fluorescence of GCaMP6s in transfected neurons caused by fluctuations in calcium activity. We estimated functional connectivity from the fluorescence traces using a cross-correlation of differenced activity for every pair of cells. We note that our analysis was carried out on continuous calcium traces rather than sparse spike trains, which can be obtained using deconvolution methods. We modeled the cell-to-cell correlation matrix as a network [23], and quantitatively characterized the network’s architecture using well-developed tools from network science [19]. Specifically, we assessed the modularity of the network structure using a commonly applied community detection technique known as modularity maximization [55-57]. Further, we assessed temporal fluctuations in this modular structure using tools for the analysis of dynamic graphs [21, 58, 59]. For further details on our methodological approach, see Materials and methods.

Networks exhibit multi-scale modular structure

One of the most important organizational principles of biological neural networks is their organization into cohesive modules [45]. These modules are thought to support specialized information processing while conferring robustness to perturbations. Moreover, converging evidence from micro- and macro-scale network analyses suggest that network modules are also organized hierarchically, with larger modules subtending broader brain function and smaller modules playing more specialized roles [24, 60, 61]. In this section, we test the hypothesis that networks reconstructed from fluorescence correlations in mouse auditory cortex exhibit hierarchically modular structure.

To address this hypothesis, we leverage recent advances in community detection methods [62]—a collection of algorithms and heuristics that use data-driven approaches to uncover the modular structure of networks. Specifically, we use an extension of the popular modularity maximization algorithm [55]. The standard version of this algorithm defines a module as a group of network nodes whose internal density of connections is maximally greater than what would be expected under a chance model. The extension of this algorithm samples modules over multiple organization scales, ranging from coarse divisions of the network into a few large modules to finer divisions of the network into many small modules. Importantly, unlike past applications, this extension also includes built-in null statistical testing capable of rejecting modular structure at different levels of the proposed hierarchy if they were consistent with a null model.

Here, we applied this approach to investigate the hierarchically modular structure of networks derived from correlated fluorescence traces (Fig 1a). The module detection method was applied separately to networks constructed from data in each recording session, which allowed us to take full advantage of all cells recorded on a given day. The algorithm resulted in a hierarchy of communities that passed statistical testing for significance (p < 0.05; Fig 1b–1f). In general, we found that the fluorescence networks exhibited hierarchical, multi-scale modular structure. Of the 43 recordings (aggregated across all mice) we observed statistically significant hierarchies in all. Across recording sessions, the average number of scales in a hierarchy was 19 (inter-quartile range of [10.25, 27.75]; Fig 1c).

Fig 1

Detection of hierarchical modular structure.

(a) Mean fluorescence of pixels, averaged over the full recording session. (b) Co-classification matrix generated using all statistically significant hierarchical levels. The dendrogram to the right depicts module splits. (c) The number of hierarchical levels aggregating data from all mice and all recording sessions. Panels (d), (e), and (f) depict module assignments at different levels of the hierarchy. The network diagrams shown in these panels are identical to one another and represent binarized matrices obtained by thresholding the jitter-adjusted correlation matrix of fluorescence traces between pairs of cells. Panels (g), (h), (i), and (j) depict distributions of z-scored mean intra-module Euclidean distance for each module and for each mouse. Panels (b), (d), (e), and (f) depict representative results from mouse 1. Here, the acronyms “IQR” and “Med.” represent interquartile range and median, respectively. Note also that nodes in panels b and d-f are ordered according to their hierarchical community assignments.

Additionally, we also computed spatial statistics for each module. Past studies have shown that communities tend to be spatially co-localized, so that other cells located near one another are more likely to belong to the same module compared to cells located far from one another [63]. To test whether this was also the case in our data, we computed the Euclidean distance from each cell to the nearest cell assigned to the same community. We then averaged this measure over all nodes in the same module. If cells were arranged in spatially dense, compact modules, then this measure would be small. Here, we calculated this measure for each module at every level of the hierarchy and compared these values against a null distribution generated by randomly and uniformly permuting the cell’s spatial locations but preserving their module assignments. For each module, we expressed the mean nearest-neighbor distance as a z-score with respect to this distribution. We found that the observed modules tended to be more spatially compact than expected by chance. For each mouse, the median z-score was less than zero and in all cases the inter-quartile range of z-scores excluded a value of zero (Fig 1g, 1h, 1i and 1j), indicating that the observed modules tended to be more spatially compact than expected by chance.

Lastly, we asked whether a module’s spatial features varied as a function of where it appeared in the hierarchy. To address this question, we aggregated z-scored nearest neighbor distance for each module along with modules’ hierarchical levels. A hierarchical level of 1 indicates large communities comprised of many nodes, while the deeper hierarchical levels refer to divisions of the network into smaller communities. Then, separately for each mouse, we computed the correlation of z-scored nearest neighbor distance with hierarchical level. We found negative association in three of the four mice, suggesting that deeper hierarchical levels, i.e. smaller communities, may be more spatially contiguous than larger communities. The Spearman correlation of hierarchical level and nearest neighbor distance was ρ = −0.22, −0.19, +0.01, and −0.27 for the four mice.

Network and module similarity decays over time

In the previous section we demonstrated that the correlation pattern of fluorescence traces exhibits modular structure across multiple scales, and that these multi-scale modules unfold as part of a hierarchy. In these analyses, the network’s modular structure was derived separately for each recording session. While this approach allowed us to characterize the modular structure on a given day, it tells us little about how those modules fluctuate over the timescales of days or weeks. Here, we address this question directly, taking advantage of the longitudinal tracking of cells across multiple recording sessions to assess the temporal consistency of the network’s overall organization, as reflected in the full correlation matrix, and in the network’s mesoscale organization, as reflected in its modular structure.

We begin by calculating the similarity between the correlation structure for any two recording sessions, u and v (Fig 2a). We first identified the set of cells from which fluorescence traces were recorded in both sessions. The average number of cells observed in any pair of recording sessions was 193.0±15.6 (Mouse 1), 132.7 ± 40.8 (Mouse 2), 121.7 ± 20.0 (Mouse 3), and 107.0 ± 35.6 (Mouse 4). We then extracted the correlation structure among those subsets of cells for each of the two recording sessions, resulting in two correlation matrices: W and W. Next, we vectorized the upper triangles of both matrices and computed the correlation of their elements, r. Finally, we expressed this correlation as a z-score, , with respect to a null distribution generated by randomly and uniformly permuting rows and columns of W and recomputing the correlation of with W (essentially the Mantel test [64]). Accordingly, large positive z-scores indicate that the correlation of r was much greater than expected in the non-parametric permutation-based null model. Aggregating z-scores across all pairs of recording sessions resulted in the z-scored similarity matrix, .

Fig 2

Reconfiguration of correlation structure over time.

(a) Analysis pipeline for comparing correlation structure. For any two correlation matrices, W and W, whose elements have been z-scored against those obtained under a “jittered” null model in which random offsets were added to timeseries (see Methods), we vectorize the upper triangular elements and compute their similarity using a Pearson correlation coefficient. We compare the observed correlation coefficient against that which we would expect under a null model in which rows and columns of W are permuted uniformly at random. In panels (b), (c), (d), and (e), we show the scatterplots of standardized similarity scores for pairs of correlation matrices with the number of days separating their respective recording sessions. In panels (f), (g), (h), and (i), each point represents the standardized similarity scores of module co-assignment matrices across pairs of recording sessions.

To assess the degree to which the similarity in correlation structure depended upon the time interval that separated the recordings, we also computed the distance matrix D = [d], which measures the distance (in number of days) between recording sessions u and v. We then compared the upper triangular elements of Z with the corresponding elements of D. In general, we observed that decayed monotonically as a function of d. Notably, this observation was consistent across all mice (mean±standard deviation Pearson correlation of ) (Fig 2b–2e). These findings indicate that the magnitude with which individual cells are correlated with one another over time varies systematically over recording sessions. Specifically, the correlation structures of recording sessions separated by a short period of time tend to be similar to one another, whereas those separated by longer periods of time tend to be dissimilar.

In addition to assessing whether cell-to-cell correlation patterns varied across recording sessions, we also aimed to assess the variability of modular structure. To address this question, we performed an analogous procedure to the one described above where we substitute the module co-assignment matrices C = [C] and C = [C] for the correlation matrices, W and W. Here, the element C indicates the fraction of all detected community partitions in which cells i and j were co-assigned to the same module. Otherwise, this procedure for relating the modular structure of networks from different recording sessions proceeded exactly as described above. We denote the z-score matrix from the module comparison as . As before, we observed that the correlation of similarity in modular structure decays with time (mean±standard deviation Pearson correlation of ), indicating the presence of marked quotidian variation (Fig 2f–2i). We also repeated this analysis after dividing communities into terciles according to their hierarchical level to assess whether there was a relationship between hierarchy and module stability. We found evidence of decaying similarity at all levels, but no clear dependence on hierarchy (S1 Fig).

Temporal core-periphery structure

In the previous section we demonstrated that, on average, hierarchical modular structure becomes increasingly dissimilar over time. However, it may be the case that some sets of brain areas maintain their modular structure despite the passage of time, forming a stable temporal core surrounded by a fluctuating and variable periphery [44]. To test this hypothesis, we focused on sequences of recording sessions and characterized the stability of modules across those sessions.

Because the number of recording sessions varied from one mouse to another, we focused on sequences of five recording sessions (the greatest number that was available for all mice). For each mouse, we modeled the thresholded connectivity data from each of these five recording sessions as the layers in a multi-layer network object [22, 65], and we used multi-layer modularity maximization [66] to track the fluctuations in modular structure across those five sessions (Fig 3a; see Materials and methods for more details). The multi-layer modularity maximization approach extends the traditional modularity maximization approach [55] by incorporating networks from all five sessions into a singular multi-layer network object, and then detecting modules in all layers simultaneously. The main advantage of this procedure is that the same set of module assignments are preserved across all layers, making it possible to track the formation and dissolution of modules over time and to seamlessly compare modules from one recording session to the next. Given such a mapping, one can then calculate measures like the local “network flexibility” [44, 67], which indicates how frequently a given node changes its module assignment across layers. Past studies have used this flexibility measure to identify temporally stable cores and variable peripheries (clusters of nodes with low and high flexibility, respectively) [44]. The multi-layer modularity maximization procedure required that we maintain a consistent set of nodes across layers, i.e. a fixed set of cells that were observed on all recording sessions. This conjunction results in a smaller set of cells than the pairwise comparisons utilized in the previous section. The number of cells retained for the multi-layer analysis was 144 (Mouse 1), 68 (Mouse 2), 67 (Mouse 3), and 54 (Mouse 4).

Fig 3

Estimation of core-periphery structure and network flexibility.

(a) Thresholded correlation matrices are separately treated as: a) layers in a multi-layer network, their communities estimated, and network flexibility estimated as the frequency with which a node changes its community assignment across layers; b) the consistency matrix is submitted to a core-periphery detection algorithm and each node’s “coreness” is estimated. Here, consistency measures the fraction of layers (recording sessions) in which a connection was present. (b) Node’s flexibility scores plotted in anatomical space. (c) Nodes’ “coreness” plotted in anatomical space. The size of nodes in panels b and c is proportional to their average weight across all five recording sessions. (d) Because flexibility is a measure of variability while “coreness” is a measure of stability, we find that the two are inversely correlated with one another (red line represents the identity line). (e) Cross-subject consistency of optimal parameters for fitting the core-periphery model. For each mouse, we calculated the difference between observed core quality and that of a null model, and we retained the top 10% of those points. These points are depicted at the level of individual mice in panel f. In panel e, we aggregate those values across all mice.

Obtaining estimates of network flexibility requires the detection of communities using multi-layer modularity maximization, which depends upon two parameters, γ and ω. These parameters control the resolution (size and number) of modules detected and their stability across layers, respectively. Here, we use a recently developed procedure that allows us to obtain a representative sample from the parameter space defined by these two variables [68]. For each such sample, we calculated a local (node-level) measure of flexibility, ranked the flexibility scores of all nodes, and subsequently averaged these ranked flexibility scores across all samples to generate an average flexibility profile for the population of cells. We note that each layer in this model represents the structure of a network estimated during different recording sessions. The multi-layer model makes it possible to aggregate these different networks into a single mathematical object rather than treating them as independent and disjoint estimates of the network [69]. This framework has been used widely in network neuroscience for modeling nervous systems whose network structure evolves over time [70, 71], differs across individuals [68], and spans multiple association modalities [72, 73].

In addition to the flexibility approach, we also used a second method to provide converging evidence of temporal core-periphery structure (Fig 3a) [74]. In this procedure, we calculated the session-averaged connectivity matrix (over the five recording sessions), and based on its organization we algorithmically assign cells (nodes) to a continuously defined core and periphery (see Materials and methods for more details). Intuitively, core nodes are nodes that maintain strong connections to one another and to the periphery across recording sessions, while peripheral nodes are those whose connections are variable (e.g., observed in only a few recording sessions or absent altogether). The size of the core and the smoothness of the transition from core to periphery are controlled by two free parameters, α and β. We systematically explored this parameter space and at each point, we fit the core-periphery model to the session-averaged network to calculated the core quality [74]. We compare the quality of cores fit to the observed session-averaged matrix against the qualities fit to random matrices generated by a permutation-based null model. This comparison allows us to identify points of interest in the parameter space: points where the observed core was of greater quality than that of the null model.

We note that, here, we focus on temporal core-periphery structure, where the core is comprised of nodes that maintain their connection to one another over time and the periphery is made up of nodes whose connections to one another are variable or infrequent [44]. This definition is distinct from core-periphery structure defined based on topological features alone [75, 76]. In the case of topological core-periphery structure, the core represents a strongly connected cluster of nodes that is weakly connected to a set of peripheral nodes. The peripheral nodes make few connections to one another.

In general, we found evidence that cortical activity in all mice exhibited temporally-stable cores of nodes that maintained community assignments and connectivity over many days. In general, the flexibility measure converged with the coreness measure, implicating roughly the same sets of nodes as temporally stable (i.e., manifesting high coreness, low flexibility) (Fig 3b–3d). Across mice, cores tended to be fairly exclusive (Fig 3e and 3f); core quality was maximally greater than the null model at points in parameter space corresponding to a small subset of cells. At these points in parameter space and across all mice, we found the average core included ≈ 20% of all nodes, an estimate that is obtained after thresholding the sigmoidal coreness measure at its midpoint. We note, however, that as a result of parametric variation we observed cores containing as few as 0.09% of all nodes and as many as 31.6% of nodes. The smoothness of the transition between core and periphery was more variable, suggesting that these networks may exhibit multiple cores with different degrees of smoothness separating the core from the periphery.

The detection of temporal core-periphery structure suggests that there exists a small subset of nodes whose modular organization is preserved across time. To test whether this was the case, we repeated the procedure from the previous section, wherein we calculated the similarity of cell-to-cell correlations across multiple days. Here, we calculated the similarity of connections among core nodes (the top 10% ranked by coreness) and connections among non-core nodes. If the core were indeed stable, we would expect that the connections within the core would be more stable across time compared to the connections within the non-core. Indeed, we found that this was the case. For mice 1–4, we found that the similarity of connections among core nodes were −0.26, +0.92, −0.16 and +0.06, compared to connections involving the non-core nodes, which were −0.29, −0.52, −0.33, −0.76. On average, the difference in correlation between core and non-core was 0.41±0.38. This observation suggests that the core nodes are more temporally stable than the periphery nodes, both in terms of their community structure and in terms of their connectivity patterns.

Discussion

Spontaneous fluctuations in neural activity at the cellular scale can modulate behavioral responses to incoming sensory stimuli [16, 17]. Yet the nature of that modulation is not well understood, in part due to the fact that such spontaneous activity does not appear to be random in nature, but instead displays heterogeneous dependencies or correlations among units. Little is known about the rules constraining the architecture of these correlations, or their variability over time. Here we sought to partially address this gap in knowledge by using recently developed techniques in network science to examine the network architecture of correlations in spontaneous activity in mouse auditory cortex as measured by two-photon microscopy and calcium imaging over the course of several weeks. We found that networks exhibited striking modular architecture, with smaller modules being located within larger modules in a multi-scale hierarchy. We also found significant temporal rearrangement of modular architecture, as indicated by the fact that the similarity in modules decreased monotonically as a function of the time interval between recording sessions, even when only considering those units that were present in both sessions. Finally, we found that the broadly observed temporal rearrangement of modules was complemented by the presence of a small number of cells whose modular allegiance remained stable throughout the 2–4 weeks of experimentation. We confirmed with additional testing that the co-existence of stable and unstable units was consistent with a temporal core-periphery model of system dynamics, where a stable core of units is accompanied by a flexibly periphery.

Multi-scale modular network structure

Biological systems generally and neural systems specifically, are frequently required to develop, adapt, and evolve in changing environments [77, 78]. This pervasive demand for adaptation is thought to be a partial explanation for the striking modular structure observed in biological systems [79, 80]. Each module is thought to have the capacity to change or adapt without adversely impacting the function of other modules. In neural systems, modules are thought to exist in order to segregate specific cognitive function or computations, allowing enhanced specialization of the organism [45, 81]. Such modular structure has also been observed in spontaneous recordings of intact zebrafish larvae, where topographically compact assemblies of functionally similar neurons reflect the tectal retinotopic map despite being independent of retinal drive [82, 83]. These data suggest that spontaneous activity displays modular structure that is a functional adaptation specifically tuned to support the system’s behavior. Similar observations have also recently been made in ferret visual cortex, where widespread modular correlation patterns in spontaneous activity accurately predict the local structure of visually evoked orientation columns several millimeters away [34].

Hierarchical modularity in biological systems is further thought to allow for a decomposability of the system’s temporal responses to the environment, with fast processes occurring in small modules at a low level of the hierarchy and slow processes occurring in large modules at a high level of the hierarchy [84]. Prior work at the large-scale has demonstrated the presence of hierarchically modular structure in neural systems specifically, and suggested that large modules support broad cognitive functions while small modules support specialized cognitive functions [24, 60, 61, 85]. Here we extend these prior observations by showing that over short time periods approximately equal to the duration of a recording session, neurons assemble into cohesive modules of varying size, ranging from large, spatially-distributed clusters of weakly coupled neurons to compact, highly correlated ensembles. In theoretical work, it is interesting to note that hierarchical modularity provides an efficient solution to the problem of evolving adaptable systems while minimizing the cost of connections [86]. This relation between hierarchical modularity and low cost yet efficient information processing in neural systems has also been supported by both theoretical work and analysis of neural data in both C. elegans and human [85]. When considering our results in this light, it is useful to note that the spatial compactedness of modules suggests that maintaining long-distance correlated activity may be metabolically costly and therefore uncommon. Overall, these findings are consistent with those observed in other micro- and macro-scale networks and suggest that the organizational principles of modular architecture and spatially-compact, low-cost clusters may be conserved across spatial scales [45, 87].

Daily variation in network architecture and module constituency

Accompanying the nascent use of tools from network science to understand interaction or connection patterns between neural units, there has been a marked interest in understanding the dynamics of interaction patterns as a function of time, and across a variety of different time scales [21, 58, 88, 89]. Particularly in the human imaging literature, efforts have begun to understand principles of dynamic network reconfiguration on the time scale of minutes or hours [70, 90], days [91, 92], weeks [93], months [46, 94], and years [95, 96]. Here we exercise that interest in the domain of network models of correlation matrices derived from spontaneous activity in mouse auditory cortex over 2 to 4 weeks of experimentation. Our findings suggest that quotidian variation in correlation structure is manifest at multiple scales: (i) at the level of cell-to-cell correlations, but also (ii) at the level of large-scale and module patterns in the network. This latter observation is particularly interesting to consider in light of findings at the macro-scale level of whole-brain networks derived from fMRI data. Specifically, at this large scale, much of the modular organization of spontaneous correlations in the human brain is conserved across the time scales of days and weeks, with notable flexibility largely present at module boundaries. One could speculate that gross temporal stability in macro-scale networks is underpinned by notable micro-scale variability. It would be interesting in future to more directly address the question of the functional role of this micro-scale network reconfiguration, and specifically to test the hypothesis that the correlation structure of fluorescence traces in mouse primary auditory cortex is reorganized over timescales of days to weeks to support cortical functional reorganization.

A stable network core accompanied by a flexible network periphery

In other natural dynamical systems, it has been noted that density tends to support temporal stability, while sparsity tends to support temporal instability [97]. In the context of networked systems, the notion can be expanded to describe the phenomenon in which a core of densely interconnected units tends to display weak or slow temporal fluctuations, while a periphery of sparsely interconnected units tends to display strong or fast temporal fluctuations [44]. In the context of the human brain, this temporal core-periphery structure has been raised as a model for the balanced constraints of task-general processes, implemented by the temporal core, and task-specific processes, implemented by the temporal periphery [98]. It is interesting to consider whether such a delineation into temporal core and periphery is also characteristic of cellular networks, and whether that separation is functionally meaningful in a similar sense. Our findings suggest that, while calcium fluorescence correlation structure changes markedly over time, there remains a relatively small set of cells whose interactions, both as single connections but also as communities, are spared and preserved. There is some evidence in theoretical studies that such core-like structures emerge early in development, and are strengthened through functional activation [99]. In analyses of macro-scale networks, core stability and peripheral flexibility have been associated with learning [44], leading us to speculate that the emergence of core-periphery structure in micro-scale networks may serve a similar role in preserving learned (auditory) relationships, while maintaining enough variability to learn and map novel stimuli. Thus, future work could be directed to investigate the functional roles of cores and peripheries during task conditions.

Our observations are not without precedent and are in line with other previous cellular-level studies. For instance, in a study of Tritonia motor programs, [100] found that in addition to a “core” set of neurons that consistently responded to stimulation, the motor network expanded to include a peripheral set of neurons whose membership varied across trials. The authors went on to show that the variable set of neurons included those that were strongly coupled to the network at rest, suggesting that spontaneous coupling patterns predispose cells to be recruited into the motor program. Similarly, [101] distinguished between sets of cells as “soloists” and “choristers,” which respectively comprise those whose activity was distinct from that of the rest of the population and those whose activity was consistent with that of the bulk of the population. Our finding of temporal core-periphery structure echoes both of these studies. In our study, however, we find analogous cores whose community affinity is preserved over days. Although the functional relevance of the core-periphery structure reported here remains unclear, large-scale studies of brain networks have linked reorganization of these types of structures to learning [44]. Outside of neuroscience, cores are thought to represent polyfunctional units in a network, positioned where modules overlap with one another [102], and to occupy positions of influence within the broader topology [103]. Future studies could better disambiguate the functional role of temporal core-periphery structure by linking those types of structures to behavior.

Methodological considerations and limitations

There are several methodological considerations and limitations that are pertinent to the interpretation and generalizability of our results. First, we note that the experimental methods allow us to sample only a subset of neurons within a specific “slice” of the auditory cortex. It is likely that most of the neurons that directly target the neurons that we image are not captured by the analysis. Therefore, the estimates for the network connectivity should not be taken as an approximation for the actual physical connectivity in the cortical circuit. Another important aspect of data collection is that we focus on a specific cortical layer: layer 2/3. Neurons in the cortex differ tremendously in their connectivity patterns across different layers [104, 105]. It would be important in future studies to sample the activity across cortical depth to better understand integration of information across cortex.

Second, we note that we have examined correlations in spontaneous activity fluorescence traces, and this approach has the strengths of computational simplicity and ease of interpretation [106]. However, we acknowledge that correlation-based approaches focus on pairwise functional interactions, and remaining agnostic to underlying structural connectivity as well as to higher-order (non-pairwise) relations between units. It would be interesting in future to consider maximum entropy models as an alternative method to estimate connections between units [107], both for its sensitivity to underlying structure [108], and for its ability to assess higher-order interactions [109]. Approaches that could then take advantage of the richer assessment of higher order interactions in these data include emerging tools from algebraic topology [110, 111], which have already proven relevant for understanding structure-function relationships at both large and small scales in neural systems [112, 113].

Finally, an additional limitation concerns the measures used to establish the presence or absence of connections between cells. Specifically, we constructed networks where nodes represented individual cells and where edges represented the correlation magnitude of fluorescence traces. Importantly, correlated activity is not a direct proxy for underlying structural connectivity [114], and thus a pair of neurons that may not be directly synaptically connected can exhibit correlated activity; rather than reflecting structural connections, functional connectivity provides information about the interactions between neurons due to their function [18]. Moreover, correlated activity also does not represent the coupling matrix that prescribes the temporal evolution of brain activity [115]. Rather, the correlation structure of neural activity represents the product of a dynamical system whose evolution is constrained by structural connections. Though correlated activity at the large-scale has proven useful for investigating the functional organization of brain networks [30, 116, 117], its utility for understanding and characterizing the structure and function of micro-scale networks remains unclear and largely untested [118]. Future studies should both investigate in greater detail the relative advantages of alternative, domain-specific measures of functional connectivity [119, 120] and the relationship of these measures to other connection modalities [121].

A related limitation concerns the properties of network connections or edges. Here, we used fully weighted and signed representations to calculate network similarity (see Fig 2) and sparse, binary networks to detect modules and core-periphery structure (See Figs 1 and 3). We adopted these particular edge definitions because they allowed us to compute similarity using the full distribution of edge weights and to avoid making future, somewhat arbitrary, decisions about the appropriate null model for weighted and signed networks in modularity maximization [122]. Nonetheless, there exists a spectrum of possible ways to define edges and their properties. Future work is needed to understand the tradeoffs between different definitions and to ultimately ground edge definition in biophysical properties of nervous systems.

Conclusion

Across many scientific disciplines from plant biology [123] to biogeodynamics [124] and the study of biodiversity [125], scientists are faced with the challenge of bridging two or more scales of investigation into the function of complex systems. For example, in evolutionary biology, a key challenge is to bridge physical scales from protein sequences to fitness of organisms and populations [126], while in the study of cancer progression a key challenge is to map genotype to phenotype [127]. Neuroscience is no exception. Ongoing efforts seek to bridge the gap between the connectome and the transcriptome [128], between brains and social groups [129], or between large-scale brain regions and small-scale cellular circuitry [130]. In each case, the development of a formal understanding will depend upon the capacity to build mathematical descriptions and theories across scales. One natural approach to this challenge is to use a formalism that is scale invariant, a characteristic that makes network science particularly appealing. Our work in this study is an example of considering tools and conceptual paradigms previously exercised at the large-scale of brain regions, and exercising them at the level of cellular circuitry. We look forward to future efforts explicitly measuring and examining the network architecture of neural systems across both of these scales simultaneously in the same animal, with the goal of better understanding and predicting behavior.

Methods

Ethics statement

All experimental procedures were performed in accordance with NIH guidelines and approved by the IACUC at the University of Pennsylvania.

Animals

All experiments were performed with equal numbers of adult male and female mice (supplier—Jackson Laboratories; age, 12-22 weeks; weight, 20-36 g; PV-Cre mice, strain: B6;129P2-Pvalbtm1(cre)Arbr/J).

Two-photon microscopy and calcium imaging

Four mice were implanted with cranial windows over auditory cortex. Briefly, the mice were anaesthetized with 1.5-3% isoflurane and a 3mm circular craniotomy was performed over auditory cortex (stereotaxic coordinates) using a 3mm biopsy punch. An adeno-associated virus (AAV) vector encoding the calcium indicator GCaMP6s (AAV1-SYN-GCAMP6s, UPENN vector core) was injected for expression in layer 2/3 neurons in left A1 within the window (750 nl, 1.89 × 10-12 genome copies per ml) [131]. After injection a glass circular 3mm coverslip (size 0, Warner Instruments) was placed in the craniotomy and fixed in place using a mix of Krazy glue and dental cement. A custom-made stainless steel head-plate (eMachine Shop) was fixed to the skull using C&B Metabond dental cement (Parkell). All imaging sessions were carried out inside a single-walled acoustic isolation booth (Industrial Acoustics) as previously described. Mice were placed in the imaging setup, and the headpost was secured to a custom base (eMachine Shop) serving to immobilize the head. Mice were gradually habituated to the apparatus over 3 days, 3-4 weeks after surgery.

Using two-photon microscopy (Ultima in vivo multiphoton microscope, Bruker) changes in fluorescence of GCaMP6s in transfected neurons caused by fluctuations in calcium activity were recorded in awake, head-fixed mice. We recorded from the same cells over many days in layer 2/3 of auditory cortex, using blood vessel architecture, depth from the surface, and the shape of cells to return to the same imaging site. Laser power at the brain surface was kept below 30 mW. Chronic imaging of the same field of view across days was carried out for the duration of the experiment.

Recordings were made at 512×512 pixels and 13-bit resolution at approximately 30 frames per second. Spontaneous activity was recorded for 10 minutes in each session. Publicly available toolboxes [132] were used to register the resulting images, select regions of interest, estimate neuropil contamination, and extract the changes in fluorescence from each cell. Upon conclusion of the imaging sessions, brains were extracted following perfusion in 0.01M phosphate buffer pH 7.4 (PBS) and 4% paraformaldehyde (PFA), post-fixed in PFA overnight and cryopreserved in 30% sucrose solution for 2 days prior to slicing. The location and spread of GCaMP6s was confirmed through fluorescent imaging. These methods are consistent with the recommendations of the American Veterinary Medical Association (AVMA) Guidelines on Euthanasia.

Cell tracking over days

To identify ROIs from different imaging sessions that correspond to the same cell we performed a multi-step routine: (1) The mean fluorescence images from each day were registered by transforming the coordinates of landmarks present in both images in MATLAB (2017a) using the fitgeotrans function. The resulting transformation was used to transform the ROIs from the second imaging session to match the first; all sessions were aligned to the first imaging session. (2) We calculated the distance between all pairs of centroids across the two sessions. For each ROI from session 2, we computed the percentage overlap of the 10 cells with the smallest centroid distances from session 1. Cells that had more than 1 ROI with higher than 20% overlap were manually inspected; the ROI that matched the current cell was selected from the overlapping ROIs or none were selected if it was unclear whether they were the same cell. A good match was determined by considering the percent overlap and the shape of the ROIs. All other cells were assigned the closest ROI as matching. (3) We manually inspected any cells that had duplicate matching ROIs; again considering the shape and the percent overlap of matching ROIs, we selected the ROI that matched that cell, or we decided that none was a good match. (4) To check for false positive matches, we manually inspected any matches where the centroid distances were greater than the mean + 1 std of all the matches, and any matches that had less than 30% ROI overlap; the match was deemed good or not depending on the match criteria. Neurons that were not matched to any ROI were counted as different or new and assigned a new cell number. This process was repeated for subsequent sessions, registering the imaging field to the first session and comparing the ROIs to the cumulative ROIs from previous sessions. A final manual inspection of all the unique ROIs was performed after all the imaging sessions were registered; ROIs that overlapped were excluded from the dataset since it was unclear whether they were the same or different cells. Examples of tracked cells and aligned ROIs are shown in S2 Fig.

Network reconstruction

We estimated functional connectivity from fluorescence traces. Let x(t) indicate the intensity of fluorescence in cell i at time t. Next, we computed the cross-correlation of differenced fluorescence traces for every pair of cells: where μ and σ are the mean and standard deviation of the differenced time series, respectively. Here, “differenced activity” refers to a preprocessing strategy sometimes used to reduce serial (auto) correlation in time series data and to make a time series approximately stationary (i.e. constant mean over time) [133]. Given a time series of length T defined as x = [x1, …, x], the differenced time series has length T − 1 and is calculated as x = [x2 − x1, …, x − x]. This difference transform is particularly appropriate when networks are constructed using Pearson correlations, where time series are assumed to be stationary and to contain uncorrelated samples.

To reduce the likelihood that the observed correlations were driven by chance fluctuations, we “jittered” the time series of cells by adding or subtracting < 1 second offsets to each cell independently. We then computed jittered cross correlations, . We repeated this procedure 1000 times. More specifically, the jittering procedure generates surrogate time series by iterating cell by cell, selecting a random starting index, p, that occurs within the first 2 seconds (the first 60 samples given the 1 second offset and 30 Hz sampling rate), and retaining the next T − 2p − 1 samples. The resulting time series have length T − 2p. The jittered time series underwent the same differencing procedure as the original time series.

We estimated for every pair of cells the probability that the jittering procedure would generate a correlation as strong as that which was observed empirically, and we made binary connections between those cells with p < 0.05. This procedure resulted in a sparse matrix, with elements A ∈ [0, 1]. We note that the fully-weighted and signed matrices were used in the analysis of the day-to-day network similarity. The thresholded, binary, and sparse matrices were used for community and core-periphery detection due to the complexity of the algorithms.

Module detection

We used modularity maximization to detect network modules based on connectivity data [55]. This method aims to divide network nodes (cells) into modules whose internal density of connections is maximally greater than what would be expected under a null model. This intuition is formalized by the modularity quality function [134]:

In this equation, k = ∑ A is the degree of node i and 2m = ∑ k is the total number of connections in the network. The term gives the expected number of connections between node i and node j given the null model in which each node’s degree is preserved but connections are formed at random. The resolution parameter, γ, scales the relative contribution of the null model. The module assignment of node i is encoded as g and δ(g, g) is the Kronecker delta, whose value is equal to unity when g = g and is zero otherwise.

In this manuscript, we used two variants of modularity maximization. First, we studied the network community structure for each recording session independently. For this analysis, we combined modularity maximization with a newly-developed multi-resolution technique that divides the network into communities of different sizes (scales) that are related to one another hierarchically [135]. This procedure allows us to examine community structure across a range of scales, from large communities to smaller communities that might support more specialized information processing.

Additionally, we used a multi-layer variant of modularity maximization that makes it possible to track the evolution, formation, and dissolution of communities across recording sessions [66]. In this procedure, the standard modularity maximization equation is modified to read: Here, the subscript s denotes network layers, u ∈ {1, …, T}. That is, each network layer corresponds to connectivity patterns estimated on a different day. So A represents the presence or absence and weight of the connection between node i and node j in layer u. Similarly, k = ∑ A is the degree of node i in layer u and g is the community to which node i is assigned in layer u. Unique to the multi-layer variant of modularity maximization is the inter-layer coupling parameter ω, which links node i to itself across layers. From the perspective of maximizing Q, non-zero values of ω make it advantageous to group node i into the same community across layers. When ω is small, the advantage is correspondingly small, and the detected communities emphasize the unique community structure of layers. On the other hand, when ω is large, the detected communities are consistent across layers and emphasize shared features of community structure.

Here, we used a recently-developed procedure to obtain estimates of community structure with the values of {γ, ω} sampled from a restricted parameter space [68]. This procedure involved first estimating the boundaries of a restricted parameter space wherein any {γ, ω} pair would result in community structure where the number of communities is > 1 and < N × T (where T is the total number of layers; T = 6, in this case), and where community structure is neither uniform across layers (flexibility of exactly 0) nor is it maximally dissimilar (flexibility of exactly 1). See [24] for more details on how these boundaries were estimated. We then sampled 10000 {γ, ω} pairs from within this parameter space and for each sample we maximized the corresponding Q(γ, ω). All subsequent analyses were carried out on these detected communities.

The principle advantage of the multi-layer formulation is that it estimates communities for all layers simultaneously and preserves nodes’ community labels across layers. This advantage makes it possible to directly compare the community assignment of a given node in layer u and in layer v ≠ u, and to identify nodes whose community assignments are flexible (varying from one layer to another) or inflexible (remaining in the same community across layers). We can quantify this intuition using the network measure flexibility [67, 70]: Intuitively, flexibility counts the fraction of times that nodes’ community assignments in layers u and u + 1 differ. Nodes that differ more frequently have flexibility values closer to 1, while nodes that differ less frequently have flexibility values closer to 0. Here, we used the flexibility measure as an index of change in network community structure across recording sessions.

Co-assignment probability

Based on the 10,000 partitions generated using single-layer modularity maximization, we were able to estimate the probability that nodes i and j were assigned to the same module. Let c denote the community assignment of node i in partition h ∈ 1, …, 10000. The co-assignment probability of nodes i and j was calculated as: The value of C ranges from 0 to 1, indicating that nodes i and j were never assigned to the same module or were always assigned to the same module, respectively.

Core-periphery detection

Separately, we also characterized the stability of network organization across recording sessions by computing a temporal core and periphery. In this context, a core refers to a group of nodes that are densely internally connected and to the periphery, which is weakly internally connected [103]. To identify temporal core-periphery structure, we first generated a connection consistency matrix, whose element G = ∑ A represented the fraction of layers (recording sessions) in which a network connection was present. In this matrix, a core refers to a group of nodes whose connections are maintained across time, while the periphery is a set of nodes whose connections are more variable.

We used a variant of a common core-periphery definition in which the transition from core to periphery varies smoothly (non-binary). We begin by defining the N × 1 vector C of non-negative elements [74]. Given this vector, we then defined the matrix C = CC subject to the constraint that ∑ C = 1. The values in the vector C are permutations of the vector: The coreness of each node is the permutation of that maximizes the core quality function: This method introduces two free parameters, α ∈ [0, 1] and β ∈ [0, 1]. The value of α determines the sharpness of the core-periphery boundary. With α = 1, the transition is binary while the transition with α = 0 is maximally fuzzy. Similarly, the value of β determines the size of the core; as β ranges from 0 to 1, the size of the core varies from N to 0. In our application, we performed a grid search of 31 linearly-spaced values of α and β, using a simulated annealing algorithm to maximize R (with 10 restarts).

Module similarity at different hierarchical levels.

(TIF)

Click here for additional data file.

Examples of cell tracking over days.

(TIF)

Click here for additional data file.