## Abstract

Mapping the flow of activity through neocortical microcircuits provides key insights into the underlying circuit architecture. Using a comparative analysis we determined the extent to which the dynamics of microcircuits in mouse primary somatosensory barrel field (S1BF) and auditory (A1) neocortex generalize. We imaged the simultaneous dynamics of up to 1126 neurons spanning multiple columns and layers using high-speed multiphoton imaging. The temporal progression and reliability of reactivation of circuit events in both regions suggested common underlying cortical design features. We used circuit activity flow to generate functional connectivity maps, or graphs, to test the microcircuit hypothesis within a functional framework. S1BF and A1 present a useful test of the postulate as both regions map sensory input anatomically, but each area appears organized according to different design principles. We projected the functional topologies into anatomical space and found benchmarks of organization that had been previously described using physiology and anatomical methods, consistent with a close mapping between anatomy and functional dynamics. By comparing graphs representing activity flow we found that each region is similarly organized as highlighted by hallmarks of small world, scale free, and hierarchical modular topologies. Models of prototypical functional circuits from each area of cortex were sufficient to recapitulate experimentally observed circuit activity. Convergence to common behavior by these models was accomplished using preferential attachment to scale from an auditory up to a somatosensory circuit. These functional data imply that the microcircuit hypothesis be framed as scalable principles of neocortical circuit design.

## Introduction

Computation in mammalian neocortex relies on specific circuits comprising individual neurons and the connections between them. Given the myriad functions that can be assigned to different regions of the brain, it is unclear whether circuitry is generalized across multiple regions of neocortex. Because all regions must perform similar basic tasks under the same biophysical constraints (Douglas et al., 1989; von Melchner et al., 2000), the cortex may use a general circuit design as described by the microcircuit hypothesis (Mountcastle, 1957; Szentágothai, 1978; Douglas et al., 1989). System level studies have provided data consistent with the postulate showing that primary sensory cortices process other modalities (Kayser et al., 2005) and are capable of taking on a primary processing role of a different modality following experimental manipulation (von Melchner et al., 2000). It is clear that microcircuitry in the neocortex is structured (Song et al., 2005; Yoshimura et al., 2005; Perin et al., 2011; Levy and Reyes, 2012); however, it is unknown how this structure manifests functionally particularly at the larger mesoscale throughout the neocortex. Elucidating the organization of functional circuitry (Gerstein et al., 1978) will provide key insights into the flow of activity through local neocortical circuitry, the underlying circuit architecture, and also has the potential to provide insight into the computational strategies used in each respective cortical region (Watts and Strogatz, 1998; Alon, 2007).

We imaged the flow of activity at the at the mesoscale level, which spans multiple columns and layers, to generate functional wiring diagrams in two areas of sensory neocortex. The lack of experimentally defined benchmarks to characterize functional microcircuitry necessitated a novel approach that would allow us to identify which statistical features of activity flow were informative. By increasing the field of view, we maximized the number of neurons imaged and the statistical power to investigate neocortical circuit dynamics. Moreover, we were able to evaluate the role, if any, of traditional anatomical boundaries in shaping the flow of activity and in turn the functional circuitry. We chose a comparative methodology (Kätzel et al., 2010; Yang and Zador, 2012) to examine the microcircuit postulate by comparing functional wiring diagrams generated from primary auditory (A1) and somatosensory barrel field (S1BF) neocortex. These two regions are an interesting test of the microcircuit hypothesis as both map sensory input anatomically and display temporally structured circuit activity (Luczak et al., 2007; Montemurro et al., 2007), but each area appears organized according to different design principles. A1 can be considered a one-dimensional tonotopic mapping of the cochlea along the rostrocaudal axis (Bandyopadhyay et al., 2010; Oviedo et al., 2010; Rothschild et al., 2010; Levy and Reyes, 2012), whereas S1BF provides a two-dimensional mapping along both the rostrocaudal and dorsoventral axes corresponding to the spatial location of the whiskers, manifested in a clear columnar organization containing barrels (Woolsey and Van der Loos, 1970; Welker, 1976; Simons, 1997; Lefort et al., 2009). Additionally, laminar cell-type composition and thalamic projections may differ slightly between these regions (Barbour and Callaway, 2008). Given that these areas map sensory information in anatomically distinct ways, similarities in emergent circuit activity would reflect common cortical organization, whereas differences would highlight the role of the distinct architecture for each region.

## Materials and Methods

##### Preparation of calcium dye-loaded slices.

C57BL/6 strain mice of either sex on postnatal day 14–17 were anesthetized by intraperitoneal injection of ketamine-xylazine, rapidly decapitated, and had their brains removed and placed in oxygenated ice-cold “cut” artificial CSF (ACSF; contents contain the following, in mm: 3 KCl, 26 NaHCO3, 1 NaH_{2}PO_{4}, 0.5 CaCl2, 3.5 MgSO_{4} 25 dextrose, 123 sucrose). Coronal slices (500 μm thick) containing the sensory region of interest was cut perpendicular to the pial surface using a vibratome (VT1000S; Leica). In a subset of experiments, alternate coronal brain slices with thalamocortical connectivity intact were cut (450 μm thick S1, 500 μm thick A1) at angles as previously described (Agmon and Connors, 1991; Cruikshank et al., 2002). Slices were placed in a 35°C oxygenated incubation fluid (Incu-ACSF; contents contain the following, in mm: 123 NaCl, 3 KCl, 26 NaHCO3, 1 NaH_{2}PO_{4}, 2 CaCl2, 6 MgSO_{4}, 25 dextrose) for 30–45 min. Calcium dye loading was then achieved by placing all slices into a small Petri dish containing ∼2 ml of Incu-ACSF, an aliquot of 50 μg Fura-2AM (Invitrogen) in 13 μl DMSO and 2 μl of Pluronic F-127 (Invitrogen) as previously described (Sadovsky et al., 2011). All procedures were performed in accordance and approved by the Institutional Animal Care and Use Committee at the University of Chicago.

