Introduction

The primate brain is a complex hierarchical network of interconnected cortical areas that perform computations essential for sensorimotor processing and cognitive functions. Among these, certain computations are repeated across cortical areas and are thus said to be “canonical computations”. Divisive normalization [1] is one such fundamental canonical principle, playing a central role in regulating neural responses across sensory and higher cortical areas.

Feedback connections are ubiquitous throughout the primate brain, yet their precise computational role remains poorly understood. A wide range of functions has been attributed to feedback, including attentional modulation [2], inter-areal communication [3], learning [4], and generative modeling [5]. While divisive normalization has been extensively studied in feedforward recurrent networks [6–9], its implementation and effects in hierarchical networks with feedback remain unexplored.

Effective communication between brain areas is essential for complex cognitive tasks, yet the precise mechanisms of communication remain an active area of investigation. One prominent hypothesis, communication through coherence (CTC) [10], proposes that synchronized oscillations (coherence) are crucial for effective inter-areal communication. Indeed, coherence between brain regions has been widely observed across cognitive, perceptual, and behavioral tasks [11, 12]. However, the causal nature of the relationship remains a key unresolved question. It is currently debated whether coherence is a fundamental mechanism enabling inter-areal communication or rather an epiphenomenon emerging from underlying interactions [13, 14]. A complementary hypothesis suggests that brain areas communicate through a low-dimensional communication subspace (CS) embedded within the high-dimensional neural activity space [15]. This concept, also identified in artificial neural networks [16], posits that only activity within these shared subspaces between the source and target areas facilitates communication. Despite strong empirical evidence for both coherence-based and subspace-based communication, a unified theoretical framework capable of jointly addressing these mechanisms and explaining how they are dynamically modulated by factors such as feedback and stimulus contrast is currently lacking.

To address these gaps, we introduce a hierarchical recurrent neural circuit theory (Fig. 1) built upon two a priori constraints: i) a hierarchy of brain areas with feedforward, recurrent, and feedback connections, and ii) recurrently implemented normalization within each brain area. Within this framework, we derive analytical expressions for inter-areal coherence and communication subspaces, offering experimentally testable predictions. Remarkably, these two constraints alone give rise to a rich set of emergent properties that explain a wide spectrum of neurophysiological phenomena, including modulation of contrast response functions, oscillatory dynamics evident in power and coherence spectra, the emergence of communication subspaces, and the dynamic control of functional connectivity. A central contribution of this work is a unified view of inter-areal communication. The theory simultaneously predicts empirical observations of both communication subspaces and inter-areal coherence, and further reveals that coherence actively shapes the communication subspace. In doing so, our framework reconciles previously competing hypotheses and provides a principled account of how coherence influences subspaces.

Figure 1:

Finally, we emphasize that throughout our analysis, we neither adjust model parameters (such as time constants, weight matrices, and noise amplitudes) to reproduce specific experimental trends nor fit data when making direct comparisons. All results presented were obtained using a single set of model parameters. Additionally, we have verified that the key results hold across a wide range of parameters; for further details, see the SI Appendix.

Results

Theory

Initially conceived to explain the responses in the primary visual cortex (V1) [1], the normalization model has emerged as a canonical computation across sensory and cognitive domains [17], including higher level vision, olfaction, audition, somatosensation, attention, multi-sensory integration, working memory, and decision-making. For citations see [6]. At its core, the normalization model proposes that the response of an individual neuron is scaled by the collective activity of its neighboring neurons (viz., the normalization pool), analogous to normalizing the length of a vector.

In addition, normalization decorrelates population responses to reduce redundancy [18, 19] and prevents excessive recurrent excitation to maintain network stability [8, 9]. It is a computational-level description of neural responses that can be implemented using various mechanisms. Experimental evidence suggests that normalization operates via recurrent amplification [20, 21], amplifying weak inputs more than strong inputs, consonant with evidence that cortical circuits are dominated by recurrent excitation [22, 23].

Oscillatory Recurrent Gated Neural Integrator Circuits (ORGaNICS) is a theoretical framework that dynamically implements divisive normalization by modulating recurrent amplification [6, 7]. Previous work on normalization models (including previous work on ORGaNICs) has focused on a single neural circuit within a brain region, such as V1, ignoring the abundant top-down feedback connections. Here, we develop a new hierarchical ORGaNICs theory with feedback connections in addition to feedforward and recurrent connections.

The hierarchical ORGaNICs theory describes neuronal population responses as dynamic processes evolving in a recurrent circuit, modeled by nonlinear differential equations (see Fig. 1 for illustration and Methods for equations). The circuit comprises distinct cell types: excitatory principal cells, excitatory modulator cells, and inhibitory modulator cells (n.b., “modulator” refers to a multiplicative computation, not necessarily implemented with neuromodulators). Output firing rates of principal cells depend on three terms: 1) feedforward input drive modulated (viz., scaled) by input gain; 2) feedback drive modulated by feedback gain; 3) the sum of recurrent drive and scaled feedback drive, collectively modulated by recurrent gain (see equation at the top of Fig. 1). Unlike generic E-I models in which all neurons perform the same computation, ORGaNICs assign unique computations to each cell type (see Discussion and Methods). This theory is uniquely well-suited to analytical treatment, e.g., to derive closed-form solutions for the power spectral density, coherence, and communication subspaces (see Methods and Appendix).

Each neuron receives input drive defined as a weighted sum of responses from neurons at lower stages of the hierarchy. For example, LGN responses provide input drive to V1 (z₁ in Fig.1), and V1 responses similarly drive V2 (solid arrow from y₁ to y₂ in Fig.1, representing principal neuron responses in each area). This weighting is defined by an embedding weight matrix composed of positive and negative synaptic weights, which determine each neuron’s stimulus preferences, such as orientation selectivity (Fig.1, V1 tuning curves), analogous to a standard receptive field model (see Methods). Input gain is determined independently by input modulator cells (β₁ and β₂ in Fig.1), whose gain may, in principle, depend on both afferent input and local population activity (e.g., responses of both LGN cells and V1 cells may contribute to the input gain in V1).

