Age-specific relationship between the modulation of brain dynamics in response to task demands and bimanual performance

Behavioral data

Effect of age group and task condition on BTT performance

A 2 x 3 Mixed Model ANOVA was performed to evaluate the effects of AGE GROUP and TASK CONDITION on BTT SCORE. The means and standard deviations for BTT SCORE are presented in Table 1 and can be visualized in Figure 1.

Figure 1. Mean BTT scores across task conditions per age group. The figure shows mean BTT scores (%) for each task condition (Line 1:1, Line 3:1, and Angle 3:1) in Younger Adults (teal, triangle) and Older Adults (orange, circle). Whiskers indicate standard deviation. Both groups achieved their highest scores in Line 1:1 and their lowest in Angle 3:1. Following a significant main effect of TASK CONDITION, post hoc Tukey Honestly Significant Difference tests revealed that Angle 3:1 differed significantly from both Line 1:1 and Line 3:1 (*** p < 0.001), while the two “Line” conditions did not differ (n.s.). Abbreviations: BTT = Bimanual Tracking Task; n.s. = non-significant.

The two-way ANOVA revealed significant main effects of AGE GROUP (F[1,129] = 29.71, p < 0.001) and TASK CONDITION (F[2,129] = 16.80, p < 0.001), indicating that BTT performance differed between younger and older adults and varied across task conditions. As expected, younger adults (mean ± SD = 50.27 ± 12.69) performed on average better than older adults (mean ± SD = 31.78 ± 14.43). Moreover, participants overall performed worse on the Angle 3:1 (mean ± SD = 29.94 ± 13.32) comparatively to Line 3:1 (mean ± SD = 46.20 ± 14.91) and Line 1:1 (mean ± SD = 47.17 ± 15.84). Post-hoc Tukey HSD comparisons confirmed that Angle 3:1 differed from “Line” conditions (q < 0.001), whereas Line 1:1 and Line 3:1 did not differ from one another (q = 0.826). Lastly, the AGE GROUP × TASK CONDITION interaction was not significant (F[2,129] = 0.019, p = 0.981), suggesting that the effect of task condition on BTT performance did not differ between age groups.

Neuroimaging data

Effect of age group and task condition on framewise displacement

We examined differences in framewise displacement (FD) across age groups and task conditions using an Aligned Rank Transform (ART) factorial ANOVA to assess the potential confounding effect of head motion on BOLD SD estimates.

A significant main effect of AGE GROUP was found (F[1,129] = 42.89, p < 0.001), indicating that older adults (mean ± SD = 0.34 ± 0.12) exhibited greater FD during BTT performance compared to younger adults (mean ± SD = 0.23 ± 0.07). However, there was no significant main effect of TASK CONDITION (F[2,129] = 0.39, p = 0.675), suggesting that movement-related artifacts did not systematically differ across task conditions. Additionally, the AGE GROUP × TASK CONDITION interaction was not significant (F[2,129] = 0.17, p = 0.843), indicating that the relationship between BTT conditions and FD did not differ across age groups. These results are depicted in Figure 2.

Figure 2. Framewise displacement (FD) across age groups. The figure shows framewise displacement (FD) across age groups (Younger Adults in teal, and Older Adults in orange) while performing the Bimanual Tracking Task (BTT). The “X” indicates mean FD for each age group. A significant main effect of AGE GROUP was observed (the corresponding F-value and p-value are reported), with older adults presenting overall more framewise displacement compared to younger adults. Individual data points reflect mean FD per subject, pooled across task conditions. Significance level: *** p < 0.001.

These findings indicated a potential bias in BOLD SD estimation in older adults compared to younger due to more head motion, reinforcing the importance of controlling for FD in subsequent analyses.

Age-group differences in BOLD variability and BOLD variability modulation during BTT performance

Effect of age group and task condition on BOLD variability

To investigate whether BOLD variability differed across age groups and task conditions in the selected brain regions, a LMM was conducted per ROI, controlling for FD. Subject ID was included as a random intercept to account for within-subject variance. Box-Cox transformations were applied in 21 ROIs to meet model assumptions (see Table 2 for the λ values used in the transformation). On average, two data points of BOLD SD were removed per ROI (i.e., 1.5 % of the total data points; the range of data points excluded was from zero to four), as they were more than three standard deviations of the ROI’s overall mean BOLD SD.

Table 2. Summary of age group, task condition and interaction effects on BOLD variability across ROIs.

ROIs		λ		Age group			Task condition			Age group x task condition
				F	q		F	q		F	q
RH	A2	0.125		6.27	0.129		-	-		-	-
LH	A37vl	-0.625		-	-		-	-		-	-
RH	A39c	-		-	-		-	-		-	-
LH	A40rd	-0.654		-	-		-	-		-	-
RH	A5l	-		-	-		-	-		-	-
LH	A6cdl	-0.033		0.87	0.986		1.63	0.628		3.11	0.279
LH	A6dl	-0.801		-	-		3.41	0.313		-	-
RH	A6dl	0.130		5.15	0.172		-	-		-	-
LH	A6vl	-0.431		-	-		-	-		-	-
RH	A7c	-1.634		-	-		-	-		-	-
RH	A7ip	-0.255		-	-		-	-		-	-
RH	A7m	-0.092		3.45	0.287		1.18	0.426		3.15	0.279
LH	A7r	-		-	-		4.43	0.278		-	-
LH	iOccG	-0.814		-	-		-	-		-	-
LH	IsOccG	-3.000		5.07	0.172		4.74	0.278		-	-
LH	mOccG	-0.860		8.12	0.079		-	-		-	-
LH	OPC	-0.547		-	-		-	-		-	-
LH	V5/MT+	-0.800		0.64	1.000		3.07	0.323		4.17	0.279
RH	V5/MT+	-1.043		0.11	1.000		2.47	0.375		3.52	0.279
LH	Crus I	0.206		1.27	0.829		0.18	1.000		3.97	0.279
RH	Crus I	-1.788		2.14	0.536		0.07	1.000		4.73	0.279
LH	VI	-		-	-		-	-		-	-
RH	VI	-1.017		0.11	1.000		2.76	0.345		4.00	0.279
LH	VIIb	0.145		-	-		-	-		-	-
LH	VIIIb	-0.330		13.52	0.014*		-	-		-	-
Summary of results from the Linear Mixed Model (LMM) conducted per ROI to investigate whether BOLD variability differed across age groups and task conditions while controlling for framewise displacement (at the significance level of q < 0.05, after False Discovery Rate [FDR] correction for multiple comparisons). Values in bold reflect significant effects before FDR correction (p < 0.05). A dash (“-”) indicates that the corresponding term was excluded from the final model following stepwise simplification, or that no Box-Cox transformation was applied to the data. Regions in which all effects are indicated with a dash, the final model solely included framewise displacement as an effect. Significance level: *q < 0.05.
LH, Left Hemisphere; RH, Right Hemisphere; ROIs, Regions of Interest.