##### Electrophysiology and calcium dye imaging.

Experimentation was performed in standard ACSF (contents contain the following, in mm: 123 NaCl, 3 KCl, 26 NaHCO3, 1 NaH_{2}PO_{4}, 2 CaCl2, 2 MgSO_{4}, and 25 dextrose, which was continuously aerated with 95% O_{2}, 5% CO_{2).} Whole-cell current-clamp recordings were made using Multiclamp 700B amplifiers (Molecular Devices). We required a stable voltage value recorded during the down state, for the duration of the experiment for inclusion of the physiology. Rapid whole-field imaging of Fura-2AM loaded neurons was achieved by taking multiple 5 min movies using the Heuristically Optimal Path Scanning technique and microscopy setup as previously detailed (Sadovsky et al., 2011), allowing us to monitor action potential generation within individual neurons. Our dwell time parameter for each experiment was fixed at a value between 16 and 20 samples/cell/frame for each experiment. Greater than or equal to four events were necessary for a field of view to be included in our dataset.

##### Laminar identification.

We used biotinylated NeuN staining along with biocytin filled neurons which acted as fiduciary markers, in combination with measures of distance from pia and bright-field, NeuN, and two-photon cell density to identify lamina.

##### Statistical analysis.

All statistical analyses were performed with MATLAB (MathWorks). Unless otherwise noted, data are presented as mean ± SD. For nonparametric distribution comparison between the two sensory regions, and unless otherwise noted, the Wilcoxon rank sum test was implemented via the MATLAB “ranksum” function. If more than two medians were being compared, as in the case of lamina, the Kruskal–Wallis test was used and noted in text. A two-tailed *t* test, noted in text, was used to compare normal distributions. A one-tailed *t* test, noted at use, was used to compare graphical topological data when we had an a priori expectation of the experimental mean being larger than that from null datasets. For tests of significance, α = 0.05 was used as the cutoff.

##### Spike and circuit event detection.

Spikes were inferred from the calcium traces of individual neurons using a modified version of fast-oopsi (Vogelstein et al., 2010; Sadovsky et al., 2011). Spikes from each cell's calcium trace were then identified and circuit events were defined as regions where the network of cells was active for at least 500 ms.

##### Spike timing precision.

We established statistical significance of temporal stereotypy for each cell in a field of view by comparing all the individual spike trains of a cell obtained from every circuit event. Cells had to be active in at least four events to be considered for analysis. To compare spike trains, we used the D^{spike} Victor metric (Victor and Purpura, 1996) with a time parameter (*q*) of 1 s. This comparison gave us a distance value corresponding to the amount each spike train needed to be modified to be identical across activations. Distance is larger when the temporal similarity between spike trains is decreased. To see whether these distance value were significant compared with what would be expected by chance of a neuron firing at a rate similar to that of the network, we created a null comparison of 5000 shuffled spike trains using an inhomogeneous Poisson process with a rate that was equal to the average network firing rate across events. We then compared our observed metric value to that of this shuffled population to obtain a *p* value (see Fig. 5*A*). Cells showing *p* < 0.05 were considered significant.

##### Circuit area.

Circuit area was determined by using a convex-hull approach analogous to the smallest convex region that can be formed that contains all active cells.

##### Fuzzy clustering.