Each neuron also receives recurrent drive arising from weighted interactions among principal cells in the same brain area (Fig. 1, curved arrows). Recurrent connectivity follows a center-surround architecture, with local excitatory and distal inhibitory connections (Fig. 1, recurrent synaptic weights). The gain of this recurrent drive is dynamically controlled by recurrent modulator cells (a₁ and a₂), which combine excitatory (u₁, u₂) and inhibitory signals (q₁, q₂) from local interneuron populations as well as feedback from higher cortical areas (solid arrow from y₂ to a₁ in Fig. 1). Divisive normalization is achieved dynamically through the activity of these modulator cells, which carry the “normalization signal”.

Finally, each area receives feedback drive, a weighted sum of responses from neurons at higher stages of the hierarchy (e.g., V2 → V1; Fig. 1). The feedback drive is specified by a feedback weight matrix, which, for simplicity, we take to be the transpose of the corresponding feedforward matrix from V1 to V2 [5]. The strength of this feedback is controlled by a feedback gain parameter (γ₁ in Fig. 1). More generally, feedback strength may be modulated either by changing the feedback gain or by changing the gain of the neural responses in the higher cortical areas providing feedback (see Discussion). Thus, the feedback gain may be a fixed property of the circuit or a dynamically regulated quantity that depends on both the inputs and outputs of the corresponding area. Importantly, even if the feedback gain is fixed, it may still be possible to manipulate the feedback strength experimentally to test predictions of the theory (elaborated in what follows).

The two-stage recurrent model with feedback implements divisive normalization exactly and simultaneously across multiple cortical areas. In particular, when the gains are set to γ = β = 1, and the recurrent weight matrix is the identity matrix (i.e., each neuron recurrently amplifies only itself), the steadystate responses of the model coincide precisely with the analytical steady-state solutions of the normalization equation (Methods, Eq. 2).

Contrast response functions

Both V1 (dashed) and V2 (solid) exhibit the characteristic saturating, sigmoidal dependence on contrast (As shown in Fig. 2a, b). The model predicts a lower semisaturation constant (the contrast at which the firing rate reaches half-saturation) and a steeper slope for cortical area V2 compared to V1. This prediction is commensurate with experimental observations in the literature that used optimal stimulus sizes (Fig. 2b) [24, 25]. Similar hierarchical trends, specifically an increase in slope and a decrease in the semisaturation constant, have also been observed across other brain areas like the LGN, V1, and MT with progression up the visual hierarchy [27].

Figure 2:

The theory offers a novel, experimentally-testable prediction that increasing either the feedback gain, γ₁, or input gain, β₁, enhances neural responses across cortical areas (Figs. 3a,b). Importantly, this response enhancement is evident primarily as a change in contrast gain, i.e., a shift of the contrast-response functions on the log contrast axis rather than a change in response amplitude.

Figure 3:

V1 power spectra

Power spectra analysis of electroencephalographic (EEG) and local field potential (LFP) signals has been used extensively to identify and quantify oscillatory activity in various frequency bands (e.g., alpha, beta, gamma, theta), hypothesized to be linked to various cognitive functions, behavioral states, and neurological conditions [28]. Since our model has a known fixed point solution, we calculate the power and coherence spectra analytically, given a model of the noise in the circuit (see Methods).

We adopted an analytically tractable noise model that adds low-pass filtered synaptic noise to the membrane potential and multiplicative Gaussian white noise to the firing rate dynamics. Furthermore, we explicitly accounted for long-range correlated noise in LFP measurements. This type of noise, propagating through brain tissue, affects LFPs with negligible impact on spiking activity and firing rates [29]. Further details are provided in the SI Appendix (“Intrinsic neural circuit noise” and “Extrinsic noise in the local field potential”).

Our theory successfully replicates key experimental observations of V1 power spectra across varying contrast levels. These theoretical V1 predictions (Fig. 2c) are compared with experimental spectra from macaque V1 (Fig. 2d) under identical contrast conditions [26]. The model accurately captures three critical experimental trends: 1) the power spectrum of the noise is shaped by recurrent normalization to exhibit a resonant peak in the gamma-frequency range (≈ 40-60 Hz); 2) the peak frequency shifts rightward with increasing contrast; and 3) the peak power decreases after reaching a maximum around 40% contrast. The oscillatory activity in the gamma-frequency range depends on normalization; removing normalization eliminates the peaks (see Suppl. Fig. S4e), commensurate with experimental observations that gamma oscillations are linked to normalization [30, 31]. The peak frequencies depend on the values of the intrinsic time constants; larger time constants would shift the peaks to lower frequencies (see Discussion).

The theory also reproduces the steep fall-off in power reported at higher frequencies [32], where power approximately scales as 1/ f ⁴ with increasing frequency (f) (Fig. 2c, inset). This occurs because the synaptic noise added to the membrane potentials is first low-pass filtered by the synapse’s own time constant, and then by the recurrent processing in the circuit that also acts as a low-pass filter (each of which contributes a factor of 1/ f ² to the power spectrum).

The theory offers a novel, experimentally-testable prediction that modulating feedback and input gains elicit distinct changes in the power spectra. Figures 3c-f illustrate these distinct effects on the normalized power spectra when feedback gain (left panels) or input gain (right panels) is varied. These comparisons are presented for both a low (3%; Figs. 3c,d) and a high (50%; Figs. 3e,f) stimulus contrast. See Supplementary Material for comprehensive results across other contrast levels (Suppl. Fig. S1).

The model predicts oscillatory activity in the alphafrequency range (≈ 10 Hz) for low stimulus contrasts and small feedback gain values (Fig. 3c). As the feedback gain to the excitatory neurons increases, this alpha peak diminishes and ultimately vanishes at high feedback gain levels, while a broad peak concurrently emerges in the gamma frequency range (Figs. 3c,e). No such oscillatory activity in the alphafrequency range is observed by manipulating the input gain, for any contrast (Figs. 3d,f).

The power spectra are dominated by oscillatory activity in the gamma-frequency range at high stimulus contrast, regardless of the input gain and feedback gain values (Figs. 3e,f). Increasing the feedback gain causes this gamma peak to shift towards higher frequencies. A similar shift of the gamma peak to higher frequencies is observed when the input gain is increased. However, a key distinction arises in their impact on the peak’s bandwidth: the gamma peak bandwidth narrows with increasing input gain but broadens with increasing feedback gain.

Inter-areal coherence

