We first identified seven consistent and representative group-wise static FBNs for both session A and session B via the DLSC approach. Figure 1 shows the representative FBNs of session A. These networks include either typically activated simple networks or complex networks. The simple networks involve the visual network (Figure 1A) and the auditory network (Figure 1B). The complex networks consist of multiple co-activated brain networks/regions, including auditory and cerebellar network (AC) (Figure 1C), the audiovisual and sensorimotor network (VAS) (Figure 1D), the partial default mode network (DMN), the salience and cerebellar network (pDSC) (Figure 1E), the DMN and cerebellar network (DC) (Figure 1F), and the dorsal attention network (DAN) (Figure 1G). Specifically, the AC network is primarily composed of auditory, cerebellar posterior crus 1,2 and vermis (Figure 1C). The VAS network is composed of visual, auditory, and sensorimotor cortex (Figure 1D). The pDSC encompasses the posterior cingulate cortex, medial prefrontal cortex, angular gyrus, anterior insula, dorsal anterior cingulate cortex, cerebellar posterior crus1,2, cerebellums 9 and vermis. Notably, the pDSC network excludes the precuneus (Figure 1E). The DC network mainly consists posterior cingulate cortex, medial prefrontal cortex angular gyrus, precuneus, cerebellar posterior crus1, 2, cerebellums 9 and vermis (Figure 1F). The DAN network includes intraparietal sulcus and the frontal eye fields (Figure 1G). A comparison between these identified FBNs and well-established resting-state templates or networks from previous studies conducted under natural stimulation is presented in Supplementary Figure S8.

The FBNs derived from session B showed a high degree of spatial consistency with those observed in session A (Supplementary Figure S1), as demonstrated by the relatively high overlap rate and Pearson correlation coefficient (PCC) values between the two sets of FBNs (Table 1). Specifically, the mean overlapping rate and the mean PCC of the seven FBNs were 0.44 ± 0.11 [Mean ± standard deviation (SD)] and 0.82 ± 0.18 (Mean ± SD), respectively, suggesting the consistency and stability of the DLSC framework in detecting FBNs across two scanning sessions.

We applied the sliding time window method with a window length of 60 repetition time (TR) units and step size of 1TR, resulting in 470 available windows. Correspondingly, 470 FBNs were obtained by applying the DLSC method, which could reflect the dynamics of time-varying large-scale networks. To provide representative visualization of these FBNs, we selected and displayed the FBNs from the first window among every 50 windows. For example, the first brain activation map in Figure 2 represents the visual network obtained during the first window (1TR to 60TR), corresponding to a time period of 1 s to 132 s.

To quantitatively explore the dynamic spatial patterns of FBNs, we assessed the dynamic temporal changes of the number of activated voxels (NAV) (Figure 2) and the intensity of activated voxels (IAV) (Figure 3) for the seven FBNs. Our results revealed that for each network, both NAV and IAV exhibited temporal variations, with relatively consistent trends between the two metrics. The IAV showed less variability compared to NAV due to the calculation method of averaging activation intensities of all voxels exceeding a predefined threshold, resulting in a relatively narrow range of variation in the overall activation strength of the whole network. The peaks of NAV and IAV curves corresponded to the FBNs that displayed more pronounced and widespread patterns of activation, whereas the troughs of these curves responded to FBNs with diminished or even absent activation patterns. These findings highlight that FBNs were dynamic and evolved temporally in response to the unfolding plot of the movie, which is also consistent with the underlying neural basis of complex perception and behavior (Calhoun et al., 2014). Additionally, the lower-order perceptual networks, including visual network, auditory network, and VAS network, exhibited relatively stable level of activation over time, whereas the higher-order networks, such as pDSC, DC, and DAN networks, showed greater fluctuations in activation curves. The AC network, specifically, comprising both lower-order network (i.e., auditory network) and higher-order networks (i.e., cerebellar network), also displayed substantial fluctuations in its activation curves (Figures 2, 3). These results suggest that different FBNs exhibit distinct temporal dynamics in response to external stimuli, which may reflect their respective roles in higher-level cognitive and attentional processes.

The results of session B were generally consistent with those of session A, as evidenced by the visual comparison of the results depicted in Supplementary Figures S2, S3. Moreover, the PCC values for both NAV and IAV curves across session A and session B were relatively high for most FBNs (Supplementary Table S1), suggesting that the identified dynamic spatial patterns of FBNs induced by the movie viewing are reproducible and consistent across two scanning sessions.

