Overview (Fig. 1)Fig. 1

Experimental flow and data processing diagram. A Experimental protocol for collecting cerebral blood flow imaging from the rat hyperacute cerebral ischemia model induced by left carotid artery occlusion. B Rat atlas template and registered CBF images. C Arterial and venous ROIs and their opposing resting signals. D Cortical grid ROIs resting state signals are obtained to construct a Pearson correlation matrix. E Graph theoretical analysis is executed on the correlation matrix to reveal altered network properties. F Spatiotemporal dynamics are examined using the coactivation patterns analysis. G The evoked cortical activation is measured by general linear model analysis and correlated with network parameters. CA: carotid artery. WS: whisker stimulation

In this study, a total of 26 aging rats (6–8 months) were used to investigate cerebral ischemia. The resting state signal was initially collected under normal physiological conditions, followed by whisker stimulation using block design methodology. Subsequently, ischemia was induced by blocking and unblocking the left carotid artery. The left carotid artery occlusion was used to facilitate result comparisons across different ischemic rodent experiments (Blaschke et al. 2021; Brunner et al. 2018; Li et al. 2022). During ischemia, resting state and task signals were acquired under both occlusive and reperfused conditions using identical parameters. The subregional signals were derived from the predefined rat grid regions of interest (ROIs) atlas following post-preprocessing steps, including vein signal regression, which was subsequently employed to construct the correlation matrix. Network properties of global and local networks were assessed using the GRETNA toolbox (Wang et al. 2015), and the connection graph was constructed using the Pajek software (Mrvar and Batagelj 2016). Lastly, the brain activation response was evaluated via a general linear model and correlated with network features.

Animal preparation and LSCI data collection

Male Sprague–Dawley rats, aged 6–8 months, were procured from the Ensiweier Biotech Company (Chengdu, China) and were housed individually with ad libitum access to food and water and kept on a 12 h/12 h light/dark cycle. The experiment operation was conducted in strict accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals by the Ethics Committee of the University of Electronic Science and Technology of China. For the cranial thinning surgery, the rats were anesthetized via 1.5% gaseous isoflurane, and their skull surfaces were meticulously polished with a cranial drill, exposing the vessels and enabling the observation of cerebral blood flow signals. Their vital signs were closely monitored using an animal monitor (Paulette multi-parameter animal monitor M6vet, Shanghai, China) to ensure a stable physiological state during surgery. Chloral hydrate (6%, 400 mg/kg) was intraperitoneally injected to maintain an appropriate level of anesthesia during CBF image collection.

The CBF images were collected using a commercial LSCI apparatus (SIM BFI HR Pro, SIM Opto-Technology Co., Ltd, Wuhan, China), which is described in detail in a prior publication (P. Hu et al. 2022a, b). Figure 1A illustrates an intact data collection protocol for the rat model of hyperacute cerebral ischemia within the left brain, which contained resting and task state signals. Here, we selected occlusion of the left carotid artery in rats, which is most commonly employed in rodent models (Blaschke et al. 2021; Brunner et al. 2018; Li et al. 2022), to simulate the pathological processes of human cerebral ischemia due to its higher prevalence in humans (Hasan et al. 2021). The resting state signal consisted of 10-min physiological resting, 30-min left carotid artery block ischemia, and 30-min reperfusion states. Task state signals were acquired after each resting period, which consisted of three repetitions of blowing whisker stimulation. Each trial encompassed a baseline (10-s) and four repetitions of whisker stimuli (5-s stimuli and 15-s interval).

Image preprocessing

To register CBF images in standard space, we initially used White's method (White et al. 2011) to create a standard rat brain template shown in Fig. 1B. Briefly, a horizontal spatial template was generated based on the spatial coordinates of each cortex region on the coronal plane. The polygon edges of brain regions were smoothed to prevent overlapping spatial positions of different regions. The template was marked with the bregma, lambda suture, and two lateral edge positions for registration between CBF images and standard space (Paxinos et al., 2004), resulting in different subregions in the rat brain. We further selected 18 corresponding ROIs within these subregions while also considering their distinct characteristics such as large vein diameter and slender multi-branch artery structure (Armitage et al. 2010). The venous signal displayed distinctive signal fluctuations compared to the arterial signal (Fig. 1C). Specifically, a total of 18 venous ROIs were selected from the central M2 and V2M, middle sensorimotor, and marginal V2L, V1, and auditory subregions, each with a diameter of 0.24 mm (28 pixels).

