Introduction

Epilepsy is a neurological disorder affecting about 70 million people worldwide. Idiopathic generalized epilepsy (IGE), a common subtype, accounts for approximately 20% of epilepsy cases, with a higher prevalence in newly diagnosed children and adolescents1,2. Early onset is often associated with a poorer prognosis, significantly impairing normal development and cognitive function in affected children and adolescents, placing a considerable burden on families3. While timely diagnosis and treatment are crucial, the etiology of IGE involves complex genetic and environmental interactions. Neurological and imaging examinations often yield normal findings, making diagnosis and treatment more challenging4,5. IGE is widely regarded as a brain network disorder. During propagation, abnormal discharges of epileptic neurons spread through the intricate neuronal connections, disrupting neurobiological networks across the brain6. Research increasingly highlights the role of interconnected brain networks in epilepsy pathogenesis and outcomes. Identifying novel neural network biomarkers linked to IGE pathogenesis may enhance its diagnosis and management.

Neuroimaging has become a pivotal tool in studying brain networks in IGE, offering valuable insights into structural and functional mechanisms. Resting-state functional magnetic resonance imaging (rs-fMRI) is widely used to detect abnormal functional connectivity and network changes in IGE by monitoring spontaneous brain activity. Several studies have linked abnormal brain network alterations to cognitive impairments in children with epilepsy7,8,9,10, advancing understanding of IGE’s complex pathological mechanisms. However, the brain operates as a dynamic system, with neuronal activity continuously evolving, even at rest11. Traditional analysis methods, which assume static brain activity during scanning, fail to capture dynamic changes and time-varying functional signal patterns among brain regions, potentially overlooking critical neural interactions. Dynamic network analysis offers a solution by assessing temporal connectivity reconfiguration, providing insights into neural activity12. These reorganizations are essential for efficient intermodular communication, flexible cognition, and rapid responses to external stimuli13,14. Recently, multilayer network analysis, a novel dynamic network model, was introduced in neuroscience and captures fluctuation in functional connectivity over time, revealing abnormal dynamic network alterations in patients with neurological disorders15,16. Although previous studies have linked abnormal dynamic temporal variations in primary and high-level networks to IGE phenotypes17,18, the dynamic brain network topology of pediatric IGE and its relationship with cognitive function remains insufficiently understood. Exploring dynamic changes in brain network topology, including module transformations and neurocognition-related patterns, can enhance understanding of the neural pathological mechanisms in children with IGE.

IGE is a highly heritable disorder. A genome-wide association meta-analysis by the International League Against Epilepsy (ILAE) involving over 29,000 patients identified 19 genetic risk loci and subtype-specific genetic architectures, reinforcing the genetic basis of IGE19. Brain gene expression profiles further link brain activity to gene variations20. Imaging transcriptomics connects micro-level genetic variations with macroscopic brain functional networks, using neuroimaging as an intermediate measure. This approach provides insights into neuropsychiatric disorders21,22,23 by elucidating correlations between genetics and neuroimaging. While traditional IGE imaging studies focus on morphology6, another key feature of IGE is its high temporal variability in brain networks caused by abnormal neuronal discharges24. Functional networks, particularly connectome dynamics, provide insights into the underlying mechanism of IGE; associations between connectome dynamics and spatial gene expression patterns remain unclear and require further investigation.

This study aims to (1) investigate changes in the dynamic topology of functional networks in children with IGE, their correlation with neurocognitive function, and associated transcriptional patterns, and (2) explore the potential of dynamic networks as predictive markers for cognitive impairment.

Methods

Participants

Ethical approval for this study was obtained from the Ethics Committee of the Affiliated Hospital of Zunyi Medical University. All procedures adhered to the principles of the Declaration of Helsinki. Written informed consent from the guardians of each participant was obtained.

Participants with IGE were recruited from the Pediatric Neurology Clinic of the Affiliated Hospital of Zunyi Medical University. Diagnoses were confirmed by two experienced pediatric neurologists using the International League Against Epilepsy (ILAE) criteria. Inclusion criteria for the IGE group were: (1) confirmed clinical diagnosis of IGE; (2) normal routine head MRI findings; (3) no history of substance abuse; and (4) age between 6 and 16 years at the time of MRI. Exclusion criteria were: (1) other neurological disorders, malignant tumors, or head trauma; (2) image artifacts, such as signal loss or distortion, that affected analysis; and (3) head motion exceeding 2 mm displacement or 2° rotation during rs-fMRI preprocessing.