Coherence measures the statistical correlation between two signals (e.g., LFPs in two different brain areas) as a function of frequency. Neuronal synchrony is hypothesized to contribute to the dynamic selection and routing of information both within and between brain areas [10, 33, 34].

Our model’s predictions for V1-V2 coherence spectra accurately capture key experimental trends observed across various stimulus contrasts. We analytically determined the coherence between the maximally firing neurons of V1 and V2 (see Methods). Figure 2e displays the theoretically predicted V1-V2 coherence spectra for different contrasts. These model predictions are compared with experimental V1-V2 coherence data from macaque monkeys, obtained under identical contrast conditions (Fig. 2f; [26]). The model accurately captures several crucial experimental observations: 1) high coherence at lower frequencies for low contrasts; 2) a shift of the peak coherence toward higher frequencies with increasing contrast; and 3) a decrease in peak coherence after reaching a maximum around 40% contrast. The coherence peaks in the gamma frequency range are due to normalization; removing normalization eliminates the peaks in coherence (see Suppl. Fig. S4f). The theory offers a novel, experimentally-testable prediction that feedback and input gain differentially modulate V1-V2 coherence across contrasts (Fig. 3g-j). At low contrast (3%), a coherence peak in the beta band (12-30 Hz) is observed with low feedback gain (Fig. 3g); this peak disappears as feedback gain is increased. Notably, at low feedback gain, a dip in V1-V2 coherence occurs in the alpha band, corresponding to the frequencies where peak power is observed in V1. This suggests that alpha band activity, despite its prominence in V1 power, may not effectively propagate to V2 or participate in inter-areal communication. No such peaks are observed with changes in the input gain for low contrast (Fig. 3h). At high contrast (50%), increasing feedback gain shifts the peak coherence towards higher frequencies and increases both its magnitude and width (Fig. 3i). A similar frequency shift in peak coherence occurs with increasing input gain, but its magnitude and width decrease (Fig. 3j).

Communication subspaces

Neural communication is hypothesized to be channeled through “communication subspaces (CS)” [15]. Although neural activity space is high-dimensional, neural communication may not utilize this vast space. CS are lower-dimensional manifolds within the broader neural activity space that preferentially engage in inter-areal communication. In contrast, within-area communication utilizes a distinct, higherdimensional “private” subspace [15].

Our theory successfully captures key characteristics of the CS (Figs. 2g,h). We derived an analytical model to examine both inter-areal (V1-V2) and within-area (V1-V1) communication (see Methods). We found that the CS is low-dimensional, much lower than the number of neurons (30 neurons for the analysis in Figs. 2g,h, and Figs. 3k,l) in each brain area. In addition, the inter-areal communication (V1-V2) subspace is lower-dimensional than the within-area (V1-V1) communication subspace. Figure 2g illustrates the model prediction performance (i.e., the ability to predict activity fluctuations of a target subpopulation from the activity fluctuations of a source subpopulation using a linear model) as a function of the dimensionality of the CS. These theoretical predictions align well with experimental data (Fig. 2h; [15]).

The theory offers a novel, experimentally-testable prediction that feedback and input gain differentially modulate communication subspaces (Figs. 3k,l). Increasing feedback gain from the higher cortical area (V2) to the lower cortical area (V1) enhanced inter-areal communication while diminishing withinarea communication (Fig. 3k), without substantially altering the subspace dimensionality (see Suppl. Figs. S3a,b). Conversely, increasing input gain reduced both the inter-area and within-area communication (Fig. 3l). These results are for 100% contrast; see Suppl. Fig. S3c for other contrasts. The observed enhancement of inter-areal communication with increased feedback suggests a key role for top-down feedback, to dynamically modulate the efficiency of information transfer between cortical areas without altering the underlying structure of the communication subspace (see below, Dynamic modulation of functional connectivity through feedback), a finding consistent with previous work on subspace dimensionality [35].

Frequency decomposition of communication subspaces