CBF data during both normal physiological conditions and ischemia were preprocessed using the same steps, including registration, noise reduction, dataset binning, linear detrending, normalization, spatial smoothing, venous signal regression, bandpass filtering, and downsampling. Before registration, a brain mask was created for each animal to exclude the environmental background (outside the brain). The CBF image was then registered with a normalized horizontal template using affine transformation at four positions. To increase the signal-to-noise ratio, we performed two smoothing steps. Firstly, smoothing was performed on RGB images. The Gaussian smooth kernel (5 × 5 pixels box with a 1.3-pixel standard deviation) was used to suppress machine noise from the LSCI image (White et al. 2011). Failure to do so may introduce the Moire pattern into subsequent images. Secondly, smoothing was applied to the gray CBF image with Gaussian smooth kernel of a full-width-at-half-maximum of 5 pixels (0.4 mm) (Bergonzi et al. 2015). Time-series denoising was executed using the nearest principle to replace outlier values with non-outlier values, avoiding missing data. Linear detrending was performed to eliminate the linear drift of data, and then Z-score normalization was applied to improve data standardization. Given the inverse effect of venous activity on the CBF dataset, nine ipsilateral venous ROI signals (Fig. 1C) were extracted and averaged to generate a hemispheric venous signal, which was regressed out from the continuous time-series. The global cerebral time trace was not regressed out due to prior investigation (Murphy and Fox 2017). The CBF signal was further bandpass filtered between 0.009–0.08 Hz (White et al. 2011), restraining the ultralow-frequency drift and high-frequency physiological noise. Finally, the signal was downsampled from 10 to 1 Hz for minimal calculations.

Network construction using functional connectivity

We selected ROIs from the rat atlas template and measured their network topological characteristics. A total of 330 ROIs were evenly spaced with 0.8 mm gap (10 pixels) on the bilateral hemispheres, which were attributed to associations of the frontal (10 ROIs), motor (66 ROIs), sensory (92 ROIs), auditory (16 ROIs), visual (94 ROIs), parietal (14 ROIs), and central regions (38 ROIs). The averaged time-series for each ROI was obtained from preprocessed CBF data by averaging the time-series within a 0.4 mm diameter circle (5 pixels) in the subregional cortex. This grid atlas facilitates the time-series extraction for each subregion, resulting in a functional connectivity matrix across different subregions. Pearson's correlation coefficients were calculated between all node pairs, yielding a 330 × 330 correlation matrix for each animal. In addition, the homologous connections were calculated between the cross-hemisphere homotopic ROIs. Fisher’s r-to-z transformation was further performed for normalization.

Network topological property assessment

Graph theoretical analysis provides insight into the topological properties of brain networks regardless of anatomical tract or functional association network. Here, we employed this technique to investigate the functional connectivity network features of resting CBF activity, which was calculated with the Graph Theoretical Network Analysis (GRETNA) (http://www.nitrc.org/projects/gretna) (Wang et al. 2015) and Pajek software (version 5.18, http://mrvar.fdv.uni-lj.si/pajek/) (Mrvar and Batagelj 2016).

Graph theoretical measures were constructed using both weighted and undirected methods, including ROI nodes and functional connectivity edges for all node pair connections of the raw measures. Topological network analysis calculated global and regional network characteristics, with the sparsity threshold set at 0.10 to 0.30 and an interval size of 0.01.

The area under the curve (AUC) (Han et al. 2020; He et al. 2009; Liang et al. 2018) served as a summary for each network topological metric reflecting the underlying brain state. The AUC of network metrics was defined as the mathematical integral of the area under a curve (Zhang et al. 2011), which was computed as.

$$^= _^\left[Y\left(_\right)+Y\left(_\right)\right]\times \Delta S/2$$

(1)

where Y denotes network metrics, \(_\) means sparsity threshold, and \(\Delta S\) represents sparsity interval. Specifically, we calculated six global network metrics and four regional network metrics.

Global metrics included characteristic path length (L), global efficiency (Eg), local efficiency (El), small-world sigma (\(}\)), clustering coefficient (C), and assortativity coefficient (r).

