Full-density multi-scale account of structure and dynamics of macaque visual cortex

We present a multi-scale spiking network model of all vision-related areas of macaque cortex that represents each area by a full-scale microcircuit with area-specific architecture. The layer- and population-resolved network connectivity integrates axonal tracing data from the CoCoMac database with recent quantitative tracing data, and is systematically refined using dynamical constraints. Simulations reveal a stable asynchronous irregular ground state with heterogeneous activity across areas, layers, and populations. Elicited by large-scale interactions, the model reproduces longer intrinsic time scales in higher compared to early visual areas. Activity propagates down the visual hierarchy, similar to experimental results associated with visual imagery. Cortico-cortical interaction patterns agree well with fMRI resting-state functional connectivity. The model bridges the gap between local and large-scale accounts of cortex, and clarifies how the detailed connectivity of cortex shapes its dynamics on multiple scales.


Introduction
Cortical activity has distinct but interdependent features on local and global scales, molded by connectivity on each scale. Globally, resting-state activity has characteristic patterns of correlations (Vincent et al., 2007;Shen et al., 2012) and propagation (Mitra et al., 2014) between areas. Locally, neurons spike with time scales that tend to increase from sensory to prefrontal areas (Murray et al., 2014) in a manner influenced by both short-range and long-range connectivity (Chaudhuri et al., 2015). We present a full-density multi-scale spiking network model in which these features arise naturally from its detailed structure.
Models of cortex have hitherto used two basic approaches. The first models each neuron explicitly in networks ranging from local microcircuits to small numbers of connected areas (Hill & Tononi, 2005;Haeusler et al., 2009).
The second represents the large-scale dynamics of cortex by simplifying the ensemble dynamics of areas or populations to few differential equations, such as Wilson-Cowan or Kuramoto oscillators (Deco et al., 2009;Cabral et al., 2011).
These models can for instance reproduce resting-state oscillations at ∼ 0.1 Hz. Chaudhuri et al. (2015) developed a mean-field multi-area model with a hierarchy of intrinsic time scales in the population firing rates, relying on a gradient of excitation across areas.
Cortical processing is not restricted to one or few areas, but results from complex interactions between many areas involving feedforward and feedback processes (Lamme et al., 1998;Rao & Ballard, 1999). At the same time, the high degree of connectivity within areas (Angelucci et al., 2002a;Markov et al., 2011) hints at the importance of local processing. Capturing both aspects requires multi-scale models that combine the detailed features of local microcircuits with realistic inter-area connectivity. Another advantage of multi-scale modeling is that it enables testing the equivalence between population models and models at cellular resolution instead of assuming it a priori.
Two main obstacles of multi-scale simulations are now gradually being overcome. First, such simulations require large resources on high-performance clusters or supercomputers and simulation technology that uses these resources efficiently. Recently, important technological progress has been achieved for the NEST simulator (Kunkel et al., 2014). Second, gaps in anatomical knowledge have prevented the consistent definition of multi-area models. Recent developments in the CoCoMac database (Bakker et al., 2012) and quantitative axonal tracing (Markov et al., 2014a,b) have systematized connectivity data for macaque cortex. However, it remains necessary to use statistical regularities such as relationships between architectural differentiation and connectivity (Barbas, 1986;Barbas & Rempel-Clower, 1997) to fully specify large cortical network models. Because of these difficulties, few large-scale spiking network models have been simulated to date, and existing ones heavily downscale the number of synapses per neuron (Izhikevich & Edelman, 2008;Preissl et al., 2012), generally affecting network dynamics .
We here use realistic numbers of synapses per neuron, building on a recent model of a 1 mm 2 cortical microcircuit with ∼ 10 5 neurons (Potjans & Diesmann, 2014). This is the smallest network size where the majority of inputs per neuron (∼ 10,000) is self-consistently represented at realistic connectivity (∼ 10%). Nonetheless, a substantial fraction of synapses originates outside the microcircuit and is replaced by stochastic input. Our model reduces random input by including all vision-related areas.
The model combines simple single-neuron dynamics with complex connectivity and thereby allows us to study the influence of the connectivity itself on the network dynamics. The connectivity map customizes that of the microcircuit model to each area based on its architecture and adds inter-areal connections. By a mean-field method (Schuecker et al., 2015), we refine the connectivity to fulfill the basic dynamical constraint of nonzero and non-saturated activity.
The ground state of cortex features asynchronous irregular spiking with low pairwise correlations (Ecker et al., 2010) and low spike rates (∼ 0.1 − 30 spikes/s) with inhibitory cells spiking faster than excitatory ones (Swadlow, 1988). Our model reproduces each of these phenomena, bridging the gap between local and global brain models, and relating the complex structure of cortex to its spiking dynamics.

Results
The model comprises 32 areas of macaque cortex involved in visual processing in the parcellation of Felleman & Van Essen (1991), henceforth referred to as FV91 (Table S1). Each area contains an excitatory and an inhibitory population in each of the layers 2/3, 4, 5 and 6 (L2/3, L4, L5, L6), except area TH, which lacks L4. The model, summarized in Table 1, represents each area by a 1 mm 2 patch.

Area-specific laminar compositions
Neuronal volume densities provided in a different parcellation scheme are mapped to the FV91 scheme and partly estimated using the average density of each layer across areas of the same architectural type ( Figure 1A). Architectural types (Table 4 of Hilgetag et al., 2015) reflect the distinctiveness of the lamination as well as L4 thickness, with agranular cortices having the lowest and V1 the highest value. Neuron density increases with architectural type.
When referring to architectural types, we also use the term 'structural hierarchy'. We call areas like V1 and V2 at the bottom of the structural (or processing) hierarchy 'early', and those near the top 'higher' areas.
We find total cortical thicknesses of 14 areas to decrease with logarithmized overall neuron densities, enabling us to estimate the total thicknesses of the other 18 areas ( Figure 1B). Quantitative data from the literature combined with our own estimates from published micrographs (Table S5) determine laminar thicknesses ( Figure 1C). L4 thickness relative to total cortical thickness increases with the logarithm of overall neuron density, which predicts relative L4 thickness for areas with missing data. Since the relative thicknesses of the other layers show no notable change with architectural type, we fill in missing values using the mean of the known data for these quantities and then normalize the sum of the relative thicknesses to 1. Layer thicknesses then follow from relative thickness times total thickness (see Table S6).
Finally, for lack of more specific data, the proportions of excitatory and inhibitory neurons in each layer are taken from cat V1 (Binzegger et al., 2004). Multiplying these with the laminar thicknesses and neuron densities yields the population sizes (see Experimental procedures).
Each neuron receives synapses of four different origins ( Figure 1D). In the following, we describe how the counts for these synapse types are computed (details in Experimental procedures).