The theory predicts that the effectiveness of communication between brain areas varies systematically with frequency and contrast (Fig. 4). We quantified the efficacy of communication by analytically calculating the prediction performance for each frequency component (see Methods). Using the Wiener-Khinchin theorem, we converted the power spectra of the source and target population into their noise correlations. To isolate noise correlations at a specific frequency, we applied an idealized narrow-band pass filter. This conceptual filter is implemented via a transfer function, which is non-zero only within a very narrow frequency band around the frequency of interest. Applying this filter to the power spectra and integrating across the frequency domain yields the frequency-specific correlations (see Methods, “Frequency decomposition of the communication subspace”. This conversion enabled us to predict the power of the target population from the source population (using a linear readout) at each frequency, thus providing a measure of communication efficacy (prediction performance) as a function of frequency.

Figure 4:

V1-V2 communication is dynamically modulated by stimulus contrast in a frequency-dependent manner. At low stimulus contrasts (e.g, approximately 10%), a distinct peak in communication is observed at ≈ 20 Hz (Fig. 4b). As stimulus contrast increased, this preferred frequency for communication shifts toward higher frequencies (Fig. 4b). Furthermore, the magnitude of this communication peak exhibits a non-monotonic relationship with contrast, reaching its maximum at around 40% contrast before declining at higher contrasts. These communication peaks are due to normalization; removing normalization eliminates the peaks in both coherence and prediction performance (see Suppl. Figs. S4f,g). Under conditions of very low contrast, no distinct peak in communication is observed; instead, communication is mostly concentrated at very low frequen-cies (< 10 Hz). Notably, these observed frequencyspecific trends in communication bear a striking resemblance to the V1-V2 coherence patterns (Fig. 4a). This parallelism indicates a relationship between the degree of synchronization of neural activity between brain areas and the effectiveness of communication between these brain regions, a prediction that can be tested experimentally (see Discussion, “Unified view of inter-areal communication”).

The dimensionality of the inter-areal (V1-V2) CS also varies with frequency and stimulus contrast (Fig. 4c). The dimensionality of the communication subspace decreases at frequencies where inter-areal coherence and communication are maximal, for each contrast (Fig. 4c).

The dimensionality of the CS is dictated by the rank of the covariance matrix between the source (V1) and target (V2) subpopulations (see Methods: Interpretation of the results, Suppl. Fig. S5). The dimensionality of the CS is influenced mainly by three factors. First, dimensionality depends on normalization. Although dimensionality is relatively low in the absence of normalization (Suppl. Fig. S4h), including normalization further reduces dimensionality at coherent frequencies (Fig. 4c). The intuition for why normalization reduces dimensionality is that when only a few neurons are strongly driven, and their responses are highly correlated across areas (V1 and V2), normalization suppresses the activity of other nearby neurons that are not strongly driven. This suppressive interaction, akin to a soft winner-take-all mechanism, diminishes the contribution of other neurons to the shared covariance, thereby reducing the overall dimensionality. Second, the dimensionality is increased for a fully interconnected (all-to-all) connectivity between V1 and V2 (Suppl. Fig. S4l). Third, delocalized inputs (e.g., noise images or naturalistic images) lead to a higher dimensionality than localized inputs (e.g., gratings) (Suppl. Fig. S4p), commensurate with experimental results [15].

Dynamic modulation of functional connectivity through feedback

Anatomical connections influence neural communication but do not completely determine it. Rather, it has been observed experimentally that functional connectivity between brain areas changes dynamically [36, 37].

Our theory predicts that functional connectivity is controlled by modulating feedback gain (Fig. 5). We expanded our model to encompass three cortical areas: V1, V4, and V5 (also called MT). In this expanded model, V1 has reciprocal connections with both V4 and V5, representing the parallel dorsal and ventral visual pathways, respectively [38]. For simplicity, we assume there are no direct connections between V4 and V5 in the model. Analysis of this three-area model provides compelling evidence for the role of top-down feedback in controlling functional connectivity. When the feedback strength from V5 → V1 is greater than the feedback from V4 → V1, we observe that V1 strengthens communication with V5 (Fig. 5a). Conversely, when the feedback from V4 to V1 is stronger, V1 shows enhanced communication with V4 (Fig. 5b). This suggests that the relative strength of feedback from the higher cortical area (V5 or V4) can dictate the information flow from the lower cortical area (V1). One possible way this could be achieved is by selectively enhancing the neural activity in either V5 or V4, which will, in turn, increase the feedback to V1. This is commensurate with experimental evidence that feature-based attention and task demands may selectively boost activity in one or another cortical area [39], and specifically that attending to motion increases the transfer of motion information (communication) between V1 and V5 [40]. Therefore, top-down signals can dynamically modulate functional connectivity from lower to higher cortical areas, simply by adjusting the gain of neural activity in one higher cortical area (e.g., V4) relative to another (e.g., V5). This mechanism enables flexible cortical information routing, which is essential for cognitive flexibility.

Figure 5:

Note that, to isolate the variable of interest, we modeled V4 and V5 with identical parameters for their recurrent, feedforward (V1 → V4/V5) and feedback (V1 ← V4/V5) connectivity. This ensures that the observed increase in communication with either V4 or V5 is a direct consequence of the stronger feedback from that specific area.

Discussion

We present an analytically tractable theory of neural circuit function that implements divisive normalization dynamically across interconnected brain areas. Emergent properties of the theory account for a diverse array of neurophysiological phenomena, including oscillatory activity (measurable as peaks in LFP power and coherence spectra) and communication subspaces. Based on the theory, we propose a unified view of inter-areal communication. Our theory yields two key, experimentally-testable predictions: first, high-coherence frequencies serve as a signature for both stronger communication and reduced subspace dimensionality, and second, that increasing feedback strength selectively enhances inter-areal communication and diminishes within-area communication, without altering the subspace dimensionality. A comprehensive list of further predictions is provided in SI Appendix section “Robustness with respect to model parameters”. Furthermore, the theory proposes a mechanism for dynamic functional connectivity: topdown feedback signals enhance communication between brain areas simply by adjusting the gain of neural activity in one or more cortical areas. This modulation of neural gain offers a mechanism for dynamically changing the functional connections between brain areas over time.

Neural circuit function

The original ORGaNICs model [6, 7], in its initial form without feedback connections, serves as a biophysically-plausible counterpart to Long Short Term Memory units (LSTMs) [41], a class of artificial recurrent neural networks (RNNs) prevalent in machine-learning applications [42]. Unlike LSTMs or any other previously proposed RNNs, ORGaNICs have built-in recurrent normalization, which offers several advantages, including stability even with exuberant recurrent connections [8, 9].

The hierarchical version of the ORGaNICs theory introduced here adds feedback connections, drawing a partial analogy to “hierarchical predictive coding” models of sensory and perceptual processing [43–48]. However, our hierarchical theory diverges from predictive coding in its fundamental mechanisms. In predictive coding, it is proposed that forward-propagating signals represent prediction errors, while feedback pathways convey predictions to lower levels of hierarchy. The goal is to minimize prediction error by matching sensory input with topdown predictions. In hierarchical ORGaNICs, instead of simply propagating prediction errors, the interplay between feedforward and feedback drives adjusts neuronal responses to a value that balances the feedforward processing of its inputs with the feedback (see [5] for details). This mechanism is more consistent with neurophysiological and psychophysical phenomena than predictive coding models. For instance, neurons maintain sustained activity in response to a predictable stimulus [49], whereas predictive coding theories would suggest that the activity, representing prediction error, should be “explained away” as the stimulus becomes predictable.

Normalization

The proposed theory admits an analytical fixed point that follows the normalization equation (Eq. 2) exactly, for arbitrary non-negative normalization weights, when the input drive is constant over time (see SI Appendix for derivation) [8, 9]. Although several circuit mechanisms have been hypothesized to underlie normalization, including shunting inhibition, synaptic depression, and inhibition-stabilized networks [21, 50, 51], these models suffer from significant limitations. They either lack recurrent amplification or do not exhibit dynamics linked to normalization (e.g., gamma oscillations), or implement weighted normalization only approximately. Critically, none of them converge exactly to the normalization equation (Eq. 2). These limitations have significant practical consequences for generating testable predictions and fitting empirical data. Thus, despite the broad empirical success of the normalization model, it remains uncertain whether prior circuit models [51, 52] offer adequate mechanistic explanations.

Our model achieves normalization dynamically via recurrent amplification, which is inversely regulated by the activity of inhibitory modulator cells. When the input is weak, the modulator cells exhibit a small response, leading to strong recurrent amplification. Conversely, a strong input generates a large response from modulator cells, which suppresses the amplifycation. This process modulates both excitatory and inhibitory recurrent signals, which is consistent with the experimental observation that normalization involves a decrease in both recurrent excitatory and inhibitory conductances [21]. Because the theory successfully implements normalization, it is commensurate with a wide range of physiological and psychophysical phenomena previously explained by the normalization model (see Introduction for references).

The theory links normalization with oscillatory activity. Removing divisive normalization eliminates oscillatory peaks in power and coherence spectra (Suppl. Figs. S4e,f) and abolishes the link between coherence, prediction performance, and subspace dimensionality (Suppl. Figs. S4g,h), demonstrating that normalization is central to both the dynamics and function of the circuit.

Role of feedback

This paper provides an experimentally testable computational model of cortical function where neural activity is shaped by a combination of feedforward drive (bottom-up) and feedback drive (top-down). Feedback in the model may serve a variety of functions. First, feedback provides a means for controlling inter-areal communication (functional connectivity) simply by changing the gain of neural activity in one or more cortical areas (see Discussion section “Attention and functional connectivity”). Second, the feedback drive may contribute to selective attention (see Section “Attention and Functional Connectivity”). Third, the relative strength of the input gain and feedback gain determines whether neural responses are driven bottom-up, top-down, or a combination of the two. Fourth, feedback contributes to the normalization signal by providing excitatory drive to inhibitory neurons (see Methods). This mechanism is essential for maintaining network stability, and is commensurate with experimental evidence that normalization depends in part on loops through higher visual cortical areas [53].

The interplay between the input gain and the feedback gain determines the primary driver of neural responses. If the input gain is strong and the feedback is weak, the system behaves like traditional feedforward models. Conversely, with strong feedback and weak input, the network is similar to generative models, creating sensory representations from higher-level abstract knowledge, akin to recalling an image from memory. When input and feedback are balanced, the framework may integrate prior beliefs with sensory evidence (see [5] for details), akin to models based on Bayesian inference [54, 55].

The theory predicts differential effects of changing input gain versus feedback gain on the power, coherence spectra, and inter-areal communication (see Fig. 3). Hence, these metrics may be used as a signature for distinguishing changes in input gain from changes in feedback.

Unified view of inter-areal communication

Our theoretical framework offers a novel perspective on inter-areal communication. Rather than viewing “communication through coherence” (CTC) and “communication subspace” (CS) as alternative hypotheses, our theory unifies them, showing that coherence and communication subspaces are complementary empirical phenomena emerging naturally from the same underlying mechanism. A similar conclusion was recently reached by Ni et al. [56] using a large-scale spiking neural circuit model. We found that the frequency decomposition of communication subspaces closely mirrors the patterns of inter-areal coherence at different contrasts. Specifically, frequencies exhibiting high coherence are also those for which the communication efficacy within the subspace is maximal, and the subspace dimensionality is minimized. Consequently, coherence might dynamically shape the communication subspaces, for instance, by reducing their dimensionality at preferred frequencies [57]. Moreover, the predicted frequency-dependence of communication subspaces might serve as a marker for testing neural circuitry theories (such as ours).

Attention and functional connectivity

The theory provides a possible explanation for the efficient routing of sensory information between cortical areas, i.e., for changing functional connectivity. Lower cortical areas, such as V1, can flexibly relay information to specific higher cortical areas (e.g., V4 vs. V5) based on the relative strength of the feedback from each of the higher cortical areas. The strength of the feedback can be controlled either by changing the feedback gain or (perhaps more simply) by changing the gain of the neural responses in each of the higher cortical areas. We hypothesize that this feedback-mediated process accounts for changes in functional connectivity [58] and feature-based attention [59]. Consistent with this hypothesis, experimental evidence demonstrates that feature-based attention is spatially global, i.e., that neural activity is boosted in feature-selective neurons with receptive fields covering the visual field [60].

The theory may also provide insights into the neural processes underlying spatially-selective attention. The theory identifies two distinct pathways through which attention modulates neural gain. First, attention may modulate the input gain of the sensory neurons, consistent with established normalization model of attention [61]. Second, attention may modulate the activity of the sensory neurons via the topdown feedback drive. We propose that these mechanisms may map onto different cognitive processes: for instance, changes in input gain might primarily drive exogenous (bottom-up) attention, while the feedback drive facilitates endogenous (top-down) attention. Crucially, these two pathways are not just conceptually distinct but produce divergent spectral signatures. Our theory predicts differential effects of changing input gain versus feedback gain on the power and coherence spectra. These spectral metrics, therefore, offer an empirical means to distinguish between the two mechanisms.

Thalamo-cortical model

Alpha oscillations are a defining characteristic of thalamo-cortical dynamics, particularly prominent within the visual system during states of “relaxed wakefulness” [62, 63]. Specifically, alpha oscillations are hypothesized to be generated and modulated by thalamo-cortical loops [63, 64], suggesting that alpha rhythms signify a state of functional inhibition or “idling” in sensory pathways, potentially gated by corticothalamic feedback [64, 65]. While the hierarchical model presented in this paper was designed for cortico-cortico interactions, it may be adapted to model the thalamo-cortical interactions that contribute to the generation of alpha waves. Our theory predicts the emergence of oscillatory activity in the alpha band (8-12 Hz) under conditions of low stimulus contrast and low feedback to excitatory neurons, and high feedback to inhibitory neurons (γ < 1, Figs. 3c,g). This aligns well with the existing literature demonstrating the role of alpha oscillations in mediating cortical inhibition [65, 66]. It is, however, unclear if the current theory will suffice, on its own, to fully capture the phenomenology of alpha waves.

Furthermore, we hypothesize that the input gain modulator responses depend in part on thalamocortical loops, hence positioning the theory as a means for modeling the interplay between bottom-up sensory input and top-down feedback within the thalamocortical circuit, the gating of visual information flow, thalamic control of functional cortical connectivity and thalamocortical loops in other neural systems [68].

Limitations

Our model assumes instantaneous feedback, whereas biological feedback is delayed due to axonal conduction, synaptic transmission, and processing in higher cortical areas [69, 70]. Examining how inter-areal synchronization depends on feedback delay is an important direction for future work.

The feedforward, feedback, and recurrent weight matrices may contain both positive and negative weights. We therefore posit inhibitory interneurons to implement the sign inversion, corresponding to negative synaptic weights. These inhibitory neurons need not be one-to-one with their excitatory inputs. Rather, each may compute a partial weighted sum of their inputs corresponding to the negative terms in the matrix multiplications. Regardless, these inhibitory interneurons will introduce an additional delay.

Although most computations in the theory are biophysically plausible (see “Biological plausibility” and [7]), some mechanisms remain unresolved. In particular, while weighted sums and gain control may be implemented via balanced excitatory–inhibitory push–pull and push–push mechanisms [71, 72], we do not specify a cellular mechanism for the elementwise product of firing rates in Eq. 1b. The intrinsic time constants used in simulations are unrealistically short for individual neurons (see Table S1). Using more realistic values would shift gamma-band oscillations to lower frequencies. These parameters may therefore be better interpreted as effective population time constants rather than as the intrinsic time constants (capacitance/conductance) of individual neurons.

Finally, the model omits neural adaptation, a ubiquitous feature of cortical responses. Incorporating adaptation into the hierarchical framework remains an important direction for future studies.

Biological plausibility

We hypothesize a correspondence between the variables in the theory and cortical microcircuits composed of identified cell types. Unlike generic excitatory-inhibitory (E-I) models (e.g., [74]), our framework assigns unique computations and non-linearities to distinct cell types (see SI Appendix). Principal cells y are hypothesized to correspond to pyramidal neurons; inhibitory modulatory cells a to a combination of parvalbumin-expressing (PV+) and somatostatin-expressing (SST+) interneurons; and inhibitory interneurons q to vasoactive intestinal polypeptide (VIP) neurons.

In the proposed circuit (Fig. 1), Principal cells y excite both inhibitory modulatory cells a (via u) and inhibitory interneurons q. VIP-like cells q inhibit modulatory cells a, which in turn inhibit principal cells y. This architecture closely mirrors known cortical motifs: pyramidal neurons excite PV+, SST+, and VIP cells; VIP cells inhibit SST+ cells; and both PV+ and SST+ cells inhibit pyramidal neurons [75– 79]. We further hypothesize that modulator cells a have large receptive fields and broad orientation-selectivity (reflecting properties of the normalization pool), consistent with the response properties of SST+ and PV+ neurons, respectively [76, 80, 81].

We also posit excitatory feedback from higher cortical areas to both principal cells y and modulatory cells a (Fig. 1). Feedback to principal cells is likely received by the distal apical dendrites of layer 1, while feedback to modulator cells a is commensurate with experimental evidence that feedback from higher cortical areas targets SST+ neurons [82–84], and that normalization depends in part on loops through higher visual cortical areas [53].

One potential discrepancy between the model and the cortical circuitry is that excitatory input from principal cells y to modulatory cells a is mediated by an intermediate excitatory population u, whereas pyramidal neurons directly excite PV+ and SST+ cells in cortex. However, the theory can be reformulated so that principal cells provide direct excitatory input to modulatory cells by substituting Eq. 1b into Eq. 1c (see Methods), without altering the model’s functional behavior.

Finally, we hypothesize that the computations performed by principal cells correspond to processes distributed across dendritic compartments of pyramidal neurons (Fig. 6). The equation governing the responses of the principal cells (see Fig. 6 and Methods, Eq. 1a) expresses the following synaptic current: (input gain) (input drive) + (recurrent gain) [recurrent drive + (feedback gain) (feedback drive)]. The input drive is computed in the distal basal dendrites, commensurate with evidence that feedforward synaptic inputs contact basal dendrites. The input gain could be controlled by inhibitory synapses near the basal trunk. The feedback drive is computed in distal apical dendrites in layer 1, which are known to receive signals from higher cortical areas. The feedback gain might be attributed to the amplification of the feedback drive within the apical trunk of pyramidal cells [85]. The recurrent drive is likely computed at more proximal apical dendrites, with recurrent gain shaped by inhibitory inputs near proximal apical trunk.

Figure 6:

Methods

The Model

We propose a multistage hierarchical recurrent circuit that incorporates feedforward connections from the primary visual cortex (V1) to higher visual areas (e.g., V2, V4, MT) and feedback connections from these higher brain areas to V1 (Fig. 1). We refer to V2 throughout this paper for concreteness, but it is merely an example of a higher brain area that is reciprocally connected with V1. This model builds on the ORGaNICs framework developed by Heeger et al. [6, 7].

The response of principal excitatory neurons in the model is a sum of three key components: i) input drive, a weighted sum of neural responses from a brain area that provides feedforward input (e.g., LGN → V1 and V1 → V2), modulated by input gain; ii) feedback drive, a weighted sum of neural responses from the higher cortical area (e.g., V1 ← V2 or V1 ← V4), modulated (viz., scaled) by feedback gain; and iii) the sum of recurrent drive, a weighted sum of neuronal responses from other neurons in the same brain area (e.g., V1 → V1 and V2 → V2) and scaled feedback drive, collectively modulated by recurrent gain.