The characteristic path length (L) is defined as

$$}=\boldsymbol\frac}(}-1)}\sum_},}\in },}\ne }}}}_},}}$$

(2)

where \(}}_},}}\) indicates the shortest path length between nodes i and j, N indicates the set of all nodes in the network, and L presents the shortest nodal connection path across the entire network. Thus, the path length (L) measures the average distance between all network nodes, with a small L indicating a well-connected, integrated network and a large L indicating a fragmented, less integrated network.

Global efficiency (Eg), sometimes referred to as the inverse of path length, is given by the equation

$$}}=\boldsymbol\frac}(}-1)}\sum_},}\in },}\ne }}\frac}}_},}}}$$

(3)

Eg measures the overall efficiency of information transfer throughout the network, which could reflect the network’s capacity for connectivity and integration. A little Eg coefficient indicates low global network communication, while a high Eg indicates an effective one.

Local efficiency (El) is defined as

$$}}=\boldsymbol\frac}}\sum_}\in }}}}\left(}}_}}\right)$$

(4)

where \(}}_}}\) signifies a set of nodes directly linked with node i. This metric measures the efficiency of a node’s interactions with its interconnected nodes within the network, focusing on the node’s immediate neighbors and thus representing the nodal ability for information integration based on local network connectivity.

The clustering coefficient (C) is defined as

$$}= \frac}}\sum_}\in }}}}_}}=\frac}}}_}}(}}_}}-1)}\sum_}\in }}\sum_},}\in },}\ne }}\sqrt[3]}}}_},}}\bullet }}}_},}}\bullet }}}_},}}}$$

(5)

where \(}}_}}\) is the clustering coefficient of node i and \(}}}_},}}\) is the normalized weight adjusted for the largest value in the network. \(}}}_},}}\) is set as

$$}}}_},}}=}}_},}}/}}}(}}_},}})$$

(6)

It is a parameter indicating the node clustering capability, which reflects the potential for network nodes to cluster together, providing insight into the local information integration of the connected network.

Small-worldness \(}\),

$$}=\mathbf/\mathbf$$

(7)

presents the comprehensive capacity of network modularity levels and global integration. The small-world network provides necessary insight for understanding effective information communication and integration throughout the entire network. This represents the capacity of a well-connected network where the majority of nodes can be reached from each other node by a small number of steps. A small-worldness measure greater than one suggests a robust small-world feature.

The assortativity coefficient, \(}\), is given by the equation

$$}=\boldsymbol\frac}}^\sum_},})\in }}}}_}}}}_}}-}}^\sum_},})\in }}\boldsymbol\frac(}}_}}+\boldsymbol}}_}})\right]}^}}}^\sum_},})\in }}\boldsymbol\frac(}}_}}^\boldsymbol+\boldsymbol}}_}}^)\boldsymbol-}}^\sum_},})\in }}\boldsymbol\frac(}}_}}+\boldsymbol}}_}})\right]}^}$$

(8)

$$}}_}}=\boldsymbol\sum_}\in }}}}_},}}$$

(9)

where \(}}_}}\) is the degree of node I and \(}}_},}}\) is the connection status of node pairs i and j. If there is a link between them, then \(}}_},}}=1\), otherwise \(}}_},}}=0\). The assortativity coefficient assesses the degree of linkage between nodes possessing a similar number of degrees, reflecting the tendency of nodes to preferentially connect to other nodes with a similar number of connections. This can further facilitate our understanding of network resilience.

The regional metrics include degree centrality (DC), nodal efficiency, shortest path length, and clustering coefficient. The nodal degree centrality DC,

$$}}}_}}= \frac}}_}}}}-1}$$

(10)

reflects the abundance of nodal connections within the entire network. It provides a simple yet powerful way to understand the importance of nodes in a network. The higher the nodal degree, the more centrally important the node’s network. The other three nodal indices have a similar significance to global network properties. More information on these network metrics can be found elsewhere (Rubinov and Sporns 2010; Wang et al. 2015).