To investigate the relationship between DFNs and the unfolding of the movie, we assessed the spearman correlation between dynamic changes of NAV/IAV and movie annotations, which include language use, changepoints, the presence of positive valence of scenes (scenes_p), the presence of negative valence of scenes (scenes_n), the presence of faces with positive (face_p), and presence of faces with negative (face_n). The results showed that two DMN-related networks, i.e., the DC and pDSC networks, were significantly correlated with movie annotations. Specifically, both NVA and IVA metrics of the DC network showed statistically significant correlations with the appearance of positive facial expressions based on permutation-based testing (p < 0.05) (see Methods) (Tables 2, 3). Notably, the value of p for the IAV metric was less than 0.01 (Table 3). Additionally, the changes in INV of the pDSC network were significantly and positively correlated with the appearance of the changepoint in the movie scenes (permutation 5,000 times, p < 0.05) (Table 3).

The neural response evoked by the naturalistic stimuli exhibit not only high consistency across individuals, but also inter-subject variability and uniqueness reflecting personal experiences and intrinsically-driven processes under natural viewing condition, which varies across different brain regions/networks (Golland et al., 2007; Ren et al., 2017b). Hence, to quantify these group consistency and individual variations in defined DFNs, we adopted their corresponding group-wise static FBNs as templates to calculate the dynamic inter-subject correlation (ISC) (see Methods). Accordingly, the group-level dynamic ISC can represent the degree of temporal consistency across subjects in different FBNs (the thick blue line in Figure 4). The average values of group-level dynamic ISC during the entire period for seven FNBs (including visual, auditory, AC, VAS, pDSC, DC, and DAN networks) were 0.33 ± 0.12, 0.51 ± 0.13, 0.58 ± 0.12, 0.39 ± 0.12, 0.25 ± 0.10, 0.31 ± 0.08, and 0.34 ± 0.08 (Mean ± SD), respectively. While relatively high ISC values were observed in networks encompassing lower-level perceptual regions, especially those related to auditory processing, such as auditory and AC networks, the higher-order networks demonstrated lower ISC values that can indicate the occurrence of intrinsically-driven processes during individual movie viewing, including pDSC and DC networks, consistent with previous research (Ren et al., 2017b). Moreover, individual-level dynamic ISC also showed inter-subject variations especially in those higher-order networks under movie stimuli (colorful thin lines in Figure 4).

The dynamic ISC of session B were largely consistent with session A (Supplementary Figure S4). Quantitatively, the PCC values for seven FBNs across two sessions were relatively high (Table 4), with an average PCC value of 0.82 ± 0.18 (Mean ± SD), thereby reaffirming the consistency of DFNs identified by our DLSC framework across two sessions.

In the preceding section, there was relatively lower consistency in individual neural responses observed in higher-order brain networks, such as the pDSC and DC networks. This variability may be indicative of individual differences and unique experiences during natural viewing conditions. To investigate this assumption further, that is, exploring the potential relationship between subjective movie viewing experiences and the dynamics of DFNs, we examined whether the individual dynamic ISCs were correlated with their personal ratings of the movie. We applied an inter-subject representational similarity analysis (IS-RSA) (see Methods) to explore whether participants with similar subjective ratings also exhibited similar dynamic neural response patterns. Specifically, we employed a multidimensional scaling approach (MDS) (Carroll and Arabie, 1998) to characterize the answers to the post-movie questionnaire. Our result revealed that participants had varying experiences while watching the movie, with some reporting high engagement characterized by low boredom, high enjoyment, high emotion, and high audio quality, while others had low engagement (Figure 5A). The distances between movie ratings were measured by the Euclidean distance matrix of questionnaire answers across all individuals (Figure 5B). We computed Pearson distance to represent the inter-subject distances of dynamic ISC values for seven representative FBNs, respectively (Supplementary Figure S5). By evaluating the correlation between the movie rating distances and the inter-subject distances of the dynamic ISC, we found significant positive correlations (permutation 5,000 times, p < 0.05) in three cerebellum-related networks, that is, AC, pDSC, and DC networks. The distance matrices of dynamic ISC for these three networks are presented in Figure 5C, and their simple linear regressions are shown in Figure 5D. However, the other DFNs did not show statistically significant associations (p > 0.05) (Table 5). We did not repeat this experiment in session B as it involved a repeated viewing of the same movie, and the post-viewing questionnaire was not conducted.