Out of the 25 ROIs analyzed, five showed a significant main effect of AGE GROUP before FDR correction. In all of these regions - namely, right A2, right A6dl, left lsOccG, left mOccG, and left VIIIb – older adults exhibited greater BOLD SD compared to younger adults (see Figure 3; right A2: F[1, 51.03] = 6.27, p = 0.016, β [Younger] = -0.10; right A6dl: F[1, 50.68] = 5.15, p = 0.028, β [Younger] = -0.08; left lsOccG: F[1, 48.70] = 5.07, p = 0.029, β [Younger] = -0.95; left mOccG: F[1, 50.90] = 8.12, p = 0.006, β [Younger] = -0.24); left VIIIb: F[1, 50.64] = 13.52, p < 0.001, β [Younger] = -0.18). After correcting for multiple comparisons (with Benjamini-Hochberg FDR method), only the effect in the left VIIIb (q = 0.014) remained significant at q < 0.05.

Figure 3. BOLD variability by age group in significant ROIs for main effect of age group before FDR correction. The figure shows transformed BOLD variability (BOLD SD) across age groups (Younger Adults in teal, and Older Adults in orange) during the Bimanual Tracking Task (BTT) for regions of interest (ROIs) that showed a significant main effect of AGE GROUP prior to FDR correction for multiple comparisons (p < 0.05). From these, only the group differences in the cerebellar lobe VIIIb survived correction (with q > 0.05). Individual data points reflect BOLD SD values pooled across all task conditions, per subject. F- and p-values for the main effect of AGE GROUP are reported for each ROI. Uncorrected significance levels: * p < 0.050, ** p < 0.010, *** p < 0.001. Corrected significance level: q < 0.05. Abbreviations: A2 = primary somatosensory cortex; A6dl = dorsolateral premotor cortex; FDR = False Discovery Rate; L = left; lsOccG = lateral superior occipital gyrus; mOccG = middle occipital gyrus; R = right; VIIIb = cerebellar lobule VIIIb.

Regarding TASK CONDITION, uncorrected effects were found in the left A6dl, left A7r, and left lsOccG, with BOLD SD generally scaling up with task difficulty (i.e., Line 1:1 < Line 3:1 < Angle 3:1) in both age groups (left A6dl: F[2, 84.08] = 3.41, p = 0.038, β [Line31] = 0.02 and β [Angle31] = 0.05; left A7r: F[2, 89.36] = 4.43, p = 0.015, β [Line31] < 0.01 and β [Angle31] = 0.49; left lsOccG: F[2, 48.70] = 5.07, p = 0.011, β [Line31] = 0.30 and β [Angle31] = 0.02). However, none of these effects survived FDR correction (with significance set at q < 0.05).

An uncorrected interaction between AGE GROUP and TASK CONDITION was found in six ROIs: right A7m, bilateral V5/MT+, bilateral Crus I, and right VI (right A7m: F[2, 91.08] = 3.14, p = 0.047, β [Younger; Line31] = -0.05 and β [Younger; Angle31] = 0.02); left V5/MT+: F[2, 87.13] = 4.17, p = 0.019, β [Younger; Line31] = -0.02 and β [Younger; Angle31] = 0.08; right V5/MT+: F[2, 90.26] = 3.53, p = 0.033, β [Younger; Line31] = -0.08 and β [Younger; Angle31] = 0.02; left Crus I: F[2, 89.70] = 3.98, p = 0.022, β [Younger; Line31] = -0.03 and β [Younger; Angle31] = 0.03; right Crus I: F[2, 88.92] = 4.73, p = 0.011, β [Younger; Line31] = -0.04 and β [Younger; Angle31] = 0.08; right VI: F[2, 89.95] = 3.99, p = 0.022, β [Younger; Line31] = -0.03 and β [Younger; Angle31] = 0.05). None of the interaction effects survived FDR correction.

A visualization of the regions showing significant uncorrected results for each effect, can be found in Figure 4. For a summary of the corresponding results after FDR correction, see Table 2. In Figure 5, the group-specific patterns of BOLD SD change across task conditions are depicted for ROIs that showed a significant uncorrected effect of TASK CONDITION or AGE GROUP X TASK CONDITION.