By flattening temporal dynamics, we created binary representations of each circuit event with 1 s indicating which cells were active at any point during the event and 0 s representing cells not active in that event. Fuzzy-c means that clustering was achieved using the MATLAB function “fcm” with specified numbers of clusters *N* ranging from 2 to the total number of events for a single slice. This function returns the membership function matrix for all events indicating how strongly each event belongs to each specified *N* cluster. For our analysis, cluster sufficiency was defined as all events having at least one cluster membership larger than 1/*N* + (1/*N*)/4. This indicated that an event was nonambiguously placed into a single fuzzy cluster. To obtain the number of clusters necessary to explain all the events in a region, we iteratively ran the fuzzy clustering method with increasing values of *N*, looking for the point in which all events fell under the definition of being sufficiently clustered. Because fuzzy clustering is dependent on initial seed, we took the average output of 100 runs of the fcm method, with each run consisting of either a maximum of 100 iterations or a clustering improvement of ≥0.00001.

##### Hierarchical clustering.

Binary representations of each circuit event were created in the same manner as used in fuzzy clustering. From these time-independent flattened events, we calculated the Jaccardian distance between every pair of these binary representations to create a pairwise distance matrix where E1 = event 1, E2 = event 2:
We then clustered the resulting pairwise distance matrix, *D*, to build a binary, hierarchical clustering tree using the unweighted pair group method with arithmetic mean (UPGMA; each cluster *A*,*B* is made up of events *x* in *A* and *y* in *B*):
Then we apply a cutoff and count clusters.

A precise reliability index was made for each cell in the form of:
where participation was the fraction of times that cell was active at least once during imaged circuit events and temporal significance refers to the Victor significance *p* value for that cell.

##### Columnar/laminar flow.

For each individual circuit event, directional flow between lamina was defined by a field of view where a statistically significant correlation existed between time frames and distance from the pial surface. Intercolumnar flow was determined similarly according to distance to an arbitrary line perpendicular to the pial surface not in the field of view. Events with significant correlations in both cases were considered to have both types of flow. Events with no significant correlation in either case were considered disperse.

##### Spatial clusters similar flow analysis.

Cluster IDs for each dataset were determined using the fuzzy clustering method with *N* equal to the grand average cluster amount. For each dataset, circuits were separated into their functionally clustered spatial groups. For each group which was activated two or more times, we evaluated what percentage of those events showed similar laminar flow. The same was done with columnar flow. Permuted clusters were then made by permuting the cluster membership between circuit events, thus maintaining the total number of clusters and distribution of similar events falling into a cluster, but permuting the actual membership of each event. The same analysis was run and resulting distributions were compared.

##### Null connectivity distribution.

Our null connectivity distribution was created by looking at all the pairwise distances between 750 uniform random distributed centroids falling within a 1 mm diameter circle.

##### Graph metrics.

Graph metrics were obtained using the Brain Connectivity Toolbox for MATLAB (http://www.brain-connectivity-toolbox.net/). Modularity was performed using the toolbox's implementation of the measure presented by Leicht and Newman (2008). Some graph figures were generated using the open source graph visualization platform Gephi (https://gephi.org/).

##### Modeling.

Modeling was implemented using the NeMo (http://www.doc.ic.ac.uk/∼akf/nemo/) neural network spiking simulator using CUDA and multicore processing. Previous work in the laboratory had demonstrated the need for spike rate adapting neurons to achieve dynamics similar to experimental measures therefore networks were comprised of Izhikevich neurons (http://www.izhikevich.org/publications/spikes.htm) with 80% being excitatory (*a* = 0.02, *b* = 0.2, *c* = −65, *d* = 6) and 20% being inhibitory (*a* = 0.02, *b* = 0.25, *c* = −60, *d* = 2).

##### General network modeling.

Weights were directly drawn from experimentally generated functional adjacency matrix weights, and synaptically scaled from 0 to 5 (0 to −5 for inhibitory cells). Topology was also directly based off adjacency matrix weights for connections, yet models were not embedded into a two-dimensional anatomical space. Adjacency matrices were scanned for cells with zero degrees (no connections) and removed from the matrix. An input layer containing 50 strongly excitatory neurons (synapse strength of 10 units) was randomly connected to the network with a 0.05 probability to other excitatory cells and 0.02 probability to inhibitory cells. Each cell of the 50 neuron input layer was activated over two time frames according to a uniform random distribution, with no network input cell firing more than once. Each simulation contained a random assignment of which cells were excitatory, inhibitory, and connected to the input layer.

##### Inferred circuit model.

An adjacency matrix from each field of view was created as described in the text: edges were weighted as the fraction of observed single frame lag correlation values between two cells spiking. For each of the adjacency matrices created from data, 1000 network simulations were performed with different random seeds.

##### Idealized circuit model.

Size distributions for each area were created using a positive kernel smoothed density (MATLAB “ksdensity”) of observed data values from imaging. The lower tail of this distribution (<200) was truncated to values of 200. The distribution inferred from each area was then sampled uniformly five times and summed to size *N*, respective to the area modeled. An adjacency matrix of size *N* was then created. An out degree per node distribution (D) was created using positive kernel smoothed density with values from data normalized to the size of respective datasets. Additionally, weight distributions (W) were formed via a kernel smoothed density of the weights of respective datasets. To fill in the adjacency matrix, for each cell in the adjacency matrix we would obtain the out degree per node value according to uniform random sampling from the above distribution D. The out degree for a cell was then calculated by multiplying with the size *N* above and rounding. Then for each out degree, an out synapse was added to a cell sampled random uniformly without replacement with a weight randomly drawn from W. Simulations were performed as described above with the caveat that delays were randomly selected to be twice the duration of inferred circuits to account for the larger simulated volume of neurons. As above, all simulations with the idealized model were performed 1000 times with different random seeds.

##### Scaling circuit model.

Adjacency matrices were created by obtaining a size from A1 and S1. The A1 size would be the base size (B) and S1 would be the size to scale to (S). Initial weights/degrees were formed as done for A1 matrices above. For S1 preferential scaling, for the degree difference S-B, new cells were added using S1 weights/degrees according to a preferential attachment model (Barabási and Albert, 1999). One-thousand circuits for each region were created and simulated to create dynamics distributions.

## Results

To determine whether A1 and S1BF exhibited common functional circuit topology, we monitored neuronal activity in slices of mouse neocortex using high speed multiphoton laser calcium imaging (Vogelstein et al., 2010; Sadovsky et al., 2011). We controlled for the potential influence of slice angle and confirmed primary sensory location by using both coronal and thalamocortical slices (see Materials and Methods). Within each sensory region, data from both slice angles were indistinguishable and therefore pooled (S1 *p* = 0.0894, A1 *p* = 0.2212). We imaged the flow of activity through large populations of neurons at the mesoscale in a two-dimensional circular imaging plane with a diameter of 1.1 mm that comprised multiple layers and columns with single-cell resolution (Fig. 1*A*; Movies 1 and 2 ; intervals between events removed). This allowed us to examine the impact of anatomical boundaries on functional circuitry and also resolve the general statistical features of functional circuit topologies with greater accuracy as they are maximized with increased sample size (Milo et al., 2002). Because temporal resolution of multiphoton microscopy is compromised at these spatial scales, we used the heuristically optimized path scan technique (Sadovsky et al., 2011; Fig. 1*B*), which allowed us to achieve fast frame rates (12.1 ± 2.4 Hz) that did not differ between regions (*p* = 0.3211). Both regions of sensory cortex showed a common capacity for emergent, multineuronal activity, classified by discrete periods of correlated action potential generation within subsets of neurons. We refer to these distinct clusters of spontaneous action potentials as individual circuit events. Imaging conducted simultaneously with whole-cell patch-clamp recordings demonstrated that imaged multineuronal firing is coincident with single-cell depolarization, termed an upstate, in active neurons (Cossart et al., 2003; MacLean et al., 2005; Fig. 1*C*,*D*; *n* = 42 upstates). As previously reported (Cossart et al., 2003; Shu et al., 2003), the occurrence of spontaneous multineuronal activity required intact excitatory synaptic transmission. Application of 20 μm AMPA blocker CNQX and 50 μm NMDA blocker APV eliminated the occurrence of all emergent circuit events (*n* = 3), indicating that these events were not intrinsic to individual cells but rather required synaptic connectivity. Each imaging plane contained up to 1126 neurons, with A1 having significantly less neurons 595 ± 101, and S1BF 704 ± 157 in the imaged field of view (Fig. 2*A*; *p* = 0.0305). The majority of active cells in both cortices had a continuously detectable signal over 1 h (Fig. 2*A*) and we found that within 10 circuit events we had detected a majority of the neurons that would be active during the experiment (Fig. 2*B*). Using this approach we captured the full extent of emergent circuit events within a single field of view, which allowed comparison of spatiotemporal activity and functional circuitry across sensory cortices.

#### Similar circuit dynamics in A1 and S1BF

We analyzed spontaneous circuit dynamics in each region of neocortex to determine whether there were common features of neocortical circuits in terms of their network size and duration. We evaluated 350 datasets, comprising a sample of 21338 individual neurons spread over the two regions of sensory cortex. We contrasted the duration of each circuit event, defined as the time from the first observed spike to the last in each observation. S1BF exhibited significantly longer imaged circuit event durations (1568 ± 885 ms, *n* = 268) than A1 (1203 ± 456 ms, *n* = 82, *p* = 0.0159). Events in S1BF involved a significantly greater number of active neurons per circuit event (207 ± 123; max 577) compared with A1 (155 ± 62; max 294, *p* = 0.0040). A positive correlation existed between the number of neurons involved and the duration of an event (Fig. 3*A*; *r* = 0.6852). The number of spikes in an event also had a positive correlation with duration (*r* = 0.7451). Thus, the longest events contained the largest number of active cells and detected action potentials. This linear relationship was conserved across both areas of sensory cortex, consistent with emergent dynamics being the result of synaptic connectivity and indicating that each area showed a range of functional circuit sizes. For each imaged dataset, our data collection was restricted to a single two-dimensional plane of focus that, despite our large sample size, we expected to be insufficient to fully sample the active three-dimensional circuit. We used patch-clamp physiology and recordings of upstates as an independent measure to confirm that circuit events were longer in duration in S1BF than in A1. Upstate duration reflects the synaptic drive from the three-dimensional active circuit impinging upon the patch-clamped neuron (Hasenstaub et al., 2005) and as a result is a description of the duration of the circuit event sampled from a region larger than a two-dimensional plane. We found that the duration of the upstate was significantly longer in S1BF compared with A1 (A1 2.14 ± 0.48 s, S1BF 3.20 ± 0.66 s, *p* = 8.8514 × 10^{−4}) consistent with the significant rank order measured from the imaging data. To compare the temporal progression of each circuit event we discretized imaging resolved action potential firing into 100 ms time bins and normalized cellular activity to the maximum number of active cells in a region. This analysis revealed that the time course of neuronal recruitment was similar in both regions of neocortex. In each circuit event, activity began in a small subset of neurons, quickly rose to a peak as more cells are recruited into the circuit event and then exponentially decayed (τ_{A1} = 641 ms_{,} τ_{S1BF} = 799 ms; Fig. 3*B*). This pattern of circuit progression was present regardless of the total number of neurons recruited. Overall circuit duration correlated with the time to peak in these events (*r* = 0.8347). The number of coactive cells firing during circuit events, i.e., the functional circuit size, had a long tailed distribution across all regions (Fig. 3*C*), consistent with scale free dynamics (Beggs and Plenz, 2003). We examined whether these results could be the product of similar activity scaled to overall cell count in each slice. Together, these data indicated that functional circuits in both areas of sensory cortex fall along a continuum of numerical scaling which in turn dictated the duration and the temporal progression of the dynamics of each functional circuit.

#### Spatial distribution of functional circuits in sensory neocortex

By bounding the area encompassed by any one functional circuit (see Materials and Methods), we found that the area encompassed by neurons that were recruited into a circuit event spanned a large portion of our field of view (0.519 ± 0.1 mm^{2}; Fig. 3*D*). We normalized these active circuit areas by the largest possible convex hull size for each field of view and found that these active regions did not significantly differ in size between cortical areas (*p* = 0.1029). We observed functional circuits that overlapped and shared neurons between events in both areas. We evaluated how many circuits were present in any given field of view by observing many individual circuit activations over time and determining the contribution of neurons to one or many circuits. We applied fuzzy and hard clustering algorithms (see Materials and Methods) to these time flattened events to isolate the most robust spatial patterns. Using the fuzzy clustering of observed events, the median number of circuits differed across regions (*p* = 3.7104 × 10^{−33}; Fig. 3*E*,*F*). A1 contained the smallest number of functional circuits (3.05 ± 0.283 circuits) and S1BF the largest (4.47 ± 0.502; Fig. 3*F*). Results were similar using a hard clustering technique (see Materials and Methods), indicating that this result was independent of the clustering algorithm used. We examined the spatial extent of each functional circuit by reprojecting circuit events into anatomical space. We identified lamina, and anatomical columns in S1BF (Lefort et al., 2009; Meyer et al., 2011), in 69% of our imaged field of views (see Materials and Methods; Fig. 4*A*). In this dataset, each individual lamina was sampled equally between the two sensory regions (L1 *p* = 0.9405, L2/3 *p* = 0.1261, L4 *p* = 0.3705, L5 *p* = 0.8813). In concordance with our convex hull analysis, the full spatial extent of functional circuits spanned columns and layers. To evaluate whether anatomical boundaries constrain the flow of spontaneous circuit dynamics we categorized the time course of each circuit event as having a dominant direction of flow in relation to the pial surface, discretized as columnar (perpendicular to pia), laminar (parallel with pia), both columnar and laminar, or disperse flow (Fig. 4*B*). Flows of all types were found in both regions, characterized by a majority of disperse activity (Fig. 4*C*). Thus dynamics were often unconstrained by anatomical boundaries, even in the case of the clear anatomically defined columns in S1BF. We determined whether there was a bias for the initiation of emergent activity to originate in specific lamina. The earliest detected spiking activity in neurons was distributed across lamina with no clear significant difference toward neurons of any layer starting events (L1, L2/3, L4, L5 Kruskal-Wallis *p* = 0.1386). Finally, we determined whether the spatially clustered functional circuits were more likely to have more similar circuit flow than expected by chance. We found that both the presence or absence of laminar flow patterns and columnar flow patterns were more often shared within spatial clustered data than permuted cluster memberships (see Materials and Methods; laminar flow within-cluster data = 80% similarity, permuted-cluster = 39%; *p* = 2.5955 × 10^{−7}; columnar flow within-cluster data = 85%, permuted-cluster = 50%; *p* = 1.6160 × 10^{−4}). These data suggested that while functional circuit activity was not dictated by gross anatomical boundaries, similar spatially coactive groups of cells display similar activity flow.

#### Temporal activity is structured in both regions of sensory cortex

After establishing spatial structure and structured circuit activity flow patterns in both regions, we tested for the presence of temporally precise neuronal firing sequences across sensory cortex. Spatiotemporal structure is consistent with spontaneous activity being dictated by underlying synaptic connectivity (Luczak and MacLean, 2012). We compared single-neuron firing activity from circuit events in each field of view using a metric where increasing dissimilarity between circuit reactivation spike patterns is associated with increased cost (Victor and Purpura, 1996; Fig. 5*A*; see Materials and Methods). In both regions of sensory cortex, circuit reactivations demonstrated precision in the spike timing of a subset of individual cells (38 ± 11% of cells in A1, 31 ± 13% in S1BF, *p* = 0.4433; Fig. 5*B*). Cells that demonstrated patterned spike timing were found in all lamina (32% of cells in L1, 33% L2/3, 37% L4, 43% L5). Further, the same cells often spiked in multiple circuit activations (27.9% ± 7.5 of activations averaged across all cells in A1, 26.0 ± 08.9% S1BF, *p* = 0.9052), with some cells being active in every event (5.2 ± 5.7% of cells grand average in A1, 2.0 ± 3.2% in S1BF *p* = 0.0214). To evaluate how often the same cells fired precisely over multiple circuit events, we combined our measure of multiple circuit participation with the above defined metric of temporal precision (see Materials and Methods), with values of 1 representing temporally reliable cells present in every event and 0 representing cells that are either nonprecise, rarely firing or both. A1 had the highest reliable precision (0.58 ± 0.28), which was significantly different from S1BF (0.47 ± 0.27; *p* = 4.415 × 10^{−22}). Overall, indicators of temporal precision, firing reliability, and nonrandom circuit activation, were found in both A1 and S1BF.

#### Graph analysis of inferred functional connectivity

We used the temporal flow of activity to evaluate the statistical features of functional circuit organization in A1 and S1BF. We determined whether the imaged flow of activity through local circuits was randomly organized or whether there the activity flow reflected functional organizing principles that might be common to both sensory cortices. Graph theory provided a mathematical framework and established metrics for describing high dimensional networks (Bullmore and Sporns, 2009). Therefore, to compare our two sensory regions, we generated graphical abstractions, i.e., circuit topologies, corresponding to observed functional activity. The sparse firing of circuit events necessitated the use of correlative methods for inference of topology (Pajevic and Plenz, 2009). Consistently active neurons across events (>4 spikes) were represented as nodes in each graph. Edges between nodes were directional and formed according to a single frame lagged correlation (frame duration 86.43 ± 18.62 ms; Fig. 6*A*). Edges were weighted according to reliability of lagged activity normalized to the number of events in that field of view; their weight represented the reliability of a connection. Although it is clear that a functional relationship between neurons increases the probability of them having an effective physical connection (Ko et al., 2011, 2013), and functional structures, such as hubs, have been shown to have corresponding structural correlates (Honey et al., 2007; Bonifazi et al., 2009), it is unlikely that each functional edge, or connection, corresponds to an effective, or literal, connection (Gerstein et al., 1978). Rather, given our method of inference, our functional connectivity measure captured the flow of activity through the network during a circuit event.

We characterized the general graph properties of each functional circuit topology. Functional connectivity was sparse in both regions of neocortex and the ratio of edges to nodes did not differ between regions (*p* = 0.4508; Fig. 6*B*). We examined the strength and arrangement of connections between neurons in graphs from each region of sensory cortex. In both areas there was a distribution of strong and weak weights, reflecting both weakly and highly reliable connections (Fig. 6*C*). In both areas, the degree distribution of high weight edges, i.e., the reliable functional connections between neurons, demonstrated linear regimes on a log-log scale and scaled similarly between regions with a bias toward small degree nodes in contrast to Erdős–Rényi randomized graphs that contained the mean number of edges and nodes from an A1(shown; ER-A1) or S1BF functional graph (Fig. 6*D*). We examined the numerical size of the functional circuits and found significant differences in the sizes of the graphs. A1 graphs contained 595 ± 101 neurons and S1BF 704 ± 157 (*p* = 0.0305; Fig. 6*E*), with graphs containing only neurons with edges of weight >0.1 containing 341 ± 90 neurons in A1 and 474 ± 200 neurons in S1BF. The large size of these graphs,compared with the relatively small sizes of individual circuit activations (A1 155 ± 62 cells, S1BF 207 ± 123 cells), indicated that the full graph was the product of the connectivity inferred across many distinct observations. A1, which contained the numerically smallest circuits, also had the lowest probability for neurons with a large number of connections. We examined in and out degrees, respectively the inward and outward connections of each node, normalized by total graph size so that degree per node ranged from 0 to 1 (Fig. 6*F*). These distributions suggested that A1's lower probability of high degree nodes was linked to the size of its graph; when the distribution was scaled by the overall size, A1 closely matched the distribution of S1BF. We determined whether both areas of sensory cortex contained highly interconnected neurons, hubs, a hallmark of scale free topologies which are considered optimal for limiting the impact of the loss of any randomly selected node (Albert et al., 2000) and for coordinating activity between neurons (Bonifazi et al., 2009). Here, we defined hubs as neurons whose degree was at least 1 SD higher than the mean degree (Honey et al., 2007; Fig. 6*G*). Hub neurons did not differ in their prevalence in both areas of sensory cortex (*p* = 0.4387), regardless of whether we defined hubs using 2 (*p* = 0.2569) or 3 SDs (*p* = 0.2362) beyond the mean. These results demonstrated that the functional connectivity of both regions of sensory cortex consisted of sparsely connected neurons and a minority of highly connected hubs. Moreover, consistent with results reported by effective connectivity and theoretical studies, circuits were marked by a small number of strong reliable links among a majority of weak connections (Teramae et al., 2012).

We measured modularity in our directed networks (Leicht and Newman, 2008) to determine whether circuits were defined by one homogenous topology, marked by low modularity, or characterized by multiple interconnected modules. All of the inferred functional graphs exhibited significantly higher modularity than expected from randomly shuffled graphs in which we preserved the degree distributions (one-tailed *t* test, *p* = 0.0300). Modularity measures did not differ between brain regions (*p* = 0.3615). We measured the mean clustering coefficient of each network (Watts and Strogatz, 1998), which is a measure of how often neighboring neurons are connected to one another's neighbors. This metric provided an indicator of the small-world architecture combined with a measure of the characteristic path length through a network. The clustering coefficient was significantly greater (one-tailed *t* test, *p* = 1.316 × 10^{−15}; Fig. 7*A*) and shortest path length significantly smaller (one-tailed *t* test, *p* = 2.482 × 10^{−13}) in the experimentally derived graphs than in degree-preserved random graphs, which indicated small world topologies in both regions of sensory cortex. We tested for hierarchical modularity, a topology comprised of modules of modules that has been observed in metabolic networks and large scale fMRI data, and has been shown to allow for increased information transfer without increasing wiring cost (Ravasz et al., 2002; Gallos et al., 2012). This was achieved by evaluating the relationship between the clustering coefficient of a neuron and its number of connections. When we isolated highly reliable edges in the network (weights ≥ 0.4; degree distribution is long tailed with mean in-degree = 9 ± 11 connections; max in-degree = 117) we find a relationship where the clustering of a node, C(*k*), with *k* links followed the scaling law C(*k*) ∼ *k*^{−1}, consistent with a hierarchically modular circuit organization in S1BF and A1 (Fig. 7*B*). We determined whether there were any additional organizational principles of connectivity. As previously reported by studies that used paired patch-clamp recordings (Song et al., 2005; Perin et al., 2011), a greater number of reciprocal, or bidirectional, functional connections exists than would be expected by chance, assuming an independent probability of connection (A1 4.51 ± 2.11-fold increase, S1BF 4.48 ± 3.70; Fig. 7*C*). Together these graph metrics indicated that the functional topology of each region consisted of hierarchically organized cells that shared neighbors. Thus, despite differences in the duration and overall size of circuit events between A1 and S1BF, each sensory cortical region was made up of circuits that had similar statistically defined emergent functional architectures.

#### Projection of inferred graph topology into anatomical space

To study the inferred functional connectivity graphs in anatomical space, we examined the distribution of connections over physical neuron distance (Fig. 7*D*). Similar to benchmarks established using paired patch-clamp recordings and anatomy (Perin et al., 2011; Song et al., 2005; Hill et al., 2012), functional connections were biased toward short pairwise physical distances. (Fig. 7*D*, red line; see Materials and Methods). Additionally, we examined whether the hub neurons that we identified were over-represented in a specific anatomical location. We found neurons with hub properties in all lamina but with an overrepresentation in L1, L2/3, and L5 (5.0%, 2.5%, and 3.6% over expected) and an under-representation in L4 (−11.15%). Our data demonstrated that functional circuits were biased toward spatially local activity with high degree connections being spread throughout the lamina.

#### Circuit models are sufficient to replicate experimental circuit dynamics

Theoretical studies have shown a direct link between the connectivity in network models and the resultant dynamics (Honey et al., 2007; Galán, 2008; Roxin et al., 2011; Litwin-Kumar and Doiron, 2012). We evaluated whether these inferred functional topologies were capable of recapitulating experimentally observed dynamics in a spiking neuronal network (SNN) model. It was unclear whether these topological descriptions of the most common and reliable circuits based on function would be able to capture any of the original emergent dynamics of actual anatomical circuits. For each experiment we used experimentally defined functional connectivity as the weighted connectivity matrix in a SNN (see Materials and Methods). These network models produced mean firing dynamics that followed a similar temporal progression of neuronal recruitment to that measured experimentally with similar decay (τ_{A1-Model} = 485 ms_{,} τ_{S1BF-Model} = 878 ms; Fig. 8*A*; see Materials and Methods). Specifically the modeled network dynamics correlated with key features of the experimentally observed dynamics in terms of total cells active (*n* = 32000 simulations, *r* = 0.9740; Fig. 8*B*), total spikes (*r* = 0.9310) and duration (*r* = 0.3934). These similarities demonstrated that functional topological connectivity, inferred from reliable time lagged circuit events, was sufficient to replicate key aspects of the imaged dynamics.

We manipulated the features of our network models to determine whether the agreement between the models and experimental data were the result of an overly parameterized model (Prinz et al., 2004). In addition these manipulations allowed us to determine what aspects of the inferred connectivity were crucial to replicate experimentally observed dynamics. In line with previous work (Teramae et al., 2012), the combination of both weak and strong connections was necessary for our model to function well. Running the model without low weight edges (weights ≤ 0.2) caused mean duration of activity to decrease to 63.7% compared with experimentally generated topologies. Correspondingly, the number of active cells decreased to 36.0% and total spike count fell to 9.5% compared with models that had been generated from experimental data, indicating that a full range of weights was necessary to support firing activity (Fig. 8*C*, purple line). We manipulated experimentally generated topologies by maintaining the respective out degree distribution of each node, but permuting the edges so that individual neurons maintained the number of outgoing connections but not which cells they connected to. This decreased all activity, as indicated by a lower value of mean active cells (61.8%) and decreased total firing (69.3%), yet duration of circuit events increased (109.6%). Thus, the specific identity of connections was necessary for the model to replicate the observed dynamics. Finally, we created a null model of circuit topology by constructing edges via a random draw of all available weights. These networks had the same total degree and weight as our inferred networks, but the balance of in/out degree distributions and total weight of each node was changed. These circuit topologies demonstrated dynamics that were increased to 109.9% of mean duration of models from experimental data, yet exhibited decreased mean total spike counts (67.3%), whereas the number of active cells dropped to 62.14% (Fig. 8*C*, red line), again changing the observed relationship between active cells and duration. In both cases, these perturbations showed increased durations at the expense of activation of the entire network and also lower overall spike rates. We hypothesize that this is due to an increase in the inefficiency of the perturbed networks, where groups of cells uniquely coactivate one another, which appears to perpetuate firing within subgroups but limits the flow of activity across the functional network as a whole. Thus, we were only able to observe dynamics that corresponded to experimental data in cases in which we maintained the observed degree distributions, weight distributions, and connectivity of observed circuits demonstrating both sufficiency and necessity of the functional topologies in these models.

#### Inverse model reveals that scaling is capable of making regions similar to one another

Having found that functional circuit topologies were sufficient to replicate experimentally measured dynamics in our SNN models, we determined which graph properties were responsible for the observed differences in the observed duration between A1 and S1BF by swapping graph metrics between SNN models of each brain region. Instead of using models built according to directly observed topologies as above, we generated probabilistic models of circuits of each neocortical area using draws from distributions of size, degree, and weight from the empirically defined data distributions obtained in each specific area (see Materials and Methods). These models were capable of recapitulating the rank order of firing duration seen in the experimental data (Fig. 8*D*, left) providing a substrate upon which we were able to swap model parameters to test hybrid circuits and determine which region specific graph metric distribution accounted for observed differences. Using the idealized model, duration was significantly different between A1 and S1BF (two-tailed *t* test, *p* = 7.0647 × 10^{−212}), whereas idealized S1BF circuits using A1 sizes had similar durations as idealized A1 datasets (two-tailed *t* test, *p* = 0.1088) and idealized A1 circuits using S1BF sizes had similar durations as idealized S1BF datasets (two-tailed *t* test, *p* = 0.6031). For both A1 and S1BF, manipulations of numerical circuit size using the idealized model were sufficient to transform the firing duration of S1BF to that of A1 and vice versa.

To better understand the role that scaling plays in setting the unique observed duration of each region, we used the smallest circuit, A1, as a starting seed to evaluate whether we were able to construct a circuit that exhibits S1BF firing duration through scaling. We increased the size of A1 using two approaches. Nodes were added using either random attachment, according to a uniform probability, or by using preferential attachment where new nodes are more likely to connect to nodes that are already highly connected (Barabási and Albert, 1999; see Materials and Methods; Fig. 8*D*). Only preferential attachment of S1BF sizes resulted in circuits that captured the firing duration of S1BF networks (Fig. 8*E*).

## Discussion

Using both experimental and theoretical approaches we elucidated key topological features of functional circuitry in two areas of sensory neocortex (Table 1). Both areas of sensory cortex demonstrate emergent dynamics marked by stereotypy in both individual cell firing times as well as the set of cells active in a single event. We find that although the flow of dynamics can be canonical, that is to say perpendicular to pia, this is not the dominant observation. These data are consistent with functional circuits that are anatomically interdigitated reflecting neuronal populations that represent information from multiple octaves in A1 (Stiebler et al., 1997) and whiskers in S1BF (Hirata and Castro-Alamancos, 2008). The relationship between the duration of circuit activity, the number of neurons active, and the progression of neuronal recruitment, shows a conserved size-dependent correlation within and across sensory cortices consistent with common architectural principles. Functional circuits are structured and are significantly different from random graphs, as they contain hallmarks of scale free, small world, and hierarchically modular topologies. We replicated experimentally observed temporal firing features using prototypical circuits from each region; consistent with the hypothesis that functional circuitry provides key insights into underlying circuit architecture. Simulations using prototypical circuits of S1BF and A1 showed that scaling from the smallest circuit architecture of A1 via preferential attachment reproduced the circuit dynamics of S1BF. We conclude that S1BF and A1 sensory neocortex show common principles of microcircuit architecture, and that scaling can account for a large portion of the similarities and differences between regions.

A principle of scaling is seen throughout vertebrate CNSs, regardless of the spatial resolution at which it is examined: many brain regions scale allometrically with brain volume (Stevens, 2001), reflected in a scaled relationship between gray and white matter (Zhang and Sejnowski, 2000). Similarly, resulting dynamics are scale invariant (Beggs and Plenz, 2003). Accordant with these anatomical and modeling studies, we find a scaled relationship within and between functional circuits in A1 and S1BF. The number of neurons active in a circuit and the resulting duration of observed spontaneous circuit events are strongly correlated in both areas. Numerical scaling also accounts for many of the differences that we observe in the functional features of circuitry and dynamics between both brain regions. Interestingly, we find that S1BF functional circuits can be grown from the smaller A1 functional circuit through preferential attachment. We suggest that there are two potential, nonmutually exclusive, driving forces which underlie the scaled relationship that we see in both areas. The first is neuronal resources. If there are fewer neurons available due to differences in the brain volume occupied by a sensory area and/or differences in the density of neurons present in that area of cortex (Herculano-Houzel et al., 2008; Collins et al., 2010), the number of neurons available to comprise a basic computational unit will be less in that area of sensory cortex. Alternatively the number of neurons interconnected, as reflected by the functional circuit, could be the result of synaptic plasticity driving different circuit sizes due to the unique statistical features of each modality (Clopath et al., 2010).

Anatomical connectivity in the brain sets bounds on the dynamics which can occur. Rules clearly govern connectivity between neocortical neurons according to location and cell class (Silberberg et al., 2002; Thomson et al., 2002; Song et al., 2005; Lefort et al., 2009; Perin et al., 2011; Hill et al., 2012) and each brain region is a layered structure (Hubel and Wiesel, 1974; Rockel et al., 1980) that contains the same general classes of neurons (Silberberg et al., 2002). Our functional data indicate that at the mesoscale that there are similar architectural principles that govern functional neocortical circuitry that are likely modified to be unique in each region (Kätzel et al., 2010; Yang and Zador, 2012). We find that the flow of activity indicates a number of statistical features which generalize across both regions of sensory cortex and that wiring is not random and can be distilled to a few basic rules: (1) Circuit architectures result in emergent dynamics that are spatiotemporally structured and show a progression of neuronal recruitment. (2) As previously reported in effective connectivity studies (Song et al., 2005; Perin et al., 2011) we find that the functional relationships between neurons are dictated in part by spatial proximity, although there is a long tail in the probability distribution of a functional connection which likely reflects the presence of long connections spanning hundreds of microns (Gilbert and Wiesel, 1979; Kätzel et al., 2010; Oviedo et al., 2010) that far exceed estimates of the lateral size of a column (Lefort et al., 2009; Kätzel et al., 2010). (3) A1 and S1BF have functional architectures which are hierarchically modular, each of which exhibit hallmarks of scale free topology marked by long tailed out degree distributions and hub neurons. This quantitative description suggests that prototypical circuits from each region strike a balance between independent computation within a module and efficient transmission of information between modules (Ravasz et al., 2002; Gallos et al., 2012). (4) The functional circuit architecture, i.e., the dynamics, of both regions is a scaled manifestation of a similar basic circuit structure. These dynamics, when considered in the context of hierarchical modularity, imply that the general principles of circuit design are optimal for robust and stable activity over a large operational range.

## Footnotes

This work was supported by the DANA Foundation (J.N.M.), the National Science Foundation CAREEER Award 0952686 (A.J.S., J.N.M.), and National Institute of General Medical Sciences Grant GM007839 (A.J.S.). We thank Nicolas Brunel, Dave McLean, S. Murray Sherman, members of the MacLean laboratory for comments on the paper, and Veronika Hanko for histological work.

The authors declare no competing financial interests.

- Correspondence should be addressed to Jason MacLean, Department of Neurobiology, University of Chicago, Chicago, Illinois 60637. jmaclean{at}uchicago.edu