It was assumed that similar sensory experiences would lead to the emergence of DFNs in a consistent and reliable manner. Therefore, we assessed the level of reliability of the DFNs that develop in response to the movie storyline across two sessions. Specifically, we first calculated the scan-wise intra-group correlation coefficient (ICC) values for seven static FBNs. The results showed that the visual, AC, and VAS networks exhibited excellent reliability, the auditory and DAN networks possessed good reliability levels, and the pDSC and DC networks had moderate reliability, indicating that the networks associated with primary perceptual processes were relatively more reliable, while the higher-level networks showed less reliability, consistent with previous studies (Choe et al., 2017) (Figure 6B).

We further analyzed the reliability of the DFNs (see Methods). The reliability of the auditory, AC, VAS, and DC networks, triggered by this touching movie, gradually increased during the mid to late period (about 300–390 window, corresponding to 600–1,000 s time period), and reached its peak in the near-end period. This is consistent with the narrative structure of the movie (Supplementary Table S2), wherein the plot also reaches its climax near the end (at around 17 min or 1,020 s) (Nguyen et al., 2017). The findings suggest that as the storyline develops, individuals may experience greater behavioral constraints and engagements, leading to an increase in the test–retest reliability of brain activities.

The experiment consisted of two scanning sessions. Following a first session (session A) conducted 3 months earlier, participants underwent a second scanning session (session B) employing the same experimental paradigm. In each session, participants freely watched the 20-min short film “Butterfly Circus.” In addition, all participants completed a questionnaire immediately after session A.

The short film, “Butterfly Circus,” depicts a touching story of a limbless man who is encouraged by the showman of a renowned circus to discover his true potential. The narrative architecture of the film follows three distinctive drama acts that feature significant developments, complications, and turning points (Supplementary Table S2). Additionally, basic movie annotations were provided, including: the use of language, changepoints, the presence of positive valence of scenes, the presence of negative valence of scenes, the presence of faces with positive, and the presence of faces with negative (Supplementary Figure S6). Further details regarding the participants can be found in the Supplementary material (1.2).

Sixteen right-handed (ages 27 ± 2.7) healthy participated in this study. FMRI images were acquired from a whole-body 3 T Siemens Trio MRI scanner with the following scanning parameters: repetition time (TR) 2,200 ms, echo time (TE) 30 ms, flip angle (FA) 79°, the field of view (FOV) 192 × 192 mm, pixel bandwidth 2,003 Hz, a 64 × 64 acquisition matrix, 44 axial slices, and 3 × 3 × 3 mm 3 voxel resolution. Functional images were preprocessed using FMRI Expert Analysis Tool (FEAT) from FMRIB’s Software Library (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki), which involved motion correction, slice timing correction, spatial smoothing with 6-mm full width at half maximum Gaussian kernel, band pass filtering (0.0085 ± 0.15 Hz), linear registration to the standard Montreal Neurological Institute space (2 mm MNI152 standard template), and masking.

To discover and characterize DFNs, we proposed a computational framework comprised of two stages: (A) using group-wise dictionary learning and sparse coding (DLSC) to represent static FBNs (Figure 7A), (B) sliding-window method applying for the representation of dynamic spatial patterns of FBNs (Figure 7B).

In stage A, first, the whole-brain fMRI signals of each subject were extracted and stacked into a 2D matrix Si(Si∈ℝt×n)), followed by spatial concatenation of the individual matrix Si into a group-wise matrix 𝑺 (Si∈ℝt×(p×n)) (top panel in Figure 7A), where t represents the time length of fMRI signals, n refers to the number of the individual brain voxels, p stands for the number of subjects. Second, we applied the dictionary learning algorithm to the group-wise matrix 𝑺 to learn a meaningful group-wise dictionary 𝑫 (D∈ℝt×K). This dictionary 𝑫 consists of K atoms that can well represent the temporal features embedded in naturalistic fMRI data and are commonly shared by all the subjects (Ren et al., 2017a; Ge et al., 2020). Hence the group-wise dictionary 𝑫 can be used to sparsely represent individual fMRI signals Si, resulting in the individual spatial patterns αi(αi∈ℝK×n) (K < <n) that exhibit correspondences across subjects (middle panel in Figure 7A). Notably, we used the online dictionary learning and sparse coding algorithm, which is an effective method for extracting instinctive information from the original signal and is suitable for group-level data operations (Ponce and Sapiro, 2010; Lv, 2013). Third, to obtain the group-wise static FBNs, we performed one-sample t-test on each element of individual loading coefficient matrix αi (middle panel of Figure 7A). Specifically, for all the subjects together, we hypothesized that each element in the loading coefficient matrix αi is group-wisely null. To evaluate this assumption, we conducted one-sample t-test on the corresponding element in the loading coefficient matrix αi for different subjects, in order to test whether this hypothesis was accepted or rejected (Ren et al., 2017a). The resulting t-value was then transformed into a z-score, forming a group-wise loading coefficient matrix A composed of z-scores (Friston et al., 1994). Since each individual coefficient matrix αi is sparse, the t-test result of the group-wise loading coefficient matrix A is also sparse. Afterwards, each row of matrix A can be mapped back to brain volume with z-scores, referred to as z-score maps. Consequently, the z-score map obtained from this analysis can depict spatially consistent activation across all subjects, thus representing the static group-wise FBNs (bottom panel of Figure 1A).

