A hierarchy of time constants and reliable signal propagation in the marmoset cerebral cortex

Domain: physics
Published: 2025-11-26
Editorial outcome: minor_revision

Abstract
A hierarchy of time scales in the cerebral cortex is functionally desirable for rapid information processing in sensory areas and slow time integration in association areas. Here, through an analysis of electrocorticography (ECoG) data, we identified a timescale hierarchy in the neocortex of marmoset, a primate species commonly used in neuroscience. Constrained by the anatomical and electrophysiological data, we developed a multi-regional model of the marmoset neocortex that captures the observed timescale hierarchy phenomenon. Furthermore, we used the model to investigate how the neocortex reconciles information integration on distinct timescales in local areas with reliable signal propagation globally across regions. We found that a near-criticality state is optimal for both localized signal integration within areas and reliable signal propagation across areas in the multi-regional neocortex. Our model also mechanically explains recent experimental observations that the structural and functional connectivities are less correlated in association areas than in sensory areas.

Introduction
The cerebral cortex, together with peripheral sensory organs (e.g., retina, cochlea, spinal cord) and subcortical relays (e.g., thalamus), engages in information processing of external inputs, local integration of signals within regions, and the transmission of these signals along functionally specific pathways. A prominent and noteworthy characteristic of the neocortex is what is referred to as “timescale hierarchy”. In essence, different cortical areas exhibit varying timescales in their response to inputs, allowing for the hierarchical processing of temporal information. Primary sensory areas manifest shorter timescales, typically measured in tens of milliseconds. The short timescales enable these areas to respond rapidly and effectively to incoming stimuli. In contrast, higher-level association areas display longer timescales, often extending to hundreds of milliseconds or longer. These extended timescales support the integration of information over a relatively protracted time window. The phenomenon of timescale hierarchy emerges from a theoretical model of the multi-regional macaque neocortex1, and is supported by experimental studies conducted on multiple species, including mice2,3, macaque monkeys4,5, and humans6–11. It remains to be actively investigated on the universality of this phenomenon across different species and the corresponding large-scale circuit mechanisms.

Timescale hierarchy is a phenomenon that spans across multiple cortical areas, strongly suggesting that its underlying mechanism is intricately linked to the broad-scale structural attributes of the neocortex. On an intuitive level, the timescale of neuronal population activity is intricately tied to the strength of recurrent excitation within that population. In simpler terms, stronger recurrent excitation tends to result in slower timescales of activity. Consequently, the concept of timescale hierarchy points to the existence of a macroscopic gradient in recurrent excitation along the axis of cortical area hierarchy. This notion finds compelling support in experimental findings, particularly in studies that examine the number of dendritic spines on pyramidal neurons within the primate neocortex. These dendritic spines serve as a direct reflection of the synaptic excitation strength per neuron and notably display an increasing gradient along the cortical hierarchy12,13. Furthermore, beyond this spatial gradient in spine numbers, the neocortex also reveals macroscopic gradients in other biological properties, including neuron density14 and synaptic receptor distribution15,16. These gradients may indeed serve as fundamental principles governing the organization of the cerebral cortex on a large scale, as discussed in a comprehensive review17.

In addition to the inherent characteristics of individual cortical areas, timescale hierarchy may also be attributed to long-range connections among these areas. Recent years have witnessed significant advancements in the quantitative measurement of inter-areal connectivity in various species, including mouse18–21, macaque22–24, and marmoset25,26. It represents a departure from earlier studies that only provided a coarse assessment of primate connectomes by categorizing connections between area pairs as absent, weak, or strong27,28. In contrast, the contemporary connectome offers a more precise and quantitative characterization. The strength of inter-areal connections is characterized using a continuous metric known as the fraction of labeled neurons (FLNs). FLNs measure the relative strength of projections from a particular source area to a specific target area in relation to all source areas projecting to that target area. Such a detailed connectome framework enables the quantitative exploration of the organizational principles governing the cerebral cortex. For instance, one notable observation is that FLNs between two areas exhibit an exponential decay with increasing distance, a phenomenon known as the exponential distance rule29,30.

By incorporating the macroscopic gradient of excitation and a quantitative weighted and directed connectome of macaque monkey, a comprehensive multi-regional model of 29 areas was developed to qualitatively generate the timescale hierarchy phenomenon1. Notably, when the macroscopic gradient of excitation was removed or the connectivity was shuffled, the significance of the timescale hierarchy diminished, or in some instances, completely vanished1. To elucidate the underlying principles governing this phenomenon, a rigorous mathematical analysis was conducted to identify the importance of macroscopic gradient of excitation and detailed excitation-inhibition (E-I) balance. The balance condition stipulated that the inter-areal excitatory inputs originating from each source area must balance with the local inhibitory inputs from the target area31. As a result, cortical areas exhibited relatively weak effective interactions, leading to the emergence of timescale hierarchy in a densely connected network. While this detailed balance condition appeared to enhance signal gating, i.e., irrelevant signals are not able to propagate in the detailed balance state, as demonstrated in a computational study32, it also introduced a challenge: the model struggled to achieve reliable signal propagation, another crucial dynamic property of the large-scale neural circuits. On the other hand, signals experienced significantly less attenuation when they operated within an alternative dynamical regime known as “global balanced amplification”33. However, in this regime, the timescale hierarchy phenomenon ceased to exist. The conundrum lies in reconciling both timescale hierarchy and reliable signal propagation within the cerebral cortex. Furthermore, despite the model being informed by anatomical data, it remains unclear to which extent the model accurately captures the timescales of individual areas and the precise properties of signal propagation in a quantitative manner. A challenge stems from the absence of neuronal activity data recorded across multiple areas in the macaque neocortex.

In this study, starting from analyzing the multi-regional electrocorticography (ECoG) recording data of another primate, the marmoset, we investigated two critical aspects within the marmoset neocortex: the timescale hierarchy phenomenon and the property of signal propagation. It is notable that there has been limited prior research reporting the existence of timescale hierarchy in the marmoset neocortex. Consequently, it becomes imperative to explore whether this phenomenon also manifests in the marmoset, thereby testing the hypothesis that timescale hierarchy may serve as a universal principle governing cortical dynamics across different species. In addition, the availability of comprehensive anatomical data25,26 along with activity data34 collected from multiple cortical areas within the marmoset neocortex, provides us with an opportunity to develop a biologically realistic large-scale model that encompasses multiple cortical areas. Importantly, and in contrast to previous research, the predictions made by our model regarding the timescales of individual cortical areas and the patterns of signal propagation across multiple areas can be directly validated using experimental activity data. This anatomically and dynamically constrained model empowers us to delve deeply into the role played by the intra-areal heterogeneity (i.e., macroscopic gradient of excitation) and inter-areal connectivity in shaping local signal integration and global signal propagation dynamics within the marmoset neocortex, in particular, leading to a near-critical state of the neocortex. In addition, the model mechanically explains recent experimental observations that the structural connectivity and functional connectivity are less similar in association areas than in sensory areas.

Results
Experimental evidence of timescale hierarchy in the marmoset neocortex
We initiated our investigation with the analysis of an ECoG dataset recorded from an anesthetized marmoset in a previous study34. This dataset includes signal measurement from 96 electrodes, broadly placed on the left hemisphere of the marmoset’s brain during a state of passive listening (Fig. 1a). To circumvent any potential bias from the auditory stimulus, we selected data recorded prior to any auditory stimulus exposure, which was presumably close to the resting state (for further details, see Section “Methods”).

We adapted a spectral approach to estimate the timescales of resting-state neural activity from the ECoG dataset11. This method focuses on the Power Spectral Density (PSD) of the data, enabling reliable estimations of the time constant of neural activity from relatively short time series. In brief, we calculated the PSD of the activity data using Welch’s method, followed by a Lorentz function fitting to determine the “knee frequency” and the corresponding decay time constant (see details in Section “Methods”). Previous studies anticipated an exponent of 2 for the Lorentz function36, which has also been observed in our study, showing that the exponent clustered around 2 with minor variance in Fig. 1b.

By applying the spectral method, a prominent timescale gradient across 33 cortical areas was revealed. As illustrated in Fig. 1c, sensory areas including V1, V2, and AuA1 exhibited smaller time constants, while association areas, including A8aD, A10, and PG exhibited longer timescales. This trend was further substantiated by computing the average timescale for each category of areas, presenting a significant increase of timescales from the low-level sensory areas to the high-level association areas, as shown in Fig. 1d.

Intriguingly, our data analysis unveiled a striking consistency between the timescale hierarchy and the structural hierarchy assessed by the spine number of pyramidal neurons in each area1,12. As demonstrated in Fig. 1e, the timescale correlates strongly with the spine count across cortical areas. Furthermore, by analyzing the first 5 s of awake resting-state recordings37, we observed a clear intrinsic timescale hierarchy across cortical areas closely matching the anesthetized state (Fig. S9a, b), with a strong correlation confirming consistency in the relative ordering of timescales (Fig. S9c).

Those findings propelled us to construct a multi-regional model of the marmoset neocortex that links the structural and dynamical property of the neocortex. With such a model, we were able to investigate how network structure shapes signal integration reflected by the neural time constant within each area and signal propagation across areas in the cerebral cortex – in particular, the role of the macroscopic gradient of spine count and the inter-areal connectivity.

A multi-regional model of the marmoset neocortex
We next developed a multi-regional model of the marmoset neocortex. The form of the model was adapted from Ref. 1 that described the macaque neocortex, considering the anatomical resemblance between the macaque and marmoset cortices. The model encapsulates 55 distinct sensory and association areas of the marmoset neocortex, including the 33 areas recorded in the previous experiment34. Each area in our model encompasses an excitatory and inhibitory neural population, both of which are locally interconnected with their activity described by the population firing rates. A detailed schematic of this multi-regional model is presented in Fig. 2a, b.

