Temporal variability of brain–behavior relationships in fine-scale dynamics of edge time series Open Access

The link between individual variation in functional connectivity (FC) and behavior is of central interest in human neuroscience, with clear developmental and clinical applications. Such brain–behavior relationships are typically evaluated using FC computed from continuous and extended brain recordings, capturing minutes of time. Longer time frames and multiple scan sessions have been shown to improve reliability of FC (Birn et al., 2013; Cho et al., 2021; Elliott et al., 2019; Noble et al., 2017) and support improved ability to distinguish individuals (e.g., “fingerprinting”) (Finn et al., 2015; Pannunzi et al., 2017). However, individual differences in functional networks have been discovered at shorter timescales within a scan (Cutts et al., 2023; Peña-Gómez et al., 2018; Van De Ville et al., 2021) raising the possibility that brain–behavior associations may be time dependent. Indeed, extensive literature on time-varying functional connectivity reveals that connectivity patterns fluctuate over time and can be related to different behavioral measures (Eichenbaum et al., 2021; Liégeois et al., 2019; Lurie et al., 2020). Such approaches typically analyze changes in time series characteristics (Betzel et al., 2016; Chang & Glover, 2010; Liégeois et al., 2019; Zalesky et al., 2014), transitions or dwell times of derived network states (Allen et al., 2014; Calhoun & Adali, 2016; Vidaurre et al., 2017), instantaneous measures of FC (Liu et al., 2018; Pedersen et al., 2018; Yaesoubi et al., 2018), or most commonly through smaller sliding windows within a scan session (Calhoun & Adali, 2016; Cohen, 2018; Leonardi & Van De Ville, 2015; Lurie et al., 2020; Preti et al., 2017; Sakoğlu et al., 2010). These fluctuations as well as changes in individual differentiation suggest the possibility that brain–behavior relationships may be stronger during different connectivity patterns only visible from specific moments in time within a scan. Here we assess whether the strength and similarity of brain–behavior relationships across groups of subjects can be improved through targeted selection of time points within a scan session.

Selecting or “filtering” the time series for specific connectivity patterns across subjects could enhance the specificity of behavioral associations. Traditional windowing approaches (Calhoun & Adali, 2016; Leonardi & Van De Ville, 2015; Lurie et al., 2020; Preti et al., 2017) compare connectivity from consecutive frame subsets of varying lengths. However, discontinuous moments can also be assessed based on shared properties of the selected frames. As such, a “filter” could be utilized to retrieve similar moments across subjects, independent of temporal alignment. It has already been shown that fluctuations within a scan exhibit differences in cortex-wide connectivity structure that could be retrieved from such filters. For instance, Sporns et al. (2021) show that a resting-state scan is composed of the transient appearance of co-fluctuations (or coordination of activity) within and between different functional systems. This suggests that apparent connection strength between a pair of brain regions changes over time in relation to variations in system-level or global coordination. As such, even strong co-fluctuations between brain regions at a given moment may not always provide information regarding how that pair of regions relates to a behavior of interest. It is unclear whether data from a full scan include patterns of co-fluctuations that dilute brain–behavior correlations or whether separate moments in time better emphasize different behavioral associations. With approaches looking at instantaneous moments within a scan, similar subsets of frames could be reliably retrieved from the time series and examined across individuals for their relationship with behavioral or phenotypic measures of interest. This can be done by isolating frame subsets from the time series based on the similarity between a template and the activity occurring at specific time points. These templates can be developed for specific behaviors or phenotypes to retrieve moments that amplify brain–behavioral relationships with the measure of interest. Targeted selection of time points based on specific behavioral or phenotypic features may provide additional information unseen using the full length of data.

Past work has found evidence of differences between analyses computed on full FC versus targeted selection of frames within a scan. Edge time series (eTS) is one such method that has been used to deconstruct FC into co-fluctuations between brain regions at single time point resolution (Faskowitz et al., 2020). With this approach, Esfahlani et al. (2020) found that only a handful of high-amplitude moments in time (“events”) dominates the signal seen in full FC and these moments carry most of the individual differences traditionally studied (Betzel et al., 2022). Although these “events” contribute the most information toward full FC, they provide only a sparse representation of the whole time series. This suggests that most brain–behavior relationships have been performed on information from those few moments of high co-fluctuation amplitude. Targeted selection of time points irrespective of co-fluctuation amplitude has already been shown to improve the relationship between specific connectivity patterns with behavior as well as reveal additional or unique information from time points with lower representation in full FC. For instance, Cutts et al. (2023) show that subject identifiability (Amico & Goñi, 2018; Finn et al., 2015) can be improved above full FC with multiple selections of discontinuous framesets. Other work has found that cognitive function is best predicted by frames distributed across the time series with varying co-fluctuation amplitudes (Wehrheim et al., 2023). Additionally, binning time points based on co-fluctuation amplitude reveals distinct relationships with behavioral prediction (Sasse et al., 2023) as well as clinical phenotype (Chumin et al., 2024) unseen in full FC. These findings suggest that developing techniques to distinguish which moments in time contribute most toward behavioral relationships could aid in improved behavioral and diagnostic prediction.

Moments with higher co-fluctuation amplitude have also been found to be more similar across subjects compared with frames with lower amplitude (Esfahlani et al., 2020; Sporns et al., 2021). It is unclear the extent to which greater variance of lower amplitude frames might explain additional, largely unexplored relationships with behavioral variance. However, comparing moments between individuals from resting state proves challenging due to fluctuations occurring independently of time-locked stimuli. High co-fluctuation amplitude moments can be systematically retrieved from the time series and compared across individuals due to their similarities in features (Betzel et al., 2022; Esfahlani et al., 2020). Moments with mid to low co-fluctuation amplitude are more diverse and, therefore, harder to compare across individuals. One solution to reliably target lower amplitude frames is to compare moments across subjects with similar patterns of co-fluctuations between brain regions. This information can be derived from the nodal level based on whether two brain regions display similar or dissimilar activity at a given moment in time (Faskowitz et al., 2020). These moments can be retrieved using a “filtering template” to find the time points with activity patterns that best match the template and then perform analyses on a reconstructed component similar to FC from the selected frames (Chumin et al., 2024; Cutts et al., 2023; Sporns et al., 2021). The structure of these templates can additionally be optimized to improve selection of framesets that improve behavioral relationships. Comparing moments across individuals is then based on connectivity derived from cortex-wide patterns of activity rather than the average across longer time periods dominated by high co-fluctuation moments.

