Interindividual Variability in Memory Performance Is Related to Corticothalamic Networks during Memory Encoding and Retrieval

RUB

We enrolled 36 healthy adults for the high-resolution multiscan fMRI study approved by the Faculty of Medicine Ethics Committee of the Ruhr University Bochum. Following the Declaration of Helsinki, all participants gave written informed consent before participating. All participants had normal or corrected-to-normal vision. Vision correction during MRI was achieved via contact lenses or MRI-compatible glasses. No individual had a history of drug or alcohol abuse within the last 6 months, head trauma with loss of consciousness, or any clinically relevant medical condition. For each participant, handedness was assessed through the Edinburgh handedness inventory (Oldfield, 1971). One ambidextrous participant was excluded to reduce heterogeneity for laterality (Ferrucci et al., 2013), five did not complete the protocol due to technical failure during data collection, and one more was excluded because of data artifacts and excessive head motion during the acquisition. Excessive motion was assessed using framewise displacement (FD), i.e., the displacement of one frame relative to the previous, with subjects excluded if mean FD > 0.5 within a given scan (Power et al., 2012, 2014). A total of 29 participants were included in the final sample for the resting-state analysis (mean age ± SD = 26 ± 3 years; range, 19–33 years; gender ratio, m:f = 10:19; mean handedness ± SD = 0.89 ± 0.11). We further excluded two participants who did not complete the task fMRI, including 27 participants in the sample for the task activity analysis (mean age ± SD = 25 ± 3 years; range, 19–32 years; gender ratio, m:f = 8:19; mean handedness ± SD = 0.9 ± 0.12).

UNIBA

We recruited an independent sample of 148 healthy adults for the multiscan fMRI study approved by the institutional Ethics Committee of the University of Bari Aldo Moro. All participants gave written informed consent before participating and underwent the same screening for inclusion as RUB participants. Five participants were excluded because of data artifacts, four were excluded for technical failure during data collection, twenty-two participants were excluded for excessive motion during the resting-state fMRI acquisition, and nine more because of low performance at the associative memory task (see below, Behavioral analysis). A total of 108 participants were included in the final sample for the resting-state analysis [mean age ± SD = 28 ± 8 years; range, 18–61 years; gender ratio, m:f = 56:52; mean handedness ± SD = 0.63 ± 0.41; mean intellective quotient ± SD = 106.6 ± 9.73 measured through the Wechsler Adult Intelligence Scale—Revised (Wechsler, 1981); mean education ± SD = 17.31 ± 2.3 measured through the Hollingshead Four-Factor Index of Socioeconomic Status (Hollingshead and Frederick, 1964)]. We further excluded four participants for data artifacts and three more because of excessive motion in at least one of the task fMRI scans (mean FD > 0.5 mm). The final UNIBA sample for the task activity analysis included 101 participants (mean age ± SD = 27 ± 8 years; range, 18–59 years; gender ratio, m:f = 49:52; mean handedness ± SD = 0.64 ± 0.4; mean intellective quotient ± SD = 106.5 ± 9.71; mean education ± SD = 17.24 ± 2.25).

RUB

During each 8 min resting-state scan, participants were instructed to remain awake, with eyes open, and fixate on the crosshair in the middle of a black screen (Greicius et al., 2003; Damoiseaux et al., 2006). The event-related memory task was structured into four runs, each composed of a picture-pair association encoding phase followed by a delay and a retrieval phase (Fig. 2A). The whole scan for the associative memory task lasted approximately 30 min, with each run lasting 7 min and 30 s, including a total of 2 min and 20 s for the encoding, 1 min and 40 s for the delay, and 3 min for the retrieval phase as described in Passiatore et al. (2021).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

The associative memory task. Diagrams show the associative memory paradigms performed during fMRI in the (A) RUB and the (B) UNIBA samples.

