http://orcid.org/0000-0003-1040-2567
Figures
Abstract
Cortical activity has distinct features across scales, from the spiking statistics of individual cells to global resting-state networks. We here describe the first full-density multi-area spiking network model of cortex, using macaque visual cortex as a test system. The model represents each area by a microcircuit with area-specific architecture and features layer- and population-resolved connectivity between areas. Simulations reveal a structured asynchronous irregular ground state. In a metastable regime, the network reproduces spiking statistics from electrophysiological recordings and cortico-cortical interaction patterns in fMRI functional connectivity under resting-state conditions. Stable inter-area propagation is supported by cortico-cortical synapses that are moderately strong onto excitatory neurons and stronger onto inhibitory neurons. Causal interactions depend on both cortical structure and the dynamical state of populations. Activity propagates mainly in the feedback direction, similar to experimental results associated with visual imagery and sleep. The model unifies local and large-scale accounts of cortex, and clarifies how the detailed connectivity of cortex shapes its dynamics on multiple scales. Based on our simulations, we hypothesize that in the spontaneous condition the brain operates in a metastable regime where cortico-cortical projections target excitatory and inhibitory populations in a balanced manner that produces substantial inter-area interactions while maintaining global stability.
Author summary
The mammalian cortex fulfills its complex tasks by operating on multiple temporal and spatial scales from single cells to entire areas comprising millions of cells. These multi-scale dynamics are supported by specific network structures at all levels of organization. Since models of cortex hitherto tend to concentrate on a single scale, little is known about how cortical structure shapes the multi-scale dynamics of the network. We here present dynamical simulations of a multi-area network model at neuronal and synaptic resolution with population-specific connectivity based on extensive experimental data which accounts for a wide range of dynamical phenomena. Our model elucidates relationships between local and global scales in cortex and provides a platform for future studies of cortical function.
Citation: Schmidt M, Bakker R, Shen K, Bezgin G, Diesmann M, van Albada SJ (2018) A multi-scale layer-resolved spiking network model of resting-state dynamics in macaque visual cortical areas. PLoS Comput Biol 14(10): e1006359. https://doi.org/10.1371/journal.pcbi.1006359
Editor: Lyle J. Graham, Université Paris Descartes, Centre National de la Recherche Scientifique, FRANCE
Received: November 8, 2017; Accepted: July 12, 2018; Published: October 18, 2018
Copyright: © 2018 Schmidt et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The digitized workflow reproducing all results and figures of this work was created in compliance with the guide suggested by [Pauli18], and is available as Python code from https://github.com/INM-6/multi-area-model. The simulation data presented in this manuscript is available from https://web.gin.g-node.org/maximilian.schmidt/multi-area-model-data.
Funding: This work was supported by the Helmholtz Portfolio Supercomputing and Modeling for the Human Brain (SMHB), the European Union 7th Framework Program (grant 269921, BrainScaleS and 604102, Human Brain Project, Ramp up phase), and European Union’s Horizon 2020 research and innovation program (grant 720270, Human Brain Project, SGA1), the Jülich Aachen Research Alliance (JARA), and computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC Partition part of the supercomputer JUQUEEN at Forschungszentrum Jülich (VSR computation time grant JINB33), and Priority Program 2041 (SPP 2041) "Computational Connectomics" of the German Research Foundation (DFG). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Cortical activity has distinct but interdependent features on local and global scales, molded by multi-scale connectivity. Data from multiple species including macaque indicate that the ground state of cortex locally features asynchronous irregular spiking with low pairwise correlations [1] and low, layer-specific spike rates [2, 3] with inhibitory rates exceeding excitatory ones [4–6], and activity fluctuations on multiple timescales [7]. Globally, resting-state activity has characteristic patterns of inter-area correlations [8, 9] and propagation [10]. These interactions are layer-specific and distinct between feedback and feedforward directions [11–13]. We present a full-density multi-scale spiking network model in which these features emerge from its detailed structure.
Most cortical models concentrate on either the local or the global scale, using two basic approaches. The first approach represents each neuron explicitly in networks ranging from local microcircuits to small numbers of areas [14, 15]. The second describes large-scale cortical dynamics by simplifying ensemble dynamics to few differential equations. These models predict resting-state oscillations in a metastable regime [16–19] and reproduce the frequency specificity of inter-area interactions [20].
Cortical processing is not restricted to one or few areas, but results from complex multi-area interactions [21, 22]. Simultaneously, dense within-area connectivity [23, 24] suggests the importance of local processing, where the population-specific connectivity underlies multidimensional functional properties [25] and supports a set of computational principles that underlie sensory processing across the cortex [26, 27]. Capturing both aspects requires combining detailed features of local microcircuits with realistic inter-area connectivity. Modeling at cellular resolution enables testing the equivalence with population models instead of assuming it a priori.
Two main obstacles of multi-scale simulations are gradually being overcome. First, recent progress in simulation technology enables the efficient use of supercomputers [28]. Second, systematic connectivity data is increasingly available [29, 30]. However, statistical predictions remain necessary to fully specify large cortical network models. Consequently, few large-scale spiking network models have been simulated to date, and existing ones heavily downscale the number of synapses per neuron [31, 32] (but see [33]), affecting network dynamics [34].
We here investigate a spiking multi-area network model of macaque visual cortex, covering the scales of single neurons, microcircuits, and cortical areas. The connectivity map, derived in [35], customizes that of the microcircuit model of [36] to each area based on its architecture and adds layer-specific inter-area connections. Each area is represented by a 1 mm2 microcircuit with the full density of neurons and synapses. A mean-field method [37] refines the connectivity to fulfill the basic dynamical constraint of nonzero and non-saturated activity. By combining simple single-neuron dynamics with complex connectivity, the model enables studying the influence of the connectivity itself on the network dynamics.
We first describe the refinement of the connectivity by dynamical constraints, leading to plausible spike rates. Next, we vary cortico-cortical synaptic strengths and find that with increased coupling, connections onto inhibitory neurons must outbalance connections onto excitatory neurons for stability at low rates. The resulting network state reproduces spiking statistics of V1 resting-state data [38] and yields population bursts reflecting a metastable regime [39–42]. Outside this metastable regime, the spiking statistics deviate considerably from the experimental data. Our findings thus extend previous works demonstrating metastability of cortical networks via modeling [16–19] by unifying microscopic and macroscopic descriptions and supporting the hypothesis that plausible spiking statistics require cortex to be poised in a metastable regime. Analyzing the order of activation of the areas reveals that the population bursts propagate mainly in the feedback direction. Subsequently we show that, for intermediate cortico-cortical synaptic strengths, inter-area correlation patterns resemble fMRI functional connectivity [43]. Finally, we observe directional differences in laminar patterns of inter-area communication that reflect both structural relationships and dynamical states. Our work provides a platform for future studies addressing spiking-level functional properties and for the development of analogous models of other cortical regions. Preliminary results have been presented in abstract form [44].
Materials and methods
The model comprises 32 areas of macaque cortex involved in visual processing in the parcellation of [45], henceforth referred to as FV91 (Table 1). 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 2, represents each area by a 1 mm2 patch. We here give a brief overview of the construction principles underlying the network definition, and refer to [35] for details of the derivation and an analysis of the network structure. Fig 1 summarizes the construction principles leading to the population sizes and the area-, layer-, and population-specific connectivity map. Table 3 lists the neuron and synapse parameters.
Top: Determination of the population sizes. The size of a population is computed as the product of the layer-specific cell density and its volume computed as thickness times the fixed surface area of 1 mm2. Population sizes decrease along the gradient of architectural types. The ratio between excitatory and inhibitory neurons is taken from layer-specific data of [51] leading to an average proportion of ∼80% excitatory and ∼20% inhibitory cells. Figure adapted from [35] and [36] (with permission). Bottom: Construction of the model connectivity. Local connectivity within areas is based on the microcircuit model of [36]. Cortico-cortical connections are first determined on the area level from binary data from the CoCoMac database [29, 61] and quantitative tracing data from [30], which is completed using the exponential fall-off of connection density with inter-area distance [65]. The resulting connectivity spans six orders of magnitude, as shown in the matrix plot of area-averaged indegrees (center bottom). Synapses between each pair of cortical areas are then distributed over source and target layers based on layer-specific tracing data from [66] and CoCoMac. Synapses in the receiving area are subsequently assigned to cells according to layer- and cell-type-specific dendritic densities from [51]. These derivations result in distinct laminar patterns between feedforward, feedback, and lateral connections. Based on a theoretical method using mean-field theory [37], the resulting connectivity matrix is refined to improve the phase space of the network. The detailed derivation is published in [35].
Population sizes in the model are derived from a collection of data on laminar and total cortical thicknesses, and area-specific laminar and overall cell densities ([50]; H. Barbas and C. Hilgetag, personal communication). Cell densities in cortex and thus the number of neurons in our model decrease along the visual hierarchy from primary visual cortex V1 (197,936 neurons in total) to area TH (73,251 neurons in total). The ratio between excitatory and inhibitory cells is layer-specific and roughly 4:1 on average; we use values from cat V1 [51] as these underlie the original microcircuit model [36], but the values are close to those from macaque V1 [52]. To estimate missing data on laminar and total thicknesses and neuronal densities, we use a categorization of the areas into architectural types. The architectural type of an area reflects its laminar distinctiveness and the thickness of L4 [53–55].
Based on experimental findings [56, 57], we assume a constant overall synaptic volume density across all areas, which leads to larger numbers of synapses onto neurons in higher areas, in line with an increase in the number of dendritic spines per pyramidal neuron [58–60], due to a hierarchical reduction in neuron density. We distinguish synapses onto a neuron into four distinct types. The first two types are simulated in the model: internal synapses originating within each 1 mm2 patch (type I), and cortico-cortical synapses formed by projections from other vision-related cortical areas (type II). Two other types are beyond the scope of the network under investigation: synapses from the same cortical area, but outside the 1 mm2 patch (type III), and synapses from other cortical non-visual and non-cortical areas (type IV). For simplicity, inputs from these connections are modeled as spike trains drawn from Poisson processes with stationary rate νext. We make this choice due to the lack of comprehensive data on medium-range connections within areas, such as patchy connectivity; on connections from outside visual cortex; and on the activity levels in resting-state activity of projecting areas. In a future extension of this work, the network would benefit from including such connections or at least diversifying the random inputs.
The connectivity in our model is spatially uniform within each 1 mm2 patch. However, to determine the ratio between the numbers of type I and type III synapses we still need a model for distance-dependent connectivity. For this we assume the probability for two neurons to form one or more synapses to decay with the distance between them according to a Gaussian profile. The 1 mm2 patch is taken as a disk in the center of a larger disk that represents the full area. The radius of the patch () determines the cut-off of the Gaussian. The connectivity inside areas is based on the microcircuit model of [36]. The local connectivity of this model is tailored to a cortical area of a given surface area and population sizes of cat V1. We adapt their population-specific connectivity matrix to the compositions of the 32 areas, thereby keeping the proportion of synapses that a neuron in a given population receives from a given other population constant. This translates into requiring that the ratio of indegrees between the connection from population j to i and from population l to k is kept constant when translating from the model of [36] (K′) to our model for each area A (K):
This results in an area-specific conversion factor that scales all the indegrees for a given area. We choose this approach because it approximately preserves a defining characteristic of the local circuit, namely the mean synaptic inputs, which are proportional to the indegrees.
The corticocortical connectivity combines axonal tracing data with statistical regularities to fill in missing values, and is defined in four steps: First, a directed connection between two areas is established if it is included in the CoComac database [29, 45, 61–64] or the retrograde tracing dataset of [30]. Second, the total number of synapses between two areas is based on fractions of labeled neurons (FLN) measured in the experiment by [30] or estimated using the exponential fall-off of FLN with inter-area distance [65]. Third, we determine laminar patterns based on fractions of supragranular labeled neurons (SLN), both measured [66] and estimated using a sigmoidal fit against differences in neuronal densities, and layer-specific data on target [45, 63, 64, 67–72] and source [45, 64, 68, 71, 73, 74] patterns from CoCoMac. Finally, we account for the possibly different laminar positions of synapses and cell bodies in the target areas by computing the conditional probability that a cortico-cortical synapse in a layer ν is formed on a dendrite of a specific cell type, using the morphological reconstructions of [51].
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 from within the 1 mm2 patch and other modeled areas, 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%. Overall, external inputs amount to approximately 35% of the total inputs to each neuron in the network.
These derivations lead to a connectivity map with similar but non-identical local circuits. The connectivity between areas shows high heterogeneity with connection densities spanning six orders of magnitude, reflecting the similarly high degree of heterogeneity in the quantitative tracing data of [30] (connectivity matrix plot in Fig 1). The laminar patterns of cortico-cortical connections are distinct between feedforward and feedback connections. They mostly follow classical schemes with the important difference that L4 neurons receive a substantial amount of feedback in our network, due to the statistical assignment of synapses to their dendrites in the supragranular layers. As a consequence of the decreasing cell densities combined with a constant synaptic volume density, the average indegrees increase along the visual hierarchy from 3950 synapses per neuron in V1 to 14,000 synapses per neuron in area TH. The overall mean indegree across the entire network is 9500.
The resulting connectivity is modified according to an analytical procedure [37] that integrates fundamental constraints on cortical dynamics into the model while maintaining global stability of the system’s attractors. Details are given in the Results section.
While the number of synapses is population-, layer- and area-specific, we make a simple choice for the synaptic weights to restrict complexity: Individual weights are drawn from Gaussian distributions with mean J and width 0.1 J for excitatory synapses, and mean −g J and width 0.1 g J for inhibitory synapses, where g is the relative inhibitory synaptic strength (Table 2). There is no comprehensive knowledge about how the strength of cortico-cortical synapses differs from the strength of intra-areal synapses. Therefore, we introduce two factors χ and to scale the weight of synapses onto excitatory neurons by χ and onto inhibitory neurons by
. For χ = 1, we also choose
, so that cortico-cortical synaptic weights are identical to those of intrinsic synapses. For χ > 1, we choose
, which balances the increased excitation in the system (see Results). Delays are also drawn from Gaussian distributions with fixed mean delays inside areas and distance-dependent cortico-cortical delays (cf. Table 3).
To simulate the network dynamics, we use leaky integrate-and-fire model neurons with exponentially-shaped postsynaptic currents. We choose single-cell parameter values equal to those used by [36] (cf. Table 3). By using a simple neuron model, we limit the complexity of the simulations to bring out the influence of the network connectivity.
Network simulations
We performed simulations on the JUQUEEN supercomputer [75] with NEST version 2.8.0 [76] with optimizations for the use on the supercomputer which were subsequently released in NEST version 2.12.0 [77]. The simulations use 1024 compute nodes (corresponding to 1 rack of JUQUEEN) with 1 MPI process per node and 64 threads per MPI process. A model instance requires about 2GB of working memory on each compute node and takes about 5 minutes for the creation of the network and approximately 12 minutes per 1 s biological time for propagation of the dynamical state. All simulations use a time step of 0.1 ms and exact integration for the subthreshold dynamics of the leaky integrate-and-fire neuron model [78]. Simulations are run for 100.5 s (χ = 1.9), 50.5 s (χ ∈ [1.8, 2.0, 2.1]), and 10.5 s (χ ∈ [1., 1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 2.5]) biological time discarding the first 500 ms. Spike times of all neurons are recorded, except for the simulations shown in Fig 2A and 2B, where 1000 neurons per population are recorded. The digitized workflow reproducing all results and figures of this work was created in compliance with [79] and is available as Python code from https://github.com/INM-6/multi-area-model. The simulation data presented in this manuscript is available from https://web.gin.g-node.org/maximilian.schmidt/multi-area-model-data.
Upper part: Qualitative sketch of the phase space of the network with two stable fixed points, with low and unrealistically high activity, respectively. Modifying the connectivity according to the stabilization procedure of [37] shifts the stable low-activity fixed point towards higher rates without decreasing its global stability. The basins of attractions of the two attractors are divided by the separatrix, a sub-manifold containing one or more unstable fixed points. Lower part, upper panels: firing rates encoded in color. 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 spikes/s in gray. The color scale is identical for all three panels. Lower panels: histograms of population-averaged firing rates for excitatory (blue) and inhibitory (red) populations. The horizontal axis is split into linear- (left) and log-scaled (right) ranges. Simulation parameters: (A) g = 16, νext = 10 spikes/s, κ = 1, (B) g = 16, νext = 10 spikes/s, κ = 1.125, (C) g = 11, νext = 10 spikes/s, κ = 1.125 with the modified connectivity.
Analysis methods
Instantaneous firing rates are determined as spike histograms with bin width 1 ms averaged over the entire population or area. In Figs 3G and 5G we convolve the histograms with Gaussian kernels of optimal width using the method of [80], implemented in the Elephant package [81]. Spike-train irregularity is quantified for each population by the revised local variation LvR [82] averaged over a subsample of 2000 neurons. The cross-correlation coefficient is computed with bin width 1 ms on single-cell spike histograms of a subsample of 2000 neurons per population with at least one emitted spike per neuron. Both measures are computed on the entire population if it contains fewer than 2000 neurons. To compare the simulated with the experimental power spectrum in Fig 6K, we use the simulated spiking data from 140 neurons (equal to the number of neurons identified in the experimental data), distributed across populations in V1 in proportion to the population sizes. We compute the power spectrogram and power spectral densities using Welch’s method (signal.welch of the Python SciPy library [83] with a ‘boxcar’ window, segment length of 1024 data points and 1000 overlapping points between segments). To make our results as comparable as possible with [38], we follow these authors and disregard neurons with an average spiking rate < 0.56 spikes/s.
(A-C) Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extreme observations within 1.5×IQR beyond the IQR. (D) Population-averaged firing rates. (E) Average pairwise correlation coefficients of spiking activity. (F) Irregularity measured by revised local variation LvR [82] averaged across neurons. (G) Area-averaged firing rates, shown as raw binned spike histograms with 1 ms bin width (gray) and convolved histograms, with a Gaussian kernel (black) of optimal width [80].
We employ analytical mean-field theory to predict the stationary population-averaged firing rates of the model. In the diffusion approximation, which is valid for large numbers of sufficiently independent inputs with small synaptic weights, the dynamics of the membrane potential V and synaptic current Is of the leaky integrate-and-fire model neurons used in our model are described by [84]
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 at the population level, 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 [84]
(1)
which holds up to linear order in
and where
with ζ denoting the Riemann zeta function [85]. We solve this equation for our high-dimensional network by finding the fixed points of the first-order differential equation [86]
(2)
for different initial conditions ν0 using the continuous-time dynamics framework of NEST [87], which uses the exponential Euler algorithm, with step size h = 0.1, where s denotes a dimensionless pseudo-time. To investigate the local stability of the fixed point, we study the evolution of a small perturbation δν around the fixed point ν* to linear order,
(3)
The perturbation decays to zero if the maximal real value of the eigenvalues of the effective connectivity matrix G, the Jacobian of Φ, is smaller than 1.
To investigate inter-area propagation, we determine the temporal order of spiking based on the location of the extremum of the correlation function for each pair of areas. This measure is chosen to characterize the relative timing of activity fluctuations across areas, as opposed to measures of causal interactions like Granger causality and the directed transfer function [88]. In an analysis of the lag structure of resting-state fMRI, Mitra et al. [10] similarly characterize temporal order using the time-delay matrix derived from the lagged cross-covariance functions. Our method also resembles the assessment of propagation using the relative timing of slow waves in EEG and LFP recordings in different areas [89–91]. In analogy to structural hierarchies based on pairwise connection patterns [92, 93], we look for a temporal hierarchy that best reflects the order of activations for all pairs of areas. This overall characterization of temporal order extracts the essence of the more complex picture provided by the pairwise delays. The hierarchy is based on the cross-covariance function computed between area-averaged firing rates and subsequently convolved with Gaussian kernels with σ = 2 ms to obtain smoother curves. We use a wavelet-smoothing algorithm (signal.find_peaks_cwt of the Python SciPy library [83] with peak width Δ = 5 ms) to detect extrema for τ ∈ [−100, 100] and take the location of the extremum with the largest absolute value as the time lag. To order areas hierarchically, we determine the peak locations τAB of the cross-correlation function for each pair of areas A, B. We then define a 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,
using the optimize.minimize function of the scipy library [83] 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.
In the context of our spiking network model we define functional connectivity (FC) as the zero-time lag cross-correlation coefficient of the area-averaged synaptic inputs, which we approximate as
with the normalized post-synaptic current PSCj(t) = exp[−t/τs], * indicating convolution, synaptic time constant τs, the population firing rate νj of source population j, mean indegree Kij, and mean synaptic weight Jij of the connection from j to target population i containing Ni neurons. The population firing rate νj is defined as the spike histogram with bin width 1 ms averaged over the entire population, thus time t is in discrete increments of 1 ms.
We compute a BOLD signal from the simulated area-averaged synaptic inputs using the Balloon model [94], implemented in the neuRosim [95] package of R. Synaptic inputs IA(t) drive the responses of cerebral blood flow (CBF) f(t) and cerebral metabolic rate of oxygen (CMRO2) m(t) by linear convolutions
These responses then feed into the Balloon model which is characterized by two dynamical variables q(t), v(t):
with τMTT = 3 s, τ = 10 s, α = 0.4, E0 = 0.4. These two variables determine the relative change of the BOLD signal S:
with a1 = 3.4, a2 = 1V0 = 0.03. The parameters are chosen as in [94].
Clusters in the FC matrices are detected by optimizing the modularity of the weighted, undirected FC graph [96]. We use 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 [97]. The clustering of the structural connectivity is performed with the map equation method [98], which can handle directed connections but no negative weights. In this clustering algorithm, an agent performs random walks between graph nodes according to their degree of connectivity and a certain probability of jumping to a random network node. We choose the probability for a certain target node to be selected to be proportional to the outdegree of the connection, and p = 0.15 as the probability of a random jump. The algorithm detects clusters in the graph by minimizing the length of a binary description of the network using a Huffman code.
To investigate causal relations in the network, we compute the conditional Granger causality [99] between pairs of populations. To reduce computational load, we restrict the set of source populations for each target population i to those that form a connection with on average more than 1 synapse per target neuron. A vector autoregressive model (VAR) describes the target firing rate νi(t) based on the firing rates of other populations with a maximal time lag of 25 ms corresponding to the rounded maximal delay between any two areas in the network. For each source population j, we perform two fits: one using the set of all source populations, yielding VAR{j′}, and one using all source populations except j, yielding VAR{j′|j′ ≠ j}. To determine the causal influence j → i, we test whether the residual variances of the two VARs are significantly different using Levene’s test [100], which is more robust against non-normally distributed residuals than the F-test.
To study dominant paths in the network, we construct the weighted and directed gain matrix G with of the network at the population level, where we evaluate the terms
at the simulated population-averaged firing rates of the model with χ = 1.9. We denote the eigenvalues of G by λ and define λmax as the eigenvalue with the largest real part. To reflect the near-criticality of the brain, we perform an element-wise division by the real part of λmax: G′ = G/[Re(λmax)], so that the maximal real part of the eigenvalues λ′ of the resulting matrix G′ is max[Re(λ′)] = 1. This scaling modulates the relative strengths of direct and indirect paths: a larger value of max[Re(λ′)] increases the relative weighting of indirect paths. Subsequently we use the same method as [35] and denote the weight of the edge from population j to i as
. The logarithm of the reciprocal of the weight, dij = log(1/wij), defines the distance between two nodes in the graph so that summing the distances corresponds to a multiplication of the corresponding weights. Next, the Bellman-Ford algorithm [101–103] finds the shortest paths between any two nodes of the graph. This algorithm determines the shortest paths emanating from vertex i on a graph with N vertices in an iterative manner: it initially assigns an infinite path length to all other nodes k of the graph. Then, the algorithm loops through all edges (j, k) of the graph, tests if the path length pij plus the distance of the edge djk is smaller than the currently stored path length pik, and, if so, assigns pik ← pij + djk. By repeating the loop over all edges N − 1 times, the algorithm considers paths of increasing length on every iteration and ultimately uncovers the shortest paths between each pair of vertices. In contrast to Dijkstra’s algorithm Bellman-Ford copes with edges with negative distance values.
Recording of spiking data from [38]
The experimental recordings are described in [38] and are publicly available [104]. The data consist of sorted spike trains from a 64-electrode array implanted into primary visual area V1 of a lightly anesthetized macaque monkey. The array has 8 electrodes, called shanks, with 8 contacts sites per shank, spanning 1.4 × 1.4 mm horizontally and in depth at 200 μm spacing, covering all cortical layers. For the analysis in Fig 6, we used the 15 minutes of spontaneous activity, where no visual stimulation was provided to the animal. To obtain single-neuron spike trains, [38] performed super-paramagnetic clustering [105] on the high-frequency component (400–5000 Hz) of the recorded signal. Details of the experimental procedures are given in [38]. In our analysis, we distinguish between low-fluctuation and high-fluctuation phases, with low vs. high activity in the frequency range up to 40 Hz. We defined these phases from the power spectrum |C(ω)|2 of the spike histogram for all neurons combined at subsequent intervals of 10 s duration and assign the interval to the low-fluctuation phase if , with an empirically determined threshold θ = 0.8 ⋅ 108. This leads to 77 intervals being classified as low-fluctuation and 15 intervals as high-fluctuation.
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 Subcommittee 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 [106]. Briefly, 10 5-min resting-state fMRI scans (TR: 2 s; voxel size: 1 mm isotropic) were acquired from each subject under light anesthesia (1.5% isoflurane). Nuisance variables (six motion parameters as well as the global white matter and CSF signals) were regressed using the AFNI software package (afni.nimh.nih.gov/afni). The global mean signal was not regressed.
The FV91 parcellation was drawn on the F99 macaque standard cortical surface template [47] 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 [43]. Correlation coefficients were Fisher z-transformed and correlation matrices were averaged within animals and then across animals before transforming back to Pearson coefficients.
Results
Refinement of connectivity by dynamical constraints
From a dynamical systems perspective, we define a state of the network as a set of mean firing rates for all populations. An attractor is a state towards which the network tends to evolve for many different initial conditions. Since the network receives stochastic external input, individual neurons fluctuate around their mean firing rate. An attractor is locally stable if all eigenvalues of the effective connectivity matrix, defined as the Jacobian of the population-level transfer function obtained from mean-field theory (Eq (1)), have real values < 1. The global stability of an attractor is assessed by the size of its basin of attraction in phase space. This volume is measured by discretizing the phase space into a grid of initial conditions and defining the global stability of an attractor A as the proportion of initial conditions leading the system to evolve to A. An analysis based on mean-field theory [37] and simulations reveals that across a wide range of configurations of the external input rate νext and the relative inhibitory synaptic strength g, the network possesses a bistable activity landscape with two coexisting locally stable fixed points (Fig 2). In view of the high dimensionality of the system with 254 populations, the bistability in the mean-field theory is found numerically from a pseudo-time integration that yields the stable fixed points [37], in which the set of firing rates for the full set of populations consistently converges to one of two possible states for each combination of νext and g. The simulation results are qualitatively consistent with these mean-field results. We identify the stable fixed points based on the fact that, after a short initial simulation phase (typically ∼100 ms) and regardless of the initial condition, the network settles in either of these states. The first attractor exhibits asynchronous, irregular activity at moderate firing rates except for populations 5E and 6E, which are nearly silent (Fig 2A), while the second features highly synchronized and regular firing with excessive rates (Fig 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 the shortcoming of vanishing infragranular firing rates, we define an additional parameter κ which increases the external drive onto 5E by a factor κ = Kext,5E/Kext compared to the external drive of the other cell types. Since the rates in population 6E are even lower, we increase the external drive onto 6E linearly with κ such that κ = 1.15 results in K6E,ext/Kext = 1.5. However, even a small increase in κ already drives the network into the undesired high-activity fixed point (Fig 2B). The stabilization procedure described by [37] uses mean-field theory to determine the population-averaged firing rates characterizing the fixed points of the system (cf. Eqs (1) and (2)). By linearizing the population dynamics around the fixed points, the technique identifies connectivity components that are most critical to the global stability of the fixed points and yields targeted modifications of the connectivity within the margins of uncertainty of the anatomical data. The resulting average relative change in total indegrees (summed over source populations) is 11.3%. This allows us to increase κ while retaining the global stability of the low-activity fixed point. In the following, we choose κ = 1.125, which gives K6E,ext/Kext = 1.417, and g = 11, νext = 10 spikes/s, yielding reasonable firing rates in populations 5E and 6E (Fig 2C) with sufficient global stability of the low-activity fixed point [37].
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 [37]. In total, the 4.13 million neurons are interconnected via 2.42 ⋅ 1010 synapses in the stabilized model.
Area- and population-specific activity in the resting state
The network displays a reasonable ground state of activity with low spiking rates between 0.05 and 11 spikes/s (Fig 3). Inhibitory populations are generally more active than excitatory ones across layers and areas despite the identical intrinsic properties of the two cell types. This behavior, first found and discussed in detail in [36], is thus caused by the network connectivity which leads to a high excitation-inhibition ratio onto inhibitory cells. Spiking activity is asynchronous irregular across populations. Population activity fluctuates around its stationary point with small amplitude. Pairwise correlations are low throughout the network (Fig 3E). Excitatory neurons are less synchronized than inhibitory cells in the same layer, except for L4. Spiking irregularity is close to that of a Poisson process across areas and populations (Fig 3F). The only exception is population 6E, which features very low firing rates, so that the measure probably suffers from insufficient spiking data in single cells.
Inter-area coupling leads to metastable regime
To control interactions between areas, we scale cortico-cortical synaptic weights onto excitatory neurons by a factor and provide balance by increasing the weights
onto inhibitory neurons by twice this factor,
. For increasing χ, we observe growing fluctuations of the population spiking rates. At χ = 2 and beyond, the network enters a high-activity state at some time point in the simulation, where most populations spike at unrealistically high rates (Fig 4A). Predictions of mean-field theory show that for increasing χ, a growing proportion of initial conditions (in Eq (2)) result in states with increased activity (Fig 4B). We explain this behavior with the global phase space of the model. At any time, there are two stable attractors with basins of attraction divided by the separatrix, a hyperplane in the phase space that contains unstable fixed points. The low-activity fixed point remains locally stable for increasing χ, as determined by the maximal real part of the eigenvalues of the effective connectivity matrix which is below one for all configurations (Fig 4C). At the same time, its global stability, determined by the proportion of initial conditions leading the system to evolve to it, decreases (Fig 4B and 4D). The effect is that fluctuations around the stationary state, which are evident in a stochastic system, let the system approach the separatrix more closely. Close to the unstable fixed points, the dynamics of the system slow down, which causes the rate fluctuations to appear. From χ = 2, the system is likely to enter the high-activity state within a short amount of simulation time.
(A) Time course of spike rates averaged across all neurons in area V1 for increasing χ (graphs from top to bottom). Cortico-cortical synaptic weights onto inhibitory cells are additionally scaled by for χ > 1. (B) Histogram of analytically computed population firing rates resulting from 104 random initial conditions of Eq (2) for same χ as in (A) averaged across all populations of the network. (C) Real part of maximal eigenvalue of the effective connectivity matrix (Eq (3)) for different values of χ. The values stay below 1 for χ ≤ 2.1, i.e., the network stays locally stable. The increasing proportion of initial conditions in panel B leading to the high-activity fixed point shows that at the same time, the global stability of the low-activity fixed point is reduced. (D) Sketch of the network phase space consisting of two stable fixed points and a submanifold with unstable fixed points separating the two basins of attraction (color indicates spike rate). The network activity (gray sphere) fluctuates around the low-activity fixed point (top). For increasing χ (bottom), the network approaches the separatrix, where the fluctuations become slow. If they are large enough, the network state can pass the separatrix and undergo a transition to the high-activity fixed point. (E) Population-resolved firing rates for χ = 1.9. Same display as in Fig 2.
In the following, we choose χ = 1.9 as the parameter configuration where slow fluctuations coincide with a sufficient global stability of the LA fixed point so that the system does not enter the HA fixed point during the simulation. The corresponding activity is irregular with plausible firing rates (Fig 5A–5C). Irregularly occurring population bursts of different lengths up to several seconds (Fig 5G) arise from the asynchronous baseline activity (Fig 5A–5C) and propagate across the network. The time scales of the population bursts arise from network interactions rather than directly reflecting axonal delays or membrane and synaptic time constants, which only cover a range of 100 ∼ 101 ms. 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 (Fig 5D). The overall average firing rate is 14.6 spikes/s, with the inhibitory populations tending to have higher rates than the excitatory populations. 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 (Fig 3E). Unlike in the model without population bursts, 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 tending to fire more irregularly than inhibitory cells (Fig 3F). Higher areas exhibit bursty spiking, as illustrated by the raster plot for area FEF (Fig 5C).
(A-C) Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extreme observations within 1.5×IQR beyond the IQR. (D) Population-averaged firing rates. (E) Average pairwise correlation coefficients of spiking activity. (F) Irregularity measured by revised local variation LvR [82] averaged across neurons. (G) Area-averaged firing rates, shown as raw binned spike histograms with 1 ms bin width (gray) and convolved histograms, with a Gaussian kernel (black) of optimal width [80].
We compare the simulated spiking activity with experimental data from [38], who recorded spiking activity in 140 neurons of macaque primary visual cortex in the spontaneous condition. The experimental activity shows activity phases differing in their low-frequency power (Fig 6A). In the early stage of the recording, the population activity exhibits only small fluctuations (Fig 6B), while in later stages, the population activity fluctuates on different time scales up to the order of a second (Fig 6C). We therefore split the recorded data into low-fluctuation and high-fluctuation phases (see Materials and methods for details), distinguished by their power at frequencies up to 40 Hz (Fig 6E). To compare simulated with recorded power spectra, we compute the spike rates of 140 cells in V1 distributed across populations in proportion to the population sizes (see Materials and methods for details). We compare three different simulations with low fluctuations (χ = 1, Fig 3), meta-stable dynamics (χ = 1.9) and unrealistically high activity (χ = 2.5) in Fig 6D. The power spectral densities (PSD) of the simulations with χ = 1, 2.5 are flat while for χ = 1.9 the PSD clearly reflects the slow oscillations in the spiking activity (cf. Fig 5). We compare the PSD of this simulation with the experimental results. Overall, the simulated activity in V1 to a good approximation reproduces both the spectrum from the entire recording period and that from the low-fluctuation phase, differing mainly in its increased power between 20 and 40 Hz. The sum of squared deviations (SSD) of the logarithmized spectrum from the logarithmized experimental spectrum for the entire recording period is SSD = 44 (χ = 1.9), compared to SSD = 793 for weak cortico-cortical synapses (χ = 1) and SSD = 2180 for strong cortico-cortical synapses (χ = 2.5), showing that this match is unique to the metastable case. For the entire frequency range, the metastable case (χ = 1.9) best matches the low-fluctuation phase (SSD = 42 vs. SSD = 89 for the high-fluctuation phase). At frequencies below 3 Hz, the power spectrum of the simulations closely matches that of the high-fluctuation phase (χ = 1.9: SSD = 3.2 vs. SSD = 6.8 for the low-fluctuation phase). The horizontal stripes in Fig 5I and 5J may to some extent be due to the mixing of spike trains from excitatory and inhibitory neurons, as the spike sorting does not distinguish between these. This interpretation is supported by the fact that the simulated activity across all layers and populations of V1 closely reproduces the broad distribution of spike rates across cells (Fig 5L). The model with weak cortico-cortical synapses has an overrepresentation of near-silent neurons, while that with strong cortico-cortical synapses overrepresents neurons with high firing rates, and both settings lead to flatter spectra than observed experimentally. Thus, the close match between simulated and experimental population spectra and firing rate distributions is specific to the metastable state.
Spike recordings from 140 neurons in primary visual cortex of macaque monkey [38, 104]. (A) Spectrogram of the spike histogram across all neurons. (B) Raster plot of spiking activity for the initial phase, t ∈ [5 s, 8 s]. (C) Raster plot of spiking activity in the later phase, t ∈ [392 s, 395 s]. Neurons are sorted according to depth of the recording electrodes with neurons closest to the surface of the brain at the top. (D) Power spectral density (PSD) of spike histograms for 140 neurons randomly sampled from all populations of V1 for three different networks simulated for T = 10 s with χ = 1(red), χ = 1.9 (gray) and χ = 2.5 (blue) and the power spectral density for χ = 1.9 for a simulation time of T = 100 s (black). (E) Comparison of simulated power spectral density with χ = 1.9 for 140 neurons randomly sampled from all populations of V1 (black) with experimental recording during low-fluctuation (green) and high-fluctuation (purple) phases and the full recording (yellow). Inset: enlargement for frequencies up to 5 Hz. All power spectra were computed using Welch’s method (see Material and methods). (F) Distribution of spike rates across single cells in experimental data (green, purple, yellow) and simulated spike trains across all layers and populations of V1.
Structural and hierarchical directionality of spontaneous activity
To investigate inter-area propagation, we determine the temporal order of spiking (Fig 7A) based on the correlation between areas. We detect the location of the extremum of the correlation function for each pair of areas (Fig 7B) and collect the corresponding time lags in a matrix (Fig 7C). In analogy to structural hierarchies based on pairwise connection patterns [92, 93], we look for a temporal hierarchy that best reflects the order of activations for all pairs of areas (see Materials and methods). The result (Fig 7D) 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. Fig 7D visualizes the consistency of the hierarchy with the pairwise lags: positive (negative) time lags are placed in the upper (lower) triangle of the reordered time lag matrix. To quantify the goodness of the hierarchy, we counted the pairs of areas for which it indicates the reverse order compared to that of the cross-correlation peaks. The number of such violations is 196 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 (Fig 7G). However, activity troughs in frontal areas follow peaks in occipital activity and thus appear last. This predominant feedback propagation occurs despite feedforward connections on average constituting a greater proportion of the connections onto the neurons in the network than feedback connections, indicating that the dynamical state to an important extent determines the effective strength of anatomical connections. In particular, the high firing rates of the excitatory populations in layer 5 compared to those in layer 2/3 enhance the influence of feedback compared to feedforward projections, as feedforward and feedback projections arise predominantly from the supragranular and infragranular layers, respectively. We analyze the eigenspectrum of the effective connectivity matrix (Fig 7E, cf. Fig 4) and find that the most critical eigenvector (whose eigenvalue has the largest real part, marked in red in Fig 7E) has the largest contributions in the areas at the bottom of the temporal hierarchy. To test whether the local structure of the areas, particularly the increased indegree in higher areas, alone predicts the relative instability in higher areas, we perform another test: We compute the maximum real part of the eigenvalues of the gain matrix in isolated areas with mean-field theory, where we replace inputs from other populations by Poisson input with a global rate of νext = 10 Hz. This does not yield a systematic correlation with the position of the areas in the temporal hierarchy. We conclude that areas at the bottom of the temporal hierarchy are the most unstable areas in the network, i.e., fluctuations in these areas are less suppressed than in temporal and occipital areas, and that this is not a local effect but rather caused by the global network structure.
(A) Area-averaged firing rates (color code) for a sample period (horizontal), with areas (vertical) ordered according to the onset of increased activity. (B) Covariance functions of the area-averaged firing rates of V1 with areas V2 (gray) and FEF (light gray), and auto-covariance function of V1 (black). Dashed lines mark time lags detected by a wavelet smoothing algorithm (see Materials and methods). (C) Matrix of time lags of the correlation function for all pairs of areas. Area MDP is neglected (gray matrix elements) because it has only outgoing connections to but no incoming connections from other visual areas presently included in CoCoMac. Areas are ordered according to their architectural type. (D) Same matrix as in C, but with areas ordered according to the temporal hierarchy. (E) Eigenvalues of the effective connectivity matrix (see Materials and methods). Red, eigenvalue with maximal real part. (F) Absolute value of the corresponding critical eigenvector projected onto the populations of the model on logarithmic scale. Areas are ordered according to their position in the temporal hierarchy. (G) Lateral (left) and medial (right) view on the left hemisphere of an inflated macaque cortical surface showing the order in which areas are preferentially activated (color code). Created with the “view/map 3d surface” tool on https://scalablebrainatlas.incf.org.
Emerging interactions resemble 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 [107]. The FC matrix exhibits a rich structure, similar to experimental resting-state fMRI (Fig 8A and 8C, see Materials and methods for details). In addition, we use the Balloon model of [94] to compute a BOLD signal from the area-averaged synaptic inputs (see Materials and methods for details). The resulting matrix displays a similar structure as the FC matrix based on synaptic input currents. Overall, the values tend to be more extreme (closer to +1 or −1). There are several possible reasons for this: Since our network model comprises the visual cortex only and does not consider neuromodulation, it is potentially in a more confined state than a larger, neuromodulated model of this type and may therefore be less noisy than real brains. Furthermore, the spatial convergence and divergence of cortico-cortical connections presumably also contribute to the variability while in the model all cortico-cortical connections emanate from and target 1 mm2 of spatially unresolved microcircuitry. Lastly, the experimental data are averaged over six monkeys and inter-individual variability decreases the average absolute FC values due to a spread between both negative and positive values. In particular, the positive and negative FC values considered separately are larger in individual monkeys than in the averaged matrix.
(A) Simulated functional connectivity (FC) for χ = 1.9 measured by the zero-time-lag correlation coefficient of synaptic input currents. (B) Functional connectivity obtained by applying the Balloon model of [94] to the simulated synaptic input currents. (C) FC of macaque resting-state fMRI (see Materials and methods). The matrix elements are sorted according to simulated clusters determined with the Louvain algorithm [112] (see Materials and methods). (D) Pearson correlation coefficient of simulated FC vs. experimentally measured FC as a function of scaling factor χ for cortico-cortical synaptic weights with (dots) and
(triangle, cf. Fig 3). Red dot, Pearson correlation coefficient between simulated FC based on BOLD signal and experimentally measured FC. Dashed line, Pearson correlation coefficient of structural connectivity vs. experimentally measured FC. Bottom row: Alluvial diagrams [113] showing (colors distinguish clusters) the differences in the clusters for the structural connectivity and the simulated FC (E), and simulated and experimentally measured FC (F). The clusters of the structural connectivity are extracted using the map equation method [98].
In the simulation, frontal areas 46 and FEF are more weakly coupled with the rest of the network compared to the experiment, but the anticorrelation with V1 as also found by [108] is captured by the model (the corresponding entries in Fig 8A are light blue, see also Figs 8B, 5A and 5C). Area MDP sends connections to, but does not receive connections from other areas according to CoCoMac, limiting its functional coupling to the network. The structural connectivity of our model shows higher correlation with the experimental FC (rPearson = 0.34) than the binary connectivity matrices from both a previous [109] and the most recent release of CoCoMac (rPearson = 0.20), indicating the importance of taking into account connection weights. For , areas interact weakly, resulting in low correlation between simulation and experiment (Fig 8D). For increasing weight factor χ (with
), the correlation between simulation and experiment improves. For intermediate cortico-cortical connection strengths, the correlation of simulation vs. experiment exceeds that between the structural connectivity and experimental FC, both for FC based on synaptic currents (rPearson = 0.44) and simulated BOLD signal for χ = 1.9 (rPearson = 0.38, red dot in Fig 8D) indicating the enhanced explanatory power of the dynamical model. To compare the agreement of our simulated FC with the inter-individual variability between the six monkeys used in the experiments, we compute the average correlation across the monkeys to be 0.31. Comparing the simulated FC (based on synaptic currents) to the FC of individual monkeys yields an average correlation of rPearson = 0.39 ± 0.04, showing that the agreement between the simulated and experimental FC reaches the upper limit determined by the inter-individual variability. From χ = 2 on, the network is increasingly prone to transitions to the highly active state (cf. Fig 4), causing the correlation to decrease. Thus, the highest correlation occurs in a state in which the model exhibits metastable dynamics and slow rate fluctuations appear. Such dynamical slowing close to instability has previously been demonstrated in models of cortical resting-state dynamics where the individual areas were described by population rate equations [110] or small numbers of spiking neurons [111].
Louvain clustering [112], an algorithm optimizing the modularity of the weighted, undirected FC graph [96], yields two modules for both the simulated and the experimental data (Fig 8E). We denote the simulated modules by 1S, 2S and the experimental ones by 1E, 2E and compare these dynamical clusters with the community structure of the structural connectivity. In [35], we describe six clusters in the connectivity exposed using the map equation method [98]. Applying the same method to the modified connectivity matrix used for the simulations in the present work yields seven clusters that are similar to the clusters of the original connectivity, and reflect known functional groupings (Fig 8E). Three clusters contain dorsal stream areas and one large cluster gathers ventral stream areas. Furthermore, early visual areas and frontal areas each form separate communities. Ventral area VOT is grouped together with early visual area VP.
The modules exposed by simulation combine these structural clusters. Cluster 1S contains early visual along with ventral and three dorsal regions. Cluster 2S merges parahippocampal with dorsal but also frontal areas. The experimental module 2E comprises early visual areas, ventral area V4, and dorsal areas, while 1E consists of all other areas including also eight dorsal areas. This large degree of cross-over is most likely caused by an underestimation of the interactions between the areas in module 1E: ventral stream, frontal, and parahippocampal areas, and a subset of dorsal areas, since these are more strongly correlated in the experimental than in the simulated data.
In conclusion, our analysis, summarized in Fig 8, shows that the interactions between areas in the network resemble resting-state fMRI data and the agreement is highest when the network is in the metastable state at intermediate cortico-cortical connection strength.
Population-specific cortico-cortical interactions
We investigate the laminar patterns of cortico-cortical interactions by computing conditional Granger causality on pairs of connected populations in different areas. Testing for significance (p < 0.05) with Levene’s test [100] yields the pairs of causally influencing populations. We distinguish area pairs into three different categories based on the relation of their structural types. These categories approximately agree with the commonly used terminology of feedforward (∼high-to-low-type), lateral (∼horizontal) and feedback (∼low-to-high-type) directions. We observe systematic differences between the three different categories (Fig 9A). High-to-low-type interactions preferentially originate in the supragranular layer and target layer 4, while low-to-high-type interactions are controlled by population 5E and target the infragranular layers. Horizontal interactions share features of both other types. We compare the patterns of dynamical interactions with patterns of shortest paths in the structural connectivity between areas (Fig 9B). High-to-low-type interactions resemble the shortest paths, but the other two types show clear differences to structural paths. Dynamical interactions in the horizontal direction are more similar to low-to-high-type interactions while horizontal shortest paths resemble high-to-low-type paths. The low-to-high-type interactions are dominated by population 5E, while the shortest paths suggest a stronger influence of population 6E. The shortest paths in the structural connectivity differ from the results of [35] because we here consider the modified connectivity after stabilization by the method of [37] as well as synaptic weights with an inter-area scaling factor of χ = 1.9 and respect the population-specific firing rates as opposed to [35], where all populations are assumed to have equal activity levels. The detected shortest paths are nonetheless similar to the ones presented in Figure 8 of [35]. The most obvious difference is that paths onto inhibitory populations are significant in the simulated network (Fig 9D), most likely due to the increased synaptic weight of corticocortical projections onto inhibitory populations. Furthermore, in low-to-high-type connections, paths onto L6 are more relevant.
(A) Significant causal interactions between all pairs of areas categorized according to the difference between their architectural types. Arrow thickness indicates the occurrence of significant interactions (p < 0.05, Levene’s test) between the particular populations (cell types and layers displayed as in Fig 1). (B) Population-specific patterns of shortest paths between all pairs of areas. (C) Proportion of significant connections for the three categories and local connections. (D) Proportion of significant impact onto excitatory and inhibitory populations, respectively, for the three different categories. Σ indicates the respective value over all connections. (E) Histogram of relative indegrees. Bars show data for all (light gray) and significant (black) cortico-cortical connections; and for all (pink) and significant (purple) local connections. The triangles indicate the mean values of the four distributions.
Overall, the number of significant connections is low, with approximately 4% of the cortico-cortical connections leading to significant causal interactions (Fig 9C). Locally, one third of the connections carries causal interactions between the populations of an area. This observation reflects the high degree of local connectivity in cortex. The indegrees of connections with significant interactions are higher compared to all connections, for both cortico-cortical and local connections (Fig 9E). Still, there is no strict dependence of causal interactions on high indegrees, since there are weak but causal connections as well as strong, non-causal connections. One reason is that the activity level of the projecting population plays a major role. The connectivity structure alone predicts for instance that 6E plays a dominant role in low-to-high-type projections, but these connections are actually not significant in the simulations. This is caused by the low activity level of 6E rendering these populations only influential in local interactions. Also the activity of the target populations modulates the effective interactions by determining the susceptibility to inputs due to the nonlinearity of the firing rate response curve [114, 115]. In this sense, network simulations are more powerful than studies merely considering structural connectivity because they add these dynamical aspects to the cortico-cortical interactions. In conclusion, causal interactions in the network, summarized in Fig 9, follow laminar patterns that depend on the relative architectural differentiation of areas. This dependence results from a combination of structural connectivity differences and dynamical states of source and target populations.
Discussion
We present simulations of a multi-scale spiking network model of macaque visual cortex relating cortical connectivity to its dynamics. The connectivity [35] is refined by a mean-field-based stabilization procedure [37] incorporating fundamental activity constraints of non-vanishing, non-saturating spiking rates. The network produces a state with single-cell spiking statistics close to those from recordings in macaque V1 along with large-scale interactions resembling inter-area correlation patterns in macaque fMRI. The model predicts that cortex operates in a metastable regime, and exposes layer-specific, hierarchically organized channels mediating inter-area interactions.
Population firing rates are layer-specific, inhibitory rates generally exceeding excitatory ones, in line with experiments [4, 116, 117]. Since excitatory and inhibitory neurons are equally parametrized and excitatory neurons receive equal or stronger external stimulation compared to inhibitory ones, we conclude that the connectivity causes these differences. The contribution of faster intrinsic dynamics of inhibitory neurons [118, 119] merits future investigation. The mean firing rate of the network is not far above the mean spontaneous rate in V1 of alert monkeys [2]. However, common recording methods may miss (near-)silent neurons that would lower overall rates [120]. Future work can address this potential discrepancy.
With sufficiently strong cortico-cortical connections, fluctuations let the system approach an instability where the dynamics drastically slow down [110, 111, 121, 122]. In this metastable regime, the network closely reproduces the spike rate distribution across V1 neurons in lightly anesthetized macaque [38]. This state features variable-length population bursts mediating inter-area interactions, resembling cortical synchronized fluctuations in both spontaneous and stimulus conditions observed in mouse [123], rat [124, 125], and macaque [7, 104]. Overall, our simulations to a good approximation match the power spectrum of [38], matching periods of higher synchrony at low frequencies. Simultaneous recordings of spontaneous activity in V1 and eye movements of a monkey sitting in darkness suggest that higher synchrony reflects drowsiness [126]. The low-frequency fluctuations may represent irregular Up-Down fluctuations produced by a bistable network with transitions between Up and Down states driven by external input to each area rather than arising locally from an adaptation mechanism, consistent with a recent experimental and modeling study [125]. Further work could distinguish such network effects from other sources of low-frequency fluctuations including NMDA and GABAB transmission, neuromodulation, and adaptation.
The pattern of simulated interactions resembles fMRI resting-state activity from the anesthetized macaque. The functional connectivity (FC) based directly on the synaptic inputs provides similar predictions to that computed using the Balloon model, indicating that the low-frequency fluctuations present already at the level of the spiking activity drive these interactions. The agreement between simulation and experiment peaks at intermediate coupling strength in the metastable regime, potentially related to the slight subcriticality of the resting brain [127]. By combining these large-scale dynamics with plausible single-cell spiking statistics and layer-specific communication, our study extends earlier models exhibiting similar behavior with simpler local circuits [18, 128]. We compute a BOLD signal based on the Balloon model of [94] from the simulated spiking activity and find that the resulting functional connectivity resembles the experimental data, but that the absolute values of the functional connectivity are larger than for the experiment and the synaptic input currents of the model. The model may underestimate the level of noise and thereby overestimate the functional connectivity due to the restriction to visual cortex and the lack of neuromodulatory effects, as well as the lack of convergence and divergence in cortico-cortical connections beyond the 1 mm2 scale. On the other hand, the experimental functional connectivity may be underestimated due to averaging of the experimental data across monkeys, which leads to a decrease in the absolute values owing to inter-individual variability.
There is an ongoing debate in the literature regarding the extent to which fluctuations of functional connectivity during the resting state (e.g. [129, 130]) are significant. Our simulations do not yield slow alternations between the clusters in our FC matrix. There are various possible reasons: experimentally observed FC dynamics in the resting state may be largely due to sampling variability and head motion [131]; the limited duration over which we can simulate and save the corresponding spiking data, which becomes gigabytes or even terabytes in size, is still too short to enable FC dynamics to be observed [132]; or FC dynamics may depend on variations in cognitive state [133] or vigilance [130]. Because we do not simulate the rest of the brain nor neuromodulation, such variations in cognitive and vigilance states are not captured.
The model construction starts with the binary connectivity, defines the weighted structural connectivity based on connection densities, and finally yields the functional connectivity from simulations. Comparing with the experimentally observed functional connectivity, we find that each level (binary structure → weighted structure → dynamical simulation) adds explanatory power. In relation to fMRI functional connectivity, the role of anesthesia is worth investigating, as anesthetics influence BOLD responses in complex ways, although low doses as used in the experiment described here may have only mild effects, as shown in rat [134] and marmoset [135].
The model predicts that population bursts propagate stably across multiple areas, predominantly in the feedback direction, because parietal areas are more unstable than temporal and occipital areas, as revealed by linear stability analysis. This occurs despite higher relative indegrees of feedforward connections, indicating the importance of the dynamical state for the effective influence of the structural connectivity. Specifically, the comparatively high firing rate of the layer 5 excitatory neurons enhances the influence of feedback compared to feedforward connections, which mainly originate in supragranular layers. The systematic activation of parietal before occipital areas in the model parallels human EEG findings on information flow during visual imagery [136] and top-down slow-wave propagation during sleep in humans [89, 90] and mice [91]. Our method for determining the order of activations resembles one recently applied to fMRI recordings [10]. It could be extended to distinguish between excitatory and inhibitory interactions like those we observe between V1 and frontal areas or to identify multiple ‘lag threads’ in the network [137].
Granger causality analysis of cortico-cortical interactions reveals layer- and population-specific communication channels that depend on the difference in the architectural types of the areas. Low architectural types correspond to low neuron density and a thin or absent layer 4 (more limbic areas); high architectural types correspond to high neuron density and a pronounced layer 4 (more sensory areas) [50, 53, 54, 138]. In terms of visual processing hierarchies, areas of high architectural type are found at the bottom of such hierarchies while areas of low type constitute the top level of visual hierarchies. Interactions from high to low types, corresponding to feedforward communication, originate in layer 2/3 and target layer 4, while interactions from high to low types, associated with feedback communication, are predominantly mediated by source neurons in layer 5 and target neurons in layers 6 and 4. Thus, layer 5 neurons are the dominant source of feedback interactions, in contrast with an equal division between layers 5 and 6 in terms of anatomical connectivity [35]. This distinction is due to the higher activity in layer 5 than in layer 6 which enters the Granger causality analysis on both source and target side, whereas it influences only the target side in the path analysis of the structural connectivity via the gain of the receiving population. These findings differ from existing theories about predictive coding in cortical microcircuits [139, 140] insofar that feedback signals preferentially reach granular and infragranular neurons rather than supragranular neurons. This suggests a more prominent role of layer 4 because in addition to feedforward signals and intrinsically produced feedback predictions (via the local L2/3→L5/6→L4 pathway), statistical mapping of synapses to target neurons according to dendritic length [51] provides it with direct feedback predictions from higher areas via apical dendrites in upper layers. Incorporating a dual counterstream organization of feedforward and feedback connections [66, 141] would allow a more refined analysis of these laminar interactions.
Comparing these results with the network structure shows that substantial effective communication arises over a range of connection strengths due to the influence of the dynamical state of the populations, but that the weakest connections do not contribute significantly to communication. Our insights open up a new perspective on the significance of weak ties [142], since our results suggest that they can gain relevance if they link highly active populations but their influence should not be overestimated by binarizing the connectivity. Incorporating local structure beyond population-specific connectivity of point neurons would enable studying how this could strengthen the influence of weak connections [143, 144]. Our analysis further stresses the dominance of local communication in cortical dynamics.
In the model, cortico-cortical indegrees are higher onto excitatory than onto inhibitory neurons, but stronger synapses onto inhibitory than onto excitatory neurons are required to maintain stability. This stabilization mechanism is similar to stabilization by inhibition locally in V1 [145], but uses a different way to balance excitation and inhibition. Locally, the activity of inhibitory neurons increases in proportion to the activity of excitatory neurons [6] via a feedback mechanism [146–149]. In the corticocortical interactions in our model, the activity of the excitatory and inhibitory neurons is increased proportionally by a common input. In principle, other features could support stable inter-area activity propagation, such as spike-frequency adaptation of excitatory neurons [150, 151], short-term synaptic plasticity [152–154], more detailed local balance [155], and fine-tuned excitatory and inhibitory pathways [156–162]. However, optogenetic photostimulation experiments in mice provide some support for our prediction: inter-area EPSCs onto parvalbumin-expressing (PV) interneurons were found to be stronger than onto pyramidal neurons in mouse visual cortex, and mean EPSCs per pixel were larger onto PV interneurons at least for feedforward connections [163]. This outbalancing of excitation by inhibition resembles the “handshake” mechanism in the microcircuit model of [36] where interlaminar projections provide stability by their inhibitory net effect. In the model, this stabilization is reflected in the slightly larger proportion of significant interactions onto inhibitory than onto excitatory populations. In contrast, [20] model feedback connections onto layer 5/6 as net excitatory and onto layer 2/3 as net inhibitory, to reflect the hypothesis that the latter convey predictions suppressing feedforward activation. However, feedback connections may be facilitatory for stimuli within the classical receptive field and suppressive outside the receptive field [139, 164–166]. Such details in feedback processing can be integrated and studied in future model versions.
Random connectivity below the population level is unlikely to suffice for the brain to perform its computations. Including higher-order structures such as neuron-level motifs [167] and patchy connections would make the connectivity more realistic and enable the network to support specific activity patterns between groups of neurons inside a population [168], such as synfire chains [169]. For instance, reciprocal connections and multiple synaptic contacts between neurons are overrepresented in cortex compared to random networks [156, 161, 162], and the existence and strength of a synaptic connection between two neurons correlates with the existence of connections with other neurons [156, 159]. Such cell-specific connectivity is likely to influence firing rate distributions, correlations, and propagation of activity; and theoretical studies have shown that they influence the computational capacity of cortical circuits [15, 170]. To determine the large-scale connectivity, we included only connections from CoCoMac and [30]. CoCoMac alone contains information on approx. 78% of all pairs of cortical areas in the scheme of [45] and approx. 66% of all pairs of visual areas. Including the data of [30], the existence or absence of a connection has been established for 76% of pairs of visual areas. It would be interesting to study the influence of additional connections predicted from graph-theoretical analyses [171]. Our model distinguishes only excitatory and inhibitory cells while cortical networks consist of many different cell types [162, 172]. In particular, integrating different types of inhibitory cells and their detailed projection patterns [173–177] would enrich the model dynamics [178, 179]. Furthermore, going beyond the simple single-neuron dynamics used in this study would enable one to study effects of intrinsic neuron dynamics, such as active calcium conductances in dendrites [180, 181], on the network state.
For tractability, the model represents each area as a 1 mm2 patch of cortex. True area sizes vary from ~3 million cells in TH to ∼300 million cells in V1 for a total of ∼8 ⋅ 108 neurons per hemisphere of macaque visual cortex, a model size that with recent advances in simulation technology [28] already fits on the most powerful supercomputers available. Approaching this size would reduce distortions imposed by downscaling [34] and enable a more realistic representation of synaptic convergence.
Overall, our model elucidates multi-scale relationships between cortical structure and dynamics, and can serve as a platform for the iterative integration of new experimental data, the creation of hypotheses, and the development of functional models of cortex.
References
- 1. Ecker AS, Berens P, Keliris GA, Bethge M, Logothetis NK. Decorrelated Neuronal Firing in Cortical Microcircuits. Science. 2010 January;327(5965):584–587. Available from: https://doi.org/10.1126/science.1179867. pmid:20110506
- 2. Gur M, Kagan I, Snodderly DM. Orientation and direction selectivity of neurons in V1 of alert monkeys: functional relationships and laminar distributions. Cereb Cortex. 2005 August;8(15):121–153.
- 3. Barth AL, Poulet JF. Experimental evidence for sparse firing in the neocortex. TINS. 2012;35(6):345–355. pmid:22579264
- 4. Swadlow HA. Efferent Neurons and Suspected Interneurons in Binocular Visual Cortex of the Awake Rabbit: Receptive Fields and Binocular Properties. J Neurophysiol. 1988 April;59(4):1162–1187. Available from: https://doi.org/10.1152/jn.1988.59.4.1162. pmid:3373273
- 5. Neske GT, Patrick SL, Connors BW. Contributions of diverse excitatory and inhibitory neurons to recurrent network activity in cerebral cortex. J Neurosci. 2015 Jan;35(3):1089–1105. Available from: http://dx.doi.org/10.1523/JNEUROSCI.2279-14.2015. pmid:25609625
- 6. Dehghani N, Peyrache A, Telenczuk B, Le Van Quyen M, Halgren E, Cash SS, et al. Dynamic Balance of Excitation and Inhibition in Human and Monkey Neocortex. Scientific Reports. 2016 March;6. Available from: https://doi.org/10.1038/srep23176.
- 7. Belitski A, Gretton A, Magri C, Murayama Y, Montemurro MA, Logothetis NK, et al. Low-frequency local field potentials and spikes in primary visual cortex convey independent visual information. J Neurosci. 2008;28(22):5696–5709. pmid:18509031
- 8. Vincent J, Patel G, Fox M, Snyder A, Baker J, Van Essen D, et al. Intrinsic functional architecture in the anaesthetized monkey brain. Nature. 2007;447(7140):83–86. pmid:17476267
- 9. Fox MD, Raichle ME. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat Rev Neurosci. 2007;8:700–711. pmid:17704812
- 10. Mitra A, Snyder AZ, Hacker CD, Raichle ME. Lag structure in resting-state fMRI. J Neurophysiol. 2014;111(11):2374–2391. pmid:24598530
- 11. Von Stein A, Chiang C, König P. Top-down processing mediated by interareal synchronization. Proc Natl Acad Sci USA. 2000;97(26):14748–14753. pmid:11121074
- 12. van Kerkoerle T, Self MW, Dagnino B, Gariel-Mathis MA, Poort J, van der Togt C, et al. Alpha and gamma oscillations characterize feedback and feedforward processing in monkey visual cortex. Proc Natl Acad Sci USA. 2014 September;111(40):14332–14341. Available from: https://doi.org/10.1073/pnas.1402773111. pmid:25205811
- 13. Bastos AM, Vezoli J, Bosman CA, Schoffelen JM, Oostenveld R, Dowdall JR, et al. Visual areas exert feedforward and feedback influences through distinct frequency channels. Neuron. 2015;85(2):390–401. pmid:25556836
- 14. Hill S, Tononi G. Modeling sleep and wakefulness in the thalamocortical system. J Neurophysiol. 2005;93(3):1671–1698. pmid:15537811
- 15. Haeusler S, Schuch K, Maass W. Motif distribution, dynamical properties, and computational performance of two data-based cortical microcircuit templates. J Physiol (Paris). 2009;103(1–2):73–87.
- 16. Deco G, Jirsa V, McIntosh AR, Sporns O, Kötter R. Key role of coupling, delay, and noise in resting brain fluctuations. Proc Natl Acad Sci USA. 2009;106(25):10302–10307. pmid:19497858
- 17. Cabral J, Hugues E, Sporns O, Deco G. Role of local network oscillations in resting-state functional connectivity. NeuroImage. 2011;57(1):130–139. pmid:21511044
- 18. Deco G, Jirsa VK. Ongoing Cortical Activity at Rest: Criticality, Multistability, and Ghost Attractors. J Neurosci. 2012;32(10):3366–3375. pmid:22399758
- 19. Cabral J, Kringelbach ML, Deco G. Exploring the network dynamics underlying brain activity during rest. Prog Neurobiol. 2014;114:102–131. pmid:24389385
- 20. Mejias JF, Murray JD, Kennedy H, Wang XJ. Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Science Advances. 2016;2(11):e1601335. pmid:28138530
- 21. Lamme VA, Super H, Spekreijse H, et al. Feedforward, Horizontal, and Feedback Processing in the Visual Cortex. Curr Opin Neurobiol. 1998;8(4):529–535. pmid:9751656
- 22. Rao RPN, Ballard DH. Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive field effects. Nat Neurosci. 1999;2(1):79–87. pmid:10195184
- 23. Angelucci A, Levitt JB, Walton EJS, Hupé JM, Bullier J, Lund JS. Circuits for Local and Global Signal Integration in Primary Visual Cortex. J Neurosci. 2002;22(19):8633–8646. pmid:12351737
- 24. Markov NT, Misery P, Falchier A, Lamy C, Vezoli J, Quilodran R, et al. Weight Consistency Specifies Regularities of Macaque Cortical Networks. Cereb Cortex. 2011;21(6):1254–1272. pmid:21045004
- 25. O’Donnell CO, Gonc JT, Portera-Cailliau C. Beyond excitation / inhibition imbalance in multidimensional models of neural circuit changes in brain disorders. eLife. 2017;6:1–28.
- 26. Douglas RJ, Martin KAC. Mapping the matrix: the ways to neocortex. Neuron. 2007;56:226–238. pmid:17964242
- 27. Douglas RJ, Martin KAC. Recurrent neuronal circuits in the neocortex. 2007;17(13):496–500.
- 28. Kunkel S, Schmidt M, Eppler JM, Masumoto G, Igarashi J, Ishii S, et al. Spiking network simulation code for petascale computers. Front Neuroinform. 2014 oct;8:78. Available from: https://doi.org/10.3389/fninf.2014.00078. pmid:25346682
- 29. Bakker R, Thomas W, Diesmann M. CoCoMac 2.0 and the future of tract-tracing databases. Front Neuroinform. 2012;6:30. pmid:23293600
- 30. Markov NT, Ercsey-Ravasz MM, Ribeiro Gomes AR, Lamy C, Magrou L, Vezoli J, et al. A Weighted and Directed Interareal Connectivity Matrix for Macaque Cerebral Cortex. Cereb Cortex. 2014;24(1):17–36. pmid:23010748
- 31. Izhikevich EM, Edelman GM. Large-scale model of mammalian thalamocortical systems. Proc Natl Acad Sci USA. 2008;105(9):3593–3598. Available from: https://doi.org/10.1073/pnas.0712231105. pmid:18292226
- 32.
Preissl R, Wong TM, Datta P, Flickner M, Singh R, Esser SK, et al. Compass: a scalable simulator for an architecture for Cognitive Computing. In: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. SC’12. Los Alamitos, CA, USA: IEEE Computer Society Press; 2012. p. 54:1–54:11. Available from: http://dl.acm.org/citation.cfm?id=2388996.2389070.
- 33.
Schumann T, Erő C, Gewaltig MO, Delalondre FJ. Towards Simulating Data-Driven Brain Models at the Point Neuron Level on Petascale Computers. In: Di Napoli E, Hermanns MA, Iliev H, Lintermann A, Peyser A, editors. High-Performance Scientific Computing: First JARA-HPC Symposium, JHPCS 2016, Aachen, Germany, October 4–5, 2016, Revised Selected Papers. vol. 10164. Springer; 2017. p. 160–169.
- 34. van Albada SJ, Helias M, Diesmann M. Scalability of asynchronous networks is limited by one-to-one mapping between effective connectivity and correlations. PLoS Comput Biol. 2015;11(9):e1004490. Available from: https://doi.org/10.1371/journal.pcbi.1004490. pmid:26325661
- 35. Schmidt M, Bakker R, Hilgetag CC, Diesmann M, van Albada SJ. Multi-scale account of the network structure of macaque visual cortex. Brain Structure & Function. 2018 April;223(3):1409–1435. Available from: https://doi.org/10.1007/s00429-017-1554-4.
- 36. Potjans TC, Diesmann M. The Cell-Type Specific Cortical Microcircuit: Relating Structure and Activity in a Full-Scale Spiking Network Model. Cereb Cortex. 2014 December;24(3):785–806. Available from: https://doi.org/10.1093/cercor/bhs358. pmid:23203991
- 37. Schuecker J, Schmidt M, van Albada SJ, Diesmann M, Helias M. Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome. PLoS Comput Biol. 2017 February;13(2):e1005179. Available from: https://doi.org/10.1371/journal.pcbi.1005179. pmid:28146554
- 38. Chu CCJ, Chien PF, Hung CP. Tuning dissimilarity explains short distance decline of spontaneous spike correlation in macaque V1. Vision Res. 2014 March;96:113–132. Available from: https://doi.org/10.1016/j.visres.2014.01.008. pmid:24486852
- 39. Friston KJ. Transients, metastability, and neuronal dynamics. NeuroImage. 1997;5(2):164–171. pmid:9345546
- 40. Kelso JS. Multistability and metastability: understanding dynamic coordination in the brain. Philos Trans R Soc Lond, B. 2012;367(1591):906–918.
- 41. Hutt A, beim Graben P. Sequences by metastable attractors: interweaving dynamical systems and experimental data. Frontiers in Applied Mathematics and Statistics. 2017;3:11.
- 42. Hudson AE. Metastability of Neuronal Dynamics during General Anesthesia: Time for a Change in Our Assumptions? Front Neural Circuits. 2017;11:58. pmid:28890688
- 43. Shen K, Bezgin G, Hutchison RM, Gati JS, Menon RS, Everling S, et al. Information Processing Architecture of Functionally Defined Clusters in the Macaque Cortex. J Neurosci. 2012;32(48):17465–17476. pmid:23197737
- 44.
Schmidt M, Bakker R, Shen K, Bezgin G, Hilgetag CC, Diesmann M, et al. A spiking network model explains multi-scale properties of cortical dynamics. BMC Neuroscience. 2016;17(Suppl 1). 25th Annual Computational Neuroscience Meeting: CNS 2016.
- 45. Felleman DJ, Van Essen DC. Distributed hierarchical processing in the primate cerebral cortex. Cereb Cortex. 1991;1:1–47. pmid:1822724
- 46. Girard P, Hupé JM, Bullier J. Feedforward and Feedback Connections Between Areas V1 and V2 of the Monkey Have Similar Rapid Conduction Velocities. J Neurophysiol. 2001;85(3):1328–1331. pmid:11248002
- 47. Van Essen DC, Drury HA, Dickson J, Harwell J, Hanlon D, Anderson CH. An integrated software suite for surface-based analyses of cerebral cortex. 2001;8(5):443–459.
- 48. Bojak I, Oostendorp TF, Reid AT, Kötter R. Towards a model-based integration of co-registered electroencephalography/functional magnetic resonance imaging data with realistic neural population meshes. Phil Trans Roy Soc Lond A. 2011;369(1952):3785–3801.
- 49. Nordlie E, Gewaltig MO, Plesser HE. Towards Reproducible Descriptions of Neuronal Network Models. PLoS Comput Biol. 2009 August;5(8):e1000456. Available from: https://doi.org/10.1371/journal.pcbi.1000456. pmid:19662159
- 50. Hilgetag CC, Medalla M, Beul SF, Barbas H. The primate connectome in context: Principles of connections of the cortical visual system. NeuroImage. 2016;134:685–702. pmid:27083526
- 51. Binzegger T, Douglas RJ, Martin KAC. A Quantitative Map of the Circuit of Cat Primary Visual Cortex. J Neurosci. 2004;39(24):8441–8453. Available from: https://doi.org/10.1523/jneurosci.1400-04.2004.
- 52. Beaulieu C, Kisvarday Z, Somogyi P, Cynader M, Cowey A. Quantitative distribution of GABA-immunopositive and-immunonegative neurons and synapses in the monkey striate cortex (area 17). Cereb Cortex. 1992;2(4):295–309. pmid:1330121
- 53. Barbas H. Pattern in the laminar origin of corticocortical connections. Journal of Comparative Neurology. 1986;252(3):415–422. pmid:3793985
- 54. Barbas H, Rempel-Clower N. Cortical structure predicts the pattern of corticocortical connections. Cereb Cortex. 1997;7(7):635–646. Available from: http://cercor.oxfordjournals.org/content/7/7/635.abstract. pmid:9373019
- 55. Dombrowski SM, Hilgetag CC, Barbas H. Quantitative Architecture Distinguishes Prefrontal Cortical Systems in the Rhesus Monkey. Cereb Cortex. 2001;11(10):975–988. pmid:11549620
- 56. Cragg BG. The density of synapses and neurones in the motor and visual areas of the cerebral cortex. 1967;101(4):639–654.
- 57. O’Kusky J, Colonnier M. A laminar analysis of the number of neurons, glia, and synapses in the visual cortex (area 17) of adult macaque monkeys. J Compar Neurol. 1982;210(3):278–290.
- 58. Elston GN, Rosa MG. Pyramidal cells, patches, and cortical columns: a comparative study of infragranular neurons in TEO, TE, and the superior temporal polysensory area of the macaque monkey. J Neurosci. 2000;20(24):RC117:1–5.
- 59. Elston GN. The pyramidal neuron in occipital, temporal and prefrontal cortex of the owl monkey (Aotus trivirgatus): regional specialization in cell structure. Eur J Neurosci. 2003;17:1313–1318. pmid:12670321
- 60. Elston GN, Benavides-Piccione R, Elston A, Manger PR, DeFelipe J. Pyramidal cells in prefrontal cortex of primates: marked differences in neuronal structure among species. Frontiers in Neuroanatomy. 2011;5:2. pmid:21347276
- 61. Stephan KE, Kamper L, Bozkurt A, Burns GAPC, Young MP, Kötter R. Advanced database methodology for the collation of connectivity data on the macaque brain (CoCoMac). Philos Trans R Soc Lond, B. 2001;356:1159–1186.
- 62. Suzuki WA, Amaral DG. Topographic organization of the reciprocal connections between the monkey entorhinal cortex and the perirhinal and parahippocampal cortices. J Neurosci. 1994;14(3):1856–1877. pmid:8126576
- 63. Rockland KS, Pandya DN. Laminar origins and terminations of cortical connections of the occipital lobe in the rhesus monkey. Brain Res. 1979;179:3–20. pmid:116716
- 64. Barnes CL, Pandya DN. Efferent cortical connections of multimodal cortex of the superior temporal sulcus in the rhesus monkey. J Compar Neurol. 1992;318(2):222–244.
- 65. Ercsey-Ravasz M, Markov NT, Lamy C, Essen DCV, Knoblauch K, Toroczkai Z, et al. A Predictive Network Model of Cerebral Cortical Connectivity Based on a Distance Rule. Neuron. 2013;80(1):184–197. pmid:24094111
- 66. Markov NT, Vezoli J, Chameau P, Falchier A, Quilodran R, Huissoud C, et al. Anatomy of hierarchy: Feedforward and feedback pathways in macaque visual cortex. J Compar Neurol. 2014;522(1):225–259.
- 67. Jones E, Coulter J, Hendry S. Intracortical connectivity of architectonic fields in the somatic sensory, motor and parietal cortex of monkeys. J Compar Neurol. 1978;181(2):291–347.
- 68. Morel A, Bullier J. Anatomical segregation of two cortical visual pathways in the macaque monkey. Vis Neurosci. 1990;4(06):555–578. pmid:2278934
- 69. Webster M, Ungerleider L, Bachevalier J. Connections of inferior temporal areas TE and TEO with medial temporal-lobe structures in infant and adult monkeys. J Neurosci. 1991;11(4):1095–1116. pmid:2010806
- 70. Distler C, Boussaoud D, Desimone R, Ungerleider LG. Cortical connections of inferior temporal area TEO in macaque monkeys. J Compar Neurol. 1993;334(1):125–150.
- 71. Suzuki WL, Amaral DG. Perirhinal and parahippocampal cortices of the macaque monkey: cortical afferents. J Compar Neurol. 1994;350(4):497–533.
- 72. Webster MJ, Bachevalier J, Ungerleider LG. Connections of inferior temporal areas TEO and TE with parietal and frontal cortex in macaque monkeys. Cereb Cortex. 1994;4(5):470–483. pmid:7530521
- 73. Perkel DJ, Bullier J, Kennedy H. Topography of the afferent connectivity of area 17 in the macaque monkey: A double-labelling study. J Compar Neurol. 1986;253(3):374–402.
- 74. Seltzer B, Pandya DN. Parietal, temporal, and occipita projections to cortex of the superior temporal sulcus in the rhesus monkey: A retrograde tracer study. J Compar Neurol. 1994;343(3):445–463.
- 75.
Jülich Supercomputing Centre. JUQUEEN: IBM Blue Gene/Q® Supercomputer System at the Jülich Supercomputing Centre; 2015. Available from: http://dx.doi.org/10.17815/jlsrf-1-18.
- 76.
Eppler JM, Pauli R, Peyser A, Ippen T, Morrison A, Senk J, et al. NEST 2.8.0; 2015. Available from: https://doi.org/10.5281/zenodo.32969.
- 77.
Kunkel S, Morrison A, Weidel P, Eppler JM, Sinha A, Schenck W, et al. NEST 2.12.0; 2017. Available from: https://doi.org/10.5281/zenodo.259534.
- 78. Plesser HE, Diesmann M. Simplicity and efficiency of integrate-and-fire neuron models. Neural Comput. 2009 February;21:353–359. Available from: https://doi.org/10.1162/neco.2008.03-08-731. pmid:19431263
- 79. Pauli R, Weidel P, Kunkel S, Morrison A. Reproducing polychronization: a guide to maximizing the reproducibility of spiking network models. Front Neuroinform. 2018;12(46). pmid:30123121
- 80. Shimazaki H, Shinomoto S. Kernel bandwidth optimization in spike rate estimation. J Comput Neurosci. 2010;29(1):171–182. pmid:19655238
- 81.
Yegenoglu A, Davison A, Holstein D, Muller E, Torre E, Hagen E, et al. Elephant 0.4.1. https://github.com/NeuralEnsemble/elephant/releases/tag/0.4.1; 2017.
- 82. Shinomoto S, Kim H, Shimokawa T, Matsuno N, Funahashi S, Shima K, et al. Relating Neuronal Firing Patterns to Functional Differentiation of Cerebral Cortex. PLoS Comput Biol. 2009;5(7):e1000433. pmid:19593378
- 83.
Jones E, Oliphant T, Peterson P, et al. SciPy: Open source scientific tools for Python; 2001. http://www.scipy.org/.
- 84. Fourcaud N, Brunel N. Dynamics of the firing probability of noisy integrate-and-fire neurons. Neural Comput. 2002;14:2057–2110. pmid:12184844
- 85.
Abramowitz M, Stegun IA. Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables. New York: Dover Publications; 1974. Available from: https://doi.org/10.1119/1.15378.
- 86. Wong KF, Wang XJ. A Recurrent Network Mechanism of Time Integration in Perceptual Decisions. J Neurosci. 2006;26(4):1314–1328. pmid:16436619
- 87. Hahne J, Dahmen D, Schuecker J, Frommer A, Bolten M, Helias M, et al. Integration of continuous-time dynamics in a spiking neural network simulator. Front Neuroinform. 2017 May;11:34. Available from: https://doi.org/10.3389/fninf.2017.00034. pmid:28596730
- 88. Kamiński M, Ding M, Truccolo WA, Bressler SL. Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of signicance. Biol Cybern. 2001;85:145–157. pmid:11508777
- 89. Massimini M, Huber R, Ferrarelli F, Hill S, Tononi G. The sleep slow oscillation as a traveling wave. J Neurosci. 2004;24(31):6862–6870. pmid:15295020
- 90. Nir Y, Staba RJ, Andrillon T, Vyazovskiy VV, Cirelli C, Fried I, et al. Regional slow waves and spindles in human sleep. Neuron. 2011;70(1):153–169. pmid:21482364
- 91. Sheroziya M, Timofeev I. Global intracellular slow-wave dynamics of the thalamocortical system. J Neurosci. 2014;34(26):8875–8893. pmid:24966387
- 92. Reid AT, Krumnack A, Wanke E, Kötter R. Optimization of cortical hierarchies with continuous scales and ranges. NeuroImage. 2009;47:611–617. pmid:19398021
- 93. Krumnack A, Reid AT, Wanke E, Bezgin G, Kötter R. Criteria for optimizing cortical hierarchies with continuous ranges. Front Neuroinform. 2010;4(7). Available from: http://www.frontiersin.org/neuroinformatics/10.3389/fninf.2010.00007/abstract. pmid:20407634
- 94. Buxton RB, Uludağ K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain activation. NeuroImage. 2004;23:S220–S233. pmid:15501093
- 95. Welvaert M, Durnez J, Moerkerke B, Verdoolaege G, Rosseel Y. neuRosim: An R package for generating fMRI data. Journal of Statistical Software. 2011;44(10):1–18.
- 96. Newman MEJ. Analysis of weighted networks. Phys Rev E. 2004 Nov;70:056131. Available from: http://link.aps.org/doi/10.1103/PhysRevE.70.056131.
- 97. Rubinov M, Sporns O. Weight-conserving characterization of complex functional brain networks. NeuroImage. 2011;56(4):2068–2079. pmid:21459148
- 98. Rosvall M, Axelsson D, Bergstrom CT. The map equation. 2009;178(1):13–23.
- 99.
Ding M, Chen Y, Bressler SL. Granger causality: basic theory and application to neuroscience. In: Schelter B, Winterhalder M, Timmer J, editors. Handbook of time series analysis: recent theoretical developments and applications. Wiley-VCH; 2006. p. 437–460.
- 100.
Levene H. Robust tests for equality of variances. In: Olkin I, Ghurye SG, Hoeffding W, Madow WG, Mann HB, editors. Contributions to probability and statistics. Essays in honor of Harold Hotelling. Stanford, CA: Stanford University Press; 1961. p. 279–292.
- 101.
Shimbel A. Structure in communication nets. In: Proceedings of the Symposium on Information Networks. Polytechnic Press of the Polytechnic Institute of Brooklyn; 1955. p. 199–203.
- 102.
Ford Jr LR. Network flow theory. DTIC Document; 1956.
- 103. Bellman R. On a routing problem. Quarterly of applied mathematics. 1958;16(1):87–90.
- 104. Chu CCJ, Chien PF, Hung CP. Multi-electrode recordings of ongoing activity and responses to parametric stimuli in macaque V1. CRCNSorg. 2014; Available from: http://crcns.org/data-sets/vc/pvc-5.
- 105. Quiroga RQ, Nadasdy Z, Ben-Shaul Y. Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering. Neural Comput. 2004 Aug;16(8):1661–87. pmid:15228749
- 106. Babapoor-Farrokhran S, Hutchison RM, Gati JS, Menon RS, Everling S. Functional connectivity patterns of medial and lateral macaque frontal eye fields reveal distinct visuomotor networks. Journal of neurophysiology. 2013;109(10):2560–2570. pmid:23446697
- 107. Logothetis NK, Pauls J, Augath M, Trinath T, Oeltermann A. Neurophysiological investigation of the basis of the fMRI signal. Nature. 2001;412:150–157. pmid:11449264
- 108. Hutchison RM, Gallivan JP, Culham JC, Gati JS, Menon RS, Everling S. Functional connectivity of the frontal eye fields in humans and macaque monkeys investigated with resting-state fMRI. J Neurophysiol. 2012;107(9):2463–2474. pmid:22298826
- 109. Shen K, Hutchison RM, Bezgin G, Everling S, McIntosh AR. Network structure shapes spontaneous functional connectivity dynamics. J Neurosci. 2015;35(14):5579–5588. pmid:25855174
- 110. Cabral J, Hugues E, Kringelbach ML, Deco G. Modeling the outcome of structural disconnection on resting-state functional connectivity. NeuroImage. 2012;62(3):1342–1353. pmid:22705375
- 111. Deco G, Ponce-Alvarez A, Mantini D, Romani GL, Hagmann P, Corbetta M. Resting-state functional connectivity emerges from structurally and dynamically shaped slow linear fluctuations. J Neurosci. 2013;33(27):11239–11252. pmid:23825427
- 112. Blondel VD, Guillaume JL, Lambiotte R, Lefebvre E. Fast unfolding of communities in large networks. 2008;2008(10):P10008.
- 113. Rosvall M, Bergstrom CT. Mapping Change in Large Networks. PLOS ONE. 2010;5(1). pmid:20111700
- 114. Grytskyy D, Tetzlaff T, Diesmann M, Helias M. A unified view on weakly correlated recurrent networks. Front Comput Neurosci. 2013;7:131. pmid:24151463
- 115. Helias M, Tetzlaff T, Diesmann M. Echoes in correlated neural systems. New J Phys. 2013;15:023002.
- 116. Fujisawa S, Amarasingham A, Harrison MT, Buzsáki G. Behavior-dependent short-term assembly dynamics in the medial prefrontal cortex. Nat Neurosci. 2008;11(7):823–833. pmid:18516033
- 117. Sakata S, Harris KD. Laminar Structure of Spontaneous and Sensory-Evoked Population Activity in Auditory Cortex. Neuron. 2009;64:404–418. pmid:19914188
- 118. McCormick DA, Connors BW, Lighthall JW, Prince DA. Comparative electrophysiology of pyramidal and sparsely spiny neurons of the neocortex. J Neurophysiol. 1985;54(4):782–806. pmid:2999347
- 119. Cruikshank SJ, Lewis TJ, Connors BW. Synaptic basis for intense thalamocortical activation of feedforward inhibitory cells in neocortex. Nat Neurosci. 2007;10(4):462–468. Available from: http://dx.doi.org/10.1038/nn1861. pmid:17334362
- 120. Shoham S, O’Connor DH, Segev R. How silent is the brain: is there a “dark matter” problem in neuroscience? J Comp Phys. 2006 July;192:777–784. Available from: https://doi.org/10.1007/s00424-002-0831-z.
- 121. Hennequin G, Vogels TP, Gerstner W. Non-normal amplification in random balanced neuronal networks. Phys Rev E. 2012;86:011909.
- 122. Bos H, Diesmann M, Helias M. Identifying Anatomical Origins of Coexisting Oscillations in the Cortical Microcircuit. PLoS Comput Biol. 2016 October;12(10):e1005132. Available from: http://dx.doi.org/10.1371%2Fjournal.pcbi.1005132. pmid:27736873
- 123. Sachidhanandam S, Sreenivasan V, Kyriakatos A, Kremer Y, Petersen CC. Membrane potential correlates of sensory perception in mouse barrel cortex. Nat Neurosci. 2013;16(11):1671–1677. pmid:24097038
- 124. Vyazovskiy VV, Olcese U, Hanlon EC, Nir Y, Cirelli C, Tononi G. Local sleep in awake rats. Nature. 2011;472(7344):443–447. pmid:21525926
- 125. Jercog D, Roxin A, Barthó P, Luczak A, Compte A, de la Rocha J. UP-DOWN cortical dynamics reflect state transitions in a bistable balanced network. eLife. 2017;p. e22425. pmid:28826485
- 126. Hahn G, Ponce-Alvarez A, Monier C, Benvenuti G, Kumar A, Chavane F, et al. Spontaneous cortical activity is transiently poised close to criticality. PLoS Comput Biol. 2017;.
- 127. Priesemann V, Wibral M, Valderrama M, Pröpper R, Le Van Quyen M, Geisel T, et al. Spike avalanches in vivo suggest a driven, slightly subcritical brain state. 2014;8(108):80–96. Available from: http://www.frontiersin.org/systems_neuroscience/10.3389/fnsys.2014.00108/abstract.
- 128. Zhou C, Zemanová L, Zamora G, Hilgetag C, Kurths J. Hierarchical Organization Unveiled by Functional Connectivity in Complex Brain Networks. Phys Rev Lett. 2006 December;97(23):238103. pmid:17280251
- 129. Chang C, Glover GH. Time-frequency dynamics of resting-state brain connectivity measured with fMRI. NeuroImage. 2010;50(1):81–98. pmid:20006716
- 130. Allen E, Damaraju E, Eichele T, Wu L, Calhoun V. EEG signatures of dynamic functional network connectivity states. Brain Topography. 2018;31(1):101–116. pmid:28229308
- 131. Laumann TO, Snyder AZ, Mitra A, Gordon EM, Gratton C, Adeyemo B, et al. On the stability of BOLD fMRI correlations. Cereb Cortex. 2016;27(10):4719–4732.
- 132. Hindriks R, Adhikari MH, Murayama Y, Ganzetti M, Mantini D, Logothetis NK, et al. Can sliding-window correlations reveal dynamic functional connectivity in resting-state fMRI? NeuroImage. 2016;127:242–256. pmid:26631813
- 133. Preti MG, Bolton TA, Van De Ville D. The dynamic functional connectome: State-of-the-art and perspectives. NeuroImage. 2017;160:41–54. pmid:28034766
- 134. Masamoto K, Kim T, Fukuda M, Wang P, Kim SG. Relationship between neural, vascular, and BOLD signals in isoflurane-anesthetized rat somatosensory cortex. Cereb Cortex. 2007;17(4):942–950. pmid:16731882
- 135. Liu JV, Hirano Y, Nascimento GC, Stefanovic B, Leopold DA, Silva AC. fMRI in the awake marmoset: Somatosensory-evoked responses, functional connectivity, and comparison with propofol anesthesia. NeuroImage. 2013;78:186–195. pmid:23571417
- 136. Dentico D, Cheung BL, Chang JY, Guokas J, Boly M, Tononi G, et al. Reversal of cortical information flow during visual imagery as compared to visual perception. NeuroImage. 2014;100:237–243. Available from: http://www.sciencedirect.com/science/article/pii/S1053811914004662. pmid:24910071
- 137. Mitra A, Snyder AZ, Blazey T, Raichle ME. Lag threads organize the brain’s intrinsic activity. Proc Natl Acad Sci USA. 2015;112(17):E2235–E2244. pmid:25825720
- 138. Beul SF, Barbas H, Hilgetag CC. A predictive structural model of the primate connectome. Scientific Reports. 2017;7(43176):1–12.
- 139. Bastos AM, Usrey WM, Adams Ra, Mangun GR, Fries P, Friston KJ. Canonical microcircuits for predictive coding. Neuron. 2012 nov;76(4):695–711. pmid:23177956
- 140. Shipp S. Neural elements for predictive coding. Frontiers in Psychology. 2016;7:1792. pmid:27917138
- 141. Markov NT, Ercsey-Ravasz M, Van Essen DC, Knoblauch K, Toroczkai Z, Kennedy H. Cortical High-Density Counterstream Architectures. Science. 2013;342 (6158). Available from: http://www.sciencemag.org/content/342/6158/1238406.abstract. pmid:24179228
- 142. Goulas A, Schaefer A, Margulies DS. The strength of weak connections in the macaque cortico-cortical network. Brain Structure and Function. 2015;220(5):2939–2951. pmid:25035063
- 143. Anderson JC, Binzegger T, Martin Ka, Rockland KS. The connection from cortical area V1 to V5: a light and electron microscopic study. J Neurosci. 1998;18(24):10525–40. pmid:9852590
- 144. da Costa NM, Martin KAC. How Thalamus Connects to Spiny Stellate Cells in the Cat’s Visual Cortex. Journal of Neuroscience. 2011;31(8):2925–2937. Available from: http://www.jneurosci.org/cgi/doi/10.1523/JNEUROSCI.5961-10.2011. pmid:21414914
- 145. Ozeki H, Finn I, Schaffer E, Miller K, Ferster D. Inhibitory Stabilization of the Cortical Network Underlies Visual Surround Suppression. Neuron. 2009;62(4):578–592. pmid:19477158
- 146. van Vreeswijk C, Sompolinsky H. Chaotic Balanced State in a Model of Cortical Circuits. Neural Comput. 1998;10:1321–1371. pmid:9698348
- 147. Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci. 2000;8(3):183–208. Available from: https://doi.org/10.1023/a:1008925309027. pmid:10809012
- 148. Tetzlaff T, Helias M, Einevoll GT, Diesmann M. Decorrelation of Neural-Network Activity by Inhibitory Feedback. PLoS Comput Biol. 2012 aug;8(8):e1002596. Available from: https://doi.org/10.1371/journal.pcbi.1002596. pmid:23133368
- 149. Helias M, Tetzlaff T, Diesmann M. The correlation structure of local cortical networks intrinsically results from recurrent dynamics. PLoS Comput Biol. 2014;10(1):e1003428. Available from: https://doi.org/10.1371/journal.pcbi.1003428. pmid:24453955
- 150. Chang YM, Luebke JI. Electrophysiological diversity of layer 5 pyramidal cells in the prefrontal cortex of the rhesus monkey: in vitro slice studies. J Neurophysiol. 2007;98(5):2622–2632. pmid:17804579
- 151. Zaitsev AV, Povysheva NV, Gonzalez-Burgos G, Lewis DA. Electrophysiological classes of layer 2/3 pyramidal cells in monkey prefrontal cortex. J Neurophysiol. 2012;108(2):595–609. pmid:22496534
- 152. Varela JA, Sen K, Gibson J, Fost J, Abbott L, Nelson SB. A quantitative description of short-term plasticity at excitatory synapses in layer 2/3 of rat primary visual cortex. J Neurosci. 1997;17(20):7926–7940. pmid:9315911
- 153. Zador AM, Dobrunz LE. Dynamic synapses in the cortex. Neuron. 1997;19(1):1–4. pmid:9247258
- 154. Beierlein M, Connors BW. Short-Term Dynamics of Thalamocortical and Intracortical Synapses Onto Layer 6 Neurons in Neocortex. J Neurophysiol. 2002;88:1924–1932. pmid:12364518
- 155. Xue M, Atallah BV, Scanziani M. Equalizing excitation-inhibition ratios across visual cortical neurons. Nature. 2014 Jul;511(7511):596–600. Available from: http://dx.doi.org/10.1038/nature13321. pmid:25043046
- 156. Song S, Sjöström P, Reigl M, Nelson S, Chklovskii D. Highly nonrandom features of synaptic connectivity in local cortical circuits. 2005;3(3):e68.
- 157. Yoshimura Y, Callaway EM. Fine-scale specificity of cortical networks depends on inhibitory cell type and connectivity. Nat Neurosci. 2005;8(11):1552–1559. pmid:16222228
- 158. Yoshimura Y, Dantzker JLM, Callaway EM. Excitatory cortical neurons form fine-scale functional networks. Nature. 2005;433(24):868–873. pmid:15729343
- 159. Perin R, Berger TK, Markram H. A synaptic organizing principle for cortical neuronal groups. Proc Natl Acad Sci USA. 2011 March;108(13):5419–5424. Available from: https://doi.org/10.1073/pnas.1016051108. pmid:21383177
- 160. Pfeffer CK, Xue M, He M, Huang ZJ, Scanziani M. Inhibition of inhibition in visual cortex: the logic of connections between molecularly distinct interneurons. Nat Neurosci. 2013;16(8):1068–76. pmid:23817549
- 161. Kasthuri N, Hayworth KJ, Berger DR, Schalek RL, Conchello JA, Knowles-Barley S, et al. Saturated reconstruction of a volume of neocortex. 2015 July;162(3):648–661. Available from: https://doi.org/10.1016/j.cell.2015.06.054.
- 162. Markram H, Muller E, Ramaswamy S, Reimann MW, Abdellah M, Sanchez CA, et al. Reconstruction and simulation of neocortical microcircuitry. Cell. 2015 October;163(2):456–492. Available from: https://doi.org/10.1016/j.cell.2015.09.029. pmid:26451489
- 163. D’Souza RD, Meier AM, Bista P, Wang Q, Burkhalter A. Recruitment of inhibition and excitation across mouse visual cortex depends on the hierarchy of interconnecting areas. eLife. 2016;5(September2016):1–19.
- 164. Hupe JM, James AC, Payne BR, Lomber SG, Girard P, J B. Cortical feedback improves discrimination between figure and background by V1, V2 and V3 neurons. Nature. 1998;394(6695):784–787. pmid:9723617
- 165. Nassi JJ, Lomber SG, Born RT. Corticocortical feedback contributes to surround suppression in V1 of the alert primate. J Neurosci. 2013;33(19):8504–8517. pmid:23658187
- 166. Nassi JJ, Gómez-Laberge C, Kreiman G, Born RT. Corticocortical feedback increases the spatial extent of normalization. 2014;8:105.
- 167. Gal E, London M, Globerson A, Ramaswamy S, Reimann MW, Muller E, et al. Rich cell-type-specific network topology in neocortical microcircuitry. Nat Neurosci. 2017;20(7):1004–1013. pmid:28581480
- 168.
Abeles M. Corticonics: Neural Circuits of the Cerebral Cortex. 1st ed. Cambridge: Cambridge University Press; 1991.
- 169. Diesmann M, Gewaltig MO, Aertsen A. Stable propagation of synchronous spiking in cortical neural networks. Nature. 1999;402(6761):529–533. pmid:10591212
- 170. Timme NM, Ito S, Myroshnychenko M, Nigam S, Shimono M, Yeh FC, et al. High-degree neurons feed cortical computations. PLoS Comput Biol. 2016;12(5):e1004858. pmid:27159884
- 171. Jouve B, Rosenstiehl P, Imbert M. A mathematical approach to the connectivity between the cortical visual areas of the macaque monkey. Cereb Cortex. 1998;8(1):28–39. pmid:9510383
- 172. Jiang X, Shen S, Cadwell CR, Berens P, Sinz F, Ecker AS, et al. Principles of connectivity among morphologically defined cell types in adult neocortex. Science. 2015 November;350(6264):aac9462–aac9462. Available from: https://doi.org/10.1126/science.aac9462.
- 173. Ascoli GA, Alonso-Nanclares L, Anderson SA, Barrionuevo G, Benavides-Piccione R, Burkhalter A, et al. Petilla terminology: nomenclature of features of GABAergic interneurons of the cerebral cortex. Nat Rev Neurosci. 2008;9:557–568. pmid:18568015
- 174. Rudy B, Fishell G, Lee S, Hjerling-Leffler J. Three groups of interneurons account for nearly 100% of neocortical GABAergic neurons. Dev Neurobiol. 2011 Jan;71(1):45–61. Available from: http://dx.doi.org/10.1002/dneu.20853. pmid:21154909
- 175. DeFelipe J, López-Cruz PL, Benavides-Piccione R, Bielza C, Larrañaga P, Anderson S, et al. New insights into the classification and nomenclature of cortical GABAergic interneurons. 2013;14(3):202–216.
- 176. Karnani MM, Agetsuma M, Yuste R. A blanket of inhibition: functional inferences from dense inhibitory connectivity. Current Opinion in Neurobiology. 2014;26:96–102. SI: Inhibition: Synapses, Neurons and Circuits. Available from: http://www.sciencedirect.com/science/article/pii/S0959438813002390. pmid:24440415
- 177. Tremblay R, Lee S, Rudy B. GABAergic interneurons in the neocortex: from cellular properties to circuits. Neuron. 2016;91(2):260–292. pmid:27477017
- 178. Yang GR, Murray JD, Wang XJ. A dendritic disinhibitory circuit mechanism for pathway-specific gating. Nature Communications. 2016;7:12815. pmid:27649374
- 179. Lee JH, Koch C, Mihalas S. A Computational Analysis of the Function of Three Inhibitory Cell Types in Contextual Visual Processing. Front Comput Neurosci. 2017;11(28).
- 180. Larkum M. A cellular mechanism for cortical associations: an organizing principle for the cerebral cortex. TINS. 2013 March;36(3). pmid:23273272
- 181. Palmer L, Shai A, Reeve J, Anderson H, Paulsen O, Larkum M. NMDA spikes enhance action potential generation during sensory input. Nat Neurosci. 2014 march;17(3):383–392. pmid:24487231