In this study, we assessed whether brain–behavior correlations could be improved through targeted selection of time points irrespective of their contribution to standard FC. First, brain–behavior correlations across 58 separate behaviors were compared across levels of co-fluctuation amplitude, where lower levels contribute less information toward full FC. Second, we developed custom templates to select time points with connectivity patterns that maximize behavioral relationships. This involved implementing multiple rounds of a simulated annealing optimization to develop templates that select time points with similar activity and improve brain–behavior correlations for each behavior. Third, we analyzed the temporal properties of the time points selected both across multiple runs of the optimization for each behavior and across behaviors to assess the locality or dispersion of behavioral information across the time series. Ultimately, co-fluctuation amplitudes displayed differences in brain–behavior correlations which varied based on the behavioral variable. Over multiple optimizations, a handful of behaviors revealed that selection of filtered time points improved correlation strength and similarity of behavioral associations across groups above full FC. Associated time points were distributed across co-fluctuation amplitudes, but most selected frames loaded onto higher amplitude moments. The few moments in time that overlapped across multiple behaviors displayed greater co-fluctuation amplitude and similarity to full FC. This suggests that full FC captures most brain–behavior relationships for some behaviors but may be improved for others if time points are selectively chosen based on the behavior of interest.

Analyses were done on data from the Human Connectome Project (HCP; Van Essen et al., 2013). From the 1,200 data release, 352 subjects (52% female, mean age = 28.44 ± 3.77, age range = 22–36 years) were selected based on low subject motion, high level of data quality/completeness, low incidence of artifact removal, and no family relationships (Cole et al., 2021). Participants were instructed to fixate on a cross during the scan session. Data were collected using a 3T Siemens Connectome Skyra with a 32-channel head coil. Four 15-min resting-state fMRI scans were collected from each participant over 2 days with a gradient-echo EPI sequence, TR = 720 ms, TE = 33.1 ms, flip angle = 52, 2-mm isotropic voxel resolution, and multiband factor = 8. After removal of frames from both the start and the end of the scan session, 1,100 time points remained for a scan duration of 14:33 min. Data acquisition and the minimal preprocessing pipeline are described in detail in Glasser et al. (2013).

A selection of 58 behavioral measures (Table 1) were analyzed from a larger set from the 1,200 HCP subject release (Barch et al., 2013; Van Essen et al., 2013) based on their varied distribution across cognitive, social, and emotional measures (Kong et al., 2019; Liégeois et al., 2019). The confounds of age, sex at birth, body mass index (BMI), framewise displacement (FD), and FreeSurfer intracranial volume were regressed separately for each cross-validation training and testing group. One-hundred cross-validations were created in 80/20 subject splits, where each set consisted of 282 randomly selected training subjects and results were analyzed in the remaining 70 testing subjects.

The HCP functional data were cleaned of signal artifacts with ICA-FIX (Salimi-Khorshidi et al., 2014). Nilearn signal.clean was used to linearly detrend, band-pass filter (0.008–0.08 Hz; Parkes et al., 2018), confound regress, and standardize all functional data. The functional data were then nuisance regressed using the global signal.

Regional time series were obtained using the Schaefer 200 functional parcellation as provided in the fs_LR surface space (Schaefer et al., 2018). Following nuisance regression, functional data were averaged within each regional node per frame to create 200 spatially distinct time series. Each of the 200 nodes was matched to a set of canonical functional networks as described by Yeo et al. (2011): visual (VIS), somatomotor (SOM), dorsal attention (DAN), ventral attention (VAN), limbic (LIM), frontoparietal (FP), and default mode (DMN).

In this study, functional connectivity (FC) is defined as the statistical dependencies between BOLD time series of every pair of brain regions. Pearson correlation is most commonly used to establish the distance between nodal time series and can be calculated as

where T defines the total number of time points $(t)$ and $zi$ represents the z-scored value of region $i$⁠. Computing this for every pair of regions produces an $[N×N]$ matrix with $N$ number of nodes. All analyses in FC and on connectivity components were vectorized to include $K=N(N−1)2$ unique edges.

FC can be deconstructed into a $[K×T]$ matrix of all edges $(K)$ by time $(T)$⁠. This method is referred to as edge time series (eTS) and produces co-fluctuation time series between each pair of nodes (Faskowitz et al., 2020). Essentially, eTS removes the final averaging step from the Pearson correlation, revealing the instantaneous relationships between the activity of node pairs as

where $rij(t)$ provides a measure of co-fluctuation of regions $i$ and $j$ at time $t$⁠. Averaging eTS across time returns standard functional connectivity (FC) and, therefore, can be used to isolate the exact contribution of temporal features that produce FC. eTS was computed for each subject to derive root sum square (RSS) amplitude values at each time point as

where each time $(t)$ is assigned a value related to strength of global co-fluctuation across all edges $(K)$⁠.

Removing amplitude from the time series and preserving whether a brain region’s activity is above or below the mean maintain nearly all information utilized in FC (Sporns et al., 2021). Accordingly, each time step comprises a bipartition of brain regions (or nodes) into respective groups displaying activity either above or below their own means. Selection of frames for analysis can then be based on corresponding framewise properties or similarity between a binarized pattern of activations across brain regions with the bipartitions at each moment in time. This provides both a more efficient means of data storage and a computationally expedient alternative to previous time-varying methods. Additionally, this allows for similar selection of time points for comparison across subjects in resting state.

For each subject and run, the z-scored BOLD time series was binarized by thresholding at zero (BOLD > 0 = 1 and BOLD < 0 = 0). This removed amplitude information while retaining the sign of the signal. An agreement matrix or co-classification matrix (Jeub et al., 2018; Lancichinetti & Fortunato, 2012) provides an account of the frequency in which two nodes are assigned to the same community. This is compared with a null model of the expected frequency of community identity given random permutations of nodes while keeping the number and size of clusters fixed. The null was computed as

where $l$ is the number of samples, $Ck$ the number of clusters in the k-th sample, $|Cks|$ is the number of nodes in a given cluster of $Ck$⁠, and $N$ is the total number of nodes. Here each time point $(t)$ was considered as a sample $(k)$ and each bipartition group as a cluster $(s)$⁠. The null was subtracted from the agreement matrix prior to further analysis. Computing an agreement matrix across the bipartitioned time series produces an $[N×N]$ matrix across the selected subset of frames that is nearly identical to standard FC (Sporns et al., 2021). Agreement components (AGc) were created by computing an agreement matrix from subsets of frames within a scan session.