During the encoding phase, 20 picture-pairs per run (80 picture-pairs in total) were presented simultaneously on the screen for 2,500 ms. Each picture was selected from the FRIDa dataset (Foroni et al., 2013), representing an object. Each pair consisted of one object on the left and one on the right, with their positions randomized across participants. Participants were instructed to study the association between the picture-pairs and indicate which object in each pair was larger or smaller in reality. Responses were made using two separate button pads (“left” or “right”). To ensure a balanced left/right orientation of the stimuli, participants identified the larger object in half of the trials and the smaller object in the other half. This approach facilitated equal attention to both objects in the pair while maintaining consistency in task demands. Interleaved with the encoding trials, we included control stimuli consisting of 6 pairs per run (24 picture-pairs in total) of two-digit numbers located at the center of visually scrambled picture-pairs with a duration of 2,500 ms. To match the motor demands of “larger/smaller” trials, participants performed an “even/odd” judgment responding with two separate button pads (“left” and “right” response). At the beginning of each run, instructions randomly indicated whether the participant should identify the position of the “even” or “odd” number. These control trials served to reset the blood oxygen level-dependent (BOLD) signal in the MTL to baseline while enhancing task-related activity (Stark and Squire, 2001).

During the delay, 10 novel picture-pairs per run (40 picture-pairs in total) were presented for 2,500 ms while participants performed the same “larger/smaller” judgment task but were explicitly instructed not to memorize the associations. The picture-pairs presented during the delay were used to interfere with the previously studied pictures (Marini et al., 2017). Interleaved with the delay trials, we included control stimuli consisting of 3 pairs per run (12 pairs in total) of two-digit numbers located at the center of visually scrambled picture-pairs with a duration of 2,500 ms. As for the encoding phase, participants were instructed to identify the position of the “even” or “odd” number.

During retrieval, participants performed an “old/new” object recognition involving 30 single pictures per run (120 pictures in total) presented for 2,000 ms each. Of these, 20 objects per run (80 in total) had been studied during the encoding phase. The remaining 10 pictures per run (40 in total) were entirely novel. Participants were presented with two text boxes labeled “old” and “new,” corresponding to “left” or “right” responses. Interleaved with the retrieval trials, we included control stimuli consisting of 6 pairs per run (24 picture-pairs in total) of two-digit numbers located at the center of visually scrambled picture-pairs with a duration of 2,500 ms. As for the encoding and the delay phase, participants were instructed to identify the position of the “even” or “odd” number.

Stimuli were presented via a backprojection system, and behavioral responses were recorded through an optic-fiber response box to measure the hits (the count of correct responses) and reaction time (RT, measured in ms) for each trial. The experiment was performed using Presentation software (Version 23.0, Neurobehavioral Systems, www.neurobs.com).

Before the scanning session, all participants completed a practice run consisting of the three phases described above, using stimuli that were not included in the actual scanning session. After the scanning session, participants completed a forced-choice association test. They were first presented, one at a time, with all the pictures they had correctly identified as “old” during the retrieval phase, i.e., cue items. For each cue item, participants rated their confidence in recalling its associated pair using a four-point Likert scale, where 1 indicated certain recall and 4 indicated guesses. Next, two pictures were displayed on the screen: one was the associated picture from the encoding phase, and the other was a distractor picture shown during the delay. Participants were instructed to select the picture paired with the cue item studied in the encoding phase.

UNIBA

During the 10 min resting-state scan, participants were instructed to remain awake with their eyes open and to fixate on the crosshair in the middle of a white screen (Greicius et al., 2003; Damoiseaux et al., 2006). We used a revised version of the relational and item-specific encoding task (RISE; Ragland et al., 2012, 2015; Passiatore et al., 2023). The event-related task consisted of a picture-pair association encoding followed by a picture-pair associative retrieval (Fig. 2B). The whole scan for the associative memory task lasted approximately 10 min, with the encoding phase lasting 4 min and 12 s and the retrieval phase lasting 7 min and 16 s as described in Passiatore et al. (2023).

During the encoding phase, 35 picture-pairs drawn from the FRIDa dataset (Foroni et al., 2013) were presented simultaneously on the screen for 2,000 ms each, with one object on the left and one on the right. Participants completed 35 relational trials. They were explicitly instructed to study the association of the picture-pairs and indicate whether one object could fit inside the other, responding with the “left” and “right” buttons. Interleaved with the encoding trials, participants were presented with 10 visually scrambled picture-pairs for 2,000 ms. Half of the pairs were sized to a standard 530 × 530 pixels, while in the other half, one picture was upsized and the other downsized by 30% to maintain consistent screen brightness. Participants were instructed to indicate whether “one scrambled picture-pair could fit inside the other,” matching the “left/right” motor demands of the encoding task. The control trials, similar to those at the RUB, served to reset the BOLD signal in the MTL to baseline while enhancing task-related activity (Stark and Squire, 2001).