The network connectivity of the model was determined by the up-to-date anatomical measurement using retrograde tracing25,26. The magnitude of the directed connection between any two areas was derived from FLN determined through tracer injection experiments. The FLNs quantify the fraction of neurons projecting to a target area across all source areas, which can be viewed as a proxy of connectivity strength between two areas in our model. A heatmap representation of the FLN matrix is shown in Fig. 2c. The dynamical interactions within each area are described following a Wilson–Cowan type model (see further details in Section “Methods”).

In addition, the model accentuates the influence of the excitation gradient in scaling connection strengths between neural populations. Following prior macaque neocortex modeling1, we integrated a composite gradient of excitation level across areas in the model, which mainly reflects the areal specific spine numbers of pyramidal neurons based on a previous anatomical study26. In our model, this composite gradient index has been additionally constrained to align with the experimentally measured timescales deduced from the ECoG data, as detailed in the preceding section. A detailed methodology for estimating the composite gradient is described in Section “Methods”. The composite gradient’s bar plot is shown in Fig. 2d, alongside the average composite gradient for each category of cortical areas in Fig. 2e, revealing a statistically significant upward trend from visual areas to the prefrontal cortex. Furthermore, Fig. 2f illustrates the correlation between the composite gradient and the spine number of pyramidal neurons for all areas26, demonstrating that the composite gradient index was predominantly determined by the spine number of pyramidal neurons, thus as a proxy of excitation level for neurons in each area.

While subcortical structures (e.g., thalamus, brainstem, basal ganglia, cerebellum) and peripheral sensory pathways (e.g., retina, cochlea, spinal cord) exert influences on cerebral cortical dynamics, the present model is restricted to cerebral cortical networks due to the absence of a marmoset-wide, high-resolution subcortical and peripheral connectome. It is worth noting, however, that this work focuses on timescales of spontaneous neural activity during the resting state, when there is no experimentally designed sensory stimulation to peripheral structures, and animals are not engaged in a particular task. As future anatomical data allow the integration of subcortical and peripheral circuits, these additional components are expected to improve our model predictions of timescales and signal propagations, with the specific patterns likely influenced by interactions involving subcortical and peripheral structures.

Emergence of timescale hierarchy in the marmoset network model
We next investigated the capability of our model in replicating the experimentally observed timescales for each area in the marmoset neocortex. In simulating the model’s resting state, we applied white-noise stimuli across all brain areas, to emulate the resting state scenario in experiment. The neural activities in various areas, illustrated in Fig. 3a, reveal that lower-level areas resonate at higher frequencies, akin to the input white noise. In contrast, higher-level areas exhibit smoother variations, reflecting their differing operational timescales. Further inspection of the auto-correlation function of activity within each area, as shown in Fig. 3b, demonstrates an evident exponential decay, though with significant variability in the decay time constants across areas, thus reinforcing the existence of different operational timescales in each area.

We then examined the timescale gradient more closely by estimating timescales from the PSD of the simulated neural activity, following the same approach used for the experimental data processing. As displayed in Fig. 3c, the timescale for each area forms a gradient ranging from 50 ms to 250 ms. This gradient rises from sensory to high-level areas, echoing the experimentally observed pattern. Next, we quantitatively compared the experimentally measured and model-predicted timescales area by area. Figure 3c shows a strong correlation between them, with a scaling factor around 5. This scaling is consistent with the findings in ref. 11 where the timescale of ECoG was compared with the timescale of neuronal firing activity. These outcomes suggest our model well captures the experimental observations, thus validating the model’s effectiveness. Moreover, the timescale gradient and its close consistency with experimental data remain robust even when biologically realistic axonal conduction delays are incorporated into the model (Fig. S8). Finally, Fig. 3d demonstrates that the consistency between the time constant and the spine number across areas is also preserved in our model.

Furthermore, we substantiated the property of timescale gradient by examining the eigenmodes in the modeled system. For each eigenmode (a column of the matrix in Fig. 3e), each entry implies that the neuronal activity in the corresponding area will possess a timescale determined by the eigenvalue of that eigenmode. Thus, a sparse eigenmode indicates that only a few areas exhibit a particular timescale determined by the corresponding eigenvalue, supporting the concept of “timescale localization”. Figure 3e illustrates the phenomenon of timescale localization as indicated by the eigenmode matrix’s resemblance to a block diagonal matrix. Such an observation is instrumental as it indicates the capacity for lower-level sensory areas to respond quickly to incoming stimuli and discard older inputs. Conversely, higher-level planning and decision-making areas retain information over extended periods, aiding in integrative processes such as working memory.

We further studied the timescale gradient in the presence of stimulus in the model. By recording the post-stimulus response for all areas after activating the V1 area, Fig. 3f demonstrates that activity in high-level areas (PG, A6DR, A6VA, etc.) persisted longer than in lower-level visual areas (V1, V2, V4, etc.), reflecting a functional facet of the timescale gradient. A further quantitative investigation of this gradient by an exponential fit to the activity data (details in Section “Methods”) was depicted in Fig. 3g. Intriguingly, the auditory areas, typically considered lower-level but not directly involved in visual stimulus transmission, showed a longer timescale. This observation can potentially be attributed to the absence of direct connections between visual and auditory cortices. Hence, when the visual cortex is stimulated, the auditory cortex may be indirectly activated through high-level cortical areas operating over longer timescales. As a result, the auditory cortex, receiving inputs over extended periods, also exhibits these longer timescales, as depicted in Fig. 3g.

Figure S1a presents a direct comparison of timescales during rest and post-V1 stimulus, revealing a marked similarity between the two timescales derived under different brain states. Notably, higher-level areas exhibit longer timescales than lower-level ones, with exceptions like the auditory cortex discussed previously. We attribute these differences to the external stimulus on specific brain areas. Further mathematical analysis and numerical simulations (Supporting Information and Fig. S1b–e) clarify that the post-stimulus timescale of area activity is influenced by the system eigenmodes’ intrinsic timescales and the stimulus-induced multi-regional activity pattern. Importantly, the external stimulus does not alter the network’s intrinsic timescale or eigenmode pattern but modifies the contribution of each intrinsic timescale (and the corresponding eigenmode) to the activity of cortical regions.

Essential conditions on timescale localization
We next studied the mechanism underlying the property of timescale localization. A previous theoretical study of the macaque neocortex31 has proposed essential conditions for achieving timescale localization: the macroscopic gradient of excitation across areas and the balance of excitatory and inhibitory inputs received by each area. The macroscopic gradient refers to the gradual varying excitation strength across different cortical areas indicated by spine number of pyramidal neurons in each area, while “E-I balance” alludes to the delicate current cancellation between global long-range excitatory inputs and local inhibitory synaptic inputs to the excitatory neural population in each area.

To examine these two conditions in the marmoset cortical network, we first introduced two metrics for eigenmodes – the Inverse Participation Ratio (IPR) and the θ index to quantify the timescale localization and its spatial localization, respectively (see both definitions in Section “Methods”, θ was initially proposed in Luis Carlos Garcia and Xiao-Jing Wang, see38). A higher IPR value (closer to 1) indicates a particular timescale localized in fewer number of areas, and a higher θ index (closer to 1) describes a timescale localized in physically more proximal areas. For the model detailed previously, the scatterplot of IPR and θ for each eigenmode of the model system is shown in Fig. 4a, while the timescales of each eigenmode is indicated by the color code. Examples of eigenmodes with different IPR and θ are presented in Fig. S2. As shown in Fig. 4a (left panel), both high IPR and θ are observed for most eigenmodes, demonstrating that the timescale localization is present in the model, and each timescale is spatially localized in physically neighboring areas.

Subsequently, we examined the conditions for timescale localization in the marmoset neocortex. By removing the macroscopic composite gradients (setting the gradient h i to zero for all areas in the model), we observed a marked degradation of timescale localization and spatial localization, as evidenced by significantly smaller IPR and θ (Fig. 4a, right panel). Notably, the range of eigenmode timescales becomes narrower and more concentrated (predominantly around 50 ms), underscoring the influence of the composite excitation gradient in extending timescales. Further exploration involved adjusting the gradient’s slope (multiplying the gradient h i with a factor γ for all areas in the model) and observing the resulting changes of timescale localization. We found a progressive increase in IPR from 0.1 to 0.65 and θ from 0.65 to 0.95 as the gradient slope γ varies from 0 to 1 (Fig. 4b), supporting the pivotal role of the macroscopic gradient of excitation in timescale localization and its spatial distribution. The effect of the gradient of excitation becomes even more apparent when comparing the eigenmodes of the network with and without the gradient. Figure 4c, d showcase the fastest eigenmode for both scenarios (edged in black in Fig. 4a). While both eigenmodes predominantly involve low-level visual areas with V1 activity dominating the eigenmode, the eigenmode in the control condition is markedly sparser, highlighting the dominant role of V1 activity. This comparison demonstrates the composite gradient’s contribution to timescale localization, even within the low-level sensory cortex.

In a parallel model simulation, we disrupted the E-I balance in our model by increasing the strength of local connections from inhibitory to excitatory neurons (w EI). This modification lead to the imbalance between the global long-range excitatory inputs and the local inhibitory synaptic inputs. Consequently, it resulted in smaller IPR and θ that reflect poor timescale localization and spatial localization (Fig. 5a), similar to the case of gradient removal. With the E-I balance disturbed, the eigenmode characterized by the shortest timescale no longer involved the primary visual cortex (Fig. 5b). This shift underscores the critical function of E-I balance, not merely as an essential condition for timescale localization but also as a determinant in the gradient of timescales from low-level sensory to high-level brain areas.