Connectivity components allow for network-based analysis on non-continuous moments in time that reflect exact temporal contributions toward FC. Frames can be selectively filtered across time based on multiple criteria of interest (refer to Cutts et al., 2023). Two filtering methods were deployed: (1) Binning of frames based on a framewise metric and (2) similarity of frames to a preselected template.

Binning can be performed using any desired criteria of interest that temporally matches the functional time series. Frames can then be sorted into bins of an experimentally chosen size based on the amplitude of the framewise metric. Here, frames were sorted into deciles of RSS amplitude to directly test the appearance of brain–behavior relationships from frames which contribute less toward full FC (Esfahlani et al., 2020).

Template filtering can be performed by computing the similarity between a template and each frame of the time series. Here, distance was computed using the mutual information (MI) between a binarized template (vector of nodal length) and the bipartitioned nodal time series for each time point. MI was computed as

where $m$ and $m′$ represent modules within the larger partitions of $M$ and $M′$ under comparison. $P(m,m′)$ gives the fraction of nodes in $m$ and $m′$ that are present in both modules (Meilă, 2007). We normalized the measure to account for entropies of both $M$ and $M′$ and rescale the unit interval. The top 10% of frames with the highest normalized MI were selected for connectivity components.

MI was selected over other distance measures based on the observation that bipartitions can reconstruct full FC independent of nodal amplitude. A given bipartition and its exact inversion both provide identical information regarding the co-fluctuation of node pairs at a given moment in time. In this case, MI retrieves frames with consistent community labels of the bipartition rather than exact matches to activations. Past work suggests consistency in filtered bipartition frames using templates (Sporns et al., 2021). Further, the use of bipartitions irrespective of amplitude mirrors the strongest instances of similarities between nodes that are the most informative in standard FC analyses. Further information on template filtering can be found (Cutts et al., 2023; Sporns et al., 2021).

Irrespective of the filtering or connectivity component approach utilized, results were assessed for cortex-wide behavioral associations separately for each behavior. Spearman correlation was utilized to relate each component edge across subjects to their corresponding regressed behavioral measures. Information was then aggregated across brain–behavior correlation maps by taking the mean absolute value of all edges.

The mean correlational value was used as the cost function in the template optimizations to reduce the complexity of the full brain–behavior matrix into a single value for performance assessment (refer to Section 2.8). This metric was chosen to be more generalizable so that all but the behavioral measures remained constant throughout the optimizations. Greater generalizability was determined based on limited prior assumptions regarding the sign of correlational coefficients or number of expected behaviorally associated edges. There are many other candidate cost functions that could be used to assess improvements in brain–behavior associations. Initial pilot optimizations using metrics such as number of significant correlations (e.g., p < 0.001, p < 0.0001, p < 0.05 Bonferroni corrected) and log of the top percentile of edge p-values performed similarly to the simplified correlational mean implemented in this study.

The mean brain–behavior correlations found for each RSS decile were clustered into groups of behaviors with similar profiles. A distance metric (here Spearman correlation) was calculated across all pairs of behaviors using their RSS correlation profiles to produce a [behavior x behavior] matrix that was then clustered. Briefly, modularity maximization is a common approach used to compute data-driven communities defined by denser internal connections than would be expected by chance and weaker connections across communities (Esfahlani et al., 2021; Newman & Girvan, 2004). The modularity quality function is computed as

where $Srs$ is the similarity between behavior profiles $r$ and $s$⁠, $Prs$ is the expected similarity by chance, and $δ(σr,σs)$ is the Kronecker delta function that is equal to 1 when community assignments of $r$ and $s$⁠, denoted as $σr$ and $σs$⁠, are identical and 0 otherwise. The resolution parameter $γ$ controls the size of communities that modularity detects. The Louvain algorithm (Blondel et al., 2008) was used to maximize the value of Q, as a representation of the quality of community partition, for each entered value of the resolution parameter $γ$⁠. Here, we used the Potts null model (Traag et al., 2011) and the full FC matrix, including negative edges. We conduct a general “sweep” over the $γ$ parameter by sampling at 1,000 different values, followed by a finer sampling within the range where two-community structure appears (10,000 samples). This creates a set of multiresolution communities (Jeub et al., 2018) which are then transformed into a probabilistic agreement matrix showing the fraction of instances in which nodes $r$ and $s$ were assigned to the same communities.

An optimization algorithm was used to create custom templates for each behavior that would isolate framesets which maximized brain–behavior correlations across functional connectivity edges. Due to the similarity between the agreement matrix of binarized BOLD time series and full amplitude FC, both binary templates and AGc were utilized for computational efficiency. The algorithm was initiated with randomized binary vectors of length $N$ (number of nodes). Normalized mutual information was used to find the similarity between these initial templates and the binarized nodal activity of each time point. The top percentage of best matching frames (in this case 10%) was used to produce AGc for all subjects and runs. Brain–behavior correlations were found for each edge of the average AGc across subjects. In each new iteration, the performance of a template was compared with the template of the previous iteration in order to maximize brain–behavior correlations. If the new template outperformed the previous, it was kept and modified for the next iteration. This process was repeated over 10,000 iterations to allow the algorithm to converge on a stable end state.

With this method, there are 2^N template permutations which is too large for exhaustive analysis. Therefore, an implementation of the Metropolis-Hastings algorithm (Metropolis et al., 1953) was used to efficiently probe the search space. An objective function D was defined as $mean(abs(BB))$⁠, where $BB$ represents the brain–behavior correlations for all edges. Improvement upon this objective function was used in conjunction with simulated annealing (Kirkpatrick et al., 1983) to determine the acceptance or rejection of each template iteration. Upon each iteration, a temporary template was created by randomly inverting the binary identities of a handful of nodes from the previous template. The number of altered nodes followed a normal distribution where 1, 2, or 3 nodes were changed with a frequency of 0.68, 0.27, and 0.04, respectively. The success of the temporary template was governed by either the improvement of the objective function or the simulated “temperature” of the system. A cooling schedule was defined as $Temp=T0Texph$⁠, where $T0$ defines the initial temperature, $Texp$ the steepness of the gradient decay, and $h$ the iteration number. An annealing criterion of $e−ΔD/Temp>R(0,1)$ was used to determine acceptance of the current template, where $Temp$ refers to the simulated “temperature” governed by the cooling schedule and $R(0,1)$ is a random number between [0,1]. This annealing criterion occasionally rejects local optima in favor of continued exploration earlier in the optimization process. As the “temperature” decays based on the cooling schedule, the algorithm increasingly accepts optimal solutions as governed by the objective function. One set of parameters was selected for all behaviors to ensure stable convergence across repeated runs of the optimization within a reasonable time.