During the retrieval phase, participants were presented with 70 picture-pairs for 2,000 ms each, including 35 “intact” picture-pairs, i.e., presented in the identical configuration of the encoding phase, and 35 “rearranged” picture-pairs. Participants were instructed to indicate whether picture-pairs were “intact” or “rearranged” using the “left” and “right” buttons, respectively. Interleaved with the retrieval trials, 20 pairs of scrambled pictures were presented for 2,000 ms as control stimuli. As for the encoding phase, participants were instructed to indicate whether “one scrambled picture-pair could fit inside the other,” matching the “left/right” motor demands of the retrieval task.

Stimuli were presented via a backprojection system, and behavioral responses were recorded through an optic-fiber response box to measure accuracy (the count of correct responses) and RT (ms) for each trial. The experiment was performed using Presentation software (Version 23.0).

Before the scanning session, all participants completed a practice version of the task that included five picture-pairs for the encoding phase, interleaved with two scrambled picture-pairs. In the retrieval phase, participants were presented with five “intact” pairs and five “rearranged” pairs. The stimuli used during training were distinct from those used in the scanning session. More details on stimuli selection are described by Passiatore et al. (2023).

MRI data acquisition and processing

Characteristics of the MRI acquisitions are shown in Table 1. We followed the same preprocessing pipeline for both samples using the Statistical Parametric Mapping (SPM) Version 12 (http://www.fil.ion.ucl.ac.uk/spm), the Computational Anatomy Toolbox (CAT12, http://dbm.neuro.uni-jena.de/cat/), and the Advanced Normalization Toolbox (ANTs; Avants et al., 2009). We visually inspected all MRI scans to retain only data unaffected by technical artifacts. Quality-based inclusion criteria were the absence of blurring, ringing, or wrapping on structural (sMRI) and fMRI raw scans (Wood and Mark Henkelman, 1985; Lu et al., 2019), the absence of excessive cropping (>75%; Fornito et al., 2011; Geerligs and Henson, 2016), and the absence of excessive noise, poor image contrast, or poor boundaries in segmented sMRI scans (Song et al., 2006).

Individual sMRI and fMRI scans were reoriented to the anterior commissure–posterior commissure (AC–PC) line with the origin in the AC. The segmentation of T1-weighted scans served to estimate gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF; Ashburner and Friston, 2005). Bias correction removed intensity nonuniformities to improve coregistration with the fMRI scans. Thalamic structural segmentation on the unbiased individual T1-weighted scans was performed through the Thalamus Optimized Multi-Atlas Segmentation (THOMAS; Su et al., 2019; Vidal et al., 2024) based on multiple scan registration steps and multiatlas label fusion. We obtained twelve thalamic segmentations in the individual native space for each participant (Fig. 3, Step 1). The same procedure was performed on the Montreal Neurological Institute (MNI) 152 Template, obtaining 12 thalamic segmentations in the MNI space for group-level parcellations. Because of the coarse resolution and high correlation of the fMRI signal across neighbor voxels, we anatomically grouped thalamic nuclei into four major subdivisions: “anterior,” including anterior ventral, ventral anterior nuclei, and the mammillothalamic tract; “posterior,” including pulvinar, lateral, and medial geniculate nuclei; “medial,” including mediodorsal and centromedian nuclei; and “ventral,” including ventral lateral posterior, ventral posterior lateral, and ventral lateral nuclei (Antonucci et al., 2019; Passiatore et al., 2021). The habenula was excluded due to its small size and proximity to the third ventricle, which may affect the accuracy of the related fMRI signal detection.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Outline of the analysis. Step 1 represents the anatomical segmentation of the thalamus on the individual T1-weighted scans through THOMAS. A summary of the memory task analysis pipeline, including the analysis at the group level on the normalized task activity contrast maps and the anatomical parceling of the thalamic clusters (Step 2a); the analysis at the individual-level task activity within thalamic parcellations’ by overlapping the individual thalamic subdivisions segmented in the individual native space on the non-normalized task activity maps at α < 0.05 (Step 3a); the association between the memory index and the BOLD signal extracted from each thalamic subdivision within the individual activity contrasts across task fMRI scans through multiple regression models (Step 4a). A summary of the resting state analysis pipeline, including the analysis of group-level networks’ spatial maps extracted through ICA on the baseline and postencoding resting-state fMRI (Step 3a); the analysis of individual-level thalamic parcellations within resting state networks’ spatial maps by overlapping the individual thalamic subdivisions segmented in the individual native space on the non-normalized resting state networks’ spatial maps at α < 0.05 (Step 3b); the association of the memory index with (1) the proportion of significant voxels and (2) FC in each thalamic subdivision extracted from the resting state networks’ spatial maps across resting-state fMRI scans through two-part regression models (Step 4b). IC, independent component; FC, functional connectivity; BOLD, blood oxygenation level-dependent.

