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 ., 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 ., 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 ., 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, τ m d V d t = - V + I s ( t ) τ s d I s d t = - I s + μ + σ τ m ξ ( t ) ,, 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 ν i = τ r + τ m π ∫ V r − μ i σ i + γ τ s τ m Θ − μ i σ i + γ τ s τ m e x 2 ( 1 + erf ( x ) ) d x ≕ 1 / Φ i ( ν ) μ i = τ m ∑ j K i j J i j ν j + τ m K ext J ext ν ext σ i 2 = τ m ∑ j K i j J i j 2 ν j + τ m K ext J ext 2 ν ext , ( 1 ), which holds up to linear order in τ s / τ m and where γ = | ζ ( 1 / 2 ) | / 2 , 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, ν ˙ ≔ d ν d s = Φ ( ν ) - ν , ( 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 ,, ν i = ν i * + δ ν ˜ i = Φ i ( ν * + δ ν ) = Φ i ( ν * ) + d Φ i d ν δ ν = ν i * + ∂ Φ i ∂ μ i ∑ j d μ i d ν j δ ν j + ∂ Φ ∂ σ i 2 ∑ j d σ i 2 d ν j δ ν j = ν i * + ∑ j G i j δ ν j ⇒ δ ν ˜ i = ∑ j G i j δ ν j δ ν ˜ = G δ ν ., ( 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 ,, f ( A , B ) = h ( A ) - h ( B ) - τ A B ., To determine the hierarchical levels , we minimize the sum of f ( A , B ) over all pairs of areas ,, S = ∑ A , B f ( A , B ) ,, using the optimize . minimize function of the scipy library 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, I A ( t ) = 1 N A ∑ i ∈ A N i | I i ( t ) | = 1 N A ∑ i ∈ A N i ∑ j K i j | J i j | ( ν j * P S C j ) ( t ) ,, 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, f ( t ) = 1 + ( f 1 - 1 ) h ( t - δ t ) * I A ( t ) m ( t ) = 1 + ( m 1 - 1 ) h ( t ) * I A ( t ) with h ( t ) = 1 k τ h ( k - 1 ) !, ( t τ h ) k e - t / τ h τ f = 4 s , τ h = 0 ., 242 τ f , f 1 = 1 ., 5 , ( m 1 - 1 ) = ( f 1 - 1 ) / 2 , δ t = 1 s ., These responses then feed into the Balloon model which is characterized by two dynamical variables q ( t ) , v ( t ) :, d q d t = 1 τ MTT  f ( t ) E ( t ) E 0 - q ( t ) v ( t ) f out ( v , t )  d v d t = 1 τ MTT  f ( t ) - f out ( v , t )  with f out ( v , t ) = v 1 / α + τ d v d t E ( t ) = E 0 m ( t ) f ( t ), with τMTT = 3 s , τ = 10 s , α = 0 . 4 , E0 = 0 . 4 ., These two variables determine the relative change of the BOLD signal S:, Δ S S = V 0  a 1 ( 1 - q ) - a 2 ( 1 - v ) , 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 G i j = K i j | J i j | d μ i d ν j + K i j J i j 2 d σ i 2 d ν j of the network at the population level , where we evaluate the terms d μ i d ν j , d σ i 2 d ν i j 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 maxRe ( λ′ )  = 1 ., This scaling modulates the relative strengths of direct and indirect paths: a larger value of maxRe ( λ′ )  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 g i j ′ ., 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 ., 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 ∫ 0 Hz 40 Hz | C ( ω ) | 2 d ω ≤ θ , with an empirically determined threshold θ = 0 . 8 ⋅ 108 ., This leads to 77 intervals being classified as low-fluctuation and 15 intervals as high-fluctuation ., 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 ., 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 ., 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 ., To control interactions between areas , we scale cortico-cortical synaptic weights onto excitatory neurons by a factor χ = J cc E / J and provide balance by increasing the weights J cc I onto inhibitory neurons by twice this factor , J cc I = χ I χ J = 2 χ J . 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 ., 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 ) ., 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
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 .
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 .
action potentials, medicine and health sciences, diagnostic radiology, functional magnetic resonance imaging, neural networks, nervous system, membrane potential, vertebrates, electrophysiology, neuroscience, animals, mammals, magnetic resonance imaging, primates, brain mapping, computational neuroscience, neuronal dendrites, neuroimaging, old world monkeys, research and analysis methods, computer and information sciences, imaging techniques, monkeys, animal cells, macaque, cellular neuroscience, radiology and imaging, eukaryota, diagnostic medicine, cell biology, anatomy, synapses, physiology, neurons, single neuron function, biology and life sciences, cellular types, computational biology, amniotes, neurophysiology, organisms