Optimizations were performed across 100 cross-validations of 80/20 (train/test) groups for each of the 58 behaviors to evaluate consistency across template outputs. The selected framesets were analyzed for temporal consistency across 7 cross-validations for each subject (based on the minimum from the randomized creation of groups). Selected framesets were additionally compared with framewise properties of the time series such as RSS, similarity to full FC, and framewise displacement (FD). Generalizability was tested by applying the optimized templates to the held-out testing subjects for each cross-validation and analyzing the similarity between brain–behavior correlation maps of the training and testing sets.

The temporal location of best matching frames for each optimized behavioral template was recorded for each subject and run. This information was used to compare the consistency of iterations of the optimization as well as determine the overlap between time points retrieved from separate behaviors. The consistency across iterations of the optimization was compared with circshifted nulls (using the MATLAB circshift function) of the selected frame indices for each behavior.

Each iteration of the optimization returned 110 best matching frames regardless of overall fit. Plotted as a binary raster, successful convergence of the optimization is expected to display consistent time frames selected per optimization run. The occurrence of these frame selections was summed across the test group results from multiple cross-validations at each moment in time and the frequency of overlap was then assessed for significance against the same number of randomly shifted framesets. The null model was developed by taking the true framesets and randomly shifting them either forward or backward in time separately for each cross-validation. Frames shifted beyond the boundaries were wrapped around to the opposite side of the time course. This procedure preserved the autocorrelation and temporal spacing between frames while varying the amount of overlap between cross-validations. Ten thousand iterations of the null were performed to establish significance in overlap across optimization runs.

The significant time points established from this method were recorded for every behavior, subject, and run. Overlap across behaviors was additionally established by comparing significant moments in time across behaviors separately for each subject and run to 10,000 rounds of circshifted frames. The number of instances of behavioral overlap was tested against the null across all subjects to establish the degree to which behaviors load onto similar or disparate time points. Considering that RSS has been found to drive the signal present in full FC (Betzel et al., 2022; Esfahlani et al., 2020), behavioral overlap was also compared with co-fluctuation amplitude.

Here, we used a filtering approach to discover how multiple behavioral relationships change across co-fluctuation amplitude and whether brain–behavior correlations could be improved through targeted selection of time points. We analyzed 352 subjects from the Human Connectome Project (HCP; Van Essen et al., 2013), selected from the “1200 Subjects” HCP release based on low subject motion and no family relationships (Cole et al., 2021). Four resting-state scans were analyzed from each subject. Fifty-eight different behaviors were selected based on their distribution across cognitive, social, and emotional measures (Table 1; Kong et al., 2019; Liégeois et al., 2019).

Connectivity components similar in structure to standard FC were reconstructed using selections of frame subsets. When assessing the full time series, computed co-fluctuations between pairs of brain regions in eTS (Fig. 1A, top left; Faskowitz et al., 2020) can be averaged across time to exactly reconstruct FC. Similarly, FC can be nearly reconstructed by binarizing the z-scored time series (“bipartitions”; Sporns et al., 2021) for each brain region and then computing the agreement of partition identity between pairs of brain regions across all bipartitioned frames (Fig. 1A, bottom left). Accordingly, disparate frames can be selected based on a property of interest and reconstructed into a connectivity component either through averaging across selected eTS frames or through the agreement of bipartitions. In the case of the highest co-fluctuation time points, this reconstruction from only 10% of the time series greatly resembles full FC (Fig. 1A, right). When this handful of high co-fluctuation amplitude frames is selected across subjects, the corresponding brain–behavior correlations across edges greatly resemble behavioral relationships from full FC (Fig. 1B). Similarity between brain–behavior relationships of full FC and components derived from co-fluctuation amplitude likewise systematically decrease with amplitude across all behaviors, suggesting that lower amplitude frames reveal differences in behavioral relationships (Fig. 1C). These differences in RSS bins from full FC can be seen in their varying associations with behavioral measures, as shown for both high- and low-amplitude frames (Fig. 1C, Inset; Supplementary Fig. S1). Components created from the agreement across bipartitions (AGc) were used throughout the following analyses for computational efficiency and direct comparison across binning and optimization results (Chumin et al., 2024; Cutts et al., 2023; Sporns et al., 2021).

First, behaviors were analyzed across co-fluctuation amplitudes to assess whether lower amplitude frames with less contribution toward full FC displayed unique behavioral associations. Co-fluctuation amplitude varies across time and is computed as the root sum square (RSS) across co-fluctuations between all pairs of brain regions at each time step (refer to Section 2). We created RSS profiles for each behavior to provide a general indication of which amplitude frames of the time series express stronger brain–behavior correlations. Clustering was then used to assess which behaviors loaded onto similar time points within a scan, based on their “preference” for different parts of the RSS amplitude spectrum. We sorted the frames into decile bins based on the RSS strength to assess their relationship with behavioral correlations for each bin. AGcs were created for each RSS decile bin for every subject and averaged across runs. Correlations were computed between each of the 58 behavioral measures and all AGc edges. Strength of brain–behavior correlation was assessed using the mean absolute correlation across all edges. Each behavior consisted of a mean behavioral correlation score for each RSS decile and was compared with 100 AGc made from randomly selected frames. This provided a general mapping of a behavior’s loading onto RSS amplitude (Fig. 2), which were found to cluster behaviors into five groups related to RSS profiles (Fig. 2C; refer to Section 2). These RSS clusters reveal that there are separate groups of behaviors for which there are stronger brain–behavior correlations at different time points. Additionally, a behavior’s loading onto RSS (Fig. 2A) relates to the strength of behavioral correlations captured in full FC. Here, behaviors that mainly load onto higher amplitude frames (found in the first and second clusters) have a slight increased presence in full FC (Fig. 2B) than behaviors that load onto intermediate- to low-amplitude frames (p = 0.0227, independent sample t-test). This suggests that behaviors that load onto lower amplitude frames may be partially obscured in analyses utilizing full FC.