In stage B, in order to obtain a series of dynamic spatial patterns that evolve over time for each subject, we slid the same time window on the individual signal matrix Si and group-wise dictionary 𝑫 simultaneously. This approach not only allows us to establish correspondence of individual-level FBNs among individuals, but also generates the corresponding dynamic FBNs. Consequently, we obtained individual signal matrices Si1,Si2,…,Sij,…,Siw, which represent the individual’s signal for each window, as well as group-wise dictionaries D1,D2,…,Dj,…,Dw, which represent the group-wise dictionary for each window (top two panels in Figure 7B). The chosen window length, denoted as L (in our study, L = 60 TRs with a step size of 1 TR), resulted in a total of w = 470 windows. Consequently, the individual signal matrix Siand group-wise dictionary 𝑫 were divided into multistage signal matrices Sij and dictionaries Dj (j∈1, 2, …, w) (top two panels in Figure 7B). Next, based on Sij and Dj of each window, we leveraged sparse representation algorithms to extract a sequence of coefficient matrices αij to represent individual dynamic spatial patterns of FBNs (bottom panel in Figure 7B). The above experiments were performed on each subject in session A and session B.

To identify the matching FBNs across two sessions, we manually selected representative FBNs from session A and match them with responding FBNs in session B based on their highest Pearson Correlation Coefficient (PCC) values and the maximum number of overlapping voxels (overlapping rate) (Benesty et al., 2009; Lv et al., 2015a). A detailed pipeline for selection of representative FBNs can be found in the Supplementary material (1.4). The PCC was defined as the correlation between the representative FBNs of session A and session B, and the overlapping rate of the FBNs in session A and session B was defined quantitatively as:

where 𝑿 is the representative FBN of session A, and Y refers to responding FBN of session B.

The DLSC algorithms rely on two key parameters: the number of dictionary atoms (K) and the sparsity penalty parameter (λ). However, there is no golden standard for determining the optimal values of these parameters. Based on previous studies that used DLSC algorithms for FBNs recognition, K was experimentally and empirically set to 400, and λ was set to a range of 0.1–0.5 (Lv et al., 2015b; Zhang et al., 2017; Ge et al., 2018). Therefore, in this study, we set K = 400, while systematically tested the setting of λ (0.1, 0.5). Through our experiments, we found that the largest number of networks could be identified with manual inspection when utilizing the parameters K = 400 and λ = 0.5. Consequently, we chose this set of parameters.

The window size is a crucial parameter that can determine the tradeoff between time resolution and estimation results. Previous related studies have empirically converged to window size values between 30 and 240 s (Hutchison et al., 2013; Preti et al., 2017). Additionally, Savva et al. (2019) suggested window size of at least 120 s to ensure the reproducibility of the result. Hence, we selected a window width of 60TRs (132 s).

To quantify the dynamic changes of spatial patterns of FBNs, we employed two methods including calculating the number of activated voxels (NAV) and the intensity of activated voxels (IAV) of each FBN across all windows, respectively. Specifically, NAV was derived by counting the number of all voxels exceeding the threshold value (z = 1.65), while IAV was obtained by averaging the intensities of all voxels above this threshold. The group-wise NAV/IAV was derived by computing the average value of the NAV/IAV across all individuals.

The sliding window method produced 470 DFNs with a window length of 60 TRs and a step size of 1. Correspondingly, the duration of the scan was 530 TRs, with each TR corresponding to an annotation, resulting in 530 movie annotations in total. To establish correspondence between the DFNs and movie annotations, we selected the movie annotations occurring at the center point of each temporal window to correspond to each DFN based on previous studies (Simony et al., 2016; Tzachor and Hoshen, 2022). Specifically, we focused on a segment of movie annotations spanning from 31 to 500 TRs and examined their correlations with corresponding 470 dynamic FBNs.