In our continued exploration, we probed the influence of long-range connections on timescale localization. When connectivity strengths (FLNs) were randomly shuffled, as illustrated in Fig. 5c, the cortical network preserved high IPR values similar to the control condition (Fig. S3), indicating robust timescale localization. However, θ index reduced significantly, suggesting a decrease in spatial localization of timescales. This effect is also evident using the example of the slowest eigenmode in both the shuffled FLN network and the control condition, as illustrated in Fig. 5d and Supplementary Fig. S2e, respectively. Notably, while both eigenmodes encompass a variety of brain regions, including notably the prefrontal cortices, the eigenmode under the shuffled FLN condition exhibits a broader distribution, extending to primary visual and auditory cortices. In contrast, the eigenmode in the control condition demonstrates greater spatial localization, predominantly engaging high-level brain areas.

Besides the reduced spatial localization, the range of timescales was also narrowed compared to the control condition. These observations were further substantiated by the results from 1000 repetitions of FLN shuffling. Disrupting these long-range connections led to a noticeable contraction in the timescale range and diminished spatial localization, as evidenced in Fig. 5e, f. These results underscored the significant role of long-range connections in timescale localization, i.e., nearby areas process signals with similar timescales reflected by high θ value, while areas with some distance possess a variety of timescales enabling hierarchical signal processing. It suggests that the organization of the inter-areal connections is not random but is rather essential for brain functionalities.

Impact of near-critical state on structural and functional connectivity
As the macroscopic gradient of excitation largely influences timescales across areas, it is expected to influence the functional connectivity (FC) that relies on the temporally co-activation of pairwise areas. On the other hand, the macroscopic gradient also adjusts the effect of structural connectivity (SC) as the FLN will be scaled by the composite gradient as a proxy of the spine number of neurons in each area. Therefore, in this section, we direct our attention to how the macroscopic composite gradient shapes the correlation between SC and FC.

We initiated this investigation by removing the gradient from our model and computing the resultant FC. Under these conditions, the simulated FC from our model significantly aligns with the SC. This increased correlation, contrasted with the dissimilar case in the presence of gradient of excitation (visualized in Fig. 6a), suggests that the composite gradient can contribute to high-order connections within FC. These findings underscore that the composite gradient, as a fundamental structural attribute, can significantly shape the functional interplay between different brain areas. This phenomenon was consistent with that observed in the macaque neocortex model1, yet the mechanism remains to be elucidated.

Here we further developed a mathematical analysis to identify a quantitative relationship between FC and SC and investigated their relationship when the macroscopic gradient of excitation exists, following the analysis in ref. 39. As shown in Fig. 4a, the composite gradient gives rise to a wide range of timescale pool, in which the slowest timescale can be as long as hundreds of milliseconds or longer. The very slow timescale attributes to an eigenvalue close to zero, which indicates that the cortical network operates close to the critical state (the edge of stability). Building upon these facts, we mathematically analyzed the effects of criticality on the SC and FC relation (see Section Supporting Information), proving that the dissimilarity between FC and SC substantially increases when the network state approaches criticality.

We next scrutinized the role of the composite gradient on network criticality by adjusting its slope, akin to the model simulations conducted when investigating timescale localization. Specifically, we manipulated the slope of the gradients in two distinct ways: by either maintaining the average gradient value constant or allowing it to vary. Figure 6b (top panel) displays examples of these gradient adjustments for various gradient slope by each method. As depicted in Fig. 6b (bottom panel), the dissimilarity between FC and SC intensifies with steeper gradient slope, confirming our hypothesis that composite gradients can contribute to the system’s criticality and, subsequently, the dissimilarity between FC and SC. We further explored the gradient’s function by maintaining the average excitability of all areas as a constant. Remarkably, even when the average excitability remains constant at the whole-cortex level, composite gradients can enhance local excitability and increase the dissimilarity between FC and SC (Fig. 6b, bottom panel).

Next, we turned to analyze the impact of the composite gradient on individual cortical areas. Through mathematical analysis, we proved that the correlation between FC and SC is significantly lower in association areas compared to early sensory areas (see Section Supporting Information), which has also been reported experimentally in a monkey study40. In the model simulation, we computed the dissimilarity between FC and SC for each cortical area and visualized these differences in a marmoset brain parcellation (Fig. 6c, left panel). The simulation results confirmed that dissimilarity is substantially higher in association areas than in early sensory areas. This trend is not evident when the composite gradient is absent (Fig. 6c, right panel).

To further quantify this observation, we identified a positive and significant correlation between FC-SC dissimilarity and the rank of each area’s composite gradient (Fig. 6d, left panel). In contrast, this correlation disappears when the gradient slope is zero (Fig. 6d, right panel), indicating that this effect is not solely due to the heterogeneity of inter-areal connectivity. We further validated this by categorizing marmoset brain areas into seven subnetworks, comparable to those reported in humans41 and macaques42. As shown in Fig. 6e, FC-SC dissimilarity is relatively small for lower-level visual and somatic motor subnetworks, but becomes more significant in higher-level association subnetworks. These findings are consistent with previous findings on human cortex42,43 and are absent when the gradient slope is removed (Fig. 6e, right panel). This gradient-facilitated feature underscore the importance of higher-level areas in establishing wide-ranging functional connections, thereby supporting the execution of complex cognitive functions. These insights reaffirm the crucial role that gradients play in brain functionality.

The critical state of the model not only depends on the composite gradient of excitation across areas but also on the inter-areal connection among areas. To study the effect of inter-areal connection on the FC-SC relation, we first removed the gradient from our model and studied the impact of increasing the inter-areal coupling strength of excitatory inputs. As the coupling strength escalates, the network state gets closer to criticality, resulting in an increased deviation between FC and SC. This trend, as depicted in Fig. S4, is also consistent with our mathematical analysis (see Section Supporting Information).

Signal propagation in the near-critical state
In this section, we first investigated how the signal propagates across areas within our model, drawing parallels with observations from a previous optogenetics experiment37. In the experiment37, an input was given to the marmoset A4ab area via an LED optogenetic stimulation, and ECoG data were recorded from 64 electrodes broadly distributed on the marmoset neocortex. We replicated this experimental stimulus protocol within our model and analyzed the peak responses, defined as the highest level of neural activity observed in each cortical area post-stimulus. As shown in Figs. 7a and S6, our model well aligns with the experimental measurements for the majority of cases without additional parameter tuning, affirming the physiological relevance of our model’s responses. Similarly, this strong agreement with experiments was conserved after incorporating biologically realistic axonal conduction delays into our model (Fig. S8).

We then examined the necessary conditions for signal propagation within the model starting from a pulse stimulus applied to sensory cortices, such as the primary visual cortex (V1), a central source of visual input. Responses from critical cortical areas were depicted in Fig. 7b, showing steady signal propagation from the visual cortex to high-level areas like the prefrontal cortex.

We next investigated the role of inter-areal connectivity for signal propagation. The connectivity can be classified into feedforward and feedback connections based on anatomy25,26. As intuitively expected, eliminating feedforward connections significantly hampers signal propagation from sensory areas (data not shown). We then shifted our focus to the role of feedback loops. These loops consist of connections from high-level areas to low-level areas identified by a smaller fraction of neurons in a projection originating from the supragranular layers of the source area (SLN < 0.5)25,26. As demonstrated in Fig. 7c (orange line), the deletion of all feedback loops within the FLN precipitates a generalized decrease in peak responses across all areas, with low-level areas bearing the brunt of this impact. Subsequent analyzes investigated the influence of the composite gradient of excitation on signal propagation. As shown in Fig. 7c (green line), the absence of the composite gradient led to a significant reduction in neural responses, particularly in higher-order areas. To explore this effect further, we systematically varied the gradient slope and examined the resulting changes in signal propagation, using peak responses from the PG area as a metric. As illustrated in Fig. 8a, increasing the gradient slope enhanced signal strength, highlighting the critical role of gradients in facilitating signal propagation.

The contrast in signal propagation effects between the removal of excitation gradient and feedback connections is notable. The composite gradient of the excitation enables the high-level brain areas to receive stronger recurrent excitatory input, substantially amplifying long-range inputs from lower-level areas. In contrast, feedback connections enable indirect interactions between different lower-level areas via top-down projections from higher-level areas, significantly benefiting signal propagation, especially between areas lacking direct connections.

Building on the central role of both feedback loops and composite gradients in driving the system toward a critical state, our study delved deeper into the interplay between system criticality and the coexistence of signal propagation and integration. We examined the effects of altering global excitatory-excitatory (E-E) connectivity (μ EE), while maintaining a constant global excitatory-inhibitory (E-I) strength (μ IE), on the signal propagation and integration. Our findings, presented in Figs. 8b, c, and S7, elucidates several key insights. Firstly, the criticality of the system, implemented by the balance between μ EE and μ IE, emerges as an important determinant of neuronal activities. An imbalance, particularly at lower μ EE levels, markedly impairs both signal propagation and integration. Furthermore, we observed a trade-off regime (grey areas in Figs. 8b, c, and S7) when μ EE approaches the critical threshold, highlighting the intricate balance between signal propagation and integration. Approaching this threshold enhances signal propagation and enlarges the range of timescales, albeit at the expense of reduced timescale localization. Additionally, our results demonstrate that under conditions of strong global interaction strength (solid lines in Figs. 8b, c, and S7), both signal propagation and integration sustain high baseline values within the critical zone. This suggests that strong global connectivity, nudging the system toward criticality, is crucial in facilitating the co-existence of global signal propagation and localized integration.