Scalable scheme of local connectivity
We assume constant synaptic volume density across areas (Harrison et al., 2002). Experimental values for the average indegree in monkey visual cortex vary between 2,300 (O' Kusky & Colonnier, 1982) and 5,600 (Cragg, 1967) synapses per neuron. We take the average (3,950) as representative for V1, resulting in a synaptic density of 8.3 · 10 8 synapses mm 3 . The microcircuit model of Potjans & Diesmann (2014) serves as a prototype for all areas. The indegrees are a defining characteristic of this local circuit, as they govern the mean synaptic currents. We thus preserve their relative values when customizing the microcircuit to area-specific neuron densities and laminar thicknesses. The connectivity between populations is spatially uniform. The connection probability averages an underlying Gaussian connection profile over a disk with the surface area of the simulated area, separating simulated local synapses (type I) formed within the disk from non-simulated local synapses (type II) from outside the disk ( Figure 1D, E). In retrograde tracing experiments, Markov et al. (2011) found the fraction of labeled neurons intrinsic to each injected area (F LN i ) to be approximately constant, with a mean of 0.79. We translate this to numbers of synapses by assuming that the proportion of synapses of type I is 0.79 for realistic area size. For the 1 mm 2 model areas, we obtain an average proportion of type I synapses of 0.504.

Layer-specific heterogeneous cortico-cortical connectivity
We treat all cortico-cortical connections as originating and terminating in the 1 mm 2 patches, ignoring their spatial divergence and convergence. Two areas are connected if the connection is in CoCoMac ( Figure 1F) or reported by Markov et al. (2014a). For the latter we assume that the average number of synapses per labeled neuron is constant across projecting areas ( Figure 1G). To estimate missing values, we exploit the exponential decay of connectivity with distance (Ercsey-Ravasz et al., 2013). We first map the data from its native parcellation scheme (M132) to the FV91 scheme (see Experimental procedures) and then perform a least-squares fit ( Figure 1H). Combining the binary information on the existence of connections with the connection densities gives the area-level connectivity matrix ( Figure 1I).
Next, we distribute synapses between the populations of each pair of areas ( Figure 1K). The pattern of source layers is based on CoCoMac, if laminar data is available. Fractions of supragranular labeled neurons (SLN ) from retrograde tracing experiments yield proportions of projecting neurons in supra-and infragranular layers (Markov et al., 2014b). To predict missing values, we exploit a sigmoidal relation between the logarithmized ratios of cell densities of the participating areas and the SLN of their connection (as suggested by Beul et al. 2015; Figure 1J). Following Markov et al. (2014b), we use a generalized linear model for the fit and assume a beta-binomial distribution of source neurons. Since Markov et al. (2014b) do not distinguish infragranular layers further into L5 and L6, we use the more detailed laminar patterns from CoCoMac for this purpose, if available. We exclude L4 from the source patterns, in line with anatomical observations (Felleman & Van Essen, 1991), and approximate cortico-cortical connections as purely excitatory (Salin & Bullier, 1995;Tomioka & Rockland, 2007).
We base termination patterns on anterograde tracing studies collected in CoCoMac, if available, or on a relationship between source and target patterns (see Experimental procedures). Since neurons can receive synapses in different layers on their dendritic branches, we use laminar profiles of reconstructed cell morphologies (Binzegger et al., 2004) to relate synapse to cell-body locations. Despite the use of a point neuron model, we thus take into account the layer specificity of synapses on the single-cell level. In contrast to laminar synapse distributions, the resulting laminar distributions of target cell bodies are not highly distinct between feedforward and feedback projections.

Brain embedding
Inputs from outside the scope of our model, i.e., white-matter inputs from non-cortical or non-visual cortical areas and gray-matter inputs from outside the 1 mm 2 patch, are represented by Poisson spike trains. Corresponding numbers of synapses are not available for all areas, and laminar patterns of external inputs differ between target areas (Felleman & Van Essen, 1991;Markov et al., 2014b). Therefore, we determine the total number of external synapses onto an area as the total number of synapses minus those of type I and III, and distribute them with equal indegree for all populations.

Refinement of connectivity by dynamical constraints
Parameter scans based on mean-field theory (Schuecker et al., 2015) and simulations reveal a bistable activity landscape with two coexisting stable fixed points. The first has reasonable firing rates except for populations 5E and 6E, which are nearly silent (Figure 2A), while the second has excessive rates ( Figure 2B) in almost all populations.
Depending on the parameter configuration, either the low-activity fixed point has a sufficiently large basin of attraction for the simulated activity to remain near it, or fluctuations drive the network to the high-activity fixed point. To counter this shortcoming, we define an additional parameter κ which increases the external drive onto 5E by a factor κ = K ext,5E /K ext compared to the external drive of the other cell types. Since the rates in population 6E are even lower, we increase the external drive to 6E by a slightly larger factor than that to 5E. When applied directly to the model, even a small increase in κ already drives the network into the undesired high-activity state ( Figure 2B). Using the stabilization procedure described in Schuecker et al. (2015), we derive targeted modifications of the connectivity within the margins of uncertainty of the anatomical data, with an average relative change in total indegrees (summed over source populations) of 11.3% ( Figure S1B). This allows us to increase κ while retaining the global stability of the low-activity state. In the following, we choose κ = 1.125, which gives K 6E,ext /K ext = 1.417 and the external inputs listed in Table S11, and g = −11, ν ext = 10 spikes/s, yielding reasonable firing rates in populations 5E and 6E ( Figure 2C). In total, the 4.13 million neurons of the model are interconnected via 2.42 · 10 10 synapses.
The stabilization renders the intrinsic connectivity of the areas more heterogeneous. Cortico-cortical connection densities similarly undergo small changes, but with a notable reduction in the mutual connectivity between areas 46 and FEF. For more details on the connectivity changes, see Schuecker et al. (2015).

Community structure of anatomy relates to functional organization
We test if the stabilized network retains known organizing principles by analyzing the community structure in the weighted and directed graph of area-level connectivity. The map equation method (Rosvall et al., 2010) reveals 6 clusters ( Figure 3). We test the significance of the corresponding modularity Q = 0.32 by comparing with 1000 surrogate networks conserving the total outdegree of each area by shuffling its targets. This yields Q = −0.02 ± 0.03, indicating the significance of our clustering. The community structure reflects anatomical and functional properties of the areas. Two large clusters comprise ventral and dorsal stream areas, respectively. Ventral area VOT is grouped with early visual area VP. Early sensory areas V1 and V2 form a separate cluster, as well as parahippocampal areas TH and TF. The two frontal areas FEF and 46 form the last cluster. Nonetheless, the clusters are heavily interconnected ( Figure 3). The basic separation into ventral and dorsal clusters matches that found in the connection matrix of   Figure 1 Construction principles of the model. (A) Laminar neuron densities for the architectural types in the model. Type 2, here corresponding only to area TH, lacks L4. In the model, L1 contains synapses but no neurons. Data provided by H. Barbas and C. Hilgetag (personal communication) and linearly scaled up to account for undersampling of cells by NeuN staining relative to Nissl staining as determined by repeat measurements of 11 areas. (B) Total thickness vs. logarithmized overall neuron density and linear least-squares fit (r = −0.7, p = 0.005). (C) Relative laminar thickness (see Table S5) vs.    , and g = −11, νext = 10 spikes/s, κ = 1.125 with the modified connectivity matrix (C). The color bar holds for all three panels. Areas are ordered according to their architectural type along the horizontal axis from V1 (type 8) to TH (type 2) and populations are stacked vertically. The two missing populations 4E and 4I of area TH are marked in black and firing rates < 10 −2 Hz in gray. Bottom row: Histogram of population-averaged firing rates for excitatory (red) and inhibitory (blue) populations. The horizontal axis is split into linear-(left) and log-scaled (right) ranges.
connectivity matrix, but there are also important differences. For instance, our clustering groups areas STPa, STPp, and 7a with the dorsal instead of the ventral stream, better matching the scheme described by Nassi & Callaway (2009).

Area-and population-specific activity in the resting state
The model with cortico-cortical synaptic weights equal to local weights displays a reasonable ground state of activity but no substantial inter-area interactions ( Figure S2). To control these interactions, we scale cortico-cortical synaptic weights w cc onto excitatory neurons by a factor λ = J E cc /J and provide balance by increasing the weights J I cc onto inhibitory neurons by twice this factor, J I cc = λ I λJ = 2λJ. In the following, we choose λ = 1.9. Simulations yield irregular activity with plausible firing rates ( Figure 4A-C). Irregularly occurring population bursts of different lengths up to several seconds arise from the asynchronous baseline activity ( Figure 4G) and propagate across the network.
The firing rates differ across areas and layers and are generally low in L2/3 and L6 and higher in L4 and L5, partly due to the cortico-cortical interactions ( Figure 4D). The overall average rate is 14.6 spikes/s. Inhibitory populations are generally more active than excitatory ones across layers and areas despite the identical intrinsic properties of the two cell types. However, the strong participation of L5E neurons in the cortico-cortical interaction bursts causes these to fire more rapidly than L5I neurons. Pairwise correlations are low throughout the network ( Figure 4E). Excitatory neurons are more synchronized than inhibitory cells in the same layer, except for L6. Spiking irregularity is close to that of a Poisson process across areas and populations, with excitatory neurons consistently firing more irregularly than inhibitory cells ( Figure 4F). Higher areas exhibit bursty spiking, as illustrated by the raster plot for area FEF ( Figure 4C).  Figure 3 Community structure of the model. Clusters in the connectivity graph, indicated by the color of the nodes: Early visual areas (green), dorsal stream areas (red), areas VP and VOT (light blue), ventral stream (dark blue), parahippocampal areas (brown), and frontal areas (purple). Black, connections within clusters; gray, connections between clusters. Line thickness encodes logarithmized outdegrees. Only edges with relative outdegree> 10 −3 are shown.   Average time scale per architectural type indicated by triangles and overall trend by black curve. Area MDP (architectural type 5) has a time scale of 2 ms because it is uncoupled from other areas due to the lack of incoming connections.

Intrinsic time scales increase with structural hierarchy
We tested whether the model accounts for the hierarchical trend in intrinsic time scales observed in macaque cortex (Murray et al., 2014). Indeed, autocorrelation width in the model increases from early visual to higher areas. In early visual areas including V1, the autocorrelation decays with τ < 2.5 ms, indicating near-Poissonian spiking ( Figure 5A).
In higher areas, autocorrelations are broader with decay times ∼ 10 2 ms. The long time scales reflect bursty spike patterns of single-neuron activity (Figure 4), caused by the low neuron density in higher areas and thus high indegrees (2014), we find the time scale of area LIP to exceed that of MT, albeit by a small amount.

Structural and hierarchical directionality of spontaneous activity
To investigate inter-area propagation, we determine the temporal order of spiking ( Figure 6A) based on the correlation between areas. We detect the location of the extremum of the correlation function for each pair of areas ( Figure   6B) and collect the corresponding time lags in a matrix ( Figure 6C). In analogy to structural hierarchies based on pairwise connection patterns (Reid et al., 2009), we look for a temporal hierarchy that best reflects the order of activations for all pairs of areas (see Experimental procedures). The result ( Figure 6D) places parietal and temporal areas at the beginning and early visual as well as frontal areas at the end. The first and second halves of the time series yield qualitatively identical results ( Figure S3). Figure 6E shows the consistency of the hierarchy with the pairwise lags. To quantify the goodness of the hierarchy, we counted the pairs of areas for which it indicates a wrong ordering. The number of such violations is 190 out of 496, well below the 230 ± 12 (SD) violations obtained for 100 surrogate matrices, created by shuffling the entries of the original matrix while preserving its antisymmetric character.
This indicates that the simulated temporal hierarchy reflects nonrandom patterns. The propagation is mostly in the feedback direction not only in terms of the structural hierarchy, but also spatially: activity starts in parietal regions, and spreads to the temporal and occipital lobes ( Figure 6F). However, activity troughs in frontal areas follow peaks in occipital activity and thus appear last.

Emerging interactions mimic experimental functional connectivity
We compute the area-level functional connectivity (FC) based on the synaptic input current to each area, which has been shown to be more comparable to the BOLD fMRI than the spiking output (Logothetis et al., 2001). The FC matrix exhibits a rich structure, similar to experimental resting-state fMRI ( Figure 7A, B, see Experimental procedures for details). In the simulation, frontal areas 46 and FEF are more weakly coupled with the rest of the network, but the anticorrelation with V1 is well captured by the model ( Figure S4). Moreover, area MDP sends connections to, but does not receive connections from other areas according to CoCoMac, limiting its functional coupling to the network.
Louvain clustering (Blondel et al., 2008), an algorithm optimizing the modularity of the weighted, undirected FC graph (Newman, 2004), yields two modules for both the simulated and the experimental data. The modules from the simulation differ from those of the structural connectivity and reflect the temporal hierarchy shown in Figure   6C. Cluster 1S merges early visual with ventral and two dorsal regions with average level in the temporal hierarchy of h = 0.47 ± 0.13 (SD). Cluster 2S contains mostly temporally earlier areas (h = 0.33 ± 0.25 (SD)) merging parahippocampal with dorsal but also frontal areas. The experimental module 2E comprises only dorsal areas, while 1E consists of all other areas including also eight dorsal areas.
The structural connectivity of our model shows higher correlation with the experimental FC (r Pearson = 0.34) than the binary connectivity matrices from both a previous (Shen et al., 2015) and the most recent release of CoCoMac (r Pearson = 0.20), further validating our weighted connectivity matrix. For increasing weight factor λ, the correlation between simulation and experiment improves ( Figure 7D). For λ = 1, areas interact weakly, resulting in low correlation between simulation and experiment ( Figure S2). For intermediate cortico-cortical connection strengths, the correlation of simulation vs. experiment exceeds that between the structural connectivity and experimental FC ( Figure 7C), indicating the enhanced explanatory power of the dynamical model. From λ = 2 on, the network is prone to switch to the high-activity state ( Figure S5). Thus, the highest correlation (r Pearson = 0.47 for λ = 1.9) occurs just below the onset of a state in which the model visits both the low-activity and high-activity attractors.    Areas are ordered according to a clustering with the Louvain algorithm (Blondel et al., 2008) applied to the simulated data (top row) and to the experimental data (bottom row), respectively (see Experimental procedures). (C) Alluvial diagram showing the differences in the clusters for the structural connectivity (left), the simulated FC (center) and the experimentally measured FC (right). (D) Pearson correlation coefficient of simulated FC vs. experimentally measured FC for varying λ with λI = 2 (triangles) and λI = 1 (dot, cf. Figure S2). Dashed line, Pearson correlation coefficient of structural connectivity vs. experimentally measured FC.

Discussion
In this work, we present a full-density spiking multi-scale network model of all vision-related areas of macaque cortex.
An updated connectivity map at the level of areas, layers, and neural populations defines its structure. Simulations of the network on a supercomputer reveal good agreement with multi-scale dynamical properties of cortex and supply testable hypotheses. Consistent with experimental results, the local structure of areas supports higher firing rates in inhibitory than in excitatory populations, and a laminar pattern with low firing rates in layers 2/3 and 6 and higher rates in layers 4 and 5. When cortico-cortical interactions are substantial, the network shows dynamic characteristics reflecting both local and global structure. Individual cells spike irregularly with increasing intrinsic time scales along the visual hierarchy and activity propagates in the feedback direction. Functional connectivity in the model agrees well with that from resting-state fMRI and yields better predictions than the structural connectivity alone. These features are direct consequences of the multi-scale structure of the network.
The structure of the model integrates a wide range of anatomical data, complemented with statistical predictions.
The cortico-cortical connectivity is based on axonal tracing data collected in a new release of CoCoMac (Bakker et al., 2012) and recent quantitative and layer-specific retrograde tracing (Markov et al., 2014b,a). We fill in missing data using relationships between laminar source and target patterns (Felleman & Van Essen, 1991;Markov et al., 2014b), and statistical dependencies of cortico-cortical connectivity on distance (Ercsey-Ravasz et al., 2013) and architectural differentiation Hilgetag et al., 2015), an approach for which Barbas (1986); Barbas & Rempel-Clower (1997) laid the groundwork. The use of axonal tracing results avoids the pitfalls of diffusion MRI data, which strongly depend on tractography parameters and are unreliable for long-range connections (Thomas et al., 2014).
Direct comparison of tracing and tractography data moreover reveals that tractography is particularly unreliable at fine spatial scales, and tends to underestimate cortical connectivity (Calabrese et al., 2015b).
Our model customizes the microcircuit of Potjans & Diesmann (2014) based on the specific architecture of each area, taking into account neuronal densities and laminar thicknesses. A stabilization procedure (Schuecker et al., 2015) further diversifies the internal circuitry of areas. Neuronal densities in the model decrease up the structural hierarchy, in line with an observed caudal-to-rostral gradient (Charvet et al., 2015). Combined with a constant synaptic volume density (O'Kusky & Colonnier, 1982;Cragg, 1967) this yields higher indegrees up the hierarchy.
This trend matches an increase in dendritic spines per pyramidal neuron (Elston & Rosa, 2000;Elston, 2000;Elston et al., 2011), also used in a recent multi-area population rate model (Chaudhuri et al., 2015). The local connectivity can be further refined using additional area-specific data.
We find total cortical thickness to decrease with logarithmized total neuron density. Similarly, total thicknesses from MR measurements decrease with architectural type (Wagstyl et al., 2015), which is known to correlate strongly with cell density . In our data set, total and layer 4 thickness are also negatively correlated with architectural type, but these trends are less significant than those with logarithmized neuron density. Laminar and total cortical thicknesses are determined from micrographs, which has the drawback that this covers only a small fraction of the surface of each cortical area. For absolute but not relative thicknesses, another caveat is potential shrinkage and obliqueness of sections. It has also been found that relative laminar thicknesses depend on the sulcal or gyral location of areas, which is not offset by a change in neuron densities (Hilgetag & Barbas, 2006). However, regressing our relative thickness data against cortical depth of the areas registered to F99 revealed no significant trends of this type ( Figure S6). Laminar thickness data are surprisingly incomplete, considering that this is a basic anatomical feature of cortex. In future, more systematic estimates from anatomical studies or MRI may become available. Total thicknesses have already recently been measured across cortex (Calabrese et al., 2015a;Wagstyl et al., 2015), and could complement the dataset used here covering 14 of the 32 areas. However, when computing numbers of neurons, using histological data may be preferable, because shrinkage effects on neuronal densities and laminar thicknesses partially cancel out.
In the model, we statistically assign synapses to target neurons based on anatomical reconstructions (Binzegger et al., 2004). On the target side, this yields similar laminar cell-body distributions for feedforward and feedback projections despite distinct laminar synapse distributions, mirroring findings in early visual cortex of mouse (De Pasquale & Sherman, 2011). Prominent experimental results on directional differences in communication patterns are based on LFP, ECoG and MEG recordings (van Kerkoerle et al., 2014;Bastos et al., 2015;Michalareas et al., 2016), which mostly reflect synaptic inputs. In future, these findings may be integrated into the stabilization procedure to better capture such differential interactions. While this is expected to enhance the distinction between average connection patterns for feedforward and feedback projections, known anatomical patterns suggest that a substantial fraction of individual pairs of areas deviate from a simple rule (Felleman & Van Essen, 1991;Krumnack et al., 2010;Bakker et al., 2012). The cortico-cortical connectivity may be further refined by incorporating the dual counterstream organization of feedforward and feedback connections (Markov et al., 2014b), or by taking into account different numbers of inter-area synapses per neuron in feedforward and feedback directions (Rockland, 2004).
In the resulting connectivity, we find multiple clusters reflecting the anatomical and functional partition of visual cortex into early visual areas, ventral and dorsal streams, parahippocampal and frontal areas, showing that the model construction yields a meaningful network structure. Moreover, the graded structural connectivity of the model agrees better with the experimentally measured resting-state activity than the binary connectivity from CoCoMac.
The network exhibits an asynchronous, irregular ground state across the network with population bursts due to inter-area interactions. Population firing rates differ across layers and inhibitory rates are generally higher than excitatory ones, in line with experimental findings (Swadlow, 1988;Fujisawa et al., 2008;Sakata & Harris, 2009). This can be attributed to the connectivity, because excitatory and inhibitory neurons are equally parametrized and excitatory neurons receive equal or stronger external stimulation compared to inhibitory ones. Laminar activity patterns vary across areas due to their customized structure and cortico-cortical connectivity.
Intrinsic single-cell time scales in the model are short in early visual areas and long in higher areas, on the same order of magnitude as found experimentally (Murray et al., 2014). The long time scales in higher areas are related to bursty firing associated with the high indegrees in these areas, but only occur in the presence of cortico-cortical interactions. Thus, the model predicts that the pattern of intrinsic time scales has a multi-scale origin. Systematic differences in synaptic composition across cortical regions and layers (Zilles et al., 2004;Hawrylycz et al., 2012) may also contribute to the experimentally observed pattern of time scales.
Inter-area interactions in the model are mainly mediated by population bursts of different lengths. The degree of synchrony accompanying inter-area interactions in the brain is as yet unclear. Obtaining substantial corticocortical interactions with low synchrony may be possible with finely structured connectivity and reduced noise input.
When neurons are to a large extent driven by a noisy external input, a smaller percentage of their activity is determined by intrinsic inputs, which can decrease their effective coupling (Aertsen & Preißl, 1990). One way of reducing the external drive while preserving the mean network activity may be for the drive to be attuned to the intrinsic connectivity (Marre et al., 2009). Stronger intrinsic coupling while maintaining stability may be achieved for instance by introducing specific network structures such as synfire chains (Diesmann et al., 1999) or other feedforward structures, subnetworks, or small-world connectivity (Jahnke et al., 2014); population-specific patterns of short-term plasticity (Sussillo et al., 2007); or fine-tuned inhibition between neuronal groups (Hennequin et al., 2014).
The synchronous population events propagate stably across multiple areas, predominantly in the feedback direction. The systematic activation of parietal before occipital areas in the model is reminiscent of EEG findings on information flow during visual imagery (Dentico et al., 2014) and the top-down propagation of slow waves during sleep (Massimini et al., 2004;Nir et al., 2011;Sheroziya & Timofeev, 2014). Our method for determining the order of activations is similar to one recently applied to fMRI recordings (Mitra et al., 2014). It could be extended to distinguish between excitatory and inhibitory interactions like those we observe between V1 and frontal areas ( Figure   S4). In the network, cortico-cortical projections target both excitatory and inhibitory populations, with the majority of synapses terminating on excitatory cells. Stronger cortico-cortical synapses to enhance inter-area interactions require increased balancing of cortico-cortical inputs to preserve network stability. This is similar to the "handshake" mechanism in the microcircuit model of Potjans & Diesmann (2014) where interlaminar projections provide network stability by their inhibitory net effect.
The pattern of simulated interactions between areas resembles fMRI resting-state activity. The agreement between simulation and experiment peaks at intermediate coupling strength, where synchronized clusters also emerged most clearly in earlier models (Zhou et al., 2006;Deco & Jirsa, 2012). Furthermore, optimal agreement occurs just below a transition to a state where the network switches between attractors, supporting evidence that the brain operates in a slightly subcritical regime (Deco & Jirsa, 2012;Priesemann et al., 2014).
Time series of spiking activity reveal broad-band transmission between areas on time scales up to several seconds.
The low-frequency part of these interactions is comparable to fMRI data, which describes coherent fluctuations on the order of seconds. The long time scales in the model activity may be caused by eigenmodes of the effective connectivity that are close to instability  or non-orthogonal (Hennequin et al., 2012). A potential future avenue for research would be to distinguish between such network effects and other sources of long time scales such as NMDA and GABA B transmission, neuromodulation, or adaptation effects.
For tractability, the model represents each area as a 1 mm 2 patch of cortex. True area sizes vary from ∼ 3 million cells in TH to ∼ 300 million cells in V1 for a total of around 8 · 10 8 neurons in one hemisphere of macaque visual cortex, a model size that with recent advances in simulation technology (Kunkel et al., 2014) already fits on the most powerful supercomputers available today. Approaching this size would reduce the negative effects of downscaling . Overall, our model elucidates multi-scale relationships between cortical structure and dynamics, and can serve as a platform for the integration of new experimental data, the creation of hypotheses, and the development of functional models of cortex. source and target neurons drawn randomly with replacement (allowing autapses and multapses) with area-and population-specific connection probabilities Weights fixed, drawn from normal distribution with mean J and standard deviation δJ = 0.1J; 4E to 2/3E increased by factor 2 (cf. Potjans & Diesmann, 2014); weights of inhibitory connections increased by factor g; excitatory weights < 0 and inhibitory weights > 0 are redrawn; cortico-cortical weights onto excitatory and inhibitory populations increased by factor λ and λ I λ, respectively Delays fixed, drawn from Gaussian distribution with mean d and standard deviation δd = 0.5d; delays of inhibitory connections factor 2 smaller; delays rounded to the nearest multiple of the simulation step size h = 0.1 ms, inter-areal delays drawn from a Gaussian distribution with mean d = s/v t , with distance s and transmission speed v t = 3.5 m/s (Girard et al., 2001); and standard deviation δd = d/2, distances determined as described in Supplemental Experimental Procedures, delays < 0.1 ms before rounding are redrawn D: Neuron and synapse model Name LIF neuron Type leaky integrate-and-fire, exponential synaptic current inputs Subthreshold dynamics  Table S3) F: Measurements Spiking activity In the following, we detail how we derive the structure of the model (summarized in Table 1), i.e., the population sizes, the local and cortico-cortical connectivity and the external drive.