The proposed neural circuit comprises four distinct cell types, the principal excitatory neurons (y), the modulatory excitatory neurons (u), the modulatory inhibitory neurons (a and q), with the following dynamics:

The above equations are for V1; V2 follows a similar set of equations (see SI Appendix). The membrane potentials of neurons in cortical areas V1 and V2 are denoted by y₁ and y₂, respectively, while their corresponding firing rates are represented by and . The ⊙ symbol represents the Hadamard product, element-wise multiplication of vectors. To estimate the firing rates and from their respective membrane potentials, we employ a modified Gaussian Rectification (GR) model [86, 87]. The variable z₁ denotes the input drive to V1 from the preceding brain area (the lateral geniculate nucleus of the thalamus, LGN), computed as a weighted (W_zx) sum of the LGN firing rates (Fig. 7a illustrates the resulting orientation tuning curves). τ_k is the intrinsic time constant for each subtype of neuron, where k = y, u, a, q. and are the recurrent and feedback drive, respectively, where W₁₁ is the recurrent weight matrix for V1. W₁₁ has a centersurround structure where the closer connections are excitatory and the more distant connections are inhibitory (Fig. 7b), same as in [7]. The recurrent weight matrix for V2, W₂₂, has a similar structure. The maximum eigenvalue for these recurrent weight matrices is set to unity to ensure stability. However, we have shown that this need not be the case because normalization ensures stability [8, 9]. W₁₂ represents the inter-areal feedback connectivity matrix from V2 to V1, its entries are all non-negative, and the matrix is diagonally dominant (Fig. 7c). The feedforward matrix from V1 to V2, W₂₁, is taken to be the transpose of W₁₂. Given that W₁₂ is symmetric, it follows that W₁₂ = W₂₁. The normalization matrix in V1 (N₁) and V2 (N₂) is taken to be all ones, meaning each neuron contributes equally to the normalization. The input gain and feedback gain for V1 are modulated by β₁ and γ₁, respectively. The recurrent and feedback amplifications in V1 are modulated dynamically by the inhibitory neuron population a₁.