The Spearman correlation coefficient between group-wise NAV/IAV and annotation vectors was computed to explore the association between group-wise dynamic changes in spatial patterns of FBNs and movie annotations, which were constructed for Language, Changepoint, Positive Scenes, Negative Scenes Positive Faces, and Negative Faces. Movie annotations were converted to vectors of 0 and 1 based on their onset and offset times (Supplementary Figure S6). To establish the statistical significance of the observed correlations, the correlation between the group-wise NAV/IAV and annotation vector was re-calculated 5,000 times by shuffling the vectors in each iteration. The observed correlation was compared with a null distribution of correlation generated by 5,000 permutations. If the observed correlation falls in the extreme tails of the distribution (i.e., the value of p is less than 0.05), we can conclude that there is a significant difference in group-wise NAV/IAV and movie annotation.

Inter-subject correlation (ISC) analysis measures the inter-subject consistency for temporal responses across participants (Hasson et al., 2004; Di and Biswal, 2020). To evaluate the ISC of dynamic FBNs, we first used the group-wise static FBNs as masks to extract the time-series signals of the corresponding FBN for each participant. Next, we averaged all the time-series signals within FBN, resulting in the FBN-level time-series signals for each FBN. Subsequently, we also applied the sliding window strategy and calculated the ISC of the FBN-level time-series signals in each time window for each subject, where the time window size was set to the same value as that in the “Dynamic sparse representation” section (i.e., 60TRs). Consequently, we derived the dynamic ISC metric for each subject for each representative FBN. To calculate the dynamic group-wise ISC metric, we performed Fisher z-transformation on the ISC values of all subjects for each window and subsequently averaged the ISC value across all individuals for each window.

The study employed a post-movie questionnaire to collect participants’ subjective appraisals of the movie, which consisted of eight questions. However, the RSA analysis excluded four questions because there was insufficient variability among participants (Supplementary Data). The remaining four questions are more focused on evaluating the movie subjectively, that is, how participants rated their feeling during the first movie session, including boredom, enjoyment, feeling happy or sad, and audio quality. Regarding question 4, specifically, the audio quality does not vary while recording, and each participant said they all had a comparable understanding of the movie’s plot. The participants’ level of engagement may have influenced how they rated the scale. All questions in the survey utilized a 1 to 5 rating scale. To represent participant differences in movie ratings, we employed a multidimensional scaling method to map responses to the questionnaire onto a two-dimensional representation.

Inter-subject representational similarity analysis (IS-RSA) is a promising approach for examining the potential relationship between inter-subject variability in brain dynamics and individual differences in behavioral disposition (Kriegeskorte et al., 2008; Finn et al., 2020; Meer et al., 2020). Thus, we conducted the IS-RSA to assess the correlation between post-hoc behavioral movie ratings and dynamic ISC distances across all subjects.

We constructed inter-subject distance matrices to represent movie impressions and dynamic ISC. Specifically, inter-subject distances for movie impressions were calculated by measuring the Euclidean distance of questionnaire ratings between each possible pair of participants, resulting in 16 (number of participants) × 15 matrices. To examine the dynamics ISC distance, we calculated the Pearson distance between the dynamics ISC matrices for every possible pair of participants, producing a dynamic ISC distance matrix of size 16 × 15.

To assess the strength of associations between the movie ratings and dynamic ISC, we calculated the Spearman correlation between the lower triangular parts of the above two distance matrices. To assess the statistical significance of the results, we performed permutation testing 5,000 times. For each iteration, we squeezed the two matrices, dynamics ISC distance and movie rating distance, into row vectors, and randomly selected a new starting point for each row vector. This procedure allowed us to generate a null distribution of correlations and determine whether the observed correlation was significant.

To assess the level of reliability of dynamic FBNs during the natural viewing conditions, we calculated the test–retest reliability of the matching dynamic FBNs across two sessions. Specifically, we measured the intra-group correlation coefficient (ICC) for each window to determine the level of consistency in the FBNs across time (Shrout and Fleiss, 1979; McGraw and Wong, 1996). For comparison, we also calculated the static test–retest reliability of FBNs by calculating ICC over the entire period. ICC can be defined by the following equation:

Here, d refers to the number of observations, which in our study was equal to 2. MSp represents the mean square variation between subjects, while MSe represents the mean square variation within subjects. The test–retest reliability was divided into five levels: excellent (ICC > 0.8), good (ICC 0.6–0.79), moderate (ICC 0.4–0.59), fair (ICC 0.2–0.39), and poor (ICC < 0.2). The test–retest reliability was assessed at the scan-wise level, and the methodology for this process was carried out in accordance with the previous study (Guo et al., 2012; Wang et al., 2017).