Second, templates were created to filter time points that improved behavioral relationships across multiple optimizations using simulated annealing. Templates consisted of a binary nodal vector that was used to filter frames based on the similarity (mutual information) between the bipartition of nodes at each time step and the template. This approach allows for the analysis of similar frames across subjects irrespective of timing during resting state. The structure of the binary nodal templates was optimized to provide the strongest behavioral correlations from the top percentage of filtered frames. Here, 10% of similar frames were selected to reduce the effects of artifactual or spurious frames on the AGc. One-hundred optimizations were performed separately for the 58 behaviors where each iteration created a filtering template from 80% of subjects (282 training set) and analyzed results from the corresponding template on the remaining 20% of individuals (70 testing set). A new randomized training and testing set was used for each of the 100 cross-validations. Each optimization was initialized with a random template that was used to filter the top 10% of similar frames from all scans of each subject (Fig. 3). Filtered frames were aggregated into an AGc for each subject across runs and used to create an edgewise brain–behavior correlation map. The strength of behavior correlations from this map was used as the cost function to create filtering templates across multiple optimizations. Each step of the optimization accepts or rejects the current template based on the strength of behavior correlations. Accepted templates are permuted at randomly selected nodes and reintroduced to start a new cycle of the optimization. The end result was an optimized template that selects moments in time for which the cost function is high. This procedure was performed separately to produce unique optimized templates for each of the 58 different behaviors.

The optimization results from a high performing example behavior (vocabulary pronunciation) are shown in Figure 4 (for an additional example behavior, see Supplementary Fig. S2). Figure 4A shows the overall strength of behavior correlations across each iteration for 100 cross-validations of the training subjects. As expected, the optimization succeeded in constructing AGcs with much higher behavioral correlations compared with full FC. Figure 4B shows that the optimized behavior templates filter greater behavioral correlations compared with the initial template as well as full FC for the training subjects. Filters applied to held-out test subjects maintained a similar average of the aggregate measure of brain–behavior correlation (Fig. 4B, bottom), but did show a significant increase in correlation strength from each cross-validation over full FC for this given behavior (p < 10^-3, paired-sample t-test; effect size, d = 0.31). The mean behavioral correlations were reported given their use as the cost function in the optimizations. However, other features of the correlation maps were modified upon optimization that were not readily captured through the average. Additional analyses of the edge correlation distributions, such as kurtosis, were also compared across conditions to better determine how behavioral correlations differed after optimization (Supplementary Figs. S3 and S4). These results suggest that a greater number of edges revealed stronger behavioral correlation following optimization, despite small differences in the averaged measure.

Notably, the brain–behavior correlation maps were also significantly more similar across training and testing subjects. Figure 4C shows the brain–behavior correlation map from the best performing cross-validation for full FC, the initial template, and optimized frames. The optimized results reveal a much stronger correspondence in behavioral relationships between training and testing groups (additional examples in Supplementary Figs. S5 and S6). This improved train–test transfer (based on Spearman correlation between brain–behavior maps) can be seen across multiple cross-validations and shows significant improvements in similarity across groups (ρ = 0.33) over full FC (ρ = 0.19; p < 10^-15, paired-sample t-test, effect size, d = 0.90; Fig. 4D, top). Positive and negative edges from the training set were compared separately with corresponding edges in the testing set. The mean correlation strength of positive edges (Fig. 4D, middle) and negative edges (Fig. 4D, bottom) was significantly greater than initial templates (p < 10^-15, paired-sample t-tests) and full FC (p < 10^-13, paired-sample t-tests). Similar results were found when the number of filtered frames varied (Supplementary Figs. S7 and S8). This shows that filtering frames using optimized templates allowed greater similarity of brain–behavior correlations across training and testing groups as well as improved mean correlation strength for positive and negative edges.

The optimizations produce consistent template structures that provide additional information regarding which nodes and functional systems are selected for specific behaviors. Figure 4E displays the optimized templates for each of the cross-validations ordered based on functional systems. Here, aggregating bipartition community identity across all nodes provides the dominant nodal affiliation for each system (Fig. 4E, right). Figure 4F displays the same information as in Figure 4E but as the frequency of nodal affiliation with each bipartition community across cross-validations. Interpretation of these templates is dependent on relating bipartitions back to their relationships with functional properties. For instance, all nodes in a template are selected for when using bipartition templates, but some nodes are consistently assigned in alignment or counter to other nodes. These results reveal that there is consistency within the template structures retrieved from these optimizations that can be related back to nodal and network-level features.

Figure 5 provides a summary of the results across all analyzed behaviors. Figure 5A shows the average mean behavioral correlation from all cross-validations improved for each behavior in the training set. Similar levels of averaged brain–behavior correlations were found across all conditions of the testing set for most behaviors. Figure 5B shows the detailed results as displayed in Figure 4, but for all behaviors. Significance values and effect sizes for comparisons between optimized and full FC as well as optimized and initial frames are listed in Supplementary Tables S1 and S2 for all behaviors. Here, all behaviors with significantly greater train–test transfer of behavioral correlations compared with full FC (p < 0.01, one-tailed t-test uncorrected) also maintained significantly greater strength of positive and negative edges in the testing set. Similar results for Bonferroni corrected as well as example AGcs and brain–behavior correlation maps across multiple behaviors are displayed in Supplementary Figures S9 and S10. A handful of behaviors also showed both greater transfer and correlation strength or only greater correlation strength than full FC. Additionally, the optimized templates were not interchangeable across behaviors. Optimized templates filtered AGc with higher mean behavioral correlations within behavior than across behaviors (Supplementary Fig. S11) and the templates were notably more similar across cross-validations within behaviors than between behaviors (Supplementary Fig. S12). This suggests that the optimizations were locking onto different features of the time series at separate time points that were not reducible to a single best pattern to describe multiple behaviors. Instead, optimizing individual behaviors repeatedly selected for template patterns specific to those behaviors.