Methods
Experimental data sources
This work analyzed two experimental datasets to estimate neuronal activity’s time constants and signal propagation across cortical areas. The resting-state time constant was estimated from an ECoG dataset, gathered during an auditory task from a marmoset with a chronic implant in the left hemisphere34. The dataset features data from 96 uniformly placed electrodes throughout the hemisphere, gathered while the marmoset was in a resting state, under the influence of ketamine-based anesthesia (30 mg/kg i.m.). For the purpose of this study, we considered only the data collected prior to any auditory stimulus—~15 s per trial—resulting in a cumulative total of 60 s of data compiled from four separate trials. We analyzed neural signals in 33 areas whose inter-areal connectivity has been measured in anatomical experiments25,26. The activity of these areas can be directly modeled, which enables us to compare experimental data with simulation results.

For the analysis of signal propagation across areas, we utilized an optogenetic ECoG dataset collected from the right hemisphere of an awake marmoset37. This dataset was gathered using a combination of 64 epidurally implanted electrodes and eight LEDs designated for optogenetic stimulation. We utilized recording data with a single LED stimulating approximately areas A4ab, A6DR, A6Va, AuCM, S2E, V3A, PE, and PG. Each recording trial incorporated a 200 ms pulse stimulus followed by a 2 s post stimulation interval, and the response was quantified by averaging the outcomes across a series of 50 trials.

For the intrinsic timescales during the awake resting state, we utilized the same optogenetic ECoG dataset for signal propagation analysis37. We specifically extracted the initial 5 s of recordings from the whole session—prior to any optogenetic stimulation – to obtain a resting–state segment free of task or laser effects. Nonetheless, due to the relatively short duration of awake-state recordings, we primarily used anesthetized-state data for subsequent analysis (e.g., comparisons with our model).

Estimation of the time constant of neural activity
We derived the time constants of neural activity during the resting state using a methodology based on PSD data11. This approach is especially beneficial when data length is short and the data contains oscillations and artifacts. The process entailed several stages. Initially, we computed the PSD using a variant of Welch’s method52 that involved 1-s long Hamming windows with 0.5-s overlap. For cases where multiple electrodes represented one area, we calculated the PSD individually for each electrode, subsequently taking the average as the PSD for that area.

Following this, we applied spectral parameterization53 to decompose the PSD into a sum of Gaussians (representing temporal periodic terms) and a Lorentz function L(f) (representing temporal aperiodic terms) given by the form

L(f) = A / (k + f^χ)

After fitting the data to determine the values of A, k, and χ, we identified the ’knee’ frequency approximately as f k = k^1/χ. This frequency signifies the location where the PSD exhibits a ‘knee’ or bend. We then computed the corresponding timescale as τ = 1 / (2π f k).

For the post-stimulus activity timescale, we fitted the response following peak activity (post-stimulus) using an exponential function

r(t) = a e^(-t/τ).

We determined the value of τ that provided the best fit as the time constant of the neural response in each cortical area.

Multi-regional cortical model architecture
We built the multi-regional marmoset cortical model using a set of ordinary differential equations, adapted from previous works of modeling macaque cortical network1,33. In the model, 55 areas from sensory to association cortex were included. Each cortical area was modeled with one excitatory and one inhibitory neuron group governed by the following dynamics:

τ E d/dt r i E = -r i E + β E [1 + η E h i] (w EE r i E - w EI r i I + Σ j=1 N [1 + η E h i] μ EE FLN ij r j E + I i ext,E )+

τ I d/dt r i I = -r i I + β I [1 + η I h i] (w IE r i E - w II r i I + Σ j=1 N [1 + η I h i] μ IE FLN ij r j E + I i ext,I )+

In these equations, r i E and r i I represent the firing rates of the excitatory and inhibitory populations in area i, respectively, while τ E and τ I are the corresponding intrinsic time constants. w XY symbolizes the coupling strength from Y population to X population within the area (X, Y ∈ {E, I}). μ X E represents the coupling strength of the excitatory input from the population in area j to the population in area i. The external input to the population is expressed as I i ext,X. We used a threshold linear function f(x) = β[x]+ for the gain function for both E and I populations, with gain factors β = β E and β I, respectively.

To complement our large-scale model that uses FLN to describe inter-areal connectivity, we additionally introduced a composite gradient of excitation, h i. Because the FLN matrix is normalized per target area (each row sums to 1), it encodes only relative inter-areal connectivity but does not capture the area-wise differences in the absolute number of synapse inputs received by each target area. Anatomical evidence indicates that neurons in higher-order association cortex receive more inputs per cell than those in early sensory areas. For instance, dendritic spine counts increase systematically from sensory to association areas26. Furthermore, a uniform ~80%/20% division of local versus long-range inputs across areas54 implies that higher-order areas inherently possess greater numbers of both local and inter-areal synapses. Thus, neurons in association cortex receive a larger overall drive than those in lower areas. To account for this gradient, we define a composite gradient of excitation h i as an area-specific scaling factor applied to all synaptic weights, effectively representing the increasing synaptic density in higher cortical areas while preserving the empirically measured FLN connectivity pattern. Within our model, h i is a factor scaling both local (intra-areal) and long-range (inter-areal) excitatory inputs and varies across different areas. The quantitative estimation of h i is detailed in the subsequent subsection. The scaling factor η X governs the slope of the composite gradient on the X population.

Adapted from a previous experimental study55, we set τ E = 20 ms, τ I = 10 ms, β E = 0.066 Hz/pA, β I = 0.351 Hz/pA, η E = 0.685, η I = 0.76, w EE = 24.4 pA/Hz, w EI = 19.7 pA/Hz, w IE = 11.66 pA/Hz, w II = 12.5 pA/Hz, μ EE = 67.4 pA/Hz, μ IE = 49.81 pA/Hz.

Incorporation of axonal conduction delays
We also extended the multi-area rate-based model to include explicit axonal conduction delays on all connections, both inter-areal and intra-areal for some studies. Inter-areal transmission delays, τ ij, were defined as proportional to the anatomical distance between the source area j and the target area i. Specifically, we set τ ij = d ij / v ax, where d ij is the white-matter tract length between areas i and j25 and v ax = 3.5 m/s is the assumed axonal conduction velocity. All intra-areal (local) recurrent projections were assigned a fixed delay of 2 ms. These parameter values align with those used in previous large-scale modeling studies of macaque cerebral cortex33, thereby ensuring biologically plausible transmission latencies consistent with existing literature.

Delays were implemented by introducing a time shift in the firing-rate equations (Eqs. (1) and (2)) for each synaptic input. Specifically, any synaptic input from population β in area j to population α in area i was applied with a time offset corresponding to the assigned delay. For inter-areal connections (i ≠ j), the presynaptic activity from population β at time t−τ ij was used when computing input to the postsynaptic population α at time t. Similarly, for local recurrent connections (i = j), population activity at t−2 ms was used in the coupling terms. This delayed-input mechanism was integrated the existing rate model dynamics without modifying any other parameters. The complete formula is as follows:

τ E d/dt r i E(t) = -r i E(t) + β E [1 + η E h i] (w EE r i E(t-2) - w EI r i I(t-2) + Σ j=1 N [1 + η E h i] μ EE FLN ij r j E(t-τ ij) + I i ext,E(t))

τ I d/dt r i I(t) = -r i I(t) + β I [1 + η I h i] (w IE r i E(t-2) - w II r i I(t-2) + Σ j=1 N [1 + η I h i] μ IE FLN ij r j E(t-τ ij) + I i ext,I(t))

Importantly, introducing these biologically realistic conduction delays did not substantially change the primary results of the model: the global inter-areal signal propagation patterns and the intrinsic timescale hierarchy remained robust and unchanged (Supplementary Fig. S8).

Estimation of macroscopic gradient of excitation
We incorporated the composite gradient of excitation h i in the model to be consistent with the experimental findings related to the gradient of timescale across areas and the structural hierarchy derived from the fraction of supragranular labeled neurons (SLN) as reported in ref. 26. As the model includes more areas (55 areas) than those recorded in the ECoG experiment (33 areas), the fitting was performed in two steps: firstly, we determined the composite gradient value of the 33 areas whose activity data and structural hierarchy data were both available. Secondly, we calculated the values of h i for the rest areas based on the structural hierarchy and the values of h i for the 33 areas obtained in the first step.

For the 33 areas whose timescale and structural hierarchy were both available, we used the following objective function that combined the similarity between the excitation composite gradient parameter h i with the structural hierarchy h i,exp determined by SLN and the timescale τ i,exp derived from neural signals, i.e.,

L(h) = Corr(τ exp, 1/(1-h)) + λ Corr(h exp, h)

where h, h exp and τ exp are vectors, and λ is a hyper-parameter that determines the preference for the gradient to be closer to the timescale or the experimental hierarchy (chosen as 1 in this work). The Corr(a, b) is the Pearson’s correlation coefficient between a and b. We assessed the similarity with the structural hierarchy h exp by measuring the correlation coefficient, and evaluated the similarity with the timescale τ exp by calculating the correlation coefficient between 1/(1-h) and the experimentally measured timescale. We made this choice based on a mathematical analysis showing that 1−h is proportional to the inverse of the time constant of neural activity when the multiple brain areas in the system are decoupled (details provided in Supporting Information). To obtain the values of the composite gradient, we solved the optimization problem presented by the objective function:

h* = argmax L(h)

For areas whose timescale data could not be obtained from experiments (no electrode is placed in these areas), we determined their composite gradient h i as a weighted average of the gradients of the 33 areas whose gradient h j have been estimated in the first step. The weight was based on the similarity between structural hierarchy of two areas:

h* i = Σ j=1 33 h* j e^(-d^2 ij / σ^2)

where d ij = |h i,exp - h j,exp| is the difference of structural hierarchy for area i, j and σ = 0.05.

