Research ArticleNEUROSCIENCE

Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex

See allHide authors and affiliations

Science Advances  16 Nov 2016:
Vol. 2, no. 11, e1601335
DOI: 10.1126/sciadv.1601335


Interactions between top-down and bottom-up processes in the cerebral cortex hold the key to understanding attentional processes, predictive coding, executive control, and a gamut of other brain functions. However, the underlying circuit mechanism remains poorly understood and represents a major challenge in neuroscience. We approached this problem using a large-scale computational model of the primate cortex constrained by new directed and weighted connectivity data. In our model, the interplay between feedforward and feedback signaling depends on the cortical laminar structure and involves complex dynamics across multiple (intralaminar, interlaminar, interareal, and whole cortex) scales. The model was tested by reproducing, as well as providing insights into, a wide range of neurophysiological findings about frequency-dependent interactions between visual cortical areas, including the observation that feedforward pathways are associated with enhanced gamma (30 to 70 Hz) oscillations, whereas feedback projections selectively modulate alpha/low-beta (8 to 15 Hz) oscillations. Furthermore, the model reproduces a functional hierarchy based on frequency-dependent Granger causality analysis of interareal signaling, as reported in recent monkey and human experiments, and suggests a mechanism for the observed context-dependent hierarchy dynamics. Together, this work highlights the necessity of multiscale approaches and provides a modeling platform for studies of large-scale brain circuit dynamics and functions.

  • Large-scale brain model
  • cortical layers
  • neural oscillations
  • Granger causality
  • visual hierarchy
  • gamma
  • alpha


Inasmuch as the primate cerebral cortex is organized hierarchically, it is essential to understand interactions between feedforward (FF) (bottom-up) information processing and feedback (FB) (top-down) signaling, which may mediate brain predictions about the sensory world, attention, behavioral context, and control. With the advance of electrode recording arrays and other techniques, experimental studies in recent years have not only yielded increasing information about dynamical interactions between pairs of cortical areas but also begun to reveal the complex FF and FB signaling flow among many cortical regions in a large-scale system. Here, motivated by new experimental observations, we built an anatomically based large-scale model of the primate cortex endowed with a laminar structure, and used it to elucidate the dynamical interplay between FF and FB signals at the global brain level.

In the macaque visual cortical system, an anatomical hierarchy is characterized by specific laminar patterns of FF and FB interareal projections (1, 2). FF projections tend to originate from supragranular layers and target layer 4 (1, 2), whereas FB projections largely stem from infragranular layers and target supra- and infragranular layers while avoiding layer 4 (1, 2). In addition, the laminar restriction of these projections increases with the hierarchical distance: If two areas are widely separated, then the supragranular origin of FF projections is more pronounced than if the areas are closer to each other in the hierarchy, and the same holds true for the infragranular origin of FB projections (2).

Hierarchical interactions between visual cortical areas occur in both FF and FB directions, with FF interactions transmitting sensory information from lower to higher brain areas and FB interactions conveying top-down modulation to early sensory areas (3). Recently, it has been observed that these interactions are characterized by distinct neural oscillatory patterns: Interactions in the FF direction are associated with enhanced gamma oscillations (30 to 70 Hz), whereas FB interactions are associated with enhanced alpha or low-beta oscillations (8 to 15 Hz) (36). This frequency profile has been observed in macaques using different recording techniques [from multicontact electrodes (3) to electrocorticography grids (5)] and more recently in humans, using magnetoencephalography (6). It is plausible that these enhanced gamma/alpha signatures in FF/FB interactions may be related to the laminar preference of the anatomical projections described above (7). However, it is unclear what the neural circuit principles are through which oscillatory dynamics originating within local microcircuits drive the spatially segregated, frequency-specific interactions within dense large-scale networks.

To address this question, we built a large-scale computational model of laminar cortical microcircuits comprising 30 cortical areas that interact via their long-range cortico-cortical FF and FB pathways. This large-scale cortical network was built on anatomical interareal connectivity data from tract-tracing studies (2, 8, 9). The connectivity data that we present here extend previous studies (2, 10) by including novel data on the parietal subnetwork [lateral intraparietal (LIP) area], yielding a large-scale network of 30 cortical areas distributed across the occipital, temporal, parietal, and frontal lobes. This rich quantitative data set contains information on the directed strength of interareal projections, their laminar origin, and the propagation latency obtained from the wiring distances between areas (11). We incorporated dynamics into the large-scale cortical circuit by modeling each local circuit with interacting excitatory and inhibitory populations so that supra- and infragranular layers differentially exhibit fast (gamma) and slow (alpha-beta) neural oscillatory dynamics, respectively.

The model spans four spatial scales: intralaminar, interlaminar, interareal, and large-scale levels (see Fig. 1 for a schematic representation). At each level, we constrained the circuit using empirical anatomical data, comparing the computed dynamics with published electrophysiological findings and providing novel and valuable insight into cortical circuit mechanisms. The parsimonious model is validated by capturing a wide range of electrophysiological observations, including (i) the predominance of weakly coherent gamma rhythms in supragranular and strongly coherent alpha in infragranular layers (12); (ii) the effects of visual contrast on gamma frequency and power in V1 (1315); (iii) the interlaminar entrainment between supragranular gamma power and infragranular alpha phase (16); (iv) the inverse correlation between alpha power and supragranular firing rates, which associates strong alpha rhythms with local inhibition (17, 18); and (v) the evidence of enhanced gamma and alpha signatures in FF and FB interareal interactions (35).

Fig. 1 Scheme of the large-scale model.

The scheme shows the four levels considered: a within-layer local microcircuit consisting of an excitatory (Excit.; in red) and an inhibitory (Inhib.; in blue) population (upper left), a laminar circuit with two laminar modules (corresponding to supra- and infragranular layers, lower left), an interareal circuit with laminar-specific projections (lower right), and a large-scale network of 30 cortical areas based on macaque anatomical connectivity (upper right). Each level is anatomically constrained, and its dynamics provide insight into different electrophysiological observations in macaques. Only the connections at each level not shown at a lower level are plotted, for clarity.

With laminar circuit models embedded within each area, and interacting through the structured interareal connectivity, our main finding is that the large-scale model captures the laminar-specific FF and FB interactions on a global scale across many areas. It quantitatively predicts the emergence of frequency-dependent functional interactions and its relationship to underlying anatomical connectivity (5, 6). It provides a mechanistic explanation for the emergence of a functional hierarchy among visual cortical areas, as recently observed in macaques (5) and humans (6), and it allows the exploration of the flexible nature of the cortical hierarchy, which can be altered by behavioral context. Our work highlights the importance of multiscale approaches in the construction of large-scale brain models.


The model spans four spatial levels of description, as depicted in Fig. 1: intralaminar, interlaminar, interareal, and large-scale network. Each level is constrained with anatomical and electrophysiological data, and forms the basis upon which the next level is built. We report results of this work across the four levels sequentially.

Intralaminar level

We first consider the local or intralaminar microcircuit, which is the lowest level of the model and may be identified as a set of neurons within a given cortical area and layer. More precisely, we assume a population of pyramidal neurons and a population of inhibitory interneurons at this level. Within and between each population are recurrent and cross connections, respectively (Fig. 2A, top). Local strongly interconnected populations of excitatory and inhibitory neurons are present in both supra- and infragranular cortical layers (19), and their dynamics have been extensively studied using diverse computational approaches (2022). For this level of description, we simulate each laminar subcircuit with a nonlinear firing rate model, of the Wilson-Cowan type (see Materials and Methods), which represents the mean activities of a population of excitatory neurons and a population of inhibitory neurons.