Figure 7:

Both the V1 and V2 neurons follow the normalization equation exactly when β and γ are set to unity (i.e. when there is identical feedback to excitatory (y₁) and inhibitory (a₁) neurons) and W₁₁ is the identity matrix (i.e., each neuron recurrently amplifies itself), but we have shown that the circuit closely approximates normalization for a wide range of recurrent weight matrices (W₁₁) [9]. For V1 neurons, normalization is expressed as:

and likewise for V2

where is the drive from V1 to V2, σ is the semi-saturation constant. Normalization is achieved dynamically via recurrent amplification (amplifying weak inputs but not strong inputs) with the three modulatory neurons, u₁, a₁ and q₁ whose dynamics are governed by Eq. 1b, Eq. 1c, and Eq. 1d respectively. The equations for modulatory neurons are derived by substituting the normalization equation (Eq. 2) in the equation for the principal neurons (Eq. 1a) (see SI Appendix for derivation). An increase in the activity of principal neurons (y₁) leads to enhanced activation in u₁ and q₁. This, in turn, drives a net increase of activity in modulatory neurons (a₁). The elevated a₁ activity subsequently exerts an inhibitory effect on y₁, thereby implementing normalization by the recurrent interactions between the neural populations (see Appendix for derivation).

The parameter γ₁ modulates the feedback gain to the principal excitatory neuron population (y₁). This determines the relative feedback to the excitatory neurons (y₁) compared to the inhibitory neurons (a₁). As a rule of thumb, γ₁ < 1 indicates greater feedback inhibition, γ₁ = 1 denotes a balance between feedback excitation and inhibition, and γ₁ > 1 signifies more feedback excitation. For simplicity, we have set the feedback to the inhibitory neurons (a₁) at a fixed value of 0.50 (hence, the factor of in equation Eq. 1c and Eq. 1a, see SI Appendix). Similarly, the parameter β₁ modulates the relative input gain to the principal excitatory neuron population (y₁) compared to the modulatory excitatory neuron population (u₁). For simplicity, we set the input gain to the modulatory neurons (u₁) to 0.50 (hence, the factor of in equation Eq. 1a and Eq. 1b). Note that when β₁ and γ₁ are both set to unity and the recurrent matrix is identity, the network achieves excitatory-inhibitory balance exactly, and the normalization (Eq. 2) is followed precisely.

