Resting-state fMRI scans were obtained from the available public dataset of the University of North Carolina samples at Greensboro (Wahlheim et al., 2022; Wu et al., 2023b). The participants were 34 healthy young adults (18–32 years old, mean: 22.21, SD: 3.65). The functional MRI was acquired with an echo-planar imaging sequence: 32 slices with 4.0 mm thickness and no skip, time of echo = 30 ms; time of repetition (TR) = 2000 ms; flip angle = 70°, field of view = 220 mm, matrix size = 74 × 74 × 32 voxels. Each fMRI scan lasted for 10 min, comprising 300 volumes. Additional details about the raw fMRI data can be found in Wahlheim et al. (2022).

To generate a steady blood oxygenation level-dependent activity signal, the first five volumes of each scan were discarded to allow for magnetic stability. Similar to the previous studies (Wu et al., 2024b, 2023a), the functional data was then processed with the following steps: (1) Realignment to correct head motion for verification details; (2) Slice time correction. (3) Outlier identification (The volumes at the time point would be regarded as outliers and removed if the signal value is three standard deviations beyond the mean global signal for the entire run or if the head motion exceeded 0.5 mm in any direction). (4) Normalization (normalize to 3 mm MNI space using a template from the SPM software package; Friston, 2003). (5) Spatial smoothing with a Gaussian kernel of 8 mm full-width at half-maximum (FWHM).

A spatial group ICA was performed on the preprocessed and denoised BOLD signal using the Group ICA of FMRI Toolbox (GIFT) infomax algorithm. Specifically, high-model order ICA with a set of 100 components was obtained, comprising resting-state brain networks spanning cerebral cortical and subcortical regions. The configuration for the group ICA algorithm was developed according to the detailed description provided by (Wu et al., 2023b; Salman et al., 2019). In particular, a two-stage Principal Component Algorithm (PCA) method was first adapted to preserve the components that account for the most variance. The top 120 principal components (PCs) of all participants obtained in the first stage were concatenated across time and then further reduced to 100 in the second stage. Then, the infomax algorithm was used with 20 repeats to find steady independent components (ICs). After back reconstruction, the participant-specific spatial maps and corresponding time courses were obtained.

Three methods were employed to detect activated potential functional networks from the IC reservoir. (1) The spatial activation maps from the ICs were visually inspected to identify if they matched the large-scale functional network locations from previous studies and to make sure they were located at gray matter volumes. (2) The multiple regression method was used to produce the weight of ICs whose spatial pattern matches with the existing functional network template. The weights were used to rank the list, and the functional network that matched the most was selected. (3) The power spectrum of the ICs was checked to see if it follows a low-frequency peak and a high-frequency steady pattern (the time courses of ICs are characterized by high dynamic range). Those ICs that located cerebrospinal fluid and ICs whose highest regression weights were significantly low and whose power-spectrum curves were different (e.g., a clear peak of the wave at high frequency) were discarded.

After removing noise-related components, the time courses of retained components were triple detrended and despiked (Wu et al., 2024a). The motion parameters and global average were regressed for postprocessing, and a band-pass filter was applied (0.023-0.1 Hz). These actions ensure artifact noise can be largely eliminated and has minimal impact on the further signal analysis.

The first step of calculating the neuron dynamics is to stimulate and construct the neurons with varying levels of dendritic integrity. It resorts to initial neurons with rich morphological shapes, which we call prototype neurons, and this process is implemented by iteratively pruning their terminal dendritic compartments. Specifically, the prototype neurons were taken from the NeuroMorpho database (Ascoli et al., 2007). This database contains thousands of neurons with detailed dendritic structure, which can be used to simulate multicompartmental neurons. To avoid the bias provided by the prototype neurons' details (i.e., the species, human or animal), location (i.e., the brain regions), and topology (i.e., spatial information), 13 high-quality prototype neurons were randomly selected across 6 species and 5 brain regions. The morphology of these prototype neurons has been presented in Figure 2, and their details, including the number of their compartments and bifurcations, can be seen in Table 1.

Previous cytoarchitectonic findings suggested that along the spatial cortical gradients, the higher the anatomical level, the smaller the dendritic size of their pyramidal neuron tends to be (Beul and Hilgetag, 2019; Hilgetag et al., 2019). Hence, we utilized the proposed pruning algorithm with the prototype neurons to simulate neurons with variable dendritic sizes. Specifically, the pruning algorithm used an iterative approach to progressively remove the most terminal compartments (Kirch and Gollo, 2021). At each iteration, the terminus was detected and then removed until the entire dendritic tree was eliminated (see Figure 1C for an illustration of the pruning process). In our experiments, the maximum iteration step explored was 80. The neuronal morphologies obtained in intervals of 10 pruning iterations were explored. Consequently, for each prototype neuron, nine different dendrites were studied, resulting in a total of 117 different neuronal structures.