Fig. 2 Local circuit model at the intralaminar level.

(A) Scheme of the local circuit (top), with the excitatory and inhibitory population in red and blue, respectively, and examples of the oscillatory activity for an excitatory-inhibitory circuit in layer 2/3 (middle, in green) and layer 5/6 (bottom, in orange). (B) Power spectrum of the firing rate of an isolated layer 2/3 as a function of input strength to the excitatory population. The spectrum of the spontaneous state (with zero input) has been subtracted in each case to highlight changes induced by the input (see main text). As the input increases (which resembles the effect of increasing the contrast of a visual stimulus), the power of gamma rhythms becomes stronger, as in observations by Henrie and Shapley (13). (C) Effect of the input to the excitatory population on the power spectrum peak (left) and frequency (right) of the oscillations, for an isolated layer 2/3 (top) and an isolated layer 5/6 (bottom).

The local circuits of the supra- and infragranular layers differ in their connectivity and physiology, leading to their distinct dynamics (19, 23). For the supragranular circuits, we constrained parameter values so that populations display a noise-driven gamma (~40 Hz) rhythm (Fig. 2A, middle), as is commonly observed in layer 2/3 (3, 12, 16, 24). For the infragranular circuits, we adjusted parameter values so that the oscillations displayed by the model fall within the alpha (~10 Hz) or low-beta (~15 to 30 Hz) frequency range (Fig. 2A, bottom) (3, 12).

A simple coupled excitatory-inhibitory system as described here is useful for studying the response of early visual neurons to incoming visual stimuli. For instance, recordings of local field potentials show that increasing the contrast of a visual grating enhances the strength of gamma activity in macaque V1 for low- and medium-contrast levels (13, 15). A higher visual contrast has also been associated with higher frequencies of the gamma activity (14), whereas other factors such as the stimulus size, the level of masking noise, and the stimulus orientation likewise have an effect on gamma power and frequency (15).