Simulation scenarios
We examined multiple simulation scenarios in the study. The scenarios encompassed a range of parameter configurations and manipulations of the network connectivity. For scenarios related to the change of gradient, we removed gradients (Figs. 4a, d, and 6a) by setting all h i = 0. Furthermore, we manipulated the slope of the composite gradients by scaling the gradient h i with a factor γ (Figs. 4b, 6b) as h i(new) = h i * γ or by the gradient h i (Fig. 6b) h i(new) = (1-γ) * h + γ * h i, 0 ≤ γ ≤ 1.

For Fig. 6c, the gradient slope γ is set to 0.

In scenarios involving the change of network connectivity FLNs, we devised a condition devoid of any feedback loops by assigning zero to all FLN entries where the connection has a smaller fraction of neurons in a projection originating from the supragranular layers of the source area (Fig. 7c), namely FLN ij = 0, for all SLN ij < 0.5. For Fig. 5c–f, we shuffled the FLN by randomly permuting its entries 1000 times to obtain the distribution of the corresponding measures.

Regarding other parameter modifications, we disrupted the E-I balance by augmenting w EI by 10% (Fig. 5a). For Fig. 8b–d, we set μ IE = 24.91, 37.36, 49.81 pA/Hz, respectively and computed the corresponding upper bound μ* EE as 34.03, 50.68, 67.40 pA/Hz such that the system remains stable. μ EE was varied as λμ* EE with λ ranging from 0.95, 0.97, 0.975 to 1 for each respective μ IE setting. For Fig. S4, we varied the value of μ EE within the range from 64 to 70.

Metrics for timescale localization
We employed two quantitative measures to describe timescale localization across the whole neocortex. The first is the IPR, which was originally introduced in quantum mechanics and solid-state physics. For a normalized eigenvector v, its IPR is defined as follows:

IPR(v) = Σ j=1 N |v j|^4

where v j is the jth element of v. The IPR for the eigenvector matrix is computed as the average IPR across each column (eigenvector) in the matrix.

Furthermore, we introduced an additional metric for spatial localization of eigenvector v (Luis Carlos Garcia and Xiao-Jing Wang, see ref. 38):

θ(v) = Σ i,j=1 N |v i|^4 |v j|^4 e^(-d ij / d*)