The filtered frames from the optimization were further matched to their corresponding RSS decile bin (shown in Fig. 2). Figure 5C relates the count of optimized frames that belong in each RSS bin to the average binning identities of 100 randomly selected bins. Across all cross-validations of each behavior, optimized templates select for frames with higher RSS and similarity to FC than the initial template without selecting for increased motion artifact based on framewise displacement (Fig. 5D). However, compared with the highest RSS decile bin, optimized frames maintain lower RSS (p = 0, independent sample t-test) and lower similarity to FC (p < 10^-15, independent sample t-test), suggesting that filtering recruits frames from lower co-fluctuation amplitudes to improve behavioral correlations. These improvements in behavioral correlations and train–test transfer relate to increased variability between subjects. Filtering time points using behavioral templates increased the differentiation of component patterns from full FC (Supplementary Figs. S13 and S14) and subject variability (Supplementary Fig. S15).

Further, the down sampled functional system structures of the templates were computed across all behaviors (Fig. 5E) as shown in Figure 4E and ordered based on the RSS clusters derived in Figure 2. Most behaviors produced templates with bipartitions separately grouping VIS, SOM, DAN, VAN and LIM, FRP, DMN. This pattern was consistent across behaviors with reliable deviations of select functional systems and largely independent of RSS classification. The similarities in the down sampled network templates likely result from the higher amplitude frames selected for most behaviors. Differences across behaviors and improvements over full FC, however, appear with the addition of lower amplitude frames with differing patterns of co-fluctuations between networks. Further assessment of node co-assignment to the same bipartition communities within optimized templates further shows that most nodes are consistently paired across functional systems rather than within systems (Supplementary Fig. S16). There is little overlap across behavior templates regarding which node pairings are consistently present within the same bipartition community (Supplementary Fig. S17), suggesting that the optimizations are selecting for varied template structures across behaviors. The one notable source of overlap between behaviors, however, is the assignment of nodes from the frontoparietal system to the same bipartition community. This functional network has been associated with higher subject variability and fingerprinting capability (Amico & Goñi, 2018; Finn et al., 2015; Peña-Gómez et al., 2018), supporting the idea that these optimized templates are selecting for moments from the scan with higher intersubject variability related to behavioral measures.

Third, the temporal properties of time points selected from behavior templates were assessed and compared across behavioral measures. So far, the results have been aggregated over the selected frames of all subjects and runs as well as multiple cross-validations, but the optimization pulls frames from specific time points whose properties can be analyzed directly. Out of the random 80/20 (train/test) splits of the 100 cross-validations, subjects were selected multiple times across both groups. Over multiple cross-validations of an individual subject, the optimizations consistently targeted a select set of specific time frames. Figure 6A (top) shows an example heatmap for the frequency of frame selection for the example behavior (vocabulary pronunciation) across seven cross-validations of testing set results for a select subject and run. Each column of the matrix depicted below (Fig. 6A, middle) consists of the corresponding heatmaps of each subject and run. Heatmaps of all subjects and runs were compared with 10,000 circular shifted nulls to determine a significance score for each frame (refer to Section 2). Frames were then selected based on a threshold of p < 0.01 (uncorrected) for significantly greater temporal consistency across optimizations than chance. A stricter threshold of p < 0.001 provides similar results (shown in Supplementary Fig. S18 and a separate behavior in Supplementary Fig. S19). The number of frames found to be significantly consistent among cross-validations was much higher for the optimized templates (Fig. 6B; mean = 29.02) than for the initial template used at the start of the optimization (Fig. 6A; mean = 3.30; p = 0, paired-sample t-test). This suggests that the optimized templates allowed for the consistent selection of time points in held-out test data despite multiple independent initializations of the optimization with different groupings of subjects.

These significant frames were then compared across behaviors. In Figure 6C, the significant frames are shown for each behavior of an example subject and run. From that, selected frames across behaviors appear to largely occur independently of one another but a handful of frames reveal overlap across behaviors. Figure 6D shows the behavioral overlap of the example subject and run shown in Figure 6C and compares the frame selection with the RSS amplitude time series. Across all subjects and runs, there were time points with greater behavioral overlap compared with 10,000 circular shifted nulls (Fig. 6E). The time points with significant behavioral overlap (p < 0.01, uncorrected) were averaged across all runs of each subject and general framewise properties were compared with results from initial and all frames. In Figure 6F, moments in time with significant behavioral overlap display increased RSS amplitude, higher similarity to full FC, and nonsignificant differences in head motion compared with initial templates and all frames.

Functional connectivity is becoming increasingly used toward the development of brain-based biomarkers of healthy cognition, demographics, and clinical diagnosis (Chen et al., 2021, 2022; Chumin et al., 2024; Dhamala et al., 2021; Dizaji et al., 2021; Dubois et al., 2018; Lake et al., 2019; Parkes et al., 2020; Sripada et al., 2020; Uddin et al., 2013; Wu et al., 2023). There are already several machine learning and predictive models that have been applied toward improving behavioral prediction and developing a deeper understanding of what networks can tell us about cognition and clinical phenotype (Cai et al., 2021; Cui & Gong, 2018; Finn et al., 2015; He et al., 2020, 2022; Kawahara et al., 2017; Kong et al., 2019, 2021; Ooi et al., 2022; Pervaiz et al., 2020; Sasse et al., 2023; Shen et al., 2017; Wehrheim et al., 2023). A large portion of FC research has been done on fMRI data, so methodological constraints and the slow hemodynamic response of the BOLD signal (Buxton, 2013; Lake et al., 2020; Logothetis, 2002; Schwalm et al., 2017) have led to the preference of conducting analyses on full scan sessions of 20–40 min of data per subject to improve reliability (Birn et al., 2013; Elliott et al., 2019; Laumann et al., 2015; Noble et al., 2017). Recent studies, however, suggest that additional behavioral information can be obtained from the analysis of temporal fluctuations (Cutts et al., 2023; Eichenbaum et al., 2021; Fong et al., 2019; Liégeois et al., 2019; Lurie et al., 2020; O’Connor et al., 2022; Peña-Gómez et al., 2018; Van De Ville et al., 2021). There are already several approaches within the time-varying literature (Allen et al., 2014; Betzel et al., 2016; Calhoun & Adali, 2016; Chang & Glover, 2010; Cohen, 2018; Eichenbaum et al., 2021; Leonardi & Van De Ville, 2015; Liégeois et al., 2019; Lindquist et al., 2014; Liu et al., 2018; Lurie et al., 2020; Pedersen et al., 2018; Preti et al., 2017; Sakoğlu et al., 2010; Vidaurre et al., 2017; Yaesoubi et al., 2018; Zalesky et al., 2014), but there are few answers regarding whether methods can be applied evenly across specific behaviors or phenotypes of interest. Here, we assess whether behavioral relationships across a suite of behaviors are visible in full FC or whether there are temporal fluctuations within a scan with stronger relationship with behavioral variance. This was done by analyzing whether brain–behavior associations change based on co-fluctuation amplitude and by filtering frames based on specific patterns of activity (bipartitions). We used optimizations to develop bipartition templates to isolate time points consistently across subjects with stronger behavioral associations and determine whether separate moments are better for different behaviors. We also show that the time points filtered for these behaviors are consistent across multiple optimizations and create transferable brain–behavior associations in held-out subjects. We found that a few moments in time improve behavioral associations across multiple behaviors, but that some behavioral information loads onto moments that are less prominent in full FC.