To test the behavior of the model at the intralaminar level, we simulated the effect of an increase of stimulus contrast on gamma rhythms for layer 2/3 neurons in V1 and compared the outcome to those observed experimentally (13). We modeled the increase of visual contrast as an increase of the input to the excitatory population, because increases in contrast of a stimulus falling within the receptive field of a neuron elicit increases in firing rate (13). As shown in Fig. 2B, higher contrast values lead to stronger gamma rhythms, characterized by a higher gamma peak of the excitatory firing rate power spectrum. Note that, as in the study by Henrie and Shapley (13), the power spectrum corresponding to spontaneous activity (that is, zero input) was subtracted from the curves shown, thereby removing the power-law appearance of the power spectra but nevertheless conserving the effects of the input in a principled fashion. The enhancement of the power spectrum is strongest on the 30- to 60-Hz gamma band, as in the study by Henrie and Shapley [(13), see their Fig. 4. The contrast-mediated enhancement of the high-gamma (>100 Hz) rhythm observed experimentally is not fully captured by our model. This is to be expected given that the broadband power at this range reflects individual spiking activity (2527), and a rate model will fail to reproduce this feature. However, this does not constitute a problem because we are not considering the effects on high-gamma rhythms (that is, larger than 80 Hz) in this study.

A second effect shown in Fig. 2B is a small but consistent shift of the gamma peak toward higher frequencies (from 30 to 40 Hz, approximately) with increasing contrast. Although not observed in the study by Henrie and Shapley (13), this effect of contrast on gamma frequency has been consistently found in more recent studies (15, 27), and the model captures this effect. Figure 2C shows in more detail the influences of increasing the input to the excitatory population on the peak power (left) and frequency (right) of the oscillations, for the case of both an isolated layer 2/3 microcircuit (gamma rhythm, top) and an isolated layer 5 microcircuit (alpha rhythm, bottom). Note that the power and frequency curves saturate at high inputs, in line with evidence of nonlinear effects on V1 gamma power for strong contrast (15, 28).

The local circuit presented here displays, therefore, the same increases of power and frequency of neural oscillations with increasing input that are observed in recent electrophysiological studies of early visual areas. This makes the circuit a good starting point to understand rhythmic interactions in larger neural systems.

Interlaminar level

Having characterized the neural dynamics of an isolated layer, we built a laminar circuit by considering several layers and adding interlaminar projections between them. To investigate the interplay between gamma and alpha/low-beta rhythms, we consider two distinct laminar modules, one for the generation of gamma and one for the generation of alpha rhythms. Gamma oscillations are most prominently found in granular and supragranular layers before propagating to other layers (3, 12, 16), and modeling work suggests that layers 2/3 and 4 may locally generate gamma rhythms (22). On the other hand, slower oscillations in the alpha and low-beta frequency range are most strongly present at infragranular layers, from where they propagate to other layers (3). We therefore assume two laminar modules: a supragranular module (layer 2/3) that displays gamma rhythms and an infragranular module (layers 5 and 6) that displays alpha rhythms.

To couple the modules of both layers for parsimony, we consider a subset of strong anatomical projections in each direction as a first approximation of the interlaminar circuit. Our hypothesis is that, in the context of gamma-alpha rhythm interactions, the dynamics of the circuit may be described reasonably well by uniquely considering these strong projections, although further details add more richness and complexity to the interlaminar dynamics.

Anatomical studies indicate that, in the supragranular-to-infragranular direction, there is a predominant excitatory projection from layer 3 pyramidal neurons to layer 5 pyramidal neurons, which is consistently stronger than most other interlaminar projections (19, 23). In the infragranular-to-supragranular direction, there is a particularly strong projection from layer 5 pyramidal neurons to layer 2/3 interneurons (29, 30). Comparative studies have confirmed these two projections as the strongest ones between layers 2/3 and 5, with strong projections from layer 6 pyramidal neurons to layer 4 inhibitory cells (which later project to layer 2/3 pyramidal cells), further supporting the inhibitory role of infragranular-to-supragranular projections (31, 32). Therefore, we set these two projections (L2/3E to L5E and L5E to L2/3I) as our main interlaminar connections in the model, as shown in Fig. 3A (left).

Fig. 3 Cortical area model at the interlaminar level.

(A) Scheme of the interlaminar circuit (left panel); self-connections within a given population are omitted in the figure for clarity. Interlaminar connections considered in the model correspond to the strongest projections between layer 2/3 and layer 5/6 as found in experimental studies. Right: Power spectrum of layer 2/3 (top) and layer 5/6 (bottom) in the case of uncoupled, isolated layers (in black, for comparison) and interconnected network (green and orange, respectively). A background input of I = 8 was fed into the excitatory population of both layers. (B) Bottom: A set of 30 traces of activity in layer 5/6 (in gray) and their average (in blue). The central peak of each trace was aligned at zero before averaging. Top: A periodogram of layer 2/3 showing the averaged power for a range of frequencies for the same temporal periods as the layer 5/6 traces. We can see the existence of a strong entrainment of gamma power to alpha phase, as in the experimental findings by Spaak et al. (16). Input was I = 6 for supragranular and I = 8 for infragranular excitatory populations. (C) Effect of injecting external current to the excitatory population of layer 5/6 on the layer 2/3 gamma power and (dimensionless) firing rate (top left and right, respectively) and on layer 5/6 alpha power (bottom left). An inverse relationship between supragranular firing rate and alpha power is observed (bottom right), which highlights a possible link of enhanced alpha rhythms with activity suppression.

As a result of the interlaminar coupling, rhythms spread across layers and the oscillatory dynamics in both supra- and infragranular layers present a more dynamically rich profile, compared to the isolated layers. As Fig. 3A (right) shows, the power spectrum of layer 2/3 activity now displays a strong alpha component as a consequence of the oscillatory input coming from layer 5/6. However, the effect of layer 2/3 activity on the dynamics of layer 5/6 is slightly different. Because of the intrinsically slower dynamics of layer 5/6, the fast fluctuations due to the weak gamma rhythm are partially filtered, and the input from layer 2/3 to layer 5/6 is a quasi-constant term that shifts the alpha peak toward higher frequencies.

The alpha-modulated input arriving at layer 2/3 from layer 5/6 also has a powerful impact, leading to the entrainment of gamma and alpha rhythms (Fig. 3B). Here, the slow oscillatory input from layer 5/6 modulates the power of gamma oscillations in layer 2/3, leading to a phase-amplitude coupling (PAC) between layer 2/3 gamma and layer 5/6 alpha rhythms. This PAC phenomenon has been observed in multielectrode recordings in macaques [(16), see their Fig. 2] and constitutes a validation of our model at the interlaminar level.

To further characterize the behavior of the multilaminar circuit, we studied the effects of external input on its dynamics. The case of an external constant input to layer 5/6 is particularly interesting (Fig. 3C). A constant input arriving at the excitatory population of layer 5/6 has two main effects on layer 2/3, due to inputs to the inhibitory neurons in layer 2/3: a decrease in gamma power and a decrease in mean firing rate of layer 2/3 pyramidal cells (top panels in Fig. 3C). In addition, the input to layer 5/6 enhances the infragranular alpha rhythm (bottom left panel), as observed for the isolated layer case in Fig. 2. This reveals a significant negative correlation between supragranular mean firing rate and alpha power (bottom right panel), which strongly supports the idea that alpha rhythms reflect local inhibition of areas not involved in a particular cognitive task (18). In summary, our laminar model of cortical area consisting of two laminar modules (for supra- and infragranular) constrained by anatomical data displays PAC between gamma and alpha rhythms, as observed experimentally, and provides insight into a possible relationship between alpha rhythms and activity suppression.

Interareal level

Extending our description to the interareal level involves characterizing the anatomical projections linking the microcircuits of two or more cortical areas. Here, we assume two areas with a clearly defined hierarchical relationship (for example, areas V1 and V4). FF projections along the visual hierarchy stem from supragranular layers and preferentially target layer 4, which, in turn, projects to layer 2/3 within that area (1, 2). In our model, this is approximated as a projection from layer 2/3 pyramidal neurons in V1 to layer 2/3 pyramidal neurons in V4. By contrast, FB projections originate from infragranular layers (predominantly from layer 6) and target supra- and infragranular layers while avoiding layer 4 (1, 2). Recent findings indicate that FB projections from these layers are more diffuse than FF projections in terms of their targets (2). In our model, we assume that the FB projection stems from layer 5/6 pyramidal neurons in V4 and targets all four populations in V1 (excitatory and inhibitory; layers 2/3 and 5/6). We make connections to pyramidal cells in layer 5/6 comparatively stronger than those targeting pyramidal cells at layer 2/3, in line with anatomical data (2) and recent experimental data suggesting that top-down signals selectively activate infragranular layers in humans (33). This interareal circuit motif is displayed in Fig. 4 (A and D).

Fig. 4 Microstimulation at the interareal level.

(A) Scheme of the interareal projections between two areas (V1 and V4); the anatomical hierarchy ascends from left to right. We inject a current of I = 15 at both supra- and infragranular excitatory populations of V1 and measure at V4. In addition to this input, a background current to excitatory populations in supragranular (I = 2) and infragranular (I = 4) layers in both V1 and V4 is injected to guarantee a minimum level of activity. (B and C) Power spectrum at V4 measured at layer 2/3 (B) and layer 5/6 (C), for resting and stimulation conditions. Insets show the peak value of the power spectrum at supragranular (B) and infragranular (C) layers for the same resting and stimulation conditions. A large increase in gamma power is found, as in the microstimulation experiments by van Kerkoerle et al. (3). (D) Same as (A), but injecting a current of I = 15 at all excitatory populations of V4 and recording in V1. In addition, a background current of I = 1 to all excitatory areas in V1 and V4 is injected to guarantee a minimum of activity. (E and F) Power spectrum at V1 layer 2/3 (E) and layer 5/6 (F) for resting and stimulation conditions. Inset shows the peak value of the power spectrum for these conditions. A large increase in alpha power and a decrease in gamma power were found, in agreement with experimental observations. For (B) and (F), the blue curve corresponds to an isolated area receiving the same input as in the stimulation case, but without its rhythmic component (see main text for details). n.s., not significant. ***P < 0.001.

Interareal interactions in the visual areas have been recently studied in the context of visual attention (35). Of particular interest for informing circuit mechanisms, van Kerkoerle et al. (3) used microstimulation to reveal an enhancement of gamma/alpha power in FF/FB interactions. Using a model of two multilaminar areas, we can study these interactions and explore the underlying mechanisms involved.

We first simulated a purely FF communication (Fig. 4A): An injection of current into excitatory cells in V1 produced a highly significant increase in gamma power in layer 2/3 of V4 (Fig. 4B) and a small, nonsignificant decrease in alpha power in layer 5/6 of V4 (Fig. 4C). These effects closely resemble the results of van Kerkoerle et al. [(3), see their Fig. 8A] and support the notion that strong gamma rhythms reflect FF communication. According to our model, this gamma increase is explained by two factors: (i) The extra input arriving at V4 layer 2/3 leads to a mean-driven increase in gamma power similar to that observed in the local circuit (see also Fig. 2B), and (ii) the gamma modulation of the arriving input further enhances the intrinsic V4 gamma rhythm via interareal synchronization. By considering V4 as an isolated multilaminar area and injecting into its layer 2/3 a constant input equivalent to the one received from V1 in the two-area case, we can compare both effects (Fig. 4B). This reveals that a constant input enhancing the noise-driven gamma rhythm alone accounts for a major part of the observed increase in gamma power, although the effect, exclusively due to interareal gamma synchronization (that is, the difference between the blue and bright green curves in Fig. 4B), makes a substantial contribution.

To simulate a purely FB communication, we injected current into excitatory cells in V4 and measured responses in V1 (Fig. 4D). This led to a drastic decrease in V1 layer 2/3 gamma power and a strong and significant increase in V1 layer 5/6 alpha power (Fig. 4, E and F). This behavior can be observed in the experimental microstimulation experiment of van Kerkoerle et al. [(3), see their Fig. 8E] and suggests that strong alpha oscillations are associated with FB interactions (along with gamma rhythm suppression). This mechanistic explanation combines several factors: The alpha-modulated input arriving from V4 increases the alpha power in V1 layer 5/6 via (i) a larger average input (as in Fig. 2C) and (ii) the synchronization of both areas in the alpha range. As for FF, we can isolate area V1 and inject the equivalent constant input. This shows that the increase of alpha power due to factor (i) is much larger than the increase due to (ii). The decrease of gamma power is also observed experimentally in the study by van Kerkoerle et al. (3) and agrees with our findings at the interlaminar level (Fig. 3C). This decrease is explained by the net inhibitory effect that V1 layer 5/6 exerts on V1 layer 2/3 plus the FB from V4 to layer 2/3 inhibitory cells, which overcomes the small excitatory contribution from V4 FB to V1 layer 2/3. Other possible scenarios, such as a strong FB connection projecting to supragranular layers, can also be considered in the model (see fig. S1), and they provide interesting insights into the context of top-down attention signals (see Discussion).

The last scenario that we consider here is bidirectional communication, that is, stimulating both areas V1 and V4 and analyzing the frequency-specific profile of the signals in each direction (Fig. 5A). The recorded signal for a given area is defined as a weighted combination of activity in layers 2/3 and 5/6, mimicking depth electrode recording (see Supplementary Methods for details). By computing the spectral coherence between activity at V1 and V4, we observe two clear peaks located at the alpha/low-beta and gamma ranges (Fig. 5B), in agreement with experimental recordings by van Kerkoerle et al. [(3), see their Fig. 7B] and by Bastos et al. [(5), see their Fig. 3D]. We can further obtain information on directionality of signal flow by computing a frequency-dependent Granger causality (GC) analysis on both FF and FB directions. The results (Fig. 5C) show excellent agreement with experimental observations by van Kerkoerle et al. [(3), see their Fig. 7D] and by Bastos et al. [(5), see their Fig. 3B], supporting the notion that enhanced gamma and alpha/low-beta rhythms are associated with FF and FB communication, respectively. This result can also be easily quantified by defining, as in the study by Bastos et al. (5), the directed asymmetry index (DAI) between two areas as the (normalized) difference between spectral GC fluxes, that isEmbedded Image(1)

Fig. 5 Frequency-specific FF and FB interactions.

(A) We now inject current and record the activity of both areas (V1 and V4). An input current of I = 8 was injected in all excitatory populations of the circuit. (B) Spectral coherence between V1 and V4 activity highlights the existence of two peaks at the alpha and gamma range, respectively. (C) Spectral pairwise GC in the V1-to-V4 (green) and V4-to-V1 (orange) directions, showing that each of the peaks found in (B) corresponds to a particular direction of influence, suggesting that frequency-dependent GC analysis could be used to deduce FF versus FB signaling directionality between areas for which the hierarchical positions are not known anatomically. (D) The DAI profile of the functional connection, which is obtained by normalizing the difference between the two GC profiles in (C), can be used to characterize a directed functional connection between two cortical areas.

Note that DAIV1→V4(f) = −DAIV4→V1(f). The DAI profile, which defines a frequency-specific directed functional connection, is shown for the example areas V1 and V4 in Fig. 5D.

In summary, our model of two interconnected cortical areas reveals frequency-dependent signatures in each direction, as observed experimentally. Because of a paucity of available data, we adjusted in the model the relative weights of connections from a FB projection onto the supra- and infragranular layers. We found that the projection should predominantly target the infragranular population in order to reproduce the recent electrophysiological observations on the frequency-dependent FF versus FB signaling. Combined with the connection from infragranular pyramidal cells to supragranular interneurons, this result suggests a net inhibitory influence of top-down signals, which can potentially be compared with a bottom-up signal via FF excitation in the supragranular circuit.

Large-scale level

To develop our interareal system into a large-scale model, we make use of anatomical connectivity data from an ongoing tract-tracing study in macaques (2, 8, 11). In brief, a retrograde tracer injected into a given target area labels neurons in multiple source areas that directly project to the target area. The proportion of labeled neurons in a given source area defines a weight index as the fraction of labeled neurons (FLN) from that source to the target area. In addition, the number of labeled neurons located in the supragranular layer of a given source area (over the total number of labeled neurons in that area) defines the lamination index SLN (for fraction of supragranular layered neurons for a given pathway from source to target area). Source areas that are lower (higher) than the target area in the anatomical hierarchy tend to have a progressively higher (lower) SLN. That is, the lower (higher) the source area relative to the target area, the higher (lower) the SLN values of the source-to-target projection (see Fig. 6A). For instance, purely FF or FB projections have SLN values of 1 and 0, respectively. By repeating the process using other anatomical areas as target areas, we obtained an anatomical connectivity data set with weighted directed connections and laminar specificity (figs. S2 and S3). Elsewhere, we have shown that SLN captures the hierarchical distance so that individual SLN values allocate individual areas in the appropriate ranking according to the Felleman and Van Essen 1991 model (1, 2).

Fig. 6 Large-scale cortical network and functional hierarchy.

(A) Illustration of the anatomical tract-tracing technique used to obtain the anatomical large-scale network and, in particular, the fraction of supragranular labeled neurons (or SLNs, see main text for details). A high (low) value of SLN for a given projection indicates that the source area is lower (higher) than the target area (the injected area) in the anatomical hierarchy. (B) 3D plot of the macaque anatomical network obtained (only projections with FLN of >0.005 are plotted, for clarity), with all 30 areas in their spatial positions. Connection strength is indicated by line width. (C and D) alpha (C) and gamma (D) power for eight selected cortical areas of interest (V1, V2, V4, DP, 8m, 8l, TEO, and 7A). (E) Correlation between SLN and DAI, as a function of frequency. The correlation is positive in the gamma range and negative in the alpha range, indicating a prevalent involvement of these rhythms in FF and FB interactions, respectively. (F) Correlation between SLN and the combined DAI across gamma (30 to 70 Hz) and alpha/low-beta (6 to 18 Hz) frequency ranges (named mDAI, see text for details). (G) Functional hierarchy emerging from the frequency-specific interactions in the network and computed using the mDAI values as in Bastos et al. (5). (H) Areas belonging to the same type (early visual, ventral, dorsal, or frontal; indicated by color box) tend to be clustered in the same way as in the experimental observations. For all panels, visual input was simulated with an input current I = 8 to the supragranular excitatory population of V1, and in addition to this input, a background current of I = 6 to all excitatory populations in the network was also considered.

Here, we use the recently acquired anatomical connectivity data set (2). In addition, we have performed additional injections to expand the connectivity data and include the parietal area LIP in the connectivity data set. The new connectivity, displayed as a three-dimensional (3D) network in Fig. 6B, has 30 cortical areas with a graph density of around 66% (10). To extend our interareal model to a large-scale one, we set a multilaminar circuit at each node of this anatomical network and use (i) FLN values as an estimate of the weight or strength of the interareal connections, (ii) the estimated wiring distances between injection sites to model axonal propagation delays between cortical areas (fig. S4), and (iii) the SLN values as a measure of hierarchical distance between areas. Note that, following previous studies (2, 34), the anatomical hierarchy of the network can be obtained by using a generalized linear model (see fig. S5).

To simulate an increased intensity of the visual stimulus, we increase the level of constant current arriving at V1 layer 2/3 neurons and investigate the resulting interactions among cortical areas and their rhythmic activity. Figure 6 (C and D) shows the alpha and gamma peak power observed in the model, for a subset of cortical areas (V1, V2, V4, TEO, 8m, 8l, DP, and 7A; their anatomical relations are shown in fig. S6), whose interactions have been analyzed in the context of visual attention (5). We observe a relative increase in the alpha rhythms in areas of the ventral stream (V4 and TEO). We also find an especially strong gamma rhythm in early visual areas (V1 and V2) due to sensory input, in agreement with observations of sensory-driven gamma rhythm enhancements being prominent in early visual cortex (6). To characterize the frequency-specific interactions between areas, we compute spectral pairwise conditional GC profiles between the subset of areas, as in the study by Bastos et al. (5). For all area pairs, we observe a relationship between FF/FB interactions and gamma/alpha rhythms (fig. S7) similar to that found earlier with our two-area model (Fig. 5C).

Although this large-scale network assumes that FB projections preferentially target excitatory neurons in layer 5/6 rather than those in layer 2/3, other configurations can also be investigated with our model. It is possible to assume, for instance, that shorter FB projections preferentially target infragranular layers, whereas longer FB projections tend to target supragranular layers (2, 35). This allows the exploration of possible gamma enhancements in visual areas due to top-down signals from the frontal eye field (FEF) (fig. S8) (3638). In the following, and for the sake of simplicity, we will restrict our study to the case in which FB projections of all lengths target infragranular layers preferentially, because distance-dependent rules will only affect a small subset of FB projections and our large-scale simulation results are not strongly affected by these considerations.

From the spectral GC for all areas as displayed in fig. S7, we obtain the spectral DAI profile for each area pair (as in our two-area model), which defines a directed spectral functional connectivity. By computing for each frequency, the correlation between these functional connections and the SLN data (which provide information about the directionality of the anatomical projections), we find a clear pattern: SLN and DAI are positively correlated in the gamma range and negatively correlated in the alpha/low-beta range (Fig. 6E). This agrees with recent experimental findings (5), and it serves as a quantitative demonstration that the strong gamma/alpha signature of FF/FB interactions holds at the level of a dense, large-scale network. To further quantify this correlation, we define the multifrequency DAI (mDAI) per area pair as the averaged DAI for both gamma and alpha ranges, that isEmbedded Image(2)where the negative sign in the second term accounts for the negative correlation of SLN and DAI in the alpha range (see Supplementary Methods). Using this quantity, we obtain a highly significant correlation coefficient (Fig. 6F) between anatomical projections and functional interactions, in accordance with experimental findings (5).

The laminar pattern of anatomical FF and FB projections is the defining feature of the global anatomical hierarchy of visual areas (1, 2). Because our model displays a strong correlation between anatomical projections and functional interactions, it is interesting to test whether the model predicts the emergence of a similar hierarchy at the functional level, as recently observed in vivo (5, 6). We follow the same procedure as in the study by Bastos et al. (5) to define the hierarchical distance between area pairs from mDAI values (see Supplementary Methods), and after simulating the full large-scale model and computing its mDAI values, we observe the emergence of a clear functional hierarchy among visual areas, as shown in Fig. 6G. As in the experimental functional hierarchy (5), early visual areas lie at the bottom of the hierarchy, followed by areas of the FEF (8l and 8m) and with extrastriate visual areas of the ventral and dorsal functional streams at the top. Areas within the same functional streams are clustered around similar hierarchical values (Fig. 6H), in agreement with the study by Bastos et al. (5). These results show that the spectral functional interactions, as well as the formation of a functional hierarchy observed experimentally, can be explained within a computational framework of locally generated rhythms propagated through a multilaminar network structure.

As an interesting example of the predictive power of our large-scale model, we use it to investigate a complex phenomenon observed during visual attention tasks. It has been reported that, contrary to the anatomical hierarchy, the functional hierarchy is not fixed. More precisely, Bastos et al. (5) found that the positions of visual areas in the functional hierarchy are highly dynamic and switch locations in a context-dependent fashion. The ranking of areas in the functionally defined hierarchy in the pre-cue period of the task differs significantly from that obtained in the post-cue period, when top-down modulations from higher areas are expected to heavily influence visual areas at the bottom of the hierarchy. These “hierarchical jumps” were especially noticeable in FEF areas, 8l and 8m, and their origin and implications are as yet unknown.

To provide a mechanistic explanation for jumps in functional hierarchy, we first analyze a simple network of a multilaminar area pair with a well-defined anatomical hierarchical relationship (Fig. 7A). In this scenario, the functional hierarchical distance between both areas is given by the mDAI value of the pair, which can be decomposed into the DAI in the gamma and alpha ranges (see Eq. 2). Considering the spectral DAI in the ascending anatomical direction, the more positive the DAI in the gamma range (and the more negative in the alpha range), the larger the hierarchical distance. Because external input to specific layers modifies the gamma and/or alpha power in a given area, it is plausible that a laminar-specific input could affect the frequency-specific interactions between areas and therefore their hierarchical distance. To test this hypothesis, we measured the frequency-specific interactions between both areas for three types of laminar input patterns to the lower area: (i) a small input to layer 2/3 and a large input to layer 5/6, (ii) an equal input to both layer 2/3 and layer 5/6, and (iii) a large input to layer 2/3 and a small input to layer 5/6. Changes in the laminar-specific inputs have a substantial impact on the DAI profile, increasing the DAI in the gamma range (and therefore the FF interactions) with increases in the supragranular-to-infragranular input ratio (Fig. 7B). These gamma enhancements are due to the combined effect of the increase in direct excitation to layer 2/3 and the decrease of interlaminar inhibition from layer 5/6. Increases of DAI in the gamma range imply a higher mDAI value of the area pair (Fig. 7C), leading to an increase in the hierarchical distance for this simple network.

Fig. 7 Mechanistic explanation for the experimentally observed hierarchical jumps.

(A) Scheme of the simple two-area circuit considered (we use the same parameters as for the two-area microstimulation protocol); the area on the left, which is lower in the anatomical hierarchy, receives laminar-specific input. (B) The gamma component of the DAI increases as the input to layer 2/3 exceeds input to layer 5/6 for the lower area (ascending curves correspond to input to L2/3E increasing from I = 4 to 6 to 8, and input to L5/6E decreasing from I = 8 to 6 to 4). In the higher area, excitatory populations receive a fixed I = 6 background current. (C) Increases in the hierarchy rank of the higher area as a consequence of the laminar-specific input in (B). The laminar specificity S, defined as the difference between the input to L2/3E and to L5/6E, goes from −4 to 0 to 4 in this example. (D) We follow the same procedure but now injecting laminar-specific current into area 8l within the full 30-area network. (E and F) Changes in alpha (E) and gamma (F) band power as a consequence of the injection of laminar-specific input. Lines in gray correspond to the same input injected at both layers (that is, no laminar-specific input), whereas colored lines correspond to a laminar specificity of S = 8. (G) The spectral DAI profile from 8l to 8m increases in the gamma range as a consequence of the laminar-specific input (gray curve, S = 0; blue curve, S = 8). (H) A hierarchical jump of area 8m is observed, as in the two-area case (gray points, S = 0; blue points, S = 8). (I) We find a robust increase of the hierarchical jump distance with the strength of the laminar specificity of the input. Other parameters and background currents are as in Fig. 6.

Extending this analysis to the full 30-area network is not a straightforward process, because cortical areas tend to interact in a nontrivial fashion, and changing the laminar specificity of the input to one area influences interactions with several areas simultaneously. However, we can test our hypothesis in a concrete case that is of particular interest. According to the FLN data (see fig. S2), areas 8l and 8m are strongly connected to each other, whereas their projections to and from other areas are weaker. This provides a useful test bed for translating our findings from the case of two isolated areas to the full cortical network. In terms of laminar origin, projections between areas 8l and 8m are approximately horizontal (that is, SLN is about 0.5 in both directions; see fig. S6). Furthermore, both areas display a high hierarchical mobility according to experimental observations (5), so the pair 8l-8m is a suitable candidate for observing hierarchical jumps in our model.

Because 8l is lower than 8m in the functional hierarchy (for both the model and the experimental findings), we simulate the large-scale network of 30 areas with a highly laminar-specific input to area 8l (strong for layer 2/3, weak for layer 5/6; see Fig. 7D) and compare the DAI profiles and functional hierarchy with the original simulation (for which the input to 8l does not have laminar specificity). This laminar-specific input elevates the gamma power and decreases alpha power in area 8l, as expected (Fig. 7, E and F). Furthermore, Fig. 7G shows an important effect in the spectral DAI between 8l and 8m in this situation: The DAI from 8l to 8m is stronger in the gamma range, as in the two-area case (although, because of the interactions of other areas, the profile is more complex in this case). The main effect in the functional hierarchy (Fig. 7H) is an elevation of area 8m, which effectively increases the hierarchical distance between 8m and 8l as expected and induces a hierarchical jump of area 8m, as experimentally observed in the study by Bastos et al. (5). To further test this effect, we repeat the process for several degrees of laminar specificity of the 8l input, and we observe in Fig. 7I a clear increase in the hierarchical distance between 8m and 8l as the laminar specificity of the input increases (stronger for layer 2/3 and weaker for layer 5/6).

These results indicate a strong prediction of our large-scale model: Jumps observed in the dynamic functional hierarchy may reflect a change in the laminar specificity of the input to cortical areas. Such a change in the laminar pattern of the input could be due, for instance, to context-dependent influences from higher cortical areas, known to be important in attention and other cognitive tasks.


The brain is characterized by interconnectivity and dynamics across multiple scales. Perhaps no recent experimental finding better highlights the challenge raised by this multiscale complexity than the layer-specific, and frequency-dependent, interplay between FF and FB signaling streams across the large-scale primate cortical circuit (36). To uncover the circuit mechanism behind these frequency-dependent processes, and to understand their implications for large-scale communication, we have built a multiscale model of the macaque brain covering both slow (alpha) and fast (gamma) neural oscillatory dynamics and four spatial levels of description constrained by the known anatomy. The spatial levels range from local homogeneous populations to laminar circuits, interareal interactions, and large-scale cortical networks based on precise anatomical macaque connectivity data. The parsimonious model can explain a wide range of macaque electrophysiological observations, including the effects of visual contrast on V1 gamma rhythms (13, 14), the interlaminar PAC between gamma and alpha rhythms (4, 16), the relationship between alpha power and local inhibition (18), gamma/alpha signatures of FF/FB interareal interactions (3, 5), and correlations between anatomical and functional networks (5). Notably, using the same analysis as in the experimental studies, our model captures the emergence of a dynamic functional hierarchy in the visual cortex (5, 6).

The model makes a number of experimentally testable predictions. At the interlaminar level, we suggest that the PAC between gamma and alpha rhythms is mediated by projections from layer 5/6 pyramidal neurons to layer 2/3 interneurons. This sets a clear infragranular-to-supragranular direction of the modulation [as tentatively discussed by Spaak et al. (16)] and highlights a key role of layer 2/3 interneurons in PAC. Physiological experiments have identified an inhibitory subgroup (including chandelier cells and irregular-spiking basket cells) as the recipient of strong projections from layer 5 pyramidal neurons (29, 30). Our model suggests that this subgroup of interneurons will play an important role in PAC between supragranular gamma and infragranular alpha rhythms, a prediction that could be experimentally tested.

The relevance of layer 2/3 interneurons on PAC is expected to have important implications at the functional level. Because enhanced gamma rhythms are associated with a strong drive to higher-order areas (3, 5), modulating the PAC could lead to a modulation of the temporal windows of elevated gamma synchrony. Therefore, top-down signals targeting layer 2/3 interneurons (which are compatible with the unspecific FB projections considered in our model) could act as a top-down gating mechanism, controlling the length and occurrence of the temporal windows for communication (39).

At the interareal level, our model provides insight into the mechanisms that subserve cortico-cortical communication (40). For FF communication, we have identified two mechanisms that may contribute to the experimentally observed enhancement of gamma power (3, 5): (i) mean-driven enhancement of oscillations and (ii) interareal synchronization. Both mechanisms contribute to the observed gamma rhythm enhancement, but due to the existence of the mean-driven enhancement mechanism, a tight synchronization between areas is not necessary to explain the experimental observations. Gamma rhythms tend to be very weakly coherent, and a communication mechanism purely based on interareal gamma synchronization could pose certain problems (14, 41) [but see the study by Fries (40)]. However, the presence of a certain level of interareal synchronization clearly enhances the gamma power associated to FF interactions, and its effect is stronger than the contributions of synchronization to FB interactions.

Although point-to-point FB projections between supragranular layers (2) are not explicitly incorporated into the model, interareal projections between supragranular layers in the FB direction will be present in the model once we incorporate the anatomical data into the framework. For instance, a projection with an SLN value of 0.2 will be mainly a FB projection, but because the SLN value is larger than zero, the projection will have a (comparatively weak) supragranular-to-supragranular component to reflect the small fraction of supragranular neurons anatomically identified for that specific projection. Because weak but nonzero values of SLN are common in our anatomical data set, this tends to be the norm rather than the exception, and the effect is the appearance of small peaks of GC at gamma ranges for projections that are generally considered as FB (fig. S7). This could also explain the existence of gamma rhythms associated with particular FB pathways, such as top-down input from FEF, which lead to an increase of gamma rhythms in early visual areas (see also further considerations in the next subsection) (3638, 42). However, a more careful consideration of point-to-point FB projections from supragranular layers is beyond the scope of this work, because it would involve extending the model of cortical area to include spatial extension and lateral connectivity within nearby neurons. This extension will be considered in future studies.

The anatomical connectivity data used at the large-scale level contains detailed information about 30 cortical areas of the macaque brain (2, 8). Other brain areas, including thalamic and other subcortical structures, are currently missing from the data set and have therefore not been explicitly considered in the present model. Although it is not uncommon to assume, as we do in this study, a stereotypical thalamic input to canonical cortical microcircuits (8), the absence of an explicit thalamic model may be seen here as an opportunity to conceptually explore the role of thalamus in cortico-cortical communication. For instance, a strong sustained input to all areas is needed in our model to guarantee that both gamma and alpha rhythms are sustained in the cortex, and such an input would most likely be provided by the thalamus. In addition, thalamic input could also be partly responsible for modulating the interactions between cortical populations, and therefore influencing the PAC between cortical layers, for instance.

One could assume that dynamic thalamo-cortical interactions, or other mechanisms such as bursting, are necessary to generate cortical rhythms, as is commonly assumed for alpha oscillations (43, 44). Although considering more detailed mechanisms of alpha rhythms is an appealing future direction in this context [and several modeling studies constitute a solid starting point in that direction; see previous studies (4345)], it is unlikely that these mechanisms would significantly alter the general conclusions of the present work. This is because, although further detail can be added into the rhythm generator, model predictions on power and frequency modulations induced by external input are consistent with the experimental observations. As a consequence, the transmission of signals between cortical areas and the modulation of these signals will accurately reflect the electrophysiological evidence available, independently of the particular mechanism for rhythm generation assumed for the model.

In relation to the nature of the slow oscillations, is it still not clear whether FB interactions are mostly mediated by alpha or (low) beta rhythms. Bastos et al. (5) identified FB interactions in the low-beta range, whereas van Kerkoerle et al. (3) and Buffalo et al. (4) consistently found FB interactions in the alpha range. Our model assumes a generic low-frequency oscillation (in the range partially overlapping that of alpha and/or low beta) associated with FB interactions, so the relative importance of each rhythm in FB interactions cannot be inferred from the model at its current state—further assumptions about the rhythm generator would be needed to provide insight into this issue. Theta rhythms, which have also been associated with FF interactions [although not as prominently as with gamma rhythms (5)], could also be considered in this context, because it is possible that multiple rhythms could contribute to communication in each direction: theta and gamma for FF interactions, and alpha and beta for FB ones.

Laminar targets of interareal projections

Synaptic connections at different levels of the model have been established according to known anatomical patterns (2, 8, 19, 23). Full mapping of the model connections with anatomy is not presently possible, in part due to simplifications needed for a computationally tractable large-scale cortical model and to the current limitations in the anatomical literature. In particular, although the laminar origin of interareal projections is quantitatively measured for each pair of areas in our network (2), the laminar target of these projections can only be grossly estimated from other anatomical studies (1, 2, 46). We have assumed, for example, that FF connections arriving to a given area exclusively target their supragranular laminar module, based on the anatomical evidence that layer 4 is the major target of FF projections (1, 2, 46). There is anatomical evidence that, in interareal pathways, supragranular layers target supragranular layers and infragranular layers target infragranular layers (2, 35), which supports this assumption. There is also evidence suggesting that FF projections could also aim at infragranular layers as secondary targets (47). This novel FF pathway would not substantially affect our modeling results, because the addition of such a secondary target for FF projections could be incorporated, as a first approximation, as a small increase in the strengths of the FF projection to layer 2/3 and the interlaminar projection from layer 2/3 to layer 5/6.

The case of the laminar target of FB projections is arguably more interesting. It is thought that FB projections preferentially target supra- and infragranular layers while avoiding layer 4 (2, 35). However, the precise weight of the FB input to infragranular layers in early areas is not known. In our model, we set connections to pyramidal cells in layer 5/6 comparatively stronger than those targeting pyramidal cells at layer 2/3, motivated by anatomical studies (2, 35). As a result, FB projections in our model are net excitatory for layer 5/6 and net inhibitory for layer 2/3 (due to FB projections to interneurons in layer 2/3 plus the FB-activated local ascending inhibition from layer 5/6). Such a connectivity pattern would be consistent with a predictive coding framework, because top-down signals would provide the strong inhibitory input to layer 2/3 necessary to suppress any predicted signals (48, 49). In such a framework, only ascending information that could not be predicted and matched by the top-down input would form prediction errors that can move forward in the hierarchy through subsequent supragranular layers. This important role of strong inhibitory top-down signals in predictive coding, which has been previously discussed (49, 50), would, in principle, be consistent with the general framework of our model. However, the current version of our model cannot explicitly study predictive coding tasks, because the intra-areal modeling level lacks the sort of selectivity dynamics needed for these tasks. Further implementations, based, for instance, on recent model implementations of dynamical causal modeling (51, 52), which is a natural language for predictive coding, could be developed and included in our large-scale framework in future studies.

Despite the general net inhibitory effect of FB projections, the existence of top-down excitatory interactions (46) may potentially play an important role in other contexts, such as in selective attention (53, 54). As discussed by Bastos et al. (49), predictive coding and selective attention are not mutually exclusive, and the potentiation of behaviorally relevant bottom-up signals through top-down attentional modulation is also an important function in a predictive coding scenario. FB projections that can provide a net excitatory contribution to supragranular layers have been anatomically identified (2, 46), and our two-area model can shed light as to its potential role in selective attention. In particular, by considering the existence of net excitatory FB projections to layer 2/3 (and leaving the FB projections to pyramidal neurons of layer 5/6 as a comparatively weaker projection), our model predicts a decrease of infragranular alpha power and an enhancement of gamma rhythms and supragranular firing rate induced by top-down signals (see fig. S1). In this case, top-down input to layer 2/3 significantly enhances gamma power by increasing the level of excitation in the circuit, with interareal synchronization playing only a minor role. These results agree with attention-related V4 activity recorded from macaque, where similar gamma/alpha signatures are observed in receptive fields covering the attended visual stimulus (53, 55).

It is not a straightforward task to introduce these net excitatory FB projections to layer 2/3 in our model, because quantitative anatomical data about the laminar targets of interareal projections are not available. However, we can hypothesize about how to incorporate these projections in a simplified manner. For instance, anatomical studies show that short-range FB projections tend to target infragranular layers, whereas long-range FB projections preferentially target supragranular layers (2, 35, 56). In a set of simulations, we have assumed a simple monotonic relationship between FB projection distance and targeted layer (fig. S8A). This implies that short-range FB projections in our model (such as the FB from V2 to V1) would mainly target infragranular layers, whereas long-range FB projections (such as 8l to V4) would preferentially target supragranular layers. With such a distribution of FB laminar targets, a strong top-down input coming from frontal areas, such as FEF area 8l, leads to an increase in the firing rate and gamma power of early visual areas (see fig. S8, B and C). This agrees with experimental evidence showing that FEF stimulation leads to stronger responses and enhanced gamma interactions in visual areas (3638, 42), and with recent studies showing that top-down attention signals from prefrontal to visual areas enhance gamma rhythms (57) and suppress alpha rhythms (4). As with the case of predictive coding, careful modelization of top-down attentional signals demands a more detailed large-scale model (with selective neural populations) than the one presented here, and therefore, this will have to be addressed in future studies.

Relation to other models

It is useful to discuss how our model compares to other laminar models of a local cortical area. Lee and collaborators (58, 59), for instance, developed a laminar spiking model with three layers and different types of interneurons, which provides interesting insight into the top-down attention mechanisms mentioned above and also into the role of cholinergic neuromodulators in interareal communication. Potjans and Diesmann (32) have also developed an anatomically and physiologically based laminar spiking model of a cortical microcolumn, with connectivity patterns in agreement with our simplified laminar circuit. Our approach of considering a parsimonious model with two laminar modules has several distinctive features and advantages at the conceptual level (for example, it provides a cleaner picture to understand the origin of PAC) and in the context of scaling up the model to large-scale networks. Additional features from the above or similar laminar models can be incorporated in our large-scale framework to further increase the multiscale predictive power of our approach in the future.

At the large-scale level, recent models of the macaque brain have highlighted the importance of heterogeneity to explain the emergence of a hierarchy of time scales (34), as recently observed experimentally (60, 61), although these models lack structure at the laminar level. Including a certain level of heterogeneity in a laminar large-scale model could reveal new mechanisms for interareal communication. It will be worth examining, in future research, the incorporation into our new model of different types of heterogeneity [in particular, leading to different gamma frequencies (14, 40)] that may contribute to efficient coordination of global brain activity. This is particularly interesting because, although cell-to-cell heterogeneity is known to have a strong effect in the dynamics of local populations (62, 63), the effect of area-to-area heterogeneity on large-scale networks is much less understood.

Anatomical and functional hierarchies

Whereas the brain hierarchy is robust, it is also flexible functionally. A strong and layer-specific input can modify the functional hierarchical distance between cortical areas and induce “jumps” in the hierarchy (see Results). This jumping phenomenon may result from layer-specific input arriving from cortical areas higher in the hierarchy, such as parietal or prefrontal areas. In particular, areas involved in higher cognitive functions, such as working memory or rule representation, should be able to send layer-specific top-down signals in a context-specific manner, leading to a reorganization of sensory and association areas in the functional hierarchy, as observed experimentally (5). The thalamocortical system may also contribute to altering the brain’s functional hierarchy, a possibility that can be analyzed in our large-scale model extended to include the thalamus.

In the context of the anatomical hierarchy, the hierarchical position of FEF areas 8l and 8m deserves mention. The original Felleman and Van Essen study (1) placed FEF above V4, based on the FF nature of the projection from V4 to FEF but in the absence of data from the FB projection. The high position of FEF with respect to V4 was also observed in later analysis (34) and in our own analysis (fig. S5). On the other hand, earlier work also suggested that FEF could be below V4 (2, 64).

This apparent contradiction disappears when we take into account that the position of FEF depends strongly on the areas that are included in the analysis used to obtain the anatomical hierarchy. Prefrontal areas are strongly connected to 8l and 8m, and therefore, a hierarchy that includes prefrontal areas [as in the study by Chaudhuri et al. (34) and in fig. S5] locates both 8l and 8m very high in the hierarchy. On the other hand, if largely visual areas are considered in the estimation of the hierarchy, then the relative position of FEF will be mainly driven by visual areas (2, 64), and therefore, 8l will be especially lower in the anatomical hierarchy. The fact that in certain tasks the functional hierarchy from Bastos et al. (5) looks similar to the functional hierarchies found by Markov et al. (2) and Barone et al. (64) could be due to a smaller recruitment of the prefrontal cortex in these tasks, which would leave the FEF-V4 relationship to be dictated by visual pathways.

On top of these anatomical considerations, it is worth noting that the position of a given area in anatomical and functional hierarchies does not have to coincide. The functional hierarchy heavily depends on the dynamical interactions between areas, and although the anatomical hierarchy can be a good estimation (as shown in Fig. 6, F and G), the functional position of a given area can be context-dependent, as observed by Bastos et al. (5) and explained by our model (Fig. 7). In the anatomical hierarchy, frontal areas 8l and 8m are higher than early visual areas (see fig. S5), but in the functional hierarchy, these frontal areas seem to be highly dynamic. This suggests that, under some circumstances, 8l and 8m will rank much higher in the functional hierarchy, in agreement with the above neuroanatomical view. This would also agree with physiological evidence from behaving monkey studies that show that FEF is a source area for attentional modulation in V4 [(65) and references therein], which suggests that FEF could rank above V4 in the hierarchy. In general terms, we expect that the observed differences between anatomical and functional hierarchies will provide insight into task-dependent computations that could be analyzed in future studies.

Finally, understanding the cortical mechanisms underlying hierarchical processing will improve future estimates of the anatomical hierarchy in humans. The measurement of laminar-specific projections via tract-tracing studies defines an anatomical hierarchy in macaques (1, 2), but these anatomical data are not available with known imaging or postmortem technique in humans. Diffusion tensor imaging does not provide directionality information about interareal connectivity. Therefore, the anatomically defined hierarchy is not known for humans at the present time. However, inspired by studies showing a strong correlation between functional hierarchy obtained from frequency-dependent GC analysis with anatomically defined hierarchy, a recent human study showed that the same analysis applied to magnetoencephalography yielded a functional hierarchy in humans (6), allowing inference of anatomical hierarchy of the human cortex at least for areas showing strong homology with macaque. Computational modeling studies, such as the one presented here, will contribute to solving this problem and will strengthen the links between functional and anatomical connectomes.


Anatomical data

The anatomical connectivity data used in this work came from an ongoing tract-tracing study in macaques (2, 8, 11). In short, retrograde tracer was injected into a given target area, and it labeled neurons in several source areas projecting to the target area. The number of labeled neurons on a given source area allowed to define the FLNs from that source to the target area. In addition, the number of labeled neurons located on the supragranular layer of a given source area (over the total number of labeled neurons on that source area) allowed us to define the supragranular layered neurons (SLN) from that source area to the target area. Source areas that were lower (higher) than the target area in the anatomical hierarchy tended to have a progressively higher (lower) proportion of labeled neurons in the supragranular layer. That is, the lower (higher) the source area relative to the target area, the higher (lower) the SLN values of the source-to-target projection. By repeating the process using other anatomical areas as target areas, we obtained an anatomical connectivity matrix with weighted directed connections and laminar specificity.

We used the anatomical connectivity matrix from Markov et al. (2), and we also performed additional injections to expand the connectivity data and include the parietal area LIP in the connectivity matrix. The 30 cortical areas, which constituted the new connectivity matrix, were as follows: V1, V2, V4, TEO, 9/46d, F5, 8m, 7A, DP, 2, 5, 7B, STPr, STPi, STPc, PBr, TEpd, 24c, F1, F2, F7, ProM, 8L, 9/46v, 46d, 8B, MT, 7m, 10, and LIP. The anatomical data used in the present study can be downloaded from (LIP_Dataset.xlsx, Neuron_2015_Table.xlsx and PNAS_2013_Distance_Matrix.xlsx).

The corresponding 30 × 30 matrices of FLN, SLN, and wiring distance are shown in figs. S2, S3, and S4, respectively. The anatomical hierarchy, which was computed using the SLN values and a generalized linear model, as in previous studies (2, 34), is shown in fig. S5. The anatomical data are also shown in fig. S6 for a particular subset of cortical areas of interest (used for the study of the functional hierarchy in the model). The subset of areas of interest was V1, V2, V4, DP, 8m, 8l, TEO, and 7A. Surgical and histology procedures were in accordance with European requirements 86/609/EEC and approved by the ethics committee of the Rhône-Alpes region.

Computational model

The large-scale model of the macaque cortex was built by assembling four different levels of description of increasing sizes: (i) local excitatory-inhibitory populations of the Wilson-Cowan type that describe the activity within a given layer, with quantitative differences for the supra- and infragranular layers so that they preferentially exhibit noisy oscillations in the gamma and alpha frequency band, respectively; (ii) interlaminar circuits coupling supra- and infragranular layers; (iii) interareal couplings that consider the layer-specific influences between two given cortical areas (such as V1 and V4); and (iv) a large-scale laminar cortical network, which extends the anatomical connectivity data from previous studies (2, 8) to include the posterior parietal area LIP and a laminar structure of local areas. The network model consisted of 30 cortical areas distributed among occipital, temporal, parietal, and frontal lobes. Further details of the model, including equations, parameter values, and the routines used for data analysis, can be found in the Supplementary Materials.


Supplementary material for this article is available at

Supplementary Methods

fig. S1. Microstimulation experiments at the interareal level, for an FB projection that is strong to L2/3E and weak to L5E (more precisely, with a supra/infra ratio of 0.8).

fig. S2. FLN connectivity matrix, after a logarithmic transformation for visualization purposes, for the 30 areas of the large-scale model.

fig. S3. SLN connectivity matrix for the 30 areas of the large-scale model.

fig. S4. Wiring distances, in millimeters, for the 30 areas of the large-scale model.

fig. S5. Anatomical hierarchy obtained from the data shown in fig. S3.

fig. S6. Different matrices for the subset of eight areas of interest (V1, V2, V4, DP, 8m, 8l, TEO, and 7A) used in the functional hierarchy study.

fig. S7. Spectral pairwise-conditioned GC profiles for all the possible pairs of interactions between the eight cortical areas of interest: V1, V2, V4, DP, 8m, 8l, TEO, and 7A, with a background input of I = 6 to all areas, plus a strong extra input of I = 6 to V1.

fig. S8. Effect of introducing a distance-dependent relationship on the target of FB projections.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Funding: This work was supported by Office of Naval Research grant N00014-13-1-0297, NIH grant R01MH062349, Science and Technology Commission of Shanghai Municipality grants 14JC1404900 and 15JC1400104 (X.-J.W.), and ANR-11-BSV4-501 and LABEX CORTEX (ANR-11-LABX-0042) of Université de Lyon, program “Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (H.K.). Author contributions: J.F.M., J.D.M., and X.-J.W. conceived the study. J.F.M. performed the simulations. J.F.M., J.D.M., and X.-J.W. analyzed the results. H.K. contributed data, and all authors contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. Anatomical data can be downloaded from The code with the computational model can be downloaded from the laboratory website of X.-J.W.

Stay Connected to Science Advances

Navigate This Article