The neuronal dynamics was simulated and the activity at the soma was recorded for subsequent analyses. Specifically, we adapted a synchronous susceptible–infected (active)–refractory–susceptive model (SIRS) to simulate the dynamics of each active dendritic compartment, and particular attention was given to the somatic spiking activity (Gollo et al., 2009, 2012, 2013; Kirch and Gollo, 2020, 2021). Each excitable compartment can be in the following 3 states: susceptible, active (spike), and refectory. Dendritic compartments become active, generating dendritic spikes, either by external synaptic input or by the propagation of activity from neighboring dendritic compartments with a probability equals p_λ. As full spatial distributions of ion channels governing the dynamics of dendritic compartments remains largely unknown, here we considered the simplest case of p_λ = 0.97, which was homogeneous across the entire dendritic tree. In our experiments, the external synaptic stimuli were modeled as a stochastic process, and the active rate was a Poisson process with rate r = 0.03. Active compartments become refractory at the next time step (δt = 1ms). The refractory compartments become susceptible after a period of δr = 8ms. Unless activated externally via synapses or by propagation from an active neighbor, dendritic compartments remain in a susceptible state. For each neuron, the simulations were run for a duration of 100,000 time steps (1000s). The soma was modeled as a single-compartment that has the same properties of other dendritic compartments.

This model dynamics representing the activity of thousands of dendritic compartments was firstly proposed by Gollo et al. (2009) to show that neurons with large active dendrites, typical of somatosensory regions (Beul and Hilgetag, 2019), can significantly increase the dynamic range of neurons, allowing them to respond to a wide range of input levels. The model was further analyzed using the excitable-wave mean-field approximation (Gollo et al., 2012), and its generalized form (Gollo et al., 2013) to gain insights into the role homogeneous and non-stereotyped dendritic spikes caused by non-homogeneous distribution of ion channels. Additionally, this model (Gollo et al., 2009, 2012) was able to reproduce not only standard sigmoid response functions but also various types of complex response functions, such as double sigmoid curves obtained from the retinal ganglion cells of the mouse (Deans et al., 2002). More recently, this neuronal model has been further refined (Kirch and Gollo, 2020, 2021) to more accurately reflect real biological processes by using the real morphology of neurons obtained from digital reconstructions available in the NeuroMorpho database (Ascoli et al., 2007).

Out of the 100 ICs identified by the group ICA, 40 ICs were identified as noise components and discarded. The remaining 60 components were assigned as RSN (an example of a clean signal component taken, and another example of a noise signal component discarded are presented in Figure 3). The 60 RSN components were assigned to six functional domains that have been widely studied for brain networks: subcortical network (SC), auditory network (AUD), visual network (VIS), sensorimotor network (SM), cognitive control network (CC), and default mode network (DMN). The spatial map of 60 RSNs and corresponding peak coordinates have been provided in the Supplementary Table 1 and Supplementary Figure 1.

We utilized a one-way ANOVA followed by post-hoc two-sample t-tests to compare the intrinsic timescales across functional networks. The results (shown in Figure 4A) indicate significant group effects on the functional networks (F = 14.612, p < 0.001) for INT_f. Post-hoc t-test shows that the subcortical network exhibited the lowest intrinsic timescales compared to other high-level cortical networks (t_SC−DMN = 3.32, p < 0.001; t_SC−VIS = 2.86 p < 0.001; t_SC−SMN = 3.76, p < 0.001; t_SC−CC = 2.78 p < 0.05, t_SC−AUD = 2.32, p < 0.05, FDR corrected). The INT_a method indicates similar results (t_SC−DMN = 4.12, p < 0.001; t_SC−VIS = 4.86 p < 0.001; t_SC−SMN = 3.76, p < 0.001; t_SC−CC = 3.68 p < 0.05, t_SC−AUD = 2.12, p < 0.05, FDR corrected; Figure 4B), sustaining the findings that functional networks involved in high-level processes tend to have longer intrinsic timescales.

In addition, the correlation between the offset of exponential fitting and INT_f was calculated via regression analysis. It was found that INT_f negatively correlates with the autocorrelation offset for each functional network (Figure 4C, see the detailed results in Table 2). The offset measures the distance between the obtained timescales and the observation window, which can reflect the contribution of timescales to neural fluctuations (Murray et al., 2014).