Our approach provides a framework to select behaviorally relevant moments in time from a typical scan session consistently across individuals. These results suggest that moments in time carry different signatures of information that are relevant toward specific brain–behavior associations. We found that there are some moments that are more informationally rich than others and that behavioral relationships differ based on co-fluctuation amplitude as well as ability to be detected from full FC. Moments that enhance behavioral associations differ across behaviors, suggesting that filtering can enhance connectivity patterns based on the phenotype or behavior of interest. This is done with resting state where normally large amounts of data are needed to ensure the researcher is analyzing stable FC across individuals. Rather than analyzing summaries over longer periods of time, finding moments that are similar across subjects in resting state could fulfill this need with less data and clearer results with the added benefit of improved temporal localization of behaviorally related signatures. Recent work suggests that thousands of subjects are needed to establish strong and reliable brain–behavior associations (DeYoung et al., 2022; Liu et al., 2023; Marek et al., 2022; Rosenberg & Finn, 2022). However, filtering time points based on specific behaviors or phenotypes provides an additional avenue for exploring intersubject variability in behavior. Longer scans would still be required to identify enough time points that better relate to specific behaviors, but by sorting frames based on their behavioral relevance, more nuanced information may be achieved over otherwise static full FC.

Across multiple filtering approaches, this study shows that behaviors often load onto similar time points as groups. For instance, clusters of behaviors exhibit similar associations within frames of RSS amplitude ranges (Fig. 2). Additionally, using optimizations, we find that most temporal overlap across behaviors occurs at high amplitude moments (Fig. 6). Although there is a dispersion of behaviorally relevant information in the time series that differs across behaviors, it remains unclear what causes instances of overlap between them. Here, we highlight that behaviors that prefer higher amplitude frames are more likely to be represented in full FC. It remains unclear, however, why behaviors load onto specific amplitude moments and how those that share similar “strategies” relate to each other behaviorally. Future work focusing on specific behaviors would further reveal which aspects of connectivity are highlighted from filtered frames, whether these features are meaningfully related across behaviors, and why certain behaviors load onto lower amplitude moments with less contribution toward the time-averaged full FC.

Linking behavioral and phenotypic associations in FC back to mechanistic or neurobiologically plausible sources will likely require improvements in temporal specificity. Along with fluctuations in time-varying FC (Allen et al., 2014; Betzel et al., 2016; Calhoun & Adali, 2016; Chang & Glover, 2010; Cohen, 2018; Eichenbaum et al., 2021; Leonardi & Van De Ville, 2015; Liégeois et al., 2019; Liu et al., 2018; Lurie et al., 2020; Pedersen et al., 2018; Preti et al., 2017; Sakoğlu et al., 2010; Vidaurre et al., 2017; Yaesoubi et al., 2018; Zalesky et al., 2014), individually distinctive signatures also vary over time (Cutts et al., 2023; Peña-Gómez et al., 2018; Sasse et al., 2023; Van De Ville et al., 2021). As such, the relationship between intersubject variability in FC and behavioral measures will change based on the connectivity pattern compared across subjects. With enough data, FC stabilizes within subjects (Gratton et al., 2020; Laumann et al., 2015; Noble et al., 2017) to provide a clear baseline for intersubject comparisons. However, moments in time that emphasize connections between specific brain regions or functional systems can also provide a scaffold for comparisons that can be related back to intersubject variability of specific behaviors. The origin of these temporal fluctuations in FC remains unclear, with some studies suggesting an origin in neuronal processes (Matsui et al., 2019; Thompson, 2018) or network dynamics (Heitmann & Breakspear, 2018; John et al., 2022), while others interpret FC states as expected stochastic variations around a central tendency that can be observed at the full scan length (Ladwig et al., 2022; Laumann et al., 2017).

Regardless of the origin, finding the moments in time that contribute most toward behavioral associations allows for closer inspection of the relationship with other neurobiological signatures or noise. For instance, out of the 15 behaviors with significant improvements over full FC, 7 were classified as task performance, 4 as self-report, and 4 as unclassified (based on Liégeois et al., 2019). This suggests a potential dichotomy where behaviors perform better either from optimized frame subsets or from full FC, and this could be due to the trait-like nature of the behavior in question or the duration it is expected to occur. Additionally, these improvements in behaviors from temporal filtering might be related to associations with physiological signatures. Although behaviorally filtered frames from this study did not show systematic differences attributable toward noise (such as head motion; Figs. 5D and 6F), and potential confounding variables such as age, sex, and BMI were regressed out, the time points selected for improvements in behavioral correlations could be related to additional physiological or behaviorally relevant time stamps. Some of these improvements in behavioral associations may be driven by variation in physiological signatures typically considered as artifact, but it is important to determine both where such improvements in behavioral associations arise as well as the relationship between potential sources of artifact and their correlations with behavioral or phenotypic traits (He et al., 2022; Uddin, 2020). Additionally, future work will involve optimizing behavioral relationships with movie and task data where the exact stimuli that drive behavior relationships can be assessed through comparison with the contents of the movie or task of interest.

Determining the temporal origin of associations between phenotypic measures and functional networks could provide greater insight into how the brain works and establish more efficient methods for processing neuroimaging data. Aging or clinical cohorts is one area where applying targeted filtering methods could enhance the quality of the data under investigation. Such cohorts are often prone to discomfort or increased movement while in the scanner (Satterthwaite et al., 2019). Using a template filter to recover frames which amplify intersubject differences in a behavioral or phenotypic measure of interest could improve data quality (if such moments are not themselves related to artifact that commonly impacts time-averaged FC). Future work, for instance, assessing differences in template filtered frames might reveal age-related differences not present in standard FC. The widespread application of these methods, however, will require development of phenotypic or clinical-specific templates from representative cohorts and through utilizing domain expertise.