Healthy controls (HCs) were aged 6–16 years, with no neurological, psychiatric, or systemic illnesses and no contraindications to MRI. Exclusion criteria for the HC group were: (1) any abnormal findings on conventional brain MRI and (2) head motion exceeding 2 mm displacement or 2° rotation during rs-fMRI preprocessing.

Neuropsychological assessment

Standardized neurocognitive assessment was conducted before MRI scanning in a quiet room, with only the participant and trained evaluator present. Performance IQ (PIQ), Verbal IQ (VIQ), and Full-Scale IQ (FIQ) were assessed using the Chinese version of the Wechsler Intelligence Scale for Children, a widely used standardized tool in China. Evaluations were conducted by a trained neuropsychologist.

Image acquisition

Following the neuropsychological assessment, participants were given a rest period before undergoing an MRI. Scanning was performed using a 3.0 T magnetic resonance imaging scanner (GE Healthcare, Chicago, IL, USA). Structural T1-weighted images were acquired using a three-dimensional brain volume sequence (3D BRAVO) with a repetition time (TR) of 1900 ms, echo time (TE) of 2.1 ms, inversion time (TI) of 900 ms, flip angle of 9°, slice thickness of 1.0 mm, and a matrix of 256 × 256. Functional images were obtained using a gradient echo planar imaging sequence with TR 2000 ms, TE 30 ms, voxel size 3.75 mm × 3.75 mm × 4 mm, a field of view (FOV) 240 mm × 240 mm, 33 slices, slice thickness 4 mm, and flip angle 90°. Prior to scanning, participants were instructed to relax, remain still, close their eyes, and stay awake.

rs-fMRI data processing

Preprocessing began by discarding the initial 10 s of data to account for signal stabilization, followed by slice timing correction, realignment, and spatial normalization to a standard template. Spatial smoothing was performed using a Gaussian kernel with a full width at half maximum (FWHM) of 6 mm. Nuisance regression was conducted using Friston’s 24 head-motion parameters and signals from cerebrospinal fluid, white matter, and global signals to minimize motion-related and physiological noise. Finally, temporal filtering was applied with a bandpass filter between 0.01 and 0.1 Hz.

Network construction

Dynamic functional connectivity based on the sliding window method

To construct a multilayer temporal network model, we employed time-based layer segmentation using the sliding window methodology12. Following established principles, the minimum window duration was set as 1/fmin, with temporal progression occurring at single time-unit intervals25. A window length of 100 s (50 TRs) and a progression step of 2 s (1 TR) were employed, resulting in 151 distinct temporal segments. Within each temporal window, Pearson’s correlation coefficients were calculated for the 100 brain regions defined by the Schaefer 100 parcellation atlas26. This process generated individual-specific dynamic network matrices (N × N × W), where N represents the 100 atlas-defined regions and W denotes the 151 temporal segments. To enhance network reliability and exclude weak connections, we applied a 20% density threshold to the resulting networks, thus generating participant-specific dynamic functional networks from rs-fMRI data.

Tracking dynamic modular structures