With the natural hierarchical and spatial gap between subcortical and cortical functional networks, we computed the cumulative distribution of INT_f of multiple brain networks. The cumulative distribution of intrinsic timescales shows significant differences across the function networks (F = 13.37, p < 0.0001). While with some specific probabilities, for example, 0.5, the VIS has the highest intrinsic timescales, DMN is the second, and then followed by the SMN and AUD. According to spatial location (the subcortical network is located more inside the brain). DMN, SMN, VIS, CC, and AUD are combined together into cortical networks. Results (Figure 5) show longer INT for cortical networks and shorter INT for the subcortical network (two-sample t-test: t_{cortical−subcortical} = 6.09, 95% CI [0.50, 096] p < 0.0001). This finding indicates a distinct temporal organization between the subcortical and cortical networks, implying that neural dynamics and information processing vary depending on their specialized functions.

To obtain neurons with different levels of dendritic-tree integrity, a pruning algorithm was used to prototype neurons obtained from the NeuroMorpho database (Ascoli et al., 2007). Figure 6A presents an example of the morphology of a prototype neuron. As the pruning algorithm proceeds, the terminal dendritic compartments are progressively removed; as a consequence, the neuron integrity is reduced. The resulting morphology of a simulated neuron at iterations 20, 40, and 80 can be seen in Figure 6A respectively.

At each iteration, considering the dynamics of each compartment as an excitable unit (Kirch and Gollo, 2021), the neuronal dynamics were simulated. The tree-shaped neuron forms an excitable network, and the spiking activity of compartments determines the temporal evolution of activity and the corresponding somatic spikes. Figure 6B illustrates snapshots of neuronal activity at different time steps for neurons with different levels of dendritic integrity.

As the pruning algorithm proceeds, the neuronal integrity is progressively affected, and the number of compartments and connections to the soma is reduced (Figure 6C). This pruning process impacts the neuronal dynamics, and causes a substantial reduction in the firing rate in most neurons (Figure 6C). To gain a deeper insight into how firing rate influences intrinsic timescales, we simulated the dynamics of neurons in the absence of dendrites, that is, at the last pruning iteration when all the neurons are identical and only the neuronal soma remains. For these punctual neuron models, we vary the rate of input (modeled as a Poisson process) that controls the firing rate of neurons and measure the intrinsic timescales. The results show that the intrinsic timescales decay linearly with the firing rate of the neuron (Figure 6D).

We then examine the association between the INT and dendritic size. The statistical test (see Figure 7A) demonstrates that the prototype and pruned neuron have significantly different INT_f (F = 5.1639, p = 0.0001). Compared to the prototype neuron, the pruned neurons exhibit significant increased INT_f for iterations greater than 30 (post-hoc t-test, t₃₀ = 2.639 p = 0.014; t₄₀ = 3.443 p = 0.002; t₅₀ = 3, 243 p = 0.003; t₆₀ = 3.927 p = 0.0006; t₇₀ = 4.472 p = 0.004; t₈₀ = 4.776 p < 0.0001, FDR-corrected). The population-level results of INT_f calculation (the best fitting curve for all neurons at the same iteration) also reveal that as the pruning iteration progresses, INT_f increases. The prototype neuron exhibits the shortest INT_f of 0.9063, while the maximally pruned neuron displays the longest INT_f of 0.9277 (Figure 7B). The pruning algorithm reduces the size and complexity of the dendrites, which causes an alteration in the INT. Our results suggest an inverse relationship: neurons possessing larger dendrites exhibit shorter INT, while those with smaller dendrites display longer INT.

Cytoarchitectonic studies demonstrate that the dendritic size decreases for increasing hierarchical levels (Hilgetag et al., 2019). Therefore, the pruning iteration reveals a gradient embodying the various hierarchical levels. Considering the spatial position of the subcortical and cortical network, the neurons generated at larger pruning iterations may represent neurons situated in brain regions higher up in the hierarchy. In the absence of a precise mapping relationship, we explore different possible neuron partitions for subcortical to cortical functional networks based on the number of iterations (as seen in Figure 8A). Under this assumption, the cumulative distribution of neuron INT demonstrates a clear hierarchical order of subcortical to cortical functional networks (Figure 8B). The cumulative distribution of neuron INT with different numbers of pruning iterations (n = 10 and n = 50) shows very similar results (Supplementary Figure 2). Statistical tests also substantiate the lower INT of subcortical networks compared to cortical networks (n = 10: t_subcortical = 3.239 p = 0.012; n = 30: t_subcortical = 3.451 p = 0.002; n = 50: t_subcortical = 4.639 p = 0.001 FDR-corrected). The results indicate that regardless of the partition used, a significant distinction between INT of neurons located at subcortical and cortical regions is observed.