This study provides proof of concept for the selective filtering approach and shows that behaviors load onto disparate time points scattered throughout the scan session. The templates created here were optimized on 58 behaviors selected to represent a larger set (Barch et al., 2013; Van Essen et al., 2013). Although this provides an overview of different characteristics of behavioral measures in FC, templates for each of these behaviors were produced from the same optimization parameters. Although behaviors were assessed equally for balanced comparison, template structure is likely influenced by this assumption. Each of these behaviors would benefit from specialized attention to improve targeted frame selection. For instance, the cost function used in this study aims to improve average brain–behavior correlations across the cortex. Certain behaviors would not be expected to involve cortex-wide associations. Therefore, we suspect this choice in cost function selects high-amplitude moments with greater potential range of edge amplitudes across the cortex to capture behavioral variation. These higher amplitude moments further consist of bipartitions with characteristic decoupling between intrinsic versus task-positive systems (Sporns et al., 2021). The notable differences in behaviors that exceed full FC, therefore, likely relate to the consistent deviations away from this template structure. Further work incorporating spatial considerations into template development could significantly improve the method’s ability to isolate behaviorally relevant moments. We also know that behaviors occur at varying timescales, therefore, optimization procedures would likely benefit from considering the number of selected frames, sequences of frames, or the use of specific nodes or functional networks. Here, analyzed components were created from 10% of the time series and it remains unclear whether behaviors would benefit from varying the quantity of considered frames (results displayed for an example behavior using 5% (Supplementary Fig. S7) and 15% frames (Supplementary Fig. S8)). Ultimately, assessment by domain experts in specific behaviors or phenotypes could further clarify the relationship between the optimized templates and their measures through increased incorporation of expected spatial and temporal features characteristic of the targets of interest.

The structure of the templates themselves is also informative of the features filtered for each behavior and should be investigated further in future work. Here, we show that the templates retrieve consistent patterns within behaviors that are not reducible to a global pattern that best describes all behaviors (Supplementary Figs. S11 and S12). This approach reveals the benefits of optimizing behaviors separately compared with applying all analyses uniformly to a single set of time points. The specificity of behavioral templates suggests potential biological relationships between the nodal structure of the templates and the behavior they were optimized on. This study found that optimized templates display different bipartitions across functional systems (Fig. 5E), but it remains unclear what extent of the template is relevant toward the behavior of interest. Further work is also needed to determine how bipartitions correspond to the literature on specific behaviors. Since bipartitions represent nodes that actively coordinate within but not across communities, exact overlap between templates and known behaviorally relevant activations have yet to be assessed at a time-resolved level. Currently, this makes direct interpretation of the templates challenging regardless of their utility in retrieving desired framesets. That said, future work assessing multiple behaviors simultaneously could lock onto separate time points best related to behavioral constructs. For instance, canonical correlation analysis (CCA) could be utilized to discover time points across multiple brain–behavior relationships simultaneously and factor analysis could be used to derive optimizable scores across behavioral measures. There are many remaining avenues to explore regarding how behavioral templates relate to nodal activations and why specific time points within the scan amplify behavioral associations.

Lastly, the optimization occasionally produces an exact inversion of the brain–behavior map found in most other cross-validations (as seen in the cross-validations with worse train–test transfer in Fig. 4D, Supplementary Fig. S6). This likely occurs when the randomly selected subjects in the testing set involve outliers of the overall trend. A similar issue with the optimization is that for most behaviors a handful of cross-validations do no better than the initialized template with low train–test transfer as well as clear lack of functionally relevant structure in edgewise correlations of the brain–behavior map. However, it should be noted that the same cross-validations that experience inversions or fail to optimize clear edgewise correlations also mirror failures or inversions in full FC (Supplementary Fig. S6), likely relating issues back to differences in the groups of subjects analyzed. Overall, successful optimizations repeatedly return the same brain–behavior correlation map as well as similarities in filtered time indices across multiple iterations. Additional work is necessary to determine how these templates lock onto and improve behavioral associations.

In summary, we show that behavioral relationships vary across time and do not always occur at the same moments for different behaviors. These behavioral relationships map onto specific activity (bipartition) templates that can be used to filter similar functional connectivity patterns across subjects even during resting state. These optimized behavioral moments occur across RSS amplitudes and, therefore, may benefit from filtering approaches to emphasize brain–behavior relationships not as visible in full FC. Future applications abound, including studies probing for individual differences that involve developmental or clinical cohorts, and that leverage emerging machine learning and A.I. technology.

Neuroimaging and behavioral data from the Human Connectome Project are available after signing a data use agreement at the link: https://db.humanconnectome.org/. Code for computing edge time series and bipartitions can be accessed at https://www.brainnetworkslab.com/coderesources and additional graph theoretic metrics are found at https://sites.google.com/site/bctnet/. Code used to perform frame filtering, create connectivity components, and produce optimized templates can be found at https://github.com/sarahacutts/temporal-behavior-optimization.

S.A.C. and O.S. conceptualized and developed methodology; S.A.C. performed data analyses; O.S. supervised project; O.S., E.J.C., and R.F.B. reviewed and edited the manuscript; and S.A.C. wrote the original draft of the paper.

The Human Connectome Project (Van Essen et al., 2013), WU-Minn Consortium, collected and made available all participant data used for the above analyses. All procedures and study protocols were approved by the Washington University Institutional Review Board and informed consent was obtained from all participants.

This investigation was supported by the National Institute of Health, T32 Grant No. HD 07475 (S.A.C.). This research was supported in part by NSF Grant No. 2023985 (R.F.B. and O.S.), the National Institute on Aging Training Grant on Alzheimer’s Disease and ADRD (AG071444; E.J.C), and the Indiana University Network Science Institute. This research was supported, in part, by the Lilly Endowment, through its support for the Indiana University Pervasive Technology Institute and, in part, by the Indiana METACyt Initiative (Stewart et al., 2017). The Indiana METACyt Initiative at Indiana University was also supported, in part, by the Lilly Endowment. Data were provided, in part, by the Human Connectome Project, WU-Minn Consortium (principal investigators: D. Van Essen and K. Ugurbil; 1U54MH091657), funded by the 16 National Institutes of Health (NIH) institutes and centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University.

The authors declare no conflict of interest.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00443.