Figure 4. Regions showing significant effects of age group, task condition and their interaction on BOLD variability prior to FDR correction. The figure shows, separately for the main effects of AGE GROUP and TASK CONDITION, as well as their interaction effect, the location of the brain regions that yield a significant result prior to FDR correction. Only the effect of AGE GROUP in the left cerebellar lobule VIIIb (circled region) remained significant after correction for multiple comparisons (with q < 0.05). Color gradients indicate uncorrected p-values for the corresponding effect (red for main effect of AGE GROUP, green for main effect of TASK CONDITION and blue for their interaction) in a given region. Brain slices and volumes depicted in this figure were created with MRIcroGL (Version 1.2.20220720b; https://www.nitrc.org/projects/mricrogl) [37]. Abbreviations: A2 = primary somatosensory cortex; A6dl = dorsolateral premotor cortex; A7m = precuneus; A7r = superior parietal lobe; Crus I = cerebellar lobule Crus I; FDR = False Discovery Rate; L = left; lsOccG = lateral superior occipital gyrus; mOccG = middle occipital gyrus; R = right; V5/MT+ = motion area V5/MT+; VI = cerebellar lobule VI; VIIIb = cerebellar lobule VIIIb.

Figure 5. BOLD variability by age group and task condition in several ROIs. Figure shows BOLD variability transformed (except for A7r L, as indicated with “⁰”) across age groups (Younger Adults in varying shades of teal, and Older Adults in different shades of orange) and task conditions (Line 1:1, Line 3:1 and Angle 3:1) for ROIs with a significant main effect of TASK CONDITION (namely, A6dl L, A7r L and lsOccG L) or AGE GROUP x TASK CONDITION interaction (i.e., A7m R, bilateral V5/MT+, bilateral Crus I and VI R), before correcting for multiple comparisons. None of these results remained significant after correction (in all, q > 0.05). Mean BOLD SD, per group, for each task condition is indicated by an “X”. Abbreviations: A6dl = dorsolateral premotor cortex; A7m = precuneus; A7r = superior parietal lobe; Crus I = cerebellar lobule Crus I; L = Left; L = left; lsOccG = lateral superior occipital gyrus; OA = Older Adults; R = Right; V5/MT+ = motion area V5/MT+; VI = cerebellar lobule VI; YA = Younger Adults.

Lastly, FRAMEWISE DISPLACMENT was consistently a robust predictor across all models, showing a significant positive association with BOLD SD in all ROIs (after FDR correction all ROIs kept a significance of q < 0.001), except in left V5/MT+ (p = 0.006), right VI (p = 0.03), and left VIIIb (p = 0.02).

For detailed results of this analysis by ROI, consult Supplementary Table 1.

Impact of age-group differences in BOLD variability and BOLD variability modulation on BTT performance

Effect of age-group differences in BOLD variability modulation on BTT performance

To investigate the relationship between BOLD variability modulation, the age group and behavioral performance, we conducted a multiple linear regression for each of the sixteen ROIs (i.e., regions highlighted in bold in Table 4). In each model, MEAN BTT SCORE was predicted by BOLD SD MODULATION, AGE GROUP (Younger Adults vs. Older Adults), and their interaction. On average, per ROI, one participant (generally an older adult; corresponding to 2.2% of the total number of participants; ranging from zero to two participants) was excluded from the analysis because BOLD SD modulation was more than three standard deviations from the region’s mean BOLD SD modulation.

Table 4. Local maxima and correspondent brain regions in the atlases.

Brain region			Atlas label	Local maxima
				Z	x	y	z
Brainnetome
LH	Superior frontal gyrus	A6dl, dorsolateral area 6	7	5.62	56	65	67
RH	Superior frontal gyrus	A6dl, dorsolateral area 6	8	5.93	35	68	63
LH	Middle frontal gyrus	A6vl, ventrolateral area 6	25	6.08	56	64	61
LH	Precentral gyrus	A6cdl, caudal dorsolateral area 6	55	5.99	62	59	61
LH	Inferior temporal gyrus	A37vl, ventrolateral area 37	97	6.88	70	31	33
LH	Superior parietal lobule	A7r, rostral area 7	125	7.47	51	30	68
				6.5	53	35	68
				6.48	51	35	72
RH	Superior parietal lobule	A7c, caudal area 7	128	6.74	37	26	65
RH	Superior parietal lobule	A5l, lateral area 5	130	5.84	28	45	57
RH	Superior parietal lobule	A7ip, intraparietal area 7	134	5.7	30	38	66
RH	Superior parietal lobule	A39c, caudal area 39	136	5.69	24	31	44
				5.66	21	25	47
				5.65	24	24	50
LH	Inferior parietal lobule	A40rd, rostrodorsal area 40	139	6.08	64	43	58
				5.98	68	45	56
RH	Precuneus	A7m, medial area 7	148	6.91	42	33	65
RH	Postcentral gyrus	A2, area 2	160	6.07	24	49	57
LH	Lateral occipital cortex	mOccG, middle occipital gyrus	199	7.67	59	17	42
				6.97	63	20	38
LH	Lateral occipital cortex	V5/MT+, area V5/MT+	201	7.5	66	27	37
RH	Lateral occipital cortex	V5/MT+, area V5/MT+	202	6.69	21	32	37
				5.77	21	26	39
LH	Lateral occipital cortex	OPC, occipital polar cortex	203	6.45	53	13	40
LH	Lateral occipital cortex	iOccG, inferior occipital gyrus	205	5.7	62	22	30
LH	Lateral occipital cortex	lsOccG, lateral superior occipital gyrus	209	6.92	55	21	57
Cerebellum SUIT
LH	Cerebellum	VI	5	5.86	61	38	18
RH	Cerebellum	VI	7	5.73	32	31	23
LH	Cerebellum	Crus I	8	5.96	65	37	19
RH	Cerebellum	Crus I	10	6.17	22	32	22
				5.67	26	37	20
LH	Cerebellum	VIIb	14	6.4	47	25	16
LH	Cerebellum	VIIIb	20	6.16	51	32	10
Coordinates and Z-values of local maxima spaced at least 8 mm apart, derived from clusters containing 10 or more voxels identified in the functional activation contrast “Angle 3:1 > Line 3:1”, since the contrast “Line 3:1 > Line 1:1” did not show significant clusters of brain activation. Activation peaks were subsequently labeled using the Brainnetome Atlas (for cortical regions) and the Cerebellum SUIT Atlas (for cerebellar regions), resulting in the 25 brain regions indicated (see Figure 9). Only the regions highlighted in bold integrated all analyses in the current study (see section 4.8.2.3).
LH, Left hemisphere; RH, Right hemisphere.

After FDR correction, significant main effects of AGE GROUP were consistently observed across all analyzed ROIs (q < 0.001), indicating that age remained a strong predictor of BTT performance regardless of BOLD SD modulation.

An uncorrected main effect of BOLD SD MODULATION was found in the left lsOccG, mOccG and Crus I (left lsOccG: p = 0.028, β = 41.87; left mOccG: p = 0.011, β = 25.66; left Crus I: p = 0.003, β = - 38.18). However, only in the left mOccG and left Crus I these effects remained significant after correction for multiple comparisons (q = 0.042 and q = 0.014, respectively).

Additionally, FDR-corrected significant AGE GROUP × BOLD SD MODULATION interactions emerged in the right A2, left A40rd, left A7m, left A7r, left V5/MT+, right V5/MT+, and right VI (right A2: q = 0.001, β [younger] = - 215.52; left A40rd: q = 0.002, β [younger] = - 123.33; left A7m: q = 0.039, β [younger] = - 81.95; left A7r: q < 0.001, β [younger] = - 172.45; left V5/MT+: q < 0.001, β [younger] = -156.75; right V5/MT+: q < 0.001, β [younger] = - 168.31; right VI: q = 0.003, β [younger] = - 191.81). These results suggested that the relationship between BOLD variability modulation and task performance differed by age in these regions (check Figure 6). A summary of the FDR-corrected results for all regions analyzed can be found in Table 5. Detailed results of this analysis for each ROI are presented in Supplementary Table 3.

Figure 6. Age-related modulation of BOLD variability and its relationship with BTT performance. For each region of interest (ROI) where BOLD SD MODULATION significantly predicted BTT MEAN SCORE after FDR correction in at least one of the age groups, the figure shows: [A] BOLD variability (BOLD SD) modulation across age groups (Younger Adults in teal [triangles], and Older Adults in orange [circles]). In the right primary somatosensory cortex (A2), the left inferior parietal lobe (A40rd), as well as the left cerebellar lobules VIIb and VIIIb (not displayed in the figure), older adults exhibited, on average (“X”), more BOLD SD modulation than younger adults, with q < 0.05 (FDR-corrected). Significance level: * q < 0.050; [B] The relationship between BOLD SD modulation and mean BTT performance, shown per age-group (Younger Adults again in teal [triangles], and Older Adults in orange [circles]). After FDR correction, younger adults consistently exhibited a significant negative association between BOLD SD modulation and performance across all ROIs analyzed, except for the left inferior parietal lobe (A40rd), while older adults showed a positive relationship in the left inferior and parietal lobes (left A40rd and A7r, respectively). Dashed regression lines identify effects that did not survive FDR correction. Abbreviations: A2 = primary somatosensory cortex; A40rd = inferior parietal lobe; A7r = superior parietal lobe; FDR = False Discovery Rate; L = Left; R = Right; V5/MT+ = motion area V5/MT+; VI = cerebellar lobule VI.

Table 5. Effect of age group, BOLD SD modulation and their interaction on mean BTT score in all ROIs analyzed.

ROIs		Age group			BOLD SD modulation			Age group X BOLD SD modulation
		F	q		F	q		F	q
RH	A2	44.36	< 0.001***		2.83	0.182		13.47	0.001**
RH	A39c	131.19	< 0.001***		-	-		-	-
LH	A40rd	71.57	< 0.001***		18.85	< 0.001***		12.60	0.002**
LH	A6dl	131.19	< 0.001***		-	-		-	-
RH	A6dl	128.53	< 0.001***		-	-		-	-
RH	A7m	39.73	< 0.001***		2.49	0.188		5.86	0.039*
LH	A7r	66.43	< 0.001***		16.42	< 0.001***		17.22	<0.001***
LH	IsOccG	138.17	< 0.001***		4.94	0.064		-	-
LH	mOccG	143.00	< 0.001***		6.74	0.042*		-	-
LH	V5/MT+	81.15	< 0.001***		6.04	0.042*		23.02	<0.001***
RH	V5/MT+	60.89	< 0.001***		6.50	0.042*		17.58	<0.001***
LH	Crus I	110.74	< 0.001***		9.35	0.014*		-	-
RH	Crus I	131.19	< 0.001***		-	-		-	-
RH	VI	33.23	< 0.001***		1.17	0.410		12.17	0.003**
LH	VIIb	123.87	< 0.001***		2.90	0.182		-	-
LH	VIIIb	132.90	< 0.001***		-	-		-	-
Summary of results from the regression analyses conducted per ROI to investigate whether BOLD SD modulation predicted mean BTT performance (Mean BTT Score) and whether this relationship differed by age group. Values in bold reflect significant effects (p < 0.05) before FDR correction for multiple comparisons. A dash (“-”) indicates that the corresponding term was excluded from the final model following stepwise simplification. Significance levels: *q < 0.05. **q < 0.01. ***q < 0.001.
LH, Left Hemisphere; RH, Right hemisphere; ROIs, Regions of Interest.

For ROIs where residuals violated normality or homoscedasticity assumptions, we employed residuals bootstrapping (20,000 samples) to obtain robust estimates and 95% confidence intervals. Predictors whose bootstrapped confidence intervals excluded zero were considered statistically significant. Importantly, after applying this residuals’ bootstrapping procedure, the estimated regression coefficients and confidence intervals remained largely consistent with those from the original models (i.e., 95 % CIs [displayed in Table 6] overlapped with the original estimates). This robustness supported the reliability of the observed effects and reinforced the conclusions drawn from the standard analysis.

Table 6. Confidence intervals from percentile residual bootstrapping analysis.

ROIs		Age group		BOLD SD modulation		Age group x BOLD SD modulation
		95% bootstrap CIs		95% bootstrap CIs		95% bootstrap CIs
RH	A2	[27.33, 49.81]		[-8.21, 111.28]		[-328.50, -100.30]
LH	A40rd	[24.33, 38.98]		[58.60, 154.30]		[-190.90, -55.00]
LH	A6dl	[15.03, 37.93]		[-55.32, 78.22]		[-224.40, 33.70]
RH	A6dl	[15.32, 40.07]		[-80.60, 96.43]		[-215.80, 27.10]
RH	A7m	[20.15, 38.13]		[-8.90, 91.48]		[-146.80, -15.50]
LH	A7r	[27.68, 44.89]		[46.70, 129.30]		[-250.80, -91.40]
LH	IsOccG	[14.74, 33.49]		[10.93, 92.55]		[-135.99, 46.08]
LH	V5/MT+	[28.99, 44.77]		[8.31, 65.49]		[-218.70, -93.50]
RH	Crus I	[9.59, 27.25]		[-41.51, 20.33]		[-55.36, 52.94]
RH	VI	[29.86, 60.37]		[-22.27, 80.11]		[-298.89, -85.69]
LH	VIIIb	[11.99, 37.37]		[-26.20, 16.57]		[-103.25, 28.67]
95 % Confidence Intervals (CIs) based on percentiles for AGE GROUP, BOLD SD MODULATION and their interaction resulting from a bootstrapping analysis of the residuals conducted for each region of interest (ROI) whose residuals violated normality and/or homoscedasticity assumptions.
CIs, Confidence Intervals; LH, Left Hemisphere; RH, Right Hemisphere; ROIs, Regions of Interest.

Lastly, a post hoc simple regression analysis was conducted per age group to determine the relationship between MEAN BTT SCORE and BOLD SD MODULATION in regions showing an FDR-corrected significant AGE GROUP X BOLD SD MODULATION interaction. These analyses revealed that younger adults exhibit a significant negative association between BTT performance and BOLD SD modulation in five out of the seven regions analyzed (namely, right A2, left A7r, bilateral V5/MT+ and right cerebellar lobule VI) after FDR correction (all with q < 0.05). Conversely, in older adults, only in the left A40rd and the left A7r, higher BOLD SD modulation predicted higher BTT mean score before and after FDR correction for multiple comparisons (q < 0.001 and q = 0.001, respectively). Consult Figure 6 for a visualization of the relationship between BOLD variability modulation and BTT performance per age group, in all the ROIs analyzed. Coefficient estimates and confidence intervals from the bootstrapping analyses validated these same conclusions. Specifically, for younger adults, the 95% bootstrapped confidence intervals were the following: right A2 [- 276.30, - 57.50], left A40rd [- 126.90, 29.33], right A7m [- 107.98, 19.42], left A7r [- 190.20, 31.40], left V5/MT+ [- 180.70, - 18.50], right V5/MT+ [- 203.38, - 23.66] and right VI [- 328.70, - 36.50]. In contrast, for older adults, the confidence intervals were: right A2 [- 45.42, 142.03], left A40rd [32.8, 165.4], right A7m [- 30.90, 126.45], left A7r [19.30, 144.40], left V5/MT+ [- 4.84, 79.76], right V5/MT+ [- 14.01, 134.65], and right VI [- 53.28, 107.19].

An overview of the main findings from the analyses focusing on BOLD SD modulation was depicted in Figure 7.

Figure 7. Significant findings on BOLD SD modulation and its relationship with bimanual task performance in aging. Summary of brain regions showing significant findings (FDR-corrected, with q < 0.05) relative to BOLD SD Modulation and its association with performance on the Bimanual Tracking Task (BTT). Cortical regions were defined based on the Brainnetome Atlas [38], and cerebellar regions using the Cerebellum SUIT Atlas [39]. Brain slices and volumes depicted in this figure in the MNI152 Template were created with MRIcroGL (Version 1.2.20220720b; https://www.nitrc.org/projects/mricrogl) [37]. Abbreviations: A2 = primary somatosensory cortex; A40rd = inferior parietal lobe; A7m = precuneus; A7r = superior parietal lobe; Crus I = cerebellar lobule Crus I; FDR = False Discovery Rate; L = left; mOccG = middle occipital gyrus; R = right; V5/MT+ = motion area V5/MT+; VI = cerebellar lobule VI; VIIb = cerebellar lobule VIIb; VIIIb = cerebellar lobule VIIIb.

Behavioral results: Lower overall BTT performance in older adults, but similar effect of task condition across age groups

Our behavioral analysis revealed robust age-related differences in bimanual performance. Across all conditions of the BTT, younger adults significantly outperformed older adults. Task difficulty also influenced performance: participants performed best on the two line-tracking conditions (Line 1:1 and Line 3:1), and significantly worse on the angle condition (Angle 3:1), indicating that increased complexity impaired task execution. Post hoc comparisons confirmed that performance in the Angle 3:1 condition was significantly lower than in both line-tracking conditions, while no difference emerged between the two “Line” conditions. Importantly, the lack of an interaction between age group and task condition suggests that both younger and older adults were similarly affected by increasing task demands, even though their overall performance levels differed, which contrasts with previous findings indicating that older adults gradually lose their capacity to adapt to varying task requirements [6, 40–42]. This discrepancy between our results and existing literature may stem from differences in the type or amount of feedback provided [43, 44], or even the amount of exposure to the task conditions during the familiarization session. For instance, in all four aforementioned studies, BTT practice during the familiarization session lasted for a fixed duration of six minutes, which may have favored faster-learning in young adults who improved more during the practice, particularly under more challenging task conditions. In contrast, in our study, training was terminated once participants reached a performance threshold of 10%, thereby minimizing initial inter-individual discrepancies in task performance during the MR session. Nevertheless, these findings support the assumption that the Angle 3:1 condition imposed greater cognitive-motor coordination challenges across the lifespan. They also point to a graded behavioral pattern aligned with task complexity, as performance gradually declined from Line 1:1 to Line 3:1 and then more sharply in the Angle 3:1 condition. Although the difference between the two “Line” conditions did not reach statistical significance, the observed trend is consistent with a complexity-driven performance gradient that may help contextualize neural variability modulation.

Higher BOLD variability in older adults, but no association with BTT performance

While previous cognitive studies reported reduced task-related BOLD variability with aging [14, 15, 34], we found higher task-related BOLD SD in older relative to younger adults in the left VIIIb. This subdivision of the cerebellar lobule VIII is located inferiorly on the posterior cerebellum, and is known to be connected with the somatosensory cerebral network [45], assisting the primary motor cortex (M1) [46] in the execution and coordination of fine limb movements [47], containing, in particular, the somatomotor representations of individual left fingers [48]. This divergence between our results showing higher BOLD SD in older relative to younger adults, and findings from existing literature reporting the opposite pattern, reinforce the notion that variability-age relationships might be task-dependent [27, 31, 32, 35].

Nevertheless, given prior work from Scarapicchia et al. [31] reporting a positive relationship between white matter hyperintensities (WMH) and BOLD variability, particularly in the cerebellum, we speculate that the increase in BOLD variability we observed with aging in the left cerebellar lobule VIIIb might have a vascular or structural origin. It might as well be driven by an age-related disruption of the dopaminergic neuromodulation systems, as lower dopamine has been previously linked to greater BOLD variability, especially in subcortical structures such as the cerebellum [28, 29]. However, the current study did not examine the role of WHM nor of dopaminergic differences in BOLD variability. Hence, future research should investigate these factors as potential mechanisms underlying BOLD variability shifts with aging.

Regarding the potential behavioral relevance of these age-related differences in BOLD variability, we did not find robust evidence supporting a brain-behavior relationship in bimanual performance. This stands in contrast to prior findings in cognitive domains linking BOLD variability to performance [14, 27].

Age-group differences in the magnitude of BOLD variability modulation predict BTT performance

Our data revealed greater BOLD variability modulation across BTT conditions in older compared to younger adults within the right primary somatosensory cortex (A2), left inferior parietal lobule (A40rd) and left cerebellar lobules VIIb and VIIIb. This pattern diverges from prior findings in cognitive domains, where aging has typically been associated with reduced intra-individual brain variability [25, 33, 34], possibly reflecting domain-specific differences in how neural variability is modulated across the lifespan. The right A2 and the left A40rd support the integration of sensory input necessary to motor planning and execution [49, 50], whereas VIIb and VIIIb are implicated in the active maintenance of attentional processes and working memory [45, 47, 51, 52]. Notably, however, these were not the regions where BOLD variability modulation directly predicted bimanual task performance in the present study.

In addition to age-related differences in BOLD variability modulation, we observed region-specific relationships between BOLD variability modulation and task performance. In the left middle occipital gyrus (mOccG), greater modulation of BOLD variability was associated with better behavioral outcomes, suggesting a beneficial role of dynamic neural adaptation in visual processing [53], possibly reflecting enhanced discrimination of visual features across task condition.

In contrast, in left Crus I - a region previously implicated in higher-order cognitive processes crucial for managing the complex demands of bimanual coordination, particularly under cognitive load [4, 42, 51], such as cognitive control and working memory [34, 36, 45, 51, 52, 54–56] - higher BOLD variability modulation predicted poorer performance. One possible interpretation is that individuals who manage a broad range of motor demands more efficiently may rely on relatively stable cognitive control processes across task conditions, potentially reflecting a comparable level of automaticity [43], and consequently, reduced condition-specific modulation of variability in this region. Alternatively, increased modulation of BOLD variability in the left Crus I may be indicative of maladaptive recruitment of cognitive control processes during task execution [57]. These interpretations remain speculative and warrant confirmation in future studies specifically designed to test these relationships.

A key finding of our study is that BOLD SD modulation differentially influenced BTT performance across age groups. In younger adults, lower BOLD SD modulation in the right primary somatosensory cortex (A2), left superior parietal lobule (A7r), bilateral motion-sensitive visual areas (V5/MT+), and right cerebellar lobule VI was associated with better task outcomes. In contrast, in older adults, higher BOLD SD modulation in the left inferior and superior parietal lobules (A40rd and A7r) predicted better performance. This indicates that the neural dynamics supporting successful task execution may be age-specific, with younger and older adults recruiting and modulating distinct networks to meet task demands.

Drawing on recent literature [42, 49, 50, 57, 58], we propose that better performance in younger adults is linked to reduced BOLD variability modulation in sensorimotor and visuospatial regions because they likely possess well-established motor representations that require minimal recalibration through sensory feedback and motor planning, possibly due to greater task automatization. Conversely, in older adults, our results suggest better task performance is associated with greater adaptation of the neural dynamics in regions involved in executive control, motor attention and visuospatial integration [44, 59]. This potentially reflects the need for more active motor control in the absence of automatization. In this context, the observed increase in BOLD variability modulation across BTT conditions may represent heightened neural adaptability, enabling the dynamic integration of sensory input and top-down control processes necessary for real-time behavioral adjustments.

Existing cognitive aging models have been developed primarily based on mean BOLD activation measures, making it challenging to directly extrapolate their principles to BOLD variability findings. Nevertheless, we tentatively propose that the aforementioned age-specific shift in the relationship between brain dynamics and behavioral performance aligns with the STAC [10] and CRUNCH models [11], despite most group differences in BOLD SD modulation in the corresponding regions did not survive FDR correction (potentially due to the severity of the correction, see Figure 6). These frameworks posit that aging is accompanied by compensatory recruitment of additional higher-order neural resources to maintain task performance. Accordingly, in older adults, heightened up-regulation of BOLD SD with task complexity, particularly within the left inferior and superior parietal lobules (A40rd and A7r), may serve as a compensatory feedback-based strategy [42], potentially expressing increased reliance on visual processing and imagery [58], as well as strategic executive processes [3, 4, 7, 8, 59], to attain the same adaptive efficiency, as achieved by younger adults, through more automatized sensorimotor routines [42, 58]. Behaviorally, both age groups were similarly affected by the task conditions, strengthening our interpretation that successful adaptation to task demands is maintained through age-dependent neural strategies, characterized by specific profiles of BOLD variability modulation.

Limitations and future directions

Several limitations warrant consideration. First, the absence of direct measurements of vascular reactivity, WMH burden, or dopaminergic function limits our ability to attribute observed effects to specific physiological mechanisms. Without such data, interpretations regarding neural compensation, structural decline, or neurotransmitter shifts remain speculative.

Second, while our sample size was comparable to similar fMRI studies, the relatively small number of participants may have reduced our power to detect subtle effects, especially after correction for multiple comparisons.

Third, the use of a ROI-based approach, while justified statistically, may have masked finer-grained or network-level effects that could have emerged from a voxel-wise or connectivity-based analysis. Moreover, the choice of these regions based on activation maps may have overlooked regions where variability plays a prominent role despite minimal relevance of the average BOLD signal, given the reported orthogonality of mean and variability measures [16]. Therefore, an alternative ROI selection method, as well as whole-brain or network-level approaches should be considered in future work.

Furthermore, parametric task complexity levels were achieved by simultaneously manipulating two task components: speed of rotation and directional shift, making it difficult to determine how each is independently affected in aging, and their corresponding neural correlates. Similarly, the currently used BTT score does not allow to fully dissociate the behavioral component of speed from accuracy. Alternative BTT scoring methods may prove valuable in future research to better establish links between age-specific performance strategies and the modulation of BOLD signal variability.

Notably, although older adults exhibited greater head motion than younger adults during task performance, which can inflate BOLD SD estimation, we mitigated this unwanted effect by accounting for any post-preprocessing residual FD-related variability in our analyses, as advised by [60]. No additional correction was applied in the analyses concerning BOLD SD modulation, as these metric measures intra-individual variability across task conditions and we did not find significant differences in FD across task conditions that justified such correction. Moreover, given that FD variability might also index other confounds such as physiological noise (e.g., respiration, heart rate) [60, 61], by including it in our statistical models, we likely also minimized this source of variability. Nonetheless, although we implemented a variety of steps to control for these potential confounders, we cannot discard with full certainty the possibility of remaining residual effects on BOLD variability estimation.

Future research should aim to replicate and extend these findings by incorporating multimodal imaging techniques - such as diffusion tensor imaging, arterial spin labeling, positron emission tomography and magnetic resonance spectroscopy to better understand the vascular, structural, and neurochemical correlates of BOLD variability. Longitudinal designs will be particularly valuable in determining whether patterns of variability modulation represent dynamic states shaped by aging, experience, or intervention. Additionally, examining both cognitive and motor domains within the same participants may help clarify whether the neural dynamics of variability modulation are domain-specific or reflect broader, cross-functional signatures of aging while simultaneously shedding light onto the neural correlates of the extensively reported age-related penetration of cognition into action.

Participants

The study included 22 younger healthy adults (YA; 10 males, mean ± SD [range]: age 31.27 ± 5.42 [25–40] years) and 23 older healthy adults (OA; 12 males, mean ± SD [range]: age 69.52 ± 4.15 [65–79] years), all recruited from Flanders, Belgium. Participants were recruited via social media, as well as posters and flyers placed in locations relevant to the target population, such as universities, senior homes, sports facilities, and hospitals. All participants were right-handed based on the Edinburgh Handedness Inventory Laterality Quotient (EHI LQ; mean ± SD [range]: 94.87 ± 11.03% [50-100]; cutoff score ≥ 50) [62], were fluent Dutch speakers, and had normal or corrected-to-normal eyesight. Cognitive screening was conducted using the Montreal Cognitive Assessment (MoCA; mean ± SD [range]: 28 ± 1.32 [26-30], cutoff score ≥ 26/30) [63]. Depressive symptoms were assessed using age-specific neuropsychological instruments: the Beck’s Depression Inventory – Second Edition (BDI-II; median ± interquartile range [IQR; range]: 3 ± 4.75 [0-12], cutoff score > 13/20) [64] for younger adults, and the 15-item Geriatric Depression Scale (GDS-15; median ± IQR [range]: 1 ± 1 [0-12], cutoff score > 5/15) [65] for older adults.

Exclusion criteria included self-reported neuropsychiatric disorders (e.g., epilepsy or stroke), history of brain injury or surgery, substance abuse (alcohol or drugs), regular use of psychoactive medication (e.g., antidepressants, sedatives), compromised upper limb mobility, or contraindications for magnetic resonance (MR; e.g., cardiac pacemaker, aneurysm clip, metal prosthesis, claustrophobia). Additionally, individuals with a high level of bimanual proficiency (i.e., engaging in specific bimanual tasks for more than five hours per week, as in [66]), were excluded to ensure a comparable level of bimanual skills across participants at baseline.

Prior to enrolling, participants provided written informed consent. The study protocol was approved by the Ethics Committee Research of UZ/KU Leuven (study number S66028) and complied with the ethical principles regarding research with humans stipulated in the latest amendment of the Declaration of Helsinki [67].

Experimental protocol

The experimental protocol consisted of two sessions, each of about two hours. The first session was divided into two parts: (1) a comprehensive eligibility screening using standard neuropsychological tests and questionnaires, and (2) a familiarization session with the task used in the current study to assess bimanual coordination, the Bimanual Tracking Task (BTT).

For the screening, participants completed the MoCA [63], EHI [62], BDI-II [64], GDS-15 [65], and a questionnaire assessing MRI exclusion criteria. A summary of the group average scores at each test can be found in Table 7.

During the familiarization session, participants practiced the BTT for eight blocks or until they reached a minimum performance score of 15% in a dummy magnetic resonance (MR) scanner. Each practice block included the same movement patterns as those later performed in the real MR scanner (i.e., Line 1:1, Line 3:1 and Angle 3:1 [left:right hand movements]; see details in section 4.4), but each condition was presented in all four possible directional combinations of dial rotation (Figure 8). This session aimed to ensure participants understood the task instructions in a simulated MRI environment and to minimize the effects of early skill learning.

Figure 8. Experimental protocol. The figure shows the study protocol consisting of two sessions, each lasting approximately two hours. The first session included a screening and BTT familiarization in a dummy MR scanner. During this session, participants were introduced to the Bimanual Tracking Task (BTT) and performed three task conditions - Line 1:1, Line 3:1, and Angle 3:1 - in all possible directions. The task setup was placed over the participants' hips while they lay supine in the scanner. In the second session, participants underwent an MR protocol comprising several sequences. Each trial of the BTT during fMRI included a 1-second preparation phase, a 6-second execution phase, and a 1-second interstimulus interval. The same three BTT conditions were used during the fMRI scan, but only in a single direction, as shown in the figure. Abbreviations: BTT = Bimanual Tracking Task; fMRI = functional Magnetic Resonance Imaging; MR = Magnetic Resonance; MRS = Magnetic Resonance Spectroscopy; rs-fMRI = resting-state fMRI; T1 = T1-weighted.

In the second experimental session, participants completed an extensive imaging protocol consisting of a 7-minute high-resolution T1-weighted structural scan, a 9-minute resting-state functional MRI (fMRI) image, two 6-minute task-related fMRI images (namely, BTT and N-back), and several magnetic resonance spectroscopy (MRS) scans (which took approximately 50 minutes). MR data acquisition took place at the University Hospital of KU Leuven. For the current paper, only data concerning the structural and the functional BTT images were examined.

Imaging acquisition

All MRI data were collected using a Philips Achieva DStream 3-Tesla scanner equipped with a 32-channel receiver coil. Functional task-related data for the BTT were acquired using a T2*-weighted Echo Planar Imaging (EPI) sequence with the following parameters: repetition time (TR) = 1400 ms, echo time (TE) = 33 ms, flip angle = 70°, voxel size = 2.5 × 2.5 × 2.5 mm³, field of view (FOV) = 240 × 240 × 150 mm³, matrix size = 96 × 95, 60 transverse slices acquired in a sequential ascending order, SENSE/Multiband SENSE factors = 1.5/3, phase-encoding direction = posterior-anterior. A total of 258 dynamic scans were acquired, preceded by 7 dummy scans to allow for magnetic field stabilization. Additionally, high-resolution anatomical images were obtained using a T1-weighted (T1w) Magnetization Prepared Rapid Gradient Echo (MPRAGE) sequence with the following parameters: TR = 5.9 ms, TE = 2.6 ms, flip angle = 8°, voxel size = 0.8 × 0.8 × 0.8 mm³, FOV = 256 × 240 × 166.4 mm³, matrix size = 320 × 300, and 208 transverse slices.

fMRI task description

A nonferromagnetic version of the Bimanual Tracking Task (BTT) was used to assess visuomotor bimanual coordination during the task-related fMRI scan [68]. This setup, consisting of two handles and two rotational dials (5 cm diameter), was designed to be placed over the participants’ hips, allowing them to perform the task while lying supine in the MR scanner. The task was projected onto a screen positioned at the back of the scanner, and made visible to participants via a mirror (~14 × 9 cm).

BTT performance required the coordinated rotation of two dials using both index fingers, each positioned in a circular groove on the respective dial and controlling either horizontal (x-axis) or vertical movement (y-axis). The objective was to track a white target dot as it moved along a predetermined white trajectory at a constant speed. Participants’ real-time movement trajectory, generated by dial rotation, was displayed as a red line, providing concurrent performance feedback.

BTT conditions differed in inter-hand frequency ratios and/or rotational direction (i.e., clockwise or counterclockwise) of each dial, thereby varying in complexity. Frequency ratios refer to the relative speed of each index finger’s movement. In a 1:1 ratio, both fingers moved at the same speed (i.e., iso-frequency movement), whereas in a 3:1 ratio, the left finger moved three times faster than the right (i.e., non-iso-frequency movement). In this study, three task conditions were used: Line 1:1, Line 3:1, and Angle 3:1 (left:right hand movements). In “Line” conditions, both dials maintained the same rotational direction throughout the entire trial duration (clockwise in this experiment). In contrast, “Angle” conditions required a mid-trial change in the rotational direction of the dial controlling vertical movement (in this study, the left dial). These BTT conditions were chosen based on prior work from our research group that examined how different age groups coped with increasing task complexity using the same task [6, 42].

The BTT was performed during approximately six minutes in a block design. Each task condition was presented for five non-consecutive blocks (each 24 seconds), with three trials per block, totaling 45 trials (3 BTT conditions x 5 blocks x 3 trials/block = 45 trials). All trials had a constant duration of eight seconds. Each began with a 1-second preparation phase during which the target trajectory appeared on a black background, followed by a 6-second movement phase where the white dot moved along the target trajectory, and ended with a 1-second inter-trial interval (see Figure 8). The BTT conditions were presented in a fixed order for each participant but varied across participants in a pseudorandomized manner.

Behavioral outcome measure

Participant performance on the BTT was expressed in a percentage, with 100 % indicating perfect accuracy. The final BTT score (S) was derived from the multiplication of two key components (S = P × D): a preliminary score (P) reflecting movement accuracy and a distance correction factor (D) accounting for deviations from the target trajectory.

The preliminary score (P), calculated as described in [41], resulted from dividing the proportion of unique target points successfully covered in the correct order during the participant’s movement trajectory, by the total number of points forming the target trajectory, all multiplied by 100. These unique target points corresponded to those with the minimal Euclidean distance from the actual movement path, considering a sampling rate of 100 Hz (i.e., every 10 ms; see Figure 1D in [41]). Hence, a higher P value indicated greater accuracy in following the intended path. However, a limitation of this method was that it did not distinguish between movements that closely follow the target trajectory and those that are parallel to it.

To penalize movement trajectories parallel to the target trajectory over closely following it, a correction factor (D) was applied as previously done in studies from our research group [66, 69]. This factor adjusted the preliminary score based on the average distance ($\overline{d}$) between the actual movement path and the target trajectory throughout a trial. The correction was implemented using the formula:

$$D=\left(1-\frac{\overline{d}}{N}\right)$$

Unlike prior studies in which the scaling factor (N) was set to 5 [66, 69], in this study it was adjusted to 8 to better account for performance variability across different age groups. This adjusted N was determined by analyzing unpublished independent BTT data from our research group with a similar set of participants, and aimed to mitigate floor and ceiling effects, ensuring a more accurate assessment of skill differences.

Imaging preprocessing

Following acquisition, all MRI data were converted to BIDS format and preprocessed at the individual level using fMRIPrep 24.0.01 [70, 71], a robust and standardized preprocessing pipeline implemented in Nipype 1.8.6 [72]. The structural T1-weighted (T1w) image was first corrected for intensity non-uniformity (INU) using ANTs 2.5.1 (Advanced Normalization Tools) [73] and subsequently used as the T1w reference. Skull-stripping was performed using ANTs, and the image was segmented into cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM) using FAST from FSL 6.0.5.2 (FMRIB Software Library) [74]. FreeSurfer 7.3.2 [75] was then used to reconstruct brain surfaces. The T1w image was normalized to MNI152NLin2009cAsym space with an isotropic resolution of 2 mm using nonlinear registration with ANTs, utilizing the brain-extracted versions of both the T1w reference and the T1w template.

Functional image preprocessing began with motion correction of BOLD volumes using MCFLIRT [76] from FSL with a reference BOLD image. Susceptibility distortion correction was performed using a fieldmap-less distortion approach (i.e., SyN correction) [77], in which deformation fields were estimated by co-registering the EPI reference to the same-subject T1w reference with its intensity inverted. Functional images were then corrected for slice timing differences with 3DTshift from AFNI 24.1.22 (Analysis of Functional NeuroImages) [78], and the structural and functional images were co-registered using FreeSurfer with boundary-based registration (BBR) [79], which was implemented with 12 degrees of freedom. To address noise and confounding signals, nuisance regressors, including framewise displacement (FD; a measure of frame-to-frame head movement determined from the realignment estimates), DVARS (i.e., spatial standard deviation of the data after temporal differing), and global signals from WM and CSF, were identified using CompCor [80] and Nilearn [81].

After implementing the standard preprocessing steps in fMRIPrep, further temporal denoising was performed using the Denoiser Toolbox (https://github.com/arielletambini/denoiser). This process of temporal detrending involved regressing out the WM and CSF signals, as well as the three translational (trans_x, trans_y, trans_z) and three rotational (rot_x, rot_y, rot_z) head-movement parameters and their temporal derivatives. High-pass filtering was applied to the functional data with a cutoff of 0.006 Hz (~167 seconds) in the BTT fMRI. Finally, spatial smoothing was performed in SPM12 [82] with MATLAB R2022a (version 9.12.0.2009381) [83], using a Gaussian kernel of 6 mm full width at half maximum (FWHM).

Brain variability measure

After preprocessing, a whole-brain detrended standard deviation map of the BOLD series (BOLD SD) was computed for each BTT condition (i.e., Line 1:1, Line 3:1, and Angle 3:1) and for the total 9-minute rest scan using the Variability Toolbox (VarTbx; https://github.com/LNDG/vartbx) in SPM12 [82]. An additional whole brain BOLD SD map containing the absolute magnitude of change in BOLD SD across BTT conditions (whereafter referred to as BOLD SD modulation) was posteriorly generated by calculating, per voxel, the sum of the absolute differences in BOLD SD between the different task conditions, as described in the following mathematical formula:

$$\begin{array}{cc}Modulation& =\left|BOLDS{D}_{Line31}-BOLDS{D}_{Line11}\right|\\ & +\left|BOLDS{D}_{Angle31}-BOLDS{D}_{Line11}\right|\end{array}$$

Hence, this composite, between task condition measure was intended to capture the overall amount of change in BOLD SD across all BTT conditions, irrespective of direction (i.e., increase or decreases in BOLD SD). An illustrative example on how to calculate BOLD SD modulation considering the provided formula can be visualized on Supplementary Figure 1.

In this study, a region-of-interest (ROI)-based approach was chosen over a voxel-wise approach for all statistical analyses given the limited number of subjects per age group (i.e., 22 young adults and 23 older adults). In total, 25 ROIs were analyzed. These regions, listed in Table 4 and displayed in Figure 9, were selected based on a preliminary analysis in FSL, which examined brain activation during BTT performance using independent data from a similar study conducted by our research group. In particular, the included ROIs correspond to local maxima at least 8 mm apart, present in clusters of 10 or more significant voxels after applying a threshold of Z > 5.5 to the contrast Z-maps “Line 3:1 > Line 1:1” and “Angle 3:1 > Line 3:1” (even though no local maxima were found for the contrast map “Line 3:1 > Line 1:1”; thus, consult Table 4 for the coordinates of each local maxima found for the contrast “Angle 3:1 > Line 3:1”). To label the brain regions corresponding to the peak activation coordinates, two widely recognized atlases were used: the Brainnetome Atlas [38] for cortical regions, and the SUIT atlas [39] for cerebellar regions. Subsequently, masks derived from these atlases were created in the same T1w template space and isotropic resolution of 2 mm³ as the preprocessed functional images. The resulting masks were applied to the BOLD SD brain maps previously generated per task condition using MATLAB R2022a (version 9.12.0.2009381). For each of these BOLD SD maps, a mean BOLD SD value was extracted per ROI by averaging the BOLD SD across all voxels within the masked regions. Similarly, a mean value of BOLD SD modulation was calculated for each ROI.

Figure 9. Location of the brain regions analyzed. Brain regions analyzed, visualized in the coronal plane, in the MNI152 Template. Cortical regions were labeled using the Brainnetome Atlas [38], and cerebellar regions with Cerebellum SUIT Atlas [39]. Brain slices and volumes depicted in this figure were created with MRIcroGL (Version 1.2.20220720b; https://www.nitrc.org/projects/mricrogl) [37].

Statistical analyses

All statistical tests were conducted in R (version 4.4.2) [84], with the significance threshold (α) set at 0.05 unless stated otherwise. Prior to analysis, normality was assessed using the Shapiro-Wilk test for the data, and by visually inspecting quantile–quantile (Q-Q) plots for the model residuals. Homoscedasticity was examined using residual-by-predicted plots.

For Linear Mixed Models (LMMs) and Multiple Linear Regressions, model simplification followed a stepwise backward elimination approach, where non-significant terms were sequentially removed, unless they were part of a significant interaction. The significance of model terms was assessed using Type III analysis of variance (ANOVA). When appropriate, post hoc Tukey Honestly Significant Difference tests were performed. Moreover, before running the statistical models, per ROI, datapoints of BOLD SD or BOLD SD modulation more than three standard deviations from the mean were excluded.

To control for Type I error in region-of-interest (ROI) analyses involving multiple comparisons, false discovery rate (FDR) correction was applied using the Benjamini -Hochberg (BH) [85] procedure, with q < 0.05 (FDR threshold; q-values are FDR-adjusted p-values).

Behavioral data

Neuroimaging data

Age-group differences in BOLD variability and BOLD variability modulation during BTT performance