Numbers of neurons
We estimate the number of neurons N (A, i) in population i of area A in three steps: 1. We translate neuronal volume densities to the FV91 scheme from the most representative area in the original scheme (Table S4). For areas not covered by the data set, we take the average laminar densities for areas of the same architectural type. Table 4 of Hilgetag et al. (2015) lists the architectural types, which we translate to the FV91 scheme according to Table S4. To the previously unclassified areas MIP and MDP we manually assign type 5 like their neighboring area PO, which is similarly involved in visual reaching Galletti et al., 2003), and was placed at the same hierarchical level by Felleman & Van Essen (1991).
2. We determine total and laminar thicknesses as detailed in Results.
3. The fraction γ(v) of excitatory neurons in layer v is taken to be identical across areas. For the laminar dependency, values from cat V1 (Binzegger et al., 2004) are used with 78% excitatory neurons in layer 2/3, 80% in L4, 82% in L5, and 83% in L6.
The resulting number of neurons in population i of area A is where v i denotes the layer of population i, S(A) the surface area of area A (cf . Table S7), D(A, v i ) the thickness of layer v i , and E, I the pool of excitatory and inhibitory populations, respectively. Table S8 gives the population sizes corresponding to the modeled 1 mm 2 area size.

Local connectivity
The connection probabilities of the microcircuit model (Potjans & Diesmann, 2014,  form the basis for the local circuit of each area. The connectivity between any pair of populations is spatially uniform. However, we take the underlying probability C for a given neuron pair to establish one or more contacts to decay with distance according to a Gaussian with standard deviation σ = 297 µm (Potjans & Diesmann, 2014). We approximate each brain area as a flat disk with (area-specific) radius R and assign polar coordinates r and θ to each neuron.
The radius determines the cut-off of the Gaussian and hence the precise connectivities. The average connection probability is obtained by integrating over all possible positions of the two neurons: π 0 exp − r 2 1 + r 2 2 − 2r 1 r 2 cos(θ 1 − θ 2 ) 2σ 2 r 1 r 2 dθ 1 dr 1 dθ 2 dr 2 , with C 0 the connection probability at zero distance. This can be reduced to a simpler form (Sheng, 1985), Averaged across population pairs, C 0 is 0.143 (computed from Eq. 8 and Table S1 in Potjans & Diesmann, 2014).
Note that Potjans & Diesmann (2014) only vary the position of one neuron, keeping the other neuron fixed in the center of the disk (Eq. 9 in that paper). Henceforth, we denote connection probabilities computed with the latter approach with the subscript PD14 and use primes for all variables referring to a network with the population sizes of the microcircuit model.
The parameters of the microcircuit model are reported for a 1 mm 2 patch of cortex, corresponding to R = 1/π mm, which we call R 0 . For each source population j and target population i, we first translate the connection probabilities of the 1 mm 2 model to area-dependent R via with c A (R) an area-specific conversion factor, which is larger for areas with smaller neuron densities because of the assumption of constant synaptic volume density. It is computed as in F99 space. On the target side we use the coordinates of the injection sites registered to the F99 atlas available via the Scalable Brain Atlas (Bakker et al., 2015) to identify the equivalent area in the FV91 parcellation (cf. Table   S9). There is data for 11 visual areas in the FV91 scheme with repeat injections in six areas, for which we take the arithmetic mean. To map data on the source side from M132 to FV91, we count the number of overlapping triangles on the F99 surface between any given pair of regions and distribute the F LN proportionally to the amount of overlap, using the F99 region overlap tool at the CoCoMac site (http://cocomac.g-node.org). To estimate values for the areas not included in the data set, we use an exponential decay of connectivities with distance (Ercsey-Ravasz et al., 2013), As a next step, we determine the distribution of connections across source and target layers. On the source side, the laminar projection pattern can be expressed as the fraction of supragranular labeled neurons (SLN ) in retrograde tracing experiments (Markov et al., 2014b). To determine the SLN entering into the model, we use the exact coordinates of the injections to determine the corresponding target area A in the FV91 parcellation, and for each pair of areas we take the mean SLN across injections. To map the data from M132 to FV91, we weight the SLN by the overlap c B,β between area β in the former and area B in the latter scheme and the F LN to take into account the overall strength of the connection, We estimate missing values using a sigmoidal fit of SLN vs. the logarithmized ratio of overall cell densities of the two areas ( Figure 1J). A relationship between laminar patterns and log ratios of neuron densities was suggested by Beul et al. (2015). Following Markov et al. (2014b), we use a generalized linear model and assume the numbers of labeled neurons in the source areas to sample from a beta-binomial distribution (e.g. Weisstein, 2005). This distribution arises as a combination of a binomial distribution with probability p of supragranular labeling in a given area, and a beta distribution of p across areas with dispersion parameter φ. With the probit link function g (e.g. McCulloch et al., 2008), the measured SLN relates to the log ratio of neuron densities for each pair of areas as where and SLN are vectors and {a 0 , a 1 } are scalar fit parameters. To fit SLN vs. log ratios of cell densities, we map the FV91 areas to the Markov et al. (2014b) scheme with the overlap tool of CoCoMac (see above) and compute the cell density of each area in the M132 scheme as a weighted average over the relevant FV91 areas. For areas with identical names in both schemes, we simply take the neuron density from the FV91 scheme. Figure 1J shows the result of the SLN fit in R (R Core Team, 2015) with the betabin function of the aod package (Lesnoff & Lancelot, 2012). In contrast to Markov et al. (2014b), who exclude certain areas when fitting SLN vs. hierarchical distances in view of ambiguous hierarchical relations, we take all data points into account to obtain a simple and uniform rule.
As a further step, we combine SLN with CoCoMac data. The data sets complement each other: SLN provides quantitative data on laminar patterns of incoming projections for about one quarter of the connected areas. CoCoMac has values for all six layers, but limited to a qualitative strength ranging from 0 (absent) to 3 (strong) which we take to represent numbers of synapses in orders of magnitude (see Supplemental Experimental Procedures). Whether or not to include a layer in source pattern P s is based on CoCoMac (Felleman & Van Essen, 1991;Barnes & Pandya, 1992;Suzuki & Amaral, 1994b;Morel & Bullier, 1990;Perkel et al., 1986;Seltzer & Pandya, 1994) if the corresponding data is available (45 % coverage); otherwise, we include L2/3, L5 and L6 and exclude L4 (Felleman & Van Essen, 1991). We model cortico-cortical connections as purely excitatory, a good approximation to experimental findings (Salin & Bullier, 1995;Tomioka & Rockland, 2007). If a layer is included in the source pattern, we assign a fraction of the total outgoing synapses to it according to the SLN . Since the SLN do not further distinguish between the infragranular layers 5 and 6, we use the rough connection densities from CoCoMac for this purpose when available, and otherwise we distribute synapses in proportion to the numbers of neurons. On the target side, we determine the pattern of target layers P t from anterograde tracer studies in CoCoMac (Jones et al., 1978;Rockland & Pandya, 1979;Morel & Bullier, 1990;Webster et al., 1991;Felleman & Van Essen, 1991;Barnes & Pandya, 1992;Distler et al., 1993;Suzuki & Amaral, 1994b;Webster et al., 1994)  , and we distribute synapses among the layers in the termination pattern in proportion to their thickness.
Since we use a point neuron model, we have to account for the possibly different laminar positions of cell bodies and synapses. The data of Binzegger et al. (2004) deliver three quantities that allow us to relate synapse to cell body locations: first, the probability P(s cc |c B s ∈ v) for a synapse in layer v on a cell of type c B (e.g., a pyramidal cell with soma in L5) to be of cortico-cortical origin; second, the relative occurrence P(c B ) of the cell type c B ; and third, the total numbers of synapses N syn (v, c B ) in layer v onto the given cell type. We map these data to our model by computing the conditional probability P(i|s cc ∈ v) for the target neuron to belong to population i if a cortico-cortical synapse s cc is in layer v. This probability equals the sum of probabilities that a synapse is established on the different Binzegger et al. subpopulations making up our populations, where The numerator gives the joint probability that a cortico-cortical synapse is formed in layer v on cell type c B , and the denominator is the probability of a cortico-cortical synapse in layer v, computed by summing over cell types, N syn,CC (v, c B ) represents the number of cortico-cortical synapses in layer v on cell type c B , which can be directly determined from the data. Combining these equations, we obtain the number of cortico-cortical (type III) synapses from excitatory population j of area B to population i of area A (cf. Figure 1K): if j ∈ I P s but no CoCoMac data available if no CoCoMac data available .
Z i is an additional factor which takes into account that cortico-cortical feedback connections preferentially target excitatory rather than inhibitory neurons (Johnson & Burkhalter, 1996;Anderson et al., 2011).  Figure S1 shows the resulting connection probabilities between all population pairs in the model.

External, random input
Since quantitative area-specific data on non-visual and subcortical inputs are highly incomplete, we use a simple scheme to determine numbers of external inputs: For each area, we compute the total number of external synapses as the difference between the total number of synapses and those of type I and III and distribute these such that all neurons in the given area have the same indegree for Poisson sources. In area TH, we compensate for the missing granular layer 4 by increasing the external drive onto populations 2/3E and 5E by 20 %. With the modified connectivity matrix yielded by the analytical procedure described in Schuecker et al. (2015), we set κ = 1.125 to increase the external indegree onto population 5E by 12.5 % and onto 6E by 42 % to elevate the firing rates in these populations. Table S11 lists the resulting external indegrees.

Network simulations
We performed simulations on the JUQUEEN supercomputer (Jülich Supercomputing Centre, 2015) with NEST version 2.8.0 (Eppler et al., 2015) with optimizations for the use on the supercomputer which will be included in a future release. All simulations use a time step of 0.1 ms and exact integration for the subthreshold dynamics of the LIF neuron model (reviewed in Plesser & Diesmann, 2009). Simulations were run for 100.5 s (λ = 1.9), 50.5 s (λ ∈ [1.8, 2.0, 2.1]), and 10.5 ms (λ ∈ [1., 1.5, 1.7, 2.5]) biological time discarding the first 500 ms. Spike times were recorded from all neurons, except for the simulations shown in Figure 2A,B, where we recorded from 1000 neurons per population.

Analysis methods
We investigate the structural properties of the model with the map equation method (Rosvall et al., 2010). In this clustering algorithm, an agent performs random walks between graph nodes with probability proportional to the outdegree of the present node and a probability (p = 0.15) of jumping to a random network node. The algorithm detects clusters in the graph by minimizing the length of a binary description of the network using a Huffman code.
To assess the quality of the clustering, we use a modularity measure which extends a measure for unweighted, directed networks (Leicht & Newman, 2008) to weighted networks, analogous to Newman 2004, Instantaneous firing rates are determined as spike histograms with bin width 1 ms averaged over the entire population or area. In Figure 4G, Figure S2G, and to determine the temporal hierarchy, we convolve the histograms with Gaussian kernels with σ = 2 ms. Spike-train irregularity is quantified for each population by the revised local The temporal hierarchy is based on the cross-covariance function between area-averaged firing rates. We use a wavelet-smoothing algorithm (signal.find peaks cwt of python scipy library (Jones et al., 2001) with peak width ∆ = 20 ms) to detect extrema for τ ∈ [−100, 100] and take the location of the extremum with the largest absolute value as the time lag.
Functional connectivity (FC) is defined as the zero-time lag cross-correlation coefficient of the area-averaged synaptic inputs with the normalized post-synaptic current P SC j (t) = exp[−t/τ s ], the population firing rate ν j of source population j, indegree K ij , and synaptic weight J ij of the connection from j to target population i containing N i neurons.
The clustering of the FC matrices was performed using the function modularity louvain und sign of the Brain Connectivity Toolbox (BCT; http://www.brain-connectivity-toolbox.net) with the Q * option, which weights positive weights more strongly than negative weights, as introduced by Rubinov & Sporns (2011).

Macaque resting-state fMRI
Data were acquired from six male macaque monkeys (4 Macaca mulatta and 2 Macaca fascicularis). All experimental protocols were approved by the Animal Use Subcomittee of the University of Western Ontario Council on Animal Care and in accordance with the guidelines of the Canadian Council on Animal Care. Data acquisition, image preprocessing and a subset of subjects (5 of 6) were previously described (Babapoor-Farrokhran et al., 2013). Briefly, 10 5-min resting-state fMRI scans (TR: 2 s; voxel size: 1 mm isotropic) were acquired from each subject under light anaesthesia (1.5 isoflurane). Additional processing for the current study included the regression of nuisance variables using the AFNI software package (afni.nimh.nih.gov/afni), which included six motion parameters as well as the global white matter and CSF signals. The global mean signal was not regressed.
The FV91 parcellation was drawn on the F99 macaque standard cortical surface template (Van Essen et al., 2001) and transformed to volumetric space with a 2 mm extrusion using the Caret software package (http://www.nitrc. org/projects/caret). The parcellation was applied to the fMRI data and functional connectivity computed as the Pearson correlation coefficients between probabilistically-weighted ROI timeseries for each scan .
Correlation coefficients were Fisher z-transformed and correlation matrices were averaged within animals and then across animals before transforming back to Pearson coefficients.         The thickness data is the same as in Figure 1. Cortical depth data obtained from F99 surface statistics available through the Caret Software (Van Essen, 2012). Values for each area are averaged across cortical surface and both hemispheres. The data is obtained using the F99 Sulcal depth tool on http://cocomac.gnode.org and can be directly accessed via these two links: http://cocomac.g-node.org/cocomac2/services/f99_ sulcal_depth.php?atlas=FV91&shape=Depth-Right&mode=avg&output=tsv&run=1 and http://cocomac.g-node.org/ cocomac2/services/f99_sulcal_depth.php?atlas=FV91&shape=Depth-Left&mode=avg&output=tsv&run=1. Area-averaged firing rates in V1 for four different settings of λ. The simulation for λ = 2.5 was run for 10 s biological time only. Fromλ = 2 on, the network spontaneously enters a high-activity state. For λ = 2.5, the network is in this state from the outset.

Supplemental Experimental Procedures
Cortical areas in the model     Table 4 of Hilgetag et al. (2015) to the modeled areas in the parcellation scheme of Felleman & Van Essen (1991). Entries marked with a star are used to translate the overall neuron density and cortical thickness which are not available in the finer of the two parcellations used by Hilgetag et al. (2015).   Derivation of the conversion factor c A (R) for the local connectivity The indegrees of the microcircuit model (Potjans & Diesmann, 2014) K ij (R) are adapted to the area-specific laminar compositions of the multi-area model with an area-specific factor c A (R), where i, j denote single populations in the 1 mm 2 patch of the cortical area. The total number of synapses local to the patch (type I) is the sum over the projections between all populations of the area: We thus obtain c A (R) by determining N syn,I . To this end, we use retrograde tracing data from Markov et al. (2011) consisting of fractions of labeled neurons (F LN ) per area as a result of injections into one area at a time. The fraction intrinsic to the injected area, F LN i , is approximately equal for all 9 areas where this fraction was determined, with a mean of 0.79. For areas modeled with reduced size, this fraction is smaller because, in that case, synapses of both type I and II contribute to the value of 0.79 ( Figure 1E). We approximate the increasing contribution of type I synapses with the modeled area size as the increase in indegrees averaged over population pairs, N syn,I (R)/N syn,tot (R) N syn,I (R full )/N syn,tot (R full ) where in the last step we use (4). Using N syn,I (R full )/N syn,tot (R full ) = F LN i , we obtain N syn,I (R) = N syn,tot (R) where N syn,tot (R) = ρ syn πR 2 D with D the total thickness of the given area. The conversion factor can thus be obtained with We substitute this into (4) for the modeled areas where R = R 0 and obtain the population-specific indegrees for type I synapses: K ij,I : =K ij (R = R 0 )

Processing of CoCoMac data
We use a new release of CoCoMac, in which mappings from brain regions in other nomenclatures were scrutinized to ensure a consistent transfer of connections into the FV91 name space. The CoCoMac database provides information on laminar patterns on the source side from retrograde tracing studies as well as on the target side from anterograde trac-ing studies. The data was extracted by using the following link, which specifies all search options: http://cocomac. g-node.org/cocomac2/services/connectivity_matrix.php?dbdate=20141022&AP=AxonalProjections_FV91&constrai &origins=&terminals=&square=1&merge=max&laminar=both&format=json&cite=1 Furthermore, we obtained the numbers of confirmative studies for each area-level connection with the following link: http://cocomac.g-node.org/cocomac2/services/connectivity_matrix.php?dbdate=20141022&AP=AxonalProjec

FV91&constraint=&origins=&terminals=&square=1&merge=count&laminar=off&format=json&cite=1
To process these data, we applied the following steps: • A connection is assumed to exist if there is at least one confirmative study reporting it.
• A connection from layer 2/3 is modeled if CoCoMac indicates a connection from either or both of layers 2 and 3.
• In the database, some layers carry an 'X' indicating a connection of unknown strength. We interpret these as '2' (corresponding to medium connection strength).
• We take connection strengths in CoCoMac to represent numbers of synapses in orders of magnitude, i.e., the relative number of synapses N ν syn in layer ν of area A with connection strength s(ν) is computed as N ν syn = 10 s(v) / v ∈A 10 s(ν ) .

Mapping of synapse to cell-body locations
Detailed calculation in section Experimental procedures. The numbers are listed in Table S10. Table S10 Conditional probabilities P(i|scc ∈ v) for the target neuron to belong to population i if a cortico-cortical synapse scc is located in layer v, computed with (8) applied to the data set of Binzegger et al. (2004). Empty cells signal zero probabilities.

Analytical mean-field theory
In Schuecker et al. (2015), analytical mean-field theory is derived describing the stationary population-averaged firing rates of the model. In the diffusion approximation, which is valid for high indegrees and small synaptic weights, the dynamics of the membrane potential V and synaptic current I s are described by (Fourcaud & Brunel, 2002) τ where the input spike trains are replaced by a current fluctuating around the mean µ with variance σ with fluctuations drawn from a random Gaussian process ξ(t) with ξ(t) = 0 and ξ(t)ξ(t ) = δ(t − t ). Going from the single-neuron level to a description of populations, we define the population-averaged firing rate ν i due to the population-specific input µ i , σ i . The stationary firing rates ν i are then given by (Fourcaud & Brunel, 2002) 1 which holds up to linear order in τ s /τ m and where γ = |ζ(1/2)|/ √ 2, with ζ denoting the Riemann zeta function (Abramowitz & Stegun, 1974).

Algorithm for the temporal hierarchy
To determine a temporal hierarchy for the onset of population bursts, we determine the peak locations τ AB of the cross-correlation function for each pair of areas A, B. We then define a scalar function for the deviation between the distance of hierarchical levels h(A), h(B) and peak locations, To determine the hierarchical levels, we minimize the sum of f (A, B) over all pairs of areas,

S =
A,B f (A, B) , using the optimize.minimize function of the scipy library (Jones et al., 2001) with random initial hierarchical levels. We verified that the initial choice of hierarchical levels does not influence the final result. We obtain hierarchical levels on an arbitrary scale, which we normalize to values h(A) ∈ [0, 1] ∀A.