All echo planar imaging (EPI) fMRI scans were corrected for the differences in acquisition time between slices within a single volume, realigned to the mean scan, coregistered with the unbiased T1-weighted scan, resampled to isotropic voxel sizes of 2 mm³ for RUB and 3 mm³ for UNIBA (reflecting their respective original acquisition resolutions) and smoothed using a 6 mm full-width at half-maximum isotropic 3D Gaussian kernel. EPI fMRI scans were examined for excessive motion exceeding the initial acquisition resolution (RUB, >2 mm in maximum translation and 1.5° in maximum rotation; UNIBA, >3 mm in maximum translation and 1.5° in maximum rotation; mean FD > 0.5 mm for both samples).

At the first level, a general linear model (GLM) was specified to estimate task-related BOLD responses. Each participant's smoothed task-related scans were modeled into a design matrix consisting of one row for each scan and one column for each explanatory variable. Explanatory variables included a stimuli onset vector representing the task events [i.e., hits, false alarms (FA), misses, no response (NR), correct rejections (CR), and control trials], convolved with a canonical hemodynamic response function implemented by SPM12. For the memory task collected at the RUB, the onset vector incorporated the classification of hits responses derived from the postscanning memory test, specifically categorized as hits with a subsequent recall of paired picture (H⁺); hits with incorrect subsequent recall, with confidence ratings between 1 and 3 (H⁻); and hits followed by a correct or incorrect recall of paired pictures with confidence ratings equal to 4 (H⁰), as detailed in Table 2. Additionally, 24 motion parameters (Friston, 2003), WM, and CSF mean signals were included as nuisance regressors to account for residual motion artifacts and non-GM signals. The design matrix was high-pass filtered with a cutoff of 128 s to remove low-frequency drifts. Global scaling was not applied to the data.

To identify event-related activations, contrasts of interest were computed at the first level (Table 3). For the task collected at the RUB, we modeled the following contrast vectors of +1 for the first condition and −1 for the second condition specified:

• trials involving the successful encoding of pictures, represented by the union of H⁺ and H⁻ responses (hereafter referred to as H⁺ ∪ H⁻), excluding trials associated with ratings equal to 4 (H⁰) at the postscanning test, were contrasted with control trials (“even/odd” judgment on scrambled picture-pairs) during the encoding phase;

• the successful retrieval of “old” single items (H⁺ ∪ H⁻) was contrasted with control trials during the retrieval phase; and

• the successful retrieval of “old” single items (H⁺ ∪ H⁻) was contrasted with the correct identification of “new” single items (CR) during the retrieval phase.

For the task collected at UNIBA, we modeled the following contrast vectors of +1 for the first condition and −1 for the second condition specified:

• trials involving the successful encoding of picture-pairs (hits; including only trials successfully retrieved during the retrieval phase) were contrasted with control trials (scrambled picture-pairs) during the encoding phase;

• the successful retrieval of “intact” picture-pairs (hits) was contrasted with control trials during the retrieval phase; and

• the successful retrieval of “intact” picture-pairs (hits) was contrasted with the correct identification of “rearranged” picture-pairs (CR).