Here, d ij represents the distance between a distance between a distance between a pair of areas i and j, which, in our context, denotes the physical distances between two areas in the marmoset neocortex; d* signifies the average distance across all pairs.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Thank you for inviting me to review this manuscript, in which the authors create a biophysical model of the cerebral cortex in an effort to explain observed variations in the timescale of electrophysiological data collected from the macaque brain.\nComments\n1. The authors use the term 'cortex' to refer to a particular cortical structure - namely, the cerebral cortex. It is important to disambiguate this from other cortical structures in the brain (such as the cerebellar cortex), which have very different organisational properties.\n2. The authors treat the cerebral cortex as though it is a disembodied structure, however neuroantomically, the structure is massively (albeit indirectly) connected to a wide range of peripheral (retina, cochlea, spinal cord) and subcortical structures that impose substantial computational constraints on the operation of the cerebral cortex. The assumption that features of the cerebral cortex *alone* are responsible for any observed variance in a measured statistic is an example of a mereological fallacy that should be corrected throughout the text.\n3. There is a substantial literature evidencing slightly-subcritical dynamics in neural recordings that is completely missing from this paper. Key examples include https://www.nature.com/articles/s41467-018-04725-4, https://pubmed.ncbi.nlm.nih.gov/25788679/, https://www.nature.com/articles/s41467-023-42465-2, https://www.nature.com/articles/ncomms3521, but this list is nowhere near exhaustive. Without the framing imposed by these manuscripts, the paper is unreasonably disconnected from a very relevant literature.\n4. Did the authors add temporal delays into their model - e.g., proportional to the distance traversed by the estimated white matter connections between regions?\n5. From the structural connectivity map in Fig 2c, it appears as though there are few (if any) local connections imposed in their model. The data shown in the source material (e.g., Ref #25) do appear to have many local connections, however their figures also appear to be plotted at a finer spatial resolution than the model used in this manuscript. This detail is important, as the density of local connections could change across the same hierarchy observed in this study, and thus may provide an alternative explanation for the observed results.\n6. The eigenmode/modelling analysis was informative and interesting.", "reviewer_2": "Li et al.'s manuscript \"A hierarchy of time constants and reliable signal propagation in marmoset cortex\" attempts to ground a multi-area cortex model in experimental ground truth to determine the connectivity between areas that can explain the hierarchy of timescales across the cortical processing hierarchy, and also enable reliable signal propagation. The model leverages marmoset EEG data that allows measurement of the timescale of cortical activity, and anatomical tracing data that includes the number of neurons projecting between the cortical areas. They found that when the network is operating in a state near criticality, both local integration (i.e. timescale hierarchies) and signal propagation are more optimal. Furthermore, the model mechanistically explains why the patterns of structural and functional connectivity are increasingly dissimilar across the cortical hierarchy.\nI think that the study is of great interest across systems and computational neuroscience, as it provides a mechanistic explanation for the prominent and systematic changes in activity dynamics (and integration) across the cortical processing hierarchy - and the anatomical basis for the model's connectivity is a major plus. I do have a few comments:\nMajor\n• In the data used to fit the model, the animals were anesthetized with ketamine. I worry that this could impact activity timescales, as well as the activity states and signal propagation through the network. Perhaps they authors could check in the awake data used for signal propagation that the timescales are still similar in that state?\n• Anatomical tracing data is used to determine FLN. Does this also somehow inform the relative interareal 'connection strength' onto E and I cell types, or were these chosen by some other means?\n• I think that the authors could make it clearer in places like the results and the methods on what basis the additional composite gradient h is needed (on top of FLN). The increasing number of spines in higher cortex seems to be the reason used. Is the implication then that a given neuron that projects to an area in higher cortex may branch its axons more and form more synapses than a neuron projecting to a 'lower' area? Or is it that more spines = more local recurrence = amplified impact of an input from another area? Further justification and explanation would be helpful in justifying the use of the composite gradient.\nMinor\n• Some typos:\no line 2 \"information process\".\no Line137\"descried\"\no Fig6:\"excitability\"ismisspelledinafewpanels\n• What about cortico-thalamo-cortical connections: can the authors comment on how these could play a role in dynamics and their timescales? Are these more or less captured in the 'mismatch' between correlations and structural connectivity?\n• It would be helpful to have more information about the structural data used in this study. Spine counts are used and correlated to the activity timescales - are these total spine counts on entire single neurons, or on specific dendrites? In which layer?\n• Fig 6d, e: It is a little confusing for the reader to have the dissimilarity axes in these two panels go from low to high and high to low."}, "author_response": "Response to Reviewers\n\nReviewer #1 (Remarks to the Author):\nThank you for inviting me to review this manuscript, in which the authors create a biophysical model of the cerebral cortex in an effort to explain observed variations in the timescale of electrophysiological data collected from the macaque brain.\nComments\n1. The authors use the term 'cortex' to refer to a particular cortical structure - namely, the cerebral cortex. It is important to disambiguate this from other cortical structures in the brain (such as the cerebellar cortex), which have very different organisational properties.\nResponse: We appreciate the reviewer’s comment on the use of the term ‘cortex’. Following the reviewer’s suggestion, we have replaced the generic term “cortex” with “cerebral cortex” or “neocortex” to avoid ambiguity throughout the entire manuscript—including the title, abstract, main text and figure captions. This change ensures that all references clearly pertain to the cerebral cortex thus it will not be confused with other cortical structures such as the cerebellar cortex.\n2. The authors treat the cerebral cortex as though it is a disembodied structure, however neuroantomically, the structure is massively (albeit indirectly) connected to a wide range of peripheral (retina, cochlea, spinal cord) and subcortical structures that impose substantial computational constraints on the operation of the cerebral cortex. The assumption that features of the cerebral cortex *alone* are responsible for any observed variance in a measured statistic is an example of a mereological fallacy that should be corrected throughout the text.\nResponse: We thank the reviewer for pointing out the influence of peripheral and subcortical circuits on cortical computation. We fully agree that the cerebral cortex does not operate in isolation and that peripheral and sub‑cortical nuclei impose important computational constraints on cortical dynamics. Our current model uses the cerebral cortical connectome because, to date, the marmoset connectome encompassing the pertinent thalamic, brain‑stem, and sensory‑peripheral pathways is not available yet at a resolution comparable with the cortical data we use. Incorporating provisional or cross‑species estimates of the complete connectome would introduce additional assumptions that cannot be validated empirically. It is worth noting, however, that this work focuses on timescales of spontaneous neural fluctuations, when there is no experimentally designed sensory stimulation to peripheral structures and animals are not engaged in a particular task. Furthermore, the present work and Chaudhuri et al. 2015 previously demonstrate that a hierarchy of time constants can arise in a model of neocortex, as a robust phenomenon, although its precise pattern is likely to also shaped by interactions with subcortical structures like the thalamus, basal ganglia and cerebellum.\nIn addition, as shown in Fig. 3c and Fig. 7a, although there are strong correlations between the experimental and modeling results, our model does not perfectly reproduce the timescales and response profiles observed in each area. This discrepancy may stem from the influence of peripheral and subcortical regions, which are not currently included but can be incorporated in future work when a more complete connectome becomes available. To avoid mereological fallacy, we have revised the manuscript throughout (Introduction, lines 2-5, Discussion, lines 530-533) to replace wordings that could be interpreted as attributing observed variance solely to cerebral cortical circuitry. For instance, the opening lines of the Introduction (lines 2–5) now read as follows: “The cerebral cortex, together with peripheral sensory organs (e.g., retina, cochlea, spinal cord) and subcortical relays (e.g., thalamus), engages in information processing of external inputs, local integration of signals within regions, and the transmission of these signals along functionally specific pathways.” Furthermore, we have revised the relevant sentences in the Discussion (lines 530-533) as follows: “This conservation across species emphasizes the robustness of the model framework as a sturdy base for future research probing the relationship between structure and function at a whole-brain level, including the interplay between cortical dynamics and subcortical or peripheral neural systems.”\nWe further include a paragraph in the Results “2.2 A multi-regional model of the marmoset neocortex” (page 7, line 172-182) to emphasize that our results describe the contribution of cerebral cortical architecture in the absence of subcortical and peripheral inputs. The paragraph is as follows: “While subcortical structures (e.g., thalamus, brainstem, basal ganglia, cerebellum) and peripheral sensory pathways (e.g., retina, cochlea, spinal cord) exert influences on cerebral cortical dynamics, the present model is restricted to cerebral cortical networks due to the absence of a marmoset-wide, high-resolution subcortical and peripheral connectome. It is worth noting, however, that this work focuses on timescales of spontaneous neural activity during the resting state, when there is no experimentally designed sensory stimulation to peripheral structures and animals are not engaged in a particular task. As future anatomical data allow the integration of subcortical and peripheral circuits, these additional components are expected to improve our model predictions of timescales and signal propagations, with the specific patterns likely influenced by interactions involving subcortical and peripheral structures.”\nWe also added a paragraph in the Discussion (page 22, line 586-597) noting the potential contribution of peripheral and subcortical circuits as follows: “While our model focuses exclusively on cortico-cortical connectivity, the cerebral cortex is not anatomically isolated -- subcortical circuits (e.g., thalamus) and peripheral sensory pathways (retina, cochlea, spinal cord) impose important computational constraints on cortical activity. However, due to the lack of a marmoset-wide subcortical connectome at comparable resolution, our model necessarily focused on the corticocortical network, which can later be augmented with subcortical circuits as new data emerge. We anticipate that incorporating cortical-thalamo-cortical (CTC) and peripheral pathways in future marmoset models will refine our present findings, particularly by enhancing model predictions related to intrinsic timescale patterns, signal propagation, and functional connectivity profiles. Incorporating such pathways represents a key direction for future work as compatible anatomical datasets become available.”\nReferences:\n1. Chaudhuri, R., Knoblauch, K., Gariel, M. A., Kennedy, H., & Wang, X. J. (2015). A large-scale circuit mechanism for hierarchical dynamical processing in the primate cortex. Neuron, 88(2), 419-431.\n3. There is a substantial literature evidencing slightly-subcritical dynamics in neural recordings that is completely missing from this paper. Key examples include https://www.nature.com/articles/s41467-018-04725-4, https://pubmed.ncbi.nlm.nih.gov/25788679/, https://www.nature.com/articles/s41467-023-42465-2, https://www.nature.com/articles/ncomms3521, but this list is nowhere near exhaustive. Without the framing imposed by these manuscripts, the paper is unreasonably disconnected from a very relevant literature.\nResponse: We thank the reviewer for pointing out these important body of works. We have now integrated them and other relevant literatures explicitly into the Discussion (page 20, line 477-495). The new paragraph is as follows:\n“Converging evidence across multiple species and brain regions indicates that cortical dynamics operate near a critical point, yet remain slightly subcritical (Petermann et al., 2009; Priesemann et al., 2014; Wilting & Priesemann, 2018). This slight subcriticality is exemplified by a consistent branching ratio of ~0.98 observed in vivo (Wilting & Priesemann, 2018) and truncated neuronal avalanche size distributions in awake cortex (Petermann et al., 2009; Priesemann et al., 2014). Notably, the precise operating point can be modulated by brain states: focused attention shifts cortical activity further into the subcritical state from the near-critical state during resting (Fagerholm et al., 2015), while changes in arousal or anesthesia can also influence the distance from criticality (Munn et al., 2023; Maschke et al., 2024). Operating just below criticality is hypothesized to confer several key functional advantages—ensuring stability against runaway excitation, preserving long-range correlations and large repertoires of activity patterns, enabling efficient information processing across a critical-like broad dynamic range, and allowing flexible tuning of network sensitivity (Moretti & Muñoz, 2013; Priesemann et al., 2014). Theoretically, the cortex’s hierarchical modular architecture could yield an extended critical-like regime on the subcritical side (a Griffiths phase) that provides these functional benefits without fine-tuning (Moretti & Muñoz, 2013). Consistent with this \"subcriticality\" framework, our modeling results for the marmoset cortical demonstrate that a near-critical (slightly subcritical) state supports robust but controlled global signal propagation and naturally extended intrinsic timescales, highlighting how slight subcriticality effectively balances integration and stability while avoiding the fragility associated with exact criticality.”\nWe hope these additions situate our study within existing literature and address the reviewer’s concern about framing.\nReferences:\n1. Petermann, T., Thiagarajan, T. C., Lebedev, M. A., Nicolelis, M. A., Chialvo, D. R., & Plenz, D. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proceedings of the National Academy of Sciences, 106(37), 15921-15926.\n2. Priesemann, V., Wibral, M., Valderrama, M., Pröpper, R., Le Van Quyen, M., Geisel, T., ... & Munk, M. H. (2014). Spike avalanches in vivo suggest a driven, slightly subcritical brain state. Frontiers in systems neuroscience, 8, 108.\n3. Wilting, J., & Priesemann, V. (2018). Inferring collective dynamical states from widely unobserved systems. Nature communications, 9(1), 2325.\n4. Fagerholm, E. D., Lorenz, R., Scott, G., Dinov, M., Hellyer, P. J., Mirzaei, N., ... & Leech, R. (2015). Cascades and cognitive state: focused attention incurs subcritical dynamics. Journal of Neuroscience, 35(11), 4626-4634.\n5. Munn, B. R., Müller, E. J., Medel, V., Naismith, S. L., Lizier, J. T., Sanders, R. D., & Shine, J. M. (2023). Neuronal connected burst cascades bridge macroscale adaptive signatures across arousal states. Nature Communications, 14(1), 6846.\n6. Maschke, C., O’Byrne, J., Colombo, M. A., Boly, M., Gosseries, O., Laureys, S., ... & Blain-Moraes, S. (2024). Critical dynamics in spontaneous EEG predict anesthetic-induced loss of consciousness and perturbational complexity. Communications Biology, 7(1), 946.\n7. Moretti, P., & Muñoz, M. A. (2013). Griffiths phases and the stretching of criticality in brain networks. Nature communications, 4(1), 2521.\n4. Did the authors add temporal delays into their model - e.g., proportional to the distance traversed by the estimated white matter connections between regions?\nResponse: We thank the reviewer for raising this point. Following the reviewer’s suggestion, in the revised manuscript, we have extended the model to include explicit axonal-conduction delays both within and between cortical areas.\nThe Implementation is as follows:\ni. Inter-areal delays were computed as , where is the measured white-matter tract length between areas and (Majka et al. 2020) and a constant axonal conduction velocity was assumed. = 3.5 m/s\nii. Intra-areal delays were set to .\niii. Both the values of axonal conduction velocity and intra-areal delays follow the ones used in our earlier macaque modeling study (Joglekar et al. 2018).\nNoticeably, introducing these delays does not alter our main conclusions:\ni. Effect on signal propagation. Figure R1a–c (Supplementary Figure 8a–c) shows that introducing these delays produces the expected tens-of-milliseconds latency of activity onset after stimulus presentation (panel a), yet the overall spatiotemporal pattern, firing-rate magnitude and decay profile — remains virtually unchanged relative to the original delay-free model (compare panel b in Fig. R1 with Fig. 7a, left, in the main text). The delayed model further reproduces the experimentally measured population response with a high correlation (panel c).\nii. Effect on resting-state temporal hierarchy, Figure R1d–f (Supplementary Figure 8d-f). During spontaneous activity, cortical areas still exhibit similar autocorrelation decay profiles (panel d) and a robust timescale hierarchy (panel e). Time constants derived from the delayed model and the original model are highly correlated (r = 0.97; panel f), indicating that the principal conclusions regarding localization of intrinsic timescales are insensitive to the addition of biologically realistic conduction lags.\nAll modelling choices and equations have been added to “Methods” (page 25-26, line 686-710) as the subsection “4.4 Incorporation of axonal conduction delays,” and the model results shown in Fig. R1 are provided in Supplementary Figure S8 (page 46). We explicitly referenced these results within the main text, specifically in the Results “2.3 Emergence of timescale hierarchy in the marmoset network model” (page 9, lines 206-208): “Moreover, the timescale gradient and its close consistency with experimental data remain robust even when biologically realistic axonal conduction delays are incorporated into the model (Figure S8).” Similarly, in Results “2.6 Signal propagation in the near-critical state” (page 16, lines 403-405), we note: “Similarly, this strong agreement with experiments was conserved after incorporating biologically realistic axonal conduction delays into our model (Figure S8).” We hope these additions have addressed the reviewer’s query about temporal delays.\nFig. R1. Results from the model incorporating axonal conduction delays. (a) Activity of representative areas following V1 stimulation – focusing the first 20 milliseconds after stimulus onset. Solid line: the activity of model incorporating delays; Dash line: the activity of model without delay. (b) Activity of representative areas in the model incorporating delays following A4ab stimulation. (c) Comparison of the delay-incorporated model's peak response with optogenetics experimental data, demonstrating high model-experimental data consistency ( ). (d) The autocorrelation function of each area's activity in the resting state. (e) Hierarchical timescale−s 6extracted from simulated resting state activity r = 0.81,p = 4.57×10 of the model incorporating delays. (f) Direct comparison of intrinsic timescales between the delay-incorporated and delay-free models, demonstrating a robust consistence ( ) and indicating minimal impact of conduction delays on the overall timescale hierarchy. −32 r = 0.97,p = 1.80×10\nReferences:\n1. Majka, P., Bai, S., Bakola, S., Bednarek, S., Chan, J. M., Jermakow, N., ... & Rosa, M. G. (2020). Open access resource for cellular-resolution analyses of corticocortical connectivity in the marmoset monkey. Nature communications, 11(1), 1133.\n2. Joglekar, M. R., Mejias, J. F., Yang, G. R., & Wang, X. J. (2018). Inter-areal balanced amplification enhances signal propagation in a large-scale circuit model of the primate cortex. Neuron, 98(1), 222-234.\n5. From the structural connectivity map in Fig 2c, it appears as though there are few (if any) local connections imposed in their model. The data shown in the source material (e.g., Ref #25) do appear to have many local connections, however their figures also appear to be plotted at a finer spatial resolution than the model used in this manuscript. This detail is important, as the density of local connections could change across the same hierarchy observed in this study, and thus may provide an alternative explanation for the observed results.\nResponse: We apologize for the confusion created by Fig 2c. For visual clarity we displayed only inter‑areal projections in Fig 2c; all intra‑areal (local) synapses were deliberately omitted in the plot but are fully implemented in the model that drives every result in the paper. Specifically, each cortical area comprises excitatory and inhibitory populations connected by four local weights within each area: (see Page 24, Equation (1)-(2), in Methods, “4.3 Multi‑regional cortical model architecture” and Fig 2b).\nConsistent with the reviewer’s comment, the density of local connections indeed changes across the same hierarchy. Because our model adopts an areal resolution that aggregates neurons within each area measured in (Ref 25 Majka et al. 2020) as a single area node, the density of local connections is effectively captured by the local connection weights. Accordingly, the same composite excitatory gradient that defines the inter‑areal hierarchy is also applied to these local weights in the model. Thus, local recurrent strength increases systematically from sensory to association areas, which do indeed contribute to the observed temporal gradients. We have added this in the revised caption of Fig 2c for the clarification: “Only inter-areal (long-range) projections are shown for clarity.” We hope this clarification, together with the explicit reference to the math formulations, addresses the reviewer’s concern.\nReferences:\n1. Majka, P., Bai, S., Bakola, S., Bednarek, S., Chan, J. M., Jermakow, N., ... & Rosa, M. G. (2020). Open access resource for cellular-resolution analyses of corticocortical connectivity in the marmoset monkey. Nature communications, 11(1), 1133.\n6. The eigenmode/modelling analysis was informative and interesting.\nResponse: We thank the reviewer for the acknowledgment.\nReviewer #1 (Remarks on code availability):\nThe code appears to be well-documented.\nResponse: We thank the reviewer for the acknowledgment.\nReviewer #2 (Remarks to the Author):\nLi et al.'s manuscript \"A hierarchy of time constants and reliable signal propagation in marmoset cortex\" attempts to ground a multi-area cortex model in experimental ground truth to determine the connectivity between areas that can explain the hierarchy of timescales across the cortical processing hierarchy, and also enable reliable signal propagation. The model leverages marmoset EEG data that allows measurement of the timescale of cortical activity, and anatomical tracing data that includes the number of neurons projecting between the cortical areas. They found that when the network is operating in a state near criticality, both local integration (i.e. timescale hierarchies) and signal propagation are more optimal. Furthermore, the model mechanistically explains why the patterns of structural and functional connectivity are increasingly dissimilar across the cortical hierarchy.\nI think that the study is of great interest across systems and computational neuroscience, as it provides a mechanistic explanation for the prominent and systematic changes in activity dynamics (and integration) across the cortical processing hierarchy - and the anatomical basis for the model's connectivity is a major plus. I do have a few comments:\nMajor\n1. In the data used to fit the model, the animals were anesthetized with ketamine. I worry that this could impact activity timescales, as well as the activity states and signal propagation through the network. Perhaps they authors could check in the awake data used for signal propagation that the timescales are still similar in that state?\nResponse: We share the reviewer’s concern and therefore repeated the timescale analysis in the awake data used for signal propagation. We extracted the first 5 seconds of the session — prior to any stimulation—to obtain a resting-state segment free of task or laser effects. Using the same spectral method as for the anesthetized data, we obtain a clear hierarchy of intrinsic timescales across areas (Fig. R2a) that closely mirrors the anesthetized data. When areas are grouped by categories, timescales again increase from low-level sensory to high-level association areas (Fig. R2b). Directly comparing the time constants of neural activity in the awake state with those in the anesthetized state yields a strong correlation of (Fig. R2c), confirming that the relative ordering of timescales is approximately preserved.\nInterestingly, we also note that the absolute values of the timescales are slightly lower in the awake state. Two factors likely contribute: (i) wakefulness reduces slow, anesthesia-related oscillations, leading to faster decorrelation of neural activity (Sellers et al. 2015), and (ii) the shorter duration of the awake recording ( vs ~ ) likely underestimates the time constants. In summary, these findings suggest that the timescale hierarchy derived from the anaesthetized dataset is robust across behavioral states.\nWe have now updated the subsection “2.1 Experimental evidence of timescale hierarchy in the marmoset neocortex” in Results (pages 4-5, line 129-133) summarizing the analysis results from the awake resting state data. The sentences are as follows: “Furthermore, by analyzing the first 5 seconds of awake resting-state recordings [36], we observed a clear intrinsic timescale hierarchy across cortical areas closely matching the anesthetized state (Figure S9a-b), with a strong correlation confirming consistency in the relative ordering of timescales (Figure S9c).” All data processing and analysis details have been added into Methods (page 23, line 627-633) under the updated subsection “4.1 Experimental data sources,” and the newly obtained awake-state timescale results (Figure R2) are included as Figure S9 in the Supplementary materials (page 47). We hope these additions have addressed the reviewer’s query regarding timescale during the awake resting state.\nFig. R2. The timescale hierarchy in the marmoset cortex during the awake resting state. (a) Timescale distribution across the marmoset cortex, derived from neural activity recorded during the awake resting state. Lower-level sensory areas such as V1 and V2 exhibited shorter timescales, in contrast to higher-level areas like A46D and A10 that exhibited longer timescales. (b) Boxplot of timescale statistics for each areal category. Box plots show the media and IQR, with whiskers extending to the minimum and maximum values excluding outliers. (c) Strong correlation between timescales estimated from anesthetized and awake resting-state data ( ), confirming that the observed timescale hierarchy is consistent across brain states. −4 r = 0.77,p = 1.54×10\nReferences:\n1. Sellers, K. K., Bennett, D. V., Hutt, A., Williams, J. H., & Fröhlich, F. (2015). Awake vs. anesthetized: layer-specific sensory processing in visual cortex and functional connectivity between cortical areas. Journal of neurophysiology, 113(10), 3798-3815.\n2. Anatomical tracing data is used to determine FLN. Does this also somehow inform the relative interareal 'connection strength' onto E and I cell types, or were these chosen by some other means?\nResponse: We thank the reviewer for the important question. The fraction-of-labelled-neurons (FLN) metric we use -- derived from retrograde tracer studies (Majka et al., 2020; Theodoni et al., 2022) represent the fraction of all neurons in source area that project to target area . These tracers do not distinguish the neurotransmitter identity of the labelled neurons, and the published datasets therefore provide no direct information about whether the postsynaptic targets in the recipient area are excitatory (E) or inhibitory (I) cells.\nIn our model we therefore treat FLN as the density of long‑range excitatory axons and assume that these axons innervate both local E and I populations with identical topography but different synaptic efficacy -- hence we apply the same FLN matrix to both cell classes. The relative efficacy of inter‑areal input onto E versus I populations is controlled by two global scaling parameters, and . To be specific, for every pair of areas , the weight derived from FLN is multiplied by a global gain when the postsynaptic target is excitatory and by when the target is inhibitory. The two gain parameters were fine‑tuned so that the resulting large‑scale circuit operates near criticality while jointly reproducing (i) the observed hierarchy of intrinsic timescales during the resting state and (ii) the amplitude of signal propagation across areas evoked by transient stimulation. Once and were fixed by this procedure, they were held constant for all simulations reported in the manuscript.\nWe acknowledge that cell-type–resolved tract-tracing data would allow a more mechanistic treatment of inter-areal E/I specificity. As such datasets become available in the future, our framework can be updated by replacing the global scalars and with projection-specific ratios. We hope this clarifies how the FLN data were used and how relative E- and I-targeted connection strengths were assigned in the model.\nWe have now added a new paragraph in Discussion (page 22, line 573-585) explicitly clarifying the rationale behind setting relative inter-areal connection strengths onto E and I cell types. The new paragraph is as follows: “Another important consideration relates to the anatomical connectivity data used: the tracer-derived FLN weights currently employed in our large-scale marmoset cortical model do not specify neurotransmitter identity of target neurons (Majka et al., 2020; Theodoni et al., 2022). We therefore applied an identical FLN connectivity matrix to both excitatory (E) and inhibitory (I) neuronal populations within each cortical region, with their relative efficacies modulated by two global parameters, and . These global scaling parameters were carefully tuned to simultaneously satisfy crucial dynamical constraints—notably, the observed hierarchy of intrinsic timescales and robust inter-areal signal propagation—and were subsequently kept fixed throughout all simulations. While this simplification enables effective modeling of large-scale cortical dynamics, it neglects potential heterogeneity in cell-type specific projections. Future availability of cell-type resolved anatomical tracing data will refine the model by introducing projection-specific E/I targeting ratios instead of global scalar parameters.”\nReferences:\n1. Majka, P., Bai, S., Bakola, S., Bednarek, S., Chan, J. M., Jermakow, N., ... & Rosa, M. G. (2020). Open access resource for cellular-resolution analyses of corticocortical connectivity in the marmoset monkey. Nature communications, 11(1), 1133.\n2. Theodoni, P., Majka, P., Reser, D. H., Wójcik, D. K., Rosa, M. G., & Wang, X. J. (2022). Structural attributes and principles of the neocortical connectome in the marmoset monkey. Cerebral Cortex, 32(1), 15-28.\n3. I think that the authors could make it clearer in places like the results and the methods on what basis the additional composite gradient h is needed (on top of FLN). The increasing number of spines in higher cortex seems to be the reason used. Is the implication then that a given neuron that projects to an area in higher cortex may branch its axons more and form more synapses than a neuron projecting to a 'lower' area? Or is it that more spines = more local recurrence = amplified impact of an input from another area? Further justification and explanation would be helpful in justifying the use of the composite gradient.\nResponse: We thank the reviewer for providing us with the opportunity to clarify the rationale behind the composite gradient . We would like to point out that the FLN matrix only describes the fraction of labelled neurons projecting from source area to target area . To be more specific, FLN is a row-normalized measure, and the row sum of FLN satisfies: for every target area . Consequently, FLN does not contain information about the absolute number of long-range synapses each target area receives. The total long-range synapse number of a neuron in each target area can be reflected in the composite gradient, which highly correlates with spine count gradient based on anatomical studies below:\n(i) The spine counts on the dendrites of pyramidal neurons, rise monotonically from lower-order sensory cortex to higher-order associative cortex (Theodoni et al. 2022), implying that single neurons in higher-order areas physically possess more excitatory synapses. Thus, along the sensory-to-association hierarchy, the total synaptic input number increases.\n(ii) Markov et al. (2011; Fig. 2A) further show that, for every area examined, ~80 % of its labelled inputs are intrinsic (local) and ~20 % are extrinsic (from other areas). Because this proportion is essentially constant across the cortex, a larger total synaptic input number at higher-order areas directly implies a larger absolute number of both local recurrent and long-range synapses—even though their fraction remains unchanged.\nTaken together these findings indicate that neurons in higher-order areas receive a larger number of total excitatory synaptic inputs — both long-range and local inputs —than neurons in lower areas. The composite gradient therefore re-introduces this biologically motivated scaling by multiplying all excitatory weights by a factor proportional to the spine-count gradient. Importantly, the gradient does not assume extra axonal branching in a single projection, nor does it selectively boost local recurrence; rather, it captures the larger pool of postsynaptic spines available to receive the inputs in higher order areas, which amplifies the functional impact of any incoming drive.\nWe have now added a new paragraph in Methods (4.3 “Multi-regional cortical model architecture”, page 24-25, line 664-680) to clarify the motivation and rationale underlying the composite gradient. The new paragraph is as follows:\n“To complement our large-scale model that uses FLN to describe inter-areal connectivity, we additionally introduced a composite gradient of excitation, . Because the FLN matrix is normalized per target area (each row sums to 1), it encodes only relative inter-areal connectivity but does not capture the area-wise differences in the absolute number of synapse inputs received by each target area. Anatomical evidence indicates that neurons in higher-order association cortex receive more inputs per cell than those in early sensory areas. For instance, dendritic spine counts increase systematically from sensory to association areas (Theodoni et al., 2022). Furthermore, a uniform ~80%/20% division of local versus long-range inputs across areas (Markov et al., 2011) implies that higher-order areas inherently possess greater numbers of both local and inter-areal synapses. Thus, neurons in association cortex receive a larger overall drive than those in lower areas. To account for this gradient, we define a composite gradient of excitation as an area-specific scaling factor applied to all synaptic weights, effectively representing the increasing synaptic density in higher cortical areas while preserving the empirically measured FLN connectivity pattern. Within our model, is a factor scaling both local (intra-areal) and long-range (inter-areal) excitatory inputs and varies across different areas.”\nReferences:\n1. Theodoni, P., Majka, P., Reser, D. H., Wójcik, D. K., Rosa, M. G., & Wang, X. J. (2022). Structural attributes and principles of the neocortical connectome in the marmoset monkey. Cerebral Cortex, 32(1), 15-28.\n2. Markov, N. T., Misery, P., Falchier, A., Lamy, C., Vezoli, J., Quilodran, R., ... & Knoblauch, K. (2011). Weight consistency specifies regularities of macaque cortical networks. Cerebral cortex, 21(6), 1254-1272.\nMinor\n4. Some typos:\no line 2 \"information process\".\no Line137\"descried\"\no Fig6:\"excitability\"is misspelled in a few panels\nResponse: We thank the reviewer for pointing out these typographical errors. We have carefully reviewed and corrected all the mentioned typos.\n5. What about cortico-thalamo-cortical connections: can the authors comment on how these could play a role in dynamics and their timescales? Are these more or less captured in the 'mismatch' between correlations and structural connectivity?\nResponse: We agree with the reviewer that cortical-thalamo-cortical (CTC) loops are likely to modulate both propagation dynamics and intrinsic timescales and could partly account for part of the “structure–function mismatch” reported in our manuscript. A potential evidence has been shown in the large-scale mouse model of Ding et al. (2024), where including CTC loops helped sustain distributed persistent activity. Nonetheless, they also showed that several cortical dynamical signatures—such as whether an area displays persistent firing and the absolute firing-rate range—remain qualitatively similar when the thalamus is omitted. This finding suggests that the cortical connectomes captured in our present model sets the baseline temporal structure, whereas CTC loops may further stabilize activity and enhance inter-areal correlation.\nWe therefore predict that incorporating marmoset CTC loops would refine, rather than overturn, our main findings: (i) it would modulate the mismatch between functional correlations and corticocortical weights, and (ii) it would slightly lengthen the intrinsic timescales of higher-order areas by providing an additional reverberatory pathway. Because a high-resolution marmoset-wide dataset that links each cortical area to its thalamic partnersis not available yet, in this work we are only able to analyze the cerebral cortex-only network first and we explicitly flagged the limitations of not including CTC in the model.\nWe therefore have added a paragraph in the Discussion (page 22, line 586-597) noting the potential contribution of CTC loops as follows: “While our model focuses exclusively on cortico-cortical connectivity, the cerebral cortex is not anatomically isolated -- subcortical circuits (e.g., thalamus) and peripheral sensory pathways (retina, cochlea, spinal cord) impose important computational constraints on cortical activity. However, due to the lack of a marmoset-wide subcortical connectome at comparable resolution, our model necessarily focused on the corticocortical network, which can later be augmented with subcortical circuits as new data emerge. We anticipate that incorporating cortical-thalamo‑cortical (CTC) and peripheral pathways in future marmoset models will refine our present findings, particularly by enhancing model predictions related to intrinsic timescale patterns, signal propagation, and functional connectivity profiles. Incorporating such pathways represents a key direction for future work as compatible anatomical datasets become available.”\nReferences:\n1. Ding, X., Froudist-Walsh, S., Jaramillo, J., Jiang, J., & Wang, X. J. (2024). Cell type-specific connectome predicts distributed working memory activity in the mouse brain. elife, 13, e85442.\n6. It would be helpful to have more information about the structural data used in this study. Spine counts are used and correlated to the activity timescales - are these total spine counts on entire single neurons, or on specific dendrites? In which layer?\nResponse: We thank the reviewer’s comments, which provides us the opportunity to clarify the structural dataset. The spine values used in Fig. 2 and 3 correspond to the mean basal-dendrite spine counts of layer-3 pyramidal neurons in each cortical area. These measurements were collected from several studies and collated in Supplementary Table 4 of “Structural Attributes and Principles of the Neocortical Connectome in the Marmoset Monkey” (Theodoni et al. 2022). To make this provenance explicit, we have added the following statement in the revised caption of Fig. 2f: “Layer-3 pyramidal-cell spine counts were taken from the cross-area compilation in Supplementary Table 4 of [Theodoni et al. 2022]; sampling details are given in the cited source articles.” We hope these additions can give readers a clear picture of how the spine-count variable was defined and measured.\nReferences:\n1. Theodoni, P., Majka, P., Reser, D. H., Wójcik, D. K., Rosa, M. G., & Wang, X. J. (2022). Structural attributes and principles of the neocortical connectome in the marmoset monkey. Cerebral Cortex, 32(1), 15-28.\n7. Fig 6d, e: It is a little confusing for the reader to have the dissimilarity axes in these two panels go from low to high and high to low.\nResponse: We thank the reviewer for this suggestion. To avoid any confusion, we have re-oriented the dissimilarity axis in Fig. 6e so that both panels now display “low → high” from left to right. The figure and caption have been updated in the revised manuscript. The updated Figure 6e is as follows:\nNew Figure 6 (e): FC-SC dissimilarity across seven functional subnetworks of the marmoset brain. Left: the control model result. Right: The model result with gradient removed.\nReviewer #2 (Remarks on code availability):\nI looked at the code - I did not try running it\nResponse: We thank the reviewer for inspecting the repository.\nReviewer #1 (Remarks to the Author):\nThe authors response to my concerns was thorough and authoritative. I have no further reservations.\nResponse: We thank the reviewer for the thoughtful feedback and are pleased that our revisions addressed all concerns.\nReviewer #2 (Remarks to the Author):\nI appreciate the thorough responses and edits that the authors made, based on my previous comments. I have no further concerns.\nResponse: We thank the reviewer for the positive assessment and helpful suggestions throughout the review process.\nReviewer #2 (Remarks on code availability):\nI reviewed the code in the last round of review, but have not re-checked it.\nResponse: We appreciate the reviewer’s previous review of the code."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors response to my concerns was thorough and authoritative. I have no further reservations.", "reviewer_2": "I appreciate the thorough responses and edits that the authors made, based on my previous comments. I have no further concerns."}, "author_response": ""}]

Decision letter