We modeled each brain area using a ring model, where principal neurons are arranged in a circle to represent continuous features like stimulus orientation. The recurrent connectivity between these neurons is determined by the distance between them (Fig. 7b). Each modeled area consists of 72 neurons for each of the four specified types (y₁,u₁,a₁,q₁). To simulate sensory input, a visual stimulus with a defined contrast and orientation is presented. The stimulus is then convolved with the tuning curves of the V1 neurons to provide the input drive to the V1 principal neurons. The tuning curves are modeled as raised cosine functions, each peaking at the neuron’s preferred orientation (Fig. 7a).

Power and coherence spectra

The time evolution of neural population activity in V1 and V2 can be described by a general stochastic differential equation (SDE) with an additive noise term LdW:

where x ∈ ℝⁿ is a vector representing the activity of neurons; f(x, u) ℝⁿ represents the deterministic part of the dynamics explained above (Eq. 1); dW is a vector of independent Gaussian increments with correlation matrix D (𝔼[dW · dW^T]), and L ∈ ℝ^n×n is the dispersion matrix defining how noise enters the system. Given the analytical solution for the fixed point of the dynamical system (Eq. 2), we can linearize the system around the fixed point (x^∗):

where J_x∗ is the Jacobian computed at the fixed point x^∗. This linearization leads to a wide-sense stationary Gaussian process, x(t), whose power spectral density matrix, 𝒮(ω) ∈ ℝ^n×n, is given by [88]:

From 𝒮(ω), we can directly compute the coherence, κ_ij, between any two i and j neurons in the two areas:

Communication subspaces

Here we present the method used to characterize the communication subspaces between two subpopulations (source and target) of neurons [15]. These subpopulations may comprise distinct neurons in the same brain area or across different brain areas. The basic idea is to determine how well the fluctuations in responses of a target subpopulation can be predicted using a linear transform (a weighted sum) of the fluctuations in responses of the source subpopulation, and to determine the dimensionality of this linear transform.

The fluctuations in the responses and their covariability are fully determined by the correlation matrix at steady-state when the strength of the noise is small. Therefore the correlation matrix, denoted as C(0), for the responses x(t) can be computed by solving the Lyapunov equation [88]:

Once C(0) is obtained, we follow the method proposed by Semedo et al. [15] to construct an analytical model of the communication subspace. To achieve this, we divide the mean-subtracted responses of the entire neuronal population, denoted as x, into two subsets: s (responses of the source neurons, e.g., in V1) and t (responses of the target neurons, e.g., in V2 or V1). We consider a linear model for predicting target responses using the source responses, , where the mean squared error (MSE) is given by:

where C₁ = 𝔼[ss^T], C₂ = 𝔼[tt^T], and C₃ = 𝔼[st^T]. These matrices can be directly obtained from C(0). To find the optimal readout matrix B_opt, we minimize ε with respect to B, leading to the Ordinary Least Squares solution:

First, we define the covariance of the target activity predicted by the optimal model (note the “hats”), which we denote as Ĉ ₂:

By substituting the expression for B_opt, we arrive at a key relationship:

The principal dimensions of this predicted covariance matrix, Ĉ ₂, define the communication subspace. These are found by computing the eigendecomposition:

Here, the columns of the orthogonal matrix V are the eigenvectors of Ĉ ₂ and Λ is a diagonal matrix of the corresponding eigenvalues, sorted in descending order. The prediction performance for a reduced-rank approximation of dimension i is defined as: A rank-i approximation, B_RRR, is constructed by projecting the optimal solution B_opt onto the subspace spanned by the top i eigenvectors of Ĉ ₂.

where V(i) is a matrix containing first i columns of V. Here P_i = V(i)V^T(i) is the projection matrix. Note that P_i is symmetric and idempotent . We define the prediction performance as:

where ε_i is the MSE for the rank-i approximation

ε₀ is the MSE for the case where the readout matrix is a null matrix, corresponding to the total variance of the target population (ε₀ = Tr(C₂)).

To simplify the expressions for ε_i. We substitute B_RRR(i) = B_optP_i into the expression:

Now, we substitute and recognize that the term is Ĉ ₂. Using the cyclic property of the trace (Tr(AB) = Tr(BA)) and the idempotent property of the projection matrix , we can simplify further:

By applying the cyclic property and expanding P_i, we can see that the second term is the sum of the top i eigenvalues (λ_j) of Ĉ ₂.

Now the expression for prediction performance simplifies to:

Which can also be interpreted as:

The maximum prediction performance (full-rank approximation) is given by:

Interpretation of the results

The prediction performance (Eq. 22) is governed by the interplay between the source (C₁), target (C₂), and cross-population (C₃) covariances. The effectiveness of communication (indicated by a large predictive performance) can be understood through three key factors: 1) Source-target correlation: Prediction performance is directly enhanced by a stronger correlation between the source and target populations (i.e., larger C₃ norm). 2) Signal reliability: Prediction performance is boosted by the reliability of the signals. Less noisy activity in source and target populations, indicated by smaller diagonal entries in C₁ and C₂, leads to better prediction performance. 3) Alignment: Prediction performance is maximized when the principal directions of source covariance (eigenvectors of C₁) align with the principal directions of source-target covariance (left singular vectors of C₃).

This framework also explains the results in the frequency domain. The key metric here is the crosspower spectrum, which is the Fourier transform of the cross-covariances, and coherence, which is the squared cross-power spectrum normalized by the product of the auto-spectra, Eq. 7. Frequencies that show high coherence, indicating high cross-covariances between source and target, the prediction performance is highest at these specific frequencies. Note that the dimensionality of the communication subspace is given by the rank of the matrix . Since the source covariance matrix C₁ is a full-rank matrix, the dimensionality of the subspace is governed by the rank of the cross-covariance matrix, C₃. In the frequency domain, this means the dimensionality at a specific frequency is determined by the rank of the cross-power spectral density matrix (see Suppl. Fig. S5a). Our results show that the extent of normalization from surrounding neurons plays a critical role in shaping the dimensionality of the communication subspace. With surround normalization, dimensionality decreases specifically at frequencies where coherence is strongest (Fig. 4c). Normalization reduces dimensionality because when just a small subset of neurons is strongly driven, and their responses are highly correlated across areas (V1 and V2), normalization suppresses the activity of surrounding neurons with weaker drive. This suppression, resembling a soft winner-take-all process, diminishes the contribution of those less active neurons to the population’s shared variability, leading to a lower effective dimensionality.

In contrast, without this normalization, no such correlation between coherence and dimensionality is observed (Fig. S4h).

Frequency decomposition of the communication subspace

To analyze how communication between different brain areas varies across frequencies, we consider the frequency decomposition of the neural responses. The Wiener-Khinchin theorem states that the autocorrelation function C(τ) is the inverse Fourier transform of the power spectral density (PSD) matrix (𝒮(ω)):

To isolate C(τ) at a particular frequency (ω^∗), we apply an idealized narrow band-pass filter centered at ω^∗. Specifically, the filter’s transfer function H(ω) is given by,

where A = [ω^∗ − ϵ, ω^∗ + ϵ] ∪ [− ω^∗ − ϵ, − ω^∗ + ϵ] and ϵ is a small positive value defining the narrow bandwidth 2ϵ around ω^∗. The spectral density matrix of this filtered signal is given by,

The resulting autocorrelation of the filtered signal is expressed as:

Using the definition of H(ω) and assuming ϵ to be sufficiently small such that 𝒮 (ω) and e^ıωτ are approximately constant over the integration intervals of width 2ϵ, the frequency-specific autocorrelation at ω^∗ is:

At zero time lag (τ = 0),

Since the PSD is Hermitian, i.e., 𝒮 (ω) = 𝒮^∗(− ω)), we have

This frequency-specific autocorrelation matrix (C(0) | _ω∗) can then be used in place of C(0) in the communication subspace analysis described above to evaluate the prediction performance at any frequency. As our analysis does not depend on the absolute scaling of this matrix, the prefactor 4ϵ can be disregarded, and we effectively use Re {𝒮(ω^∗)} as the frequency-specific correlation matrix.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code is uploaded as Source code. The code to generate all the figures can be found in the following GitHub repository: https://github.com/asit-pal/Hierarchical-ORGaNICs.

Acknowledgements

We thank Pascal Fries and Adam Kohn for generously providing the experimental data in Fig 2. We acknowledge Michael Halassa, Hillel Adesnik, Adam Kohn and Flaviano Morone for their valuable feedback on the manuscript, and Ruben Coen-Cagli, John Maunsell, and Jason MacLean for insightful discussions. This work was initiated in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. This work was supported by the National Institute of Health under award numbers R01MH137669 and R01EY035242. S.M. also acknowledges support from the Simons Center for Computational Physical Chemistry (Simons Foundation grant 839534, MT). This work was supported in part through the NYU IT High-Performance Computing resources, services, and staff expertise.

Additional files

Supplementary Information.

Additional information

Funding

National Science Foundation (NSF) (PHY-2210452)

• Stefano Martiniani

National Institutes of Health (R01MH137669)

• Stefano Martiniani

National Institutes of Health (R01EY035242)

• Stefano Martiniani

Simons Foundation (SF) (839534 MT)

• Stefano Martiniani