For the group-level analysis, individual fMRI contrasts were warped into the standard MNI space using affine and nonlinear transformations via the Symmetric Normalization function in ANTs (Avants et al., 2009, 2011).

Behavioral analysis

The response classification in the two experimental paradigms is summarized in Table 2.

Group-level task-related activity

We evaluated task-related activity at the group-level by entering the individual normalized contrast maps for each memory phase defined in Table 3 into separate one-sample t tests using SPM12 (Fig. 3A, Step 2a). All statistics were nonparametrically corrected for multiple comparisons through the threshold-free cluster enhancement (TFCE) approach (Smith and Nichols, 2009). Significance was set at α < 0.05, using the false discovery rate (FDR) for multiple-comparisons correction (Benjamini and Hochberg, 1995). We extracted the individual BOLD signal from the thalamic significant cluster through Marsbar (Brett et al., 2002). We related thalamic BOLD signal and memory performance through separate multiple regression grouped by memory phases and hemisphere (dependent variable, memory index; independent variable, averaged BOLD signal across the thalamic cluster; nuisance covariates, age, sex, and handedness) in the R environment (https://www.r-project.org). Results were corrected for multiple-comparisons p_FDR < 0.05.

To localize the significant thalamic activity into specific thalamic subdivisions, we performed a descriptive anatomical parceling of the group-level activity maps (Antonucci et al., 2019). We used Marsbar to compute the percentage of spatial overlap between the significant thalamic clusters and the thalamic subdivisions in the MNI space obtained with THOMAS (Fig. 3A, Step 2a).

Individual-level task-related activity and the anatomical parceling of the thalamus

To further investigate interindividual differences in the recruitment of thalamic subdivisions related to memory, we conducted the same analysis in the individual native space to enhance interindividual differences and avoid information loss caused by averaged group-level analysis. We overlapped the individual thalamic subdivisions segmented in the native space and the standardized contrast maps in the individual native space (Fig. 3A, Step 3a) to minimize registration bias and maximize sensitivity to detect individual regional effects that can be impacted by registration errors (Friston, 2003; Avants et al., 2009). Then, we extracted the averaged BOLD signal from the significant voxels (α < 0.05) within a given subdivision.

After scaling the distributions of individual BOLD signals extracted from each thalamic subdivision grouped by memory phases and hemisphere, we explored their relationships with the memory index through separate multiple regression models (dependent variable, memory index; independent variable, BOLD signal in the thalamic subdivision; nuisance covariates, age, sex, and handedness; Fig. 3A, Step 4a), setting significance at p_FDR < 0.05. Analyses were conducted separately for the RUB and UNIBA datasets.

Resting-state functional connectivity analysis

To explore the degree of variability in corticothalamic FC during the resting state, we estimated individual-level FC maps through independent component analysis (ICA) in the baseline resting-state scans in both samples. For the data acquired at the RUB, we performed the same analysis on the postencoding resting-state scan. We used the multiobjective optimization ICA algorithm (MOO-ICAR; Du et al., 2020) implemented in the Group ICA Toolbox (https://trendscenter.org/software/gift/; Calhoun et al., 2001; Calhoun and Adali, 2012). The MOO-ICAR algorithm extracts individual-level FC spatial maps based on a spatial reference (Fig. 3B, Step 2b). We employed the spatial reference from Iraji et al. (2019) to extract reliable FC patterns in high-order functional networks, including attentive, auditory, cerebellar, default mode (DMN), language, bilateral frontoparietal (FPN), sensorimotor, salience, and primary and secondary visual networks topologically determined based on their spatial and temporal properties from previous studies (Beckmann et al., 2005; Damoiseaux et al., 2006, 2008; Smith et al., 2009; Allen et al., 2011; Iraji et al., 2016). We inversely warped the spatial reference into each individual native space to match the preprocessed non-normalized resting-state scans used as the input. Additionally, 24 motion parameters (Friston, 2003) were regressed out from the individual FC network spatial maps to mitigate the impact of motion on connectivity estimates. For the group-level analysis, the obtained individual FC network spatial maps were warped into the standard MNI space using affine and nonlinear transformations via the symmetric normalization function in ANTs.

We evaluated the thalamic recruitment within each cortical network at the group level by entering each network's normalized individual FC maps into separate one-sample t tests using SPM12. Significant clusters were defined at α < 0.05 and corrected for multiple comparisons using the familywise error rate (FWE). We employed the mask of the whole thalamus in the MNI space segmented through THOMAS to localize the effects. Then, we extracted the IC loadings from the thalamic significant cluster through Marsbar. We associated thalamic IC loadings and memory performance through separate multiple regression grouped by network and hemisphere (dependent variable, memory index; independent variable, averaged IC loadings across the thalamic cluster; nuisance covariates, age, sex, and handedness), setting significance at p_FDR < 0.05.

Association between individual-level functional connectivity of thalamic parcellations during resting-state and memory performance

To further investigate interindividual differences in terms of spatial recruitment of corticothalamic FC related to memory, we calculated the degree of spatial overlap between the thalamic subdivisions and the significant voxels included in each non-normalized spatial networks’ map at the individual level (Fig. 3B, Step 3b), defining voxel significance at α < 0.05. Networks’ spatial maps include both positive and negative FC values. Positive values denoted correlated voxels to the individual time course, while negative values denoted anticorrelated voxels to the individual time course (Calhoun et al., 2001; Calhoun and Adali, 2012). We separately extracted the number of voxels significantly correlated and anticorrelated within each network's spatial map, along with the averaged FC from each thalamic subdivision.

To investigate the relationship between corticothalamic recruitment and memory performance across individuals, we analyzed the associations of memory performance with both (1) the proportion of significant voxels in each thalamic subdivision over the number of voxels in the entire thalamus grouped by hemisphere (Fig. 3B, Step 4b) and (2) the averaged FC extracted bilaterally from each thalamic subdivision.

For the analysis of the proportion of significant voxels in each thalamic subdivision, we quantified the number of individuals showing nonzero values in thalamic subdivision recruitment, since the spatial localization of thalamic signals may vary at the individual level, possibly not covering all the thalamic subdivisions. Thus, we restricted the analyses to thalamic subdivisions, including at least N = 10 nonzero observations (Table 4). We used two-part regression models for zero-inflated data using the twopartm R package (Belotti et al., 2015). A preliminary natural logarithmic transformation was applied to the proportion of recruited voxels in each thalamic subdivision to mitigate the right skewness of the distributions. The obtained values were used as the predictor in the second part of the model, while a binary variable indicating whether thalamic recruitment was equal to zero (0) or not (1) was used as the dependent variable in the first part of the model. In the first part of the model, a binomial regression was fitted with a logit link function to measure the probability of a nonzero outcome. The second part of the model is a GLM (dependent variable, memory index; independent variable, proportion of significant voxels in the thalamic subdivision; nuisance covariates, age, sex, and handedness). To facilitate the interpretation of results, we also calculated the averaged marginal effects (AME; Belotti et al., 2015), i.e., the estimated impact of the proportions of thalamic voxels recruited on the memory index calculated by considering the combined effect of the probability of the thalamic recruitment occurring (first part) and the magnitude of the thalamic recruitment association with the memory index (second part) and the confidence intervals (CI) at 0.95.

Similarly, for the analysis of FC in each thalamic subdivision, GLMs were used to test the effect of memory performance on the averaged FC, since FC did not include zero values (dependent variable, memory index; independent variable, averaged FC in the thalamic subdivision; nuisance covariates, age, sex, and handedness). All analyses were conducted separately for the RUB and UNIBA datasets.

In the RUB sample, we also computed separate two-part regressions and GLMs for the postencoding resting-state spatial maps to assess changes in corticothalamic recruitment after the performance of the associative memory task. In addition, we analyzed the association of memory performance with the variation in terms of thalamic recruitment from the baseline to the postencoding resting state. We computed a change score representing the difference in the proportion of significant voxels between the post- and the baseline resting-state scans for each thalamic subdivision, as well as the difference of the FC between postencoding and baseline resting state. Positive values indicate an increase in the proportion of significant voxels or FC from the baseline to the postencoding resting state, while negative values indicate a decrease in the proportion of significant voxels or FC from the baseline to the postencoding resting state. Zero values indicate no change in the proportion of significant voxels across resting-state scans. We tested the relationship between memory performance and the change score for each thalamic subdivision within each network and grouped by hemisphere through separate GLM (dependent variable, memory index; independent variable, change score; nuisance covariates, age, sex, and handedness), setting significance at p_FDR < 0.05.

Comparison across connectivity metrics: FC versus the proportion of significant voxels across subdivisions

With the hypothesis that FC could flatten the interindividual variability that is highlighted when using the spatial informed metric (H₀ σ²A = σ²B vs H₁ σ²A > σ²B), we compare the connectivity metrics previously modeled, i.e., the proportion of significant voxels and the FC across the thalamic subdivisions during resting state. Since the two metrics are on different scales and thus not directly comparable, we used the min–max scaling (de Amorim et al., 2023), which scales each distribution individually by its minimum and maximum absolute value, not centering the data, thereby preserving any sparsity. To test this hypothesis, we used Levene's tests to compare the variance between the distributions of FC and the proportion of significant voxels grouped by cohort, network, and hemisphere. Additionally, we examined the correlation between the two metrics across thalamic subdivisions to assess the discriminatory power of the measurement used through the Kendall–Tau correlation. Significance was set at α < 0.05.

Path analysis on the individual-level thalamic recruitment across resting-state and memory phases

To investigate the relationships between the individual associative memory performance and the individual-level corticothalamic recruitment across baseline resting-state FC, the activity during the encoding and the retrieval phases, serial path analyses within a structural equation modeling (SEM) framework were implemented using the lavaan R package (Rosseel, 2012). We included the thalamic recruitment as the percentage of significant voxels extracted during the baseline resting-state IC spatial maps as the predictor (X); the individual task-related activity in the encoding (M₁) and retrieval (M₂) during the associative memory task as mediators; and the index of the associative memory performance as the outcome (Y). The path analysis considered the thalamic subdivisions significantly associated with associative memory performance in the previous analysis. Age, sex, and handedness were used as covariates, and a dichotomic variable was used to control for nonzero observations. Significance was set as α < 0.05 (bootstrap, 5,000 runs). We used the effectsize R package (Ben-Shachar et al., 2020) to calculate the magnitude of the effects reported as R² (Liu et al., 2023).

k-mean clustering on the individual-level thalamic recruitment during the associative memory task

To test our hypothesis that corticothalamic circuits recruited by different individuals may underlie variability in memory performance, we performed a k-mean clustering using the kmeans R function (MacQueen, 1967; Hartigan and Wong, 1979) on the individual task-related activity during (1) the encoding and (2) the retrieval phases during the associative memory task, previously included in the indirect paths of the path analysis. Also, the cluster analysis considered the thalamic subdivisions significantly associated with memory performance in the previous analysis. We performed the clustering only on the UNIBA dataset, which included enough observations to capture variability across individuals. We determined the optimal number of clusters for the left and right hemispheres through the NbClust R package (Charrad et al., 2014). We employed the majority rule to select the optimal number of clusters. The cluster stability was calculated through the bootcluster R package considering S > 0.8 (bootstrap, 5,000 runs; Yu et al., 2019). We separately compared the memory index, the proportion of significant voxels within the baseline resting-state spatial maps, the individual BOLD estimates during the encoding and retrieval phases during task performance, age, education, and intellective quotient across clusters through the Wilcoxon rank-sum test, and the sex was compared through the chi-square test. We set the test significance at p_FDR < 0.05.

Data and code accessibility

Data from RUB and UNIBA cannot be shared at the individual level because of ethical restrictions based on the protocol approved by the relative institutional ethics committees to protect the privacy of the participants. All R code, along with detailed group statistics needed to evaluate the conclusions in this work, is publicly available at https://github.com/robertapassiatore/thalamus-memory.git. The THOMAS toolbox is publicly available at https://github.com/thalamicseg/thomas_new.git. The Group ICA toolbox is publicly available at https://trendscenter.org/software/gift. The Presentation scenarios used for the associative memory tasks are available at https://doi.org/10.5281/zenodo.14374018.