Furthermore, each animal model’s correlation matrix was plotted (Wang et al. 2023, 2012) using sparsity thresholding. The functional connectivity between a pair of ROIs was considered reliably correlated when it was statistically significant (z > 0.4, p < 0.01) (Cha et al. 2016) across all node pairs with the use of false discovery rate correction. This resulted in undirected links between different subregional ROIs, which were defined as effective edges within the connection networks. The topological graph was then created using the Kamada-Kawai algorithm, which positioned strongly connected ROIs closer to each other and placed weakly connected ROIs farther apart. The node with the most links was located in the center of the network, while the node with the least links resided at the network’s boundary. This intuitive topological graph provided a comprehensive display of the functional subregional networks. We also calculated the degree centrality for each ROI to measure their connectivity importance within the functional network. ROIs that had higher degrees, within the top 10%, were considered network hubs.

Coactivation patterns state transition analysis

Coactivation pattern (CAP) analysis can identify recurring coactivation patterns in brain functional networks, thereby revealing the spatiotemporal dynamics of CBF functional networks. Previous studies indicated that the activated and deactivated network states defined in the CAP analysis strongly align with traditional brain functional networks (Liu and Duyn 2013; Yang et al. 2021). As a data-driven analysis method, it can identify activated and deactivated configurations of network regions by analyzing global CBF signal time series.

CAP analysis has been described in detail in our previous studies (Yang et al. 2021, 2022). The time series of all ROI signals are first standardized by Z-scores to represent their relative strength. Next, all rat baseline time series were merged together to generate a group-level dataset. K-means clustering was then used to identify different CAP states in the baseline group, assigning each CBF time point to the corresponding coactivation pattern. Group-level CAP states were calculated by averaging identical group CAP clusters and dividing by the within-cluster standard deviations to yield standard Z-value CAP maps. Subsequently, we used the Pearson correlation to calculate the spatial similarity between the ischemic CBF signal and the CAP states, and assigned the ischemic CBF signal to the CAP states with the highest similarity. In our current study, we examined different cluster numbers K, ranging from 2 to 21, with a step size of 1. The clustering results under different K were evaluated using the silhouette score, and the optimal K number was determined using the elbow method. In this study, K = 6 was selected in the subsequent analysis.

Moreover, we measured temporal CAP properties for the individual rat using six metrics: fraction of time, persistence, counters, resilience, in-degree, and out-degree. The state fraction of time represents the proportion of each CAP state's cumulative time to the total time, indicating the importance of each CAP state. State persistence denotes the average dwell time of each CAP state before and after the state transition, reflecting the stability of each CAP state. State resilience indicates the possibility of maintaining the same state between two consecutive time points. State counts represent the numbers of each state over the entire length of the data. In-degree refers to the number of transitions from one state to others, and out-degree refers to the number of transitions from others to a specific state.

Evoked somatosensory functional activity assessment

To assess the cerebral evoked hyperemia response in different states, we examined the activity of the barrel cortex by stimulating non-ischemic-sided whiskers with air puffs and acquired the corresponding CBF images. The general linear model (GLM) method was employed to assess the evoked hyperemia intensity via whisker stimulation. Mathematically, this model is represented as:

$$}=}}+} + \epsilon$$

(11)

where \(}\) represents the actual pixel signal time series, \(}\) signifies the averaged barrel activation curve, \(}\) represents the intercept vector, and ε represents random noise. Subsequently, the estimated \(}\) and \(}\) intensity of the barrel cortex was extracted for further comparison and correlation analysis with the network topological metrics.

Statistical analysis

Statistical analysis was conducted utilizing MATLAB (https://ww2.mathworks.cn/) and GraphPad Prism 8 (https://www.graphpad.com/), with significance levels set at *p < 0.05, **p < 0.01, ***p < 0.001. A one-way ANOVA was implemented when comparing multiple groups. Additionally, a two-tailed paired t-test (n = 26) was utilized to compare differences between the groups. False discovery rate corrections were applied to prevent false positive results when analyzing graph theory and GLM analysis parameters.