Dynamic features of connectome organization were investigated using a multilayer network model, which maintains temporal coherence by incorporating connectivity information from adjacent time windows. We employed a multilayer-variant Louvain algorithm (http://netwiki.amath.unc.edu/GenLouvain) to optimize the modularity index Q (0–1), which quantifies network module segregation. Modularity parameters (ω and γ) were set to their default value of 1. The temporal dynamics of module membership were assessed by calculating modular variability (MV) for each brain region, where higher values reflected greater flexibility in module switching (details are provided in the Supplementary Material). To account for algorithmic variability, analyses were repeated 50 times for each participant. Statistical comparisons of MV between the IGE and HC groups were conducted using two-sample t-tests at three levels: global network efficiency, network-level analysis based on Yeo 7-network parcellation27, and regional analysis of all 100 brain areas. FDR correction was applied to address multiple comparisons (Fig. 1A).

Fig. 1
figure 1

Schematic diagram of the study design. (A) Multilayer network analysis. The Schaefer 100 Atlas was used to extract regional mean signals from preprocessed rs-fMRI data, constructing dynamic functional connectivity matrices for each subject. Temporal network reorganization was analyzed using a sliding window approach with an iterative ordinal Louvain algorithm. Network switching rate were computed, compared across groups, and correlated with clinical and neuropsychological measures. (B) Machine learning. A relevance vector regression (RVR) prediction model was established based on dynamic network property variations as features to predict cognitive function in children with IGE. (C) Gene expression analysis. Gene expression values for each region were extracted from the Allen Human Brain Atlas database, forming a gene expression matrix. (D) Partial least squares regression. Partial least squares regression linked changes in modular variability in IGE children to gene expression data, followed by Gene annotation on significant genes from the third component of PLS (PLS3). (E) Enrichment Analysis. Enrichment analysis of gene ontology, biological processes, and cell types was performed for genes showing significant correlated with MV variations. PLS partial least squares, IGE idiopathic generalized epilepsy, rs-fMRI resting state function magnetic resonance imaging, RVR relevance vector regression. Brain maps in this figure were generated using the pysurfer package in python (PySurfer v0.10.0 https://pysurfer.github.io/) and BrainNet Viewer (www.nitrc.org/projects/bnv/).

Full size image

Relationships between clinical measures and alterations in modular variability

To investigate the association between MV changes and IQ scores, disease duration, and age of onset in children with IGE, Pearson’s correlation analysis was performed. The analysis examined the mean MV of the brain regions and networks that showed significant differences between the IGE and control groups, correlating these values with scores from the clinical scales. Statistical significance was set at P < 0.05, FDR correction was applied to address multiple comparisons (Fig. 1A).

Machine learning predictive model for cognitive function based on modular variability

The Relevance Vector Regression (RVR) algorithm was employed to assess the predictive capacity of MV for cognitive function scores in children with IGE. The MV of each participant served as a predictive feature, with prediction accuracy estimated using leave-one-out cross-validation, calculated as Pearson’s correlation coefficient between observed and predicted IQ. Nonparametric permutation tests were conducted to assess the model’s predictive accuracy. During this process, labels of participants’ cognitive functions were randomly shuffled prior to applying RVR and cross-validation. This randomization was repeated 5000 times, and accuracy rates were computed for each shuffled dataset. The frequency at which the accuracy rates of the randomized data exceeded the correlation coefficients obtained from the original unshuffled dataset was determined. If fewer than 250 of the 5000 randomizations produced higher correlation coefficients than the original dataset, the model’s predictive accuracy was considered statistically significant (correspond to a P value of < 0.05). This approach ensured that the observed results were unlikely to occur by chance, providing robust evidence for the model’s validity. Prediction analysis was implemented using libsvm (www.csie.ntu.edu.tw/~cjlin/libsvm/) and modified codes from Cui and Gong (https://github.com/ZaixuCui/Pattern_Regression)28 (Fig. 1B).

Association between alterations in connectome dynamics and gene expression profiles

Gene expression data preprocessing

The AHBA database (http://human.brain-map.org) provided gene expression data from six post-mortem brains. We utilized the abagen toolbox (https://www.github.com/netneurolab/abagen)29 to analyze and spatially map transcriptomic data across 100 segmented brain regions. Gene expression data preprocessing involved several steps: (1) updating probe-to-gene annotations; (2) applying an intensity-based filter; (3) selecting probes; (4) matching samples to regions; (5) addressing missing data; (6) normalizing samples; (7) normalizing genes; (8) calculating a sample-to-region combination metric; and (9) selecting stable genes. After excluding two brain regions lacking gene expression following preprocessing, a matrix of 98 regions and 12,791 genes was used for subsequent analyses (Fig. 1C).

Spatial correlation with gene expression profiles

To investigate the relationship between gene expression and brain modularity dynamics, we conducted a modularity dynamics-transcriptome association analysis using a group difference map (t-values, Fig. S1). Partial least squares (PLS)30 regression was employed to identify weighted linear combinations (components) of gene expression patterns from 12,791 genes associated with IGE-related changes in modularity dynamics. The PLS components were ranked based on the variance they explained in the relationship between gene expression and modularity dynamics. The third component (PLS3) accounted for the largest proportion of variance, capturing key patterns linking gene expression to modularity dynamics (Fig. 1D). To evaluate the significance of PLS3, we performed 10,000 spatial autocorrelation-preserved permutation tests31. Bootstrapping was used to assess the contributions of individual genes by generating PLS weights through 10,000 resampling iterations with replacement across 100 cortical regions. Standard errors from bootstrapping results were used to calculate Z-scores (weights divided by standard errors), quantifying the significance of each gene’s contribution.

Enrichment and cell-type analysis

Enrichment analysis was conducted using Metascape (https://metascape.org/gp/index.html) on significantly expressed genes identified from PLS regression (P < 0.05, FDR), categorized into up-regulated (PLS+) and down-regulated (PLS−) groups32. Biological functions and pathways were identified using the Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases33,34 (www.kegg.jp/kegg/kegg1.html), with the human genome serving as the background set. Statistical significance was determined using a cumulative hypergeometric distribution and adjusted with the Benjamini–Hochberg correction (P < 0.05), and Chiplot (www.chiplot.online) was used to draw the bubble diagram of the pathway. Cell type-specific analysis utilized data from five single-cell studies, with genes classified into seven canonical cell types, as described by Seidlitz et al.35. The significance of gene overlap between PLS+/PLS− genes and cell types was evaluated through permutation testing36, followed by pathway enrichment analysis for each cell type (P < 0.05, FDR) (Fig. 1E).

Validation analyses

To assess the validity of our data, we evaluated the potential influence of age, gender, head movement parameter (frame displacement, FD), and processing parameters, focusing particularly on the sliding window parameters (window size and step length) and multilayer network model parameters (ω and γ). Validation was performed with different window lengths (40 TRs, 70 TRs) while maintaining a step length of 1 TR. In addition to the primary study settings of ω = γ = 1, we repeated the analysis with alternate parameter combinations, such as ω = 0.7/0.8/0.9, γ = 1, and ω = 1, γ = 0.9/0.8/0.7, to explore their effects on the results35. We further assessed the influence of age, gender, and head movement on the primary outcomes. Age and gender were included as covariates in validation analyses to control for potential confounding effects. The head movement parameter was also incorporated as a covariate to adequately assess its influence on primary outcomes. Collectively, these steps ensured the authenticity and robustness of our findings (Supplementary Methods).

Results

Demographic characteristics

No significant differences in age, gender, or years of education were observed between the IGE and HC groups. Table 1 summarizes the demographic and clinical characteristics of the participants.

Table 1 Demographic characteristics of the sample.
Full size table

Switching rates vary at the global, subnet, and node levels

At the nodal level, patients with IGE exhibited significant MV alterations across brain regions compared to HCs. Increased MV was predominantly observed in brain regions of the default mode network (DMN), including the bilateral prefrontal cortex, precuneus, posterior cingulate cortex, and right superior temporal gyrus, while decreased MV was identified in the brain regions of somatomotor, attention, frontoparietal control, and DMN networks, particularly in bilateral parietal and prefrontal areas (Fig. 2A, Table S1). The network-level analysis confirmed these findings, revealing significantly increased MV in the DMN of children with IGE (Fig. 2B, Table S2), with no differences in other subnetworks.

Fig. 2
figure 2

Case–control differences in connectome dynamics. (A) Mean regional module variability values for IGE and HC group. Higher module variability primarily in the bilateral prefrontal cortex, precuneus/posterior cingulate cortex, and lower module variability primarily in the bilateral parietal regions and prefrontal areas. Case–control differences in module dynamics at the nodal level and the last panel shows regions with significant case–control differences in modular variability (FDR-corrected P < 0.05). (B) Network-level module variability comparisons. The radar chart shows the gradient score (with respect to healthy controls) of the two groups. The DMN was increased (FDR-corrected P < 0.05). (C) The MV changes in the left lateral prefrontal cortex was negatively correlated with VIQ and FIQ, MV changes in the left posterior cingulate cortex were negatively correlated with PIQ. In contrast, MV changes in the right orbital frontal cortex showed a significant positive correlation with VIQ and FIQ. Brain maps in this figure were generated using the pysurfer package in python (PySurfer v0.10.0 https://pysurfer.github.io/). IGE idiopathic generalized epilepsy, HC healthy control, VIS visual network, SMN sensorimotor network, DAN dorsal attention network, VAN ventral attention network, SUB subcortical regions, LIM limbic network, FPN frontoparietal network, DMN default mode network.

Full size image

Relationships between neurocognitive function and modular variability

We observed a significant correlation between MV changes and clinical features, including IQ scores, disease duration, and age of onset in children with IGE. MV changes in the left lateral prefrontal cortex were negatively correlated with VIQ (Pearson’s r = − 0.400, P = 0.042, q = 0.049) and FIQ (Pearson’s r = − 0.419, P = 0.032, q = 0.049), while changes in the left posterior cingulate cortex were negatively correlated with PIQ (Pearson’s r = − 0.408, P = 0.038, q = 0.049). In contrast, MV changes in the right orbitofrontal cortex showed a positive correlation with VIQ (Pearson’s r = 0.488, P = 0.011, q = 0.038) and FIQ (Pearson’s r = 0.442, P = 0.023, q = 0.049) (Fig. 2C).

Machine learning model performance on predicting cognitive function based on modular variability

MV was identified as a significant predictor of VIQ (r = 0.433, P = 0.027) (Fig. 3A). This predictive relationship was confirmed using a nonparametric permutation test (P = 0.041; Fig. 3B). The analysis highlighted that the most influential predictive features were within the DMN and frontoparietal control network (Fig. 3C). While we examined predictions for both PIQ (P = 0.06) and FIQ (P = 0.57) using identical methodology, these analyses did not yield statistically significant results (P > 0.05).

Fig. 3
figure 3

(A) Scatter plot presents the correlation between the observed VIQ and the predicted VIQ change derived from the RVR analysis (P = 0.027), Each dot represents the data from one patient, and the dashes indicate the 95% prediction error bounds. (B) Nonparametric permutation tests (P = 0.041). (C) The predictive weight of each feature was obtained and mapped onto the brain surface. Brain maps in this figure were generated using the pysurfer package in python (PySurfer v0.10.0 https://pysurfer.github.io/).

Full size image

Association between alterations in modular variability and gene expression profiles

Partial least squares regression analysis

We used PLS regression analysis to identify gene expression patterns associated with IGE-related dynamics alterations. PLS3 explained the largest spatial variance (17.42%) in modular dynamics based on the case–control difference map (P < 0.05, corrected for spatial autocorrelation). The spatial distribution of PLS3 weights was matched to the case–control difference map (Fig. 4A), Then, the PLS3 score map was spatially correlated with the group difference t map (r = 0.418, P < 0.005; Fig. 4B). After normalizing and ranking PLS3 weights, 349 significantly upregulated genes (PLS+) and 636 downregulated genes (PLS−) were identified (P < 0.05; FDR) (Fig. 4C). The complete lists of differentially expressed genes are provided in Supplementary Data File 1.

Fig. 4
figure 4

Gene expression profiles related to connectome dynamics differences. (A) The case–control t map of the regionally module variability. (B) A scatterplot of regional PLS3 scores (a weighted sum of 12,791 gene expression scores) and regional changes in connectome dynamics (r = 0.418, P < 0.005). (C) Ranked PLS3 loadings. Functional enrichment analysis of top ranked genes (from C) with PLS3+ (D) and PLS3− (E) weights (q < 0.05), most of the genes were annotated using the GO database, followed by KEGG database (www.kegg.jp/kegg/kegg1.html). Left Panel: Bubble plot of the top 20 significant enrichment terms of PLS3 ± genes and visualization of enriched ontology terms of PLS3 ± gene. Right Panel: Metascape enrichment network visualization showing intra-cluster and inter-cluster similarities of enriched terms. Each term is represented by a circle node, with size proportional to the number of input genes included, and color indicating cluster identity (nodes of the same color belong to the same cluster). Brain maps in this figure were generated using the pysurfer package in python (PySurfer v0.10.0 https://pysurfer.github.io/).

Full size image

GO, KEGG, and cell-type enrichment analysis

Metascape analysis of PLS+ and PLS− gene sets revealed that the PLS+ genes were enriched in GO biological processes, including synaptic signaling, glutamatergic synapses, and regulation of kinase activity, and in the KEGG neuroactive ligand-receptor interaction pathway (Fig. 4D). Conversely, PLS− genes were enriched in processes such as phosphorylation, kinase binding, and mRNA metabolism but showed no KEGG pathway enrichment (Fig. 4E). Mapping PLS+ and PLS− genes to seven classical cell types revealed a significant association with high oligodendrocyte precursor cells (OPCs) expression (count = 6, P = 0.009, FDR, Fig. 5A; count = 6, P = 0.01, FDR, Fig. 5B).

Fig. 5
figure 5

Cell type specificity of PLS+/− genes associated with IGE MV changes. (A) The number of overlapping PLS+ genes identified for each cell type. (B) The number of overlapping PLS− genes identified for each cell type. The x-axis represents count of overlapping PLS+/− genes, and the y-axis indicates the cell types.

Full size image

Validation

To ensure repeatability, we evaluated the effects of age, gender, head movement, and parameter variations, including window sizes (30 TRs, 40 TRs, and 70 TRs with a step of 1 TR), topological parameters (γ = 1, 0.9, 0.8, 0.7), and temporal coupling values (ω = 1, 0.9, 0.8, 0.7). The findings remained consistent, with spatial correlations ranging from 0.34 to 0.93, and all spatial permutation tests showed P < 0.001 (Supplementary Figs. S2–S6).

Discussion

Using a multilayer brain network model, our study provides neuroimaging evidence of altered connectome dynamics in children with IGE. Increased dynamics were primarily localized in the default mode network, while decreases occurred in the somatomotor and frontoparietal control networks. MV alterations in the left lateral prefrontal cortex, precuneus/posterior cingulate cortex, and right orbitofrontal cortex correlated with cognitive performance and significantly predicted VIQ. Additionally, transcription-neuroimaging spatial correlation analysis revealed that IGE-related connectome dynamics were spatially linked to the expression of 985 genes enriched in biological pathways related to neuroactive ligand-receptor interaction and synaptic signaling, bridging neuroimaging alterations with genetic components and offering novel insights into the neuropathology and cognitive impairment of childhood IGE.

Previous studies based on rs-fMRI often assumed static connectivity, focusing on time-averaged patterns37. However, dynamic reconfiguration of brain connectivity occurs across various time scales, with regions exhibiting transient stability and rapid transitions. The assumption that connectivity patterns reflect only time-averaged features between regions overlooks these dynamic changes, limiting a comprehensive understanding of brain activity38,39. Time-varying multi-brain network analysis captures these temporal fluctuations, revealing the spatiotemporal characteristics of brain networks40. In IGE, synchronous abnormal neuronal discharges result in dynamic network alterations24, with studies showing disrupted functional connectivity and excitation-inhibition imbalance, which may trigger seizures and affect brain function17,41. Elucidating these dynamic features is essential to understanding the pathophysiological mechanisms of IGE. The multilayer network model addresses this by linking corresponding regions across time points, enabling the analysis of network configurations42. In this study, at the subnetwork level, we observed abnormally elevated module switching in the DMN of children with IGE, suggesting increased transitions between the DMN and other brain networks over time, resulting in excessive functional integration, this represents an unstable network state and a failure of normal information filtering mechanisms, where the DMN becomes overly receptive to inputs from typically segregated networks, potentially contributing to seizure propagation and cognitive impairments in IGE. These findings align with a previous study43, and demonstrate the utility of multilayer brain network analysis in exploring complex neural dynamics44. Our research highlights dynamic abnormalities in childhood IGE at the system level, providing deeper insights into its pathophysiological mechanisms and emphasizing the value of dynamic connectome topology analysis.

At the nodal level, increased nodal MV was primarily observed in the DMN, including the bilateral precuneus/posterior cingulate cortex and prefrontal cortex, while decreased nodal MV occurred predominantly in the SMN, FPN, and DMN. These findings indicate abnormal alterations in the dynamic community structure of these brain regions in children with IGE. Regions with altered nodal MV align with findings from our previous studies using different analysis methods9, suggesting more severe functional impairments in these regions. Normal cognitive function relies on synchronized interactions between tightly connected functional brain networks. Disruptions in this balance may lead to deficits across multiple networks, resulting in fragmented information processing within isolated modules and impairing cognition45. The DMN and FPN are critical for cognitive function and affective processing46,47, and abnormalities in their dynamic functional networks may affect cognitive-executive control. Correlation analysis revealed significant associations between MV changes in the left lateral prefrontal cortex, precuneus/posterior cingulate cortex, and right orbital frontal cortex and IQ scores in children with IGE. This suggests a relationship between cognitive impairment and functional plasticity in specific brain regions, highlighting dynamic network parameters as potential biomarkers of brain functional impairments. Furthermore, while regional variations in nodal MV changes, the predominant elevation in the dynamic network properties of the DMN may represent a compensatory mechanism in children with IGE, enhancing network flexibility to maintain cognitive function.

A prior study utilized multilayer network connectivity patterns as predictors of higher cognitive abilities48. To explore this potential, we employed the RVR algorithm to construct a predictive model using the temporal correlation coefficient to estimate VIQ in children with IGE. This model integrates information from all brain voxels in a high-dimensional space, improving sensitivity to subtle spatial correlations and detecting voxel distributions that most effectively predict VIQ. Our results highlighted specific cortical regions with higher weights, including the right cingulate, parietal, and left prefrontal cortices. Specific regions of the parietal cortex are associated with speech perception, while the frontal cortex is with speech production, and the cingulate cortex, a crucial DMN component, is closely involved in the visuospatial contextual memory modulation49. The high feature contribution weights of these brain regions may reflect their pivotal roles in speech-memory-related functions50. Identifying brain regions with the highest feature contribution weights can enhance our understanding of neural networks underlying cognitive processes in daily behavior and brain-cognition relationships. But, as a preliminary exploratory analysis, these observations should be interpreted cautiously as they represent initial findings that require replication in independent cohorts. The model’s performance, while promising for this exploratory investigation, will need further validation to establish its reliability and clinical utility.

Previous studies suggest that imaging transcriptome analysis improves understanding of disease pathology51,52. To explore the complex pathogenesis of IGE, we integrated connectome dynamics with human brain gene expression data from the AHBA20 and used PLS regression analysis to identify genes associated with dynamic functional network changes in children with IGE. The observed correlations between MV alterations and cortical gene expression patterns may reflect potential multi-genetic, multi-pathway mechanisms through which IGE could influence brain function. The identified gene modules, which encompass diverse cell types and biological pathways, are consistent with the complex pathophysiology of IGE, though further validation is needed to establish specific causal relationships. Among these genes, GABRE has been previously implicated in genetic epilepsies through its potential role in synaptic inhibition and GABA receptor function, as suggested by studies53,54. Our findings showing spatial correlations between GABRE expression patterns and network alterations are consistent with this proposed mechanism, though direct causal relationships remain to be established. Similarly, BCL11A expression covaried with network measures, aligning with prior reports linking this gene to neuronal development55. GO and KEGG functional enrichment analyses identified pathways linked to altered brain module dynamics in IGE, including neuroactive ligand-receptor interactions, kinase activity regulation, synaptic signaling, and phosphorylation. Cell-type enrichment analysis revealed significant involvement of OPC in both upregulated and downregulated genes. Under normal conditions, OPCs proliferate and differentiate into myelinated oligodendrocytes in the central nervous system. Chronic recurrent seizures may cause myelin damage, which can exacerbate epileptic activity. Additionally, OPCs may contribute to myelin repair and aid seizure control56,57.

This study had several limitations. First, this study was conducted at a single center with a small sample size, although we employed rigorous cross-validation and permutation testing to ensure robustness, replication in larger, multi-center cohorts is essential to verify our results. Second, due to the limited sample size, we could not fully account for the potential effects of medications used in children with IGE, necessitating further analysis in a larger sample size and future studies employing medication-naïve IGE cases would provide crucial complementary evidence. Third, our study relies on spatial correlations between in vivo functional connectivity patterns and postmortem gene expression data from different individuals. While this approach has been widely adopted in imaging-transcriptomic studies, the inherent methodological differences between these data modalities require cautious interpretation of the observed associations, future studies combining in vivo imaging with subject-specific transcriptomic profiling could help validate these findings. Finally, while we developed predictive models for cognitive function, these models underwent only internal cross-validation, highlighting the need for external validation with a larger sample size.

Conclusion

In conclusion, our study demonstrated altered brain modular dynamics in children with IGE, with increased MV in the DMN and decreased MV in the SMN and FPN. MV alterations in specific brain regions correlated with cognitive function, and dynamic network-based prediction models significantly predicted cognitive performance. Integration of neuroimaging and transcriptional data identified 985 genes associated with MV changes, highlighting diverse functional characteristics. These findings provide evidence of dynamic network alterations and link neuroimaging changes to genetic mechanisms, deepening our understanding of childhood IGE.