For each trial, we assessed two measures of performance: reaction time (RT), the interval between the extinction of the fixation stimulus at the end of the delay interval and the initiation of the MGS, as well as saccadic error (SE), the signed visual angle separating the horizontal location of the target and the end point of the MGS. For the four task conditions during a session, we computed average RT and the standard deviation of RT. As additional measures of mean behavioral performance and behavioral variability, we defined saccade inaccuracy as the absolute value of the average SE for a given target and saccade imprecision as the standard deviation of SE for each target.

Using a linear mixed-effects model to account for the longitudinal nature of the data, we found, consistent with prior findings (Luna et al., 2004), that average RT and inaccuracy decreases during development following an age⁻¹ trajectory (Figure 1c,e). We estimated how much each behavioral measure changed at the group level between the ages of 8 and 33 and observed a reduction in average RT amounting to approximately 239 ms, a 39% decline (t(1346)=9.30; p=6.03e-20 Inaccuracy decreased by approximately 0.25°, or 47% (t(1346)=3.11; p=0.0019).

Importantly, both measures of behavioral variability also decrease with development following an age⁻¹ trajectory (Figure 1d,f). The developmental change in the standard deviation of RT amounts to approximately 156 ms, a 67% reduction (t(1346)=7.84; p=9.4e-15). The imprecision of the MGS decreases with age, resulting in an estimated reduction of 0.87° of visual angle, a change, which amounts to roughly 35% (t(1346)=4.03; p=5.83e-5).

In order to determine whether behavioral variability provides additional information beyond mean behavioral measures about the developmental status of a subject, we compared the performance of linear models that predicted subject age from either mean RT or inaccuracy (null models) to matched linear models containing the corresponding behavioral variability factor (full models). We found that a null model predicting age from only mean RT was significantly improved by including the standard deviation of RT (null model: DF=4 AIC=−2212.7, Log-Likelihood=1110.4; full model: DF=5 AIC=−2220.7, Log-Likelihood=1115.6; p=0.002)(Burnham and Anderson, 2004). Including imprecision as a term improved model performance, but the difference did not achieve statistical significance (null model: DF=4, AIC=−2108.4, Log-Likelihood=1058.2; full model: DF=5 AIC=−2110.1, Log-Likelihood=1060.0, p=0.052). Thus, of the two measures of behavioral variability, reaction time variability, but not imprecision, appears to measurably reflect a unique aspect of cognitive development.

Performance of the MGS task requires coordinated activation of many brain regions involved in visuomotor, WM maintenance, and retrieval processes. Each of these processes is associated with a distinct whole-brain pattern of activity, or ‘brain state’, which is expressed at the appropriate times during a trial. Developmental changes in average global gain would result in changes in the average amplitude of expression of task-related brain states across trials; subjects with greater mean global gain would express the task-related brain states to a greater extent, while subjects with reduced global gain would express the brain states with reduced amplitude. We therefore sought to estimate the canonical brain state patterns associated with the visuomotor, maintenance, and retrieval processes involved in the MGS task and to determine whether adolescent development is associated with changes in the mean amplitude of their expression.

Using a simple dimensionality reduction technique based on linear regression (see Materials and methods), we constructed whole brain patterns of activity, that is, brain states, associated with visuomotor/encoding (VME), WM maintenance, and retrieval processes. These task-related brain state patterns were extracted from idealized time courses of BOLD activity (see Materials and methods) observed during the long delay trials (Figure 2a).

Figure 2

Download asset Open asset

To represent VME processes we extracted the average patterns of BOLD signal occurring 6 s after the initial ‘encoding’ saccade. To represent retrieval processes, we again extracted patterns of activity present 6 s after the MGS. This delay allowed the BOLD signals associated with these processes to reach their peaks (Handwerker et al., 2004). The pattern of activity associated with working memory maintenance was extracted from the TR immediately prior to the execution of the MGS, allowing us the purest estimate of delay period activation (furthest from the preceding visually-guided saccade, without intruding into the subsequent MGS). We orthogonalized each of the brain state patterns to ensure that they captured unique aspects of task activity by regressing the VME-related patterns from maintenance-related patterns, and regressing both VME- and maintenance-related patterns from the retrieval related patterns. This process removed remaining components of VME-related activity from the maintenance activity and importantly, allowed us to remove the pattern of activity associated with visuomotor responses from the retrieval-related pattern occurring during the MGS. Implicit in this procedure is the assumption that VME, maintenance, and retrieval processes are associated with distinct and consistent patterns of whole-brain BOLD activation that are expressed with the time course of a hemodynamic response.

Some neuronal gain modulators, particularly those acting through cholinergic pathways, alter gain with hemispheric specificity, similar to the effects of directed spatial attention (Sarter and Bruno, 1997; Furey et al., 2000; Bentley et al., 2004). To account for such potential hemispheric differences in gain signaling, we decomposed each brain state into two component patterns: a target hemifield specific ‘spatial’ component and a target hemifield non-specific, ‘mean’ component. Mean brain state components correspond to the average whole-brain patterns of activity associated with VME, maintenance, or retrieval, regardless of which visual hemifield the target was presented in (Figure 2b upper panels). Spatial components were constructed by computing the difference —right minus left— between brain state patterns determined separately for right and left side targets; they reflect the differences in activity during each task epoch resulting from the changes in the target’s visual hemifield (Figure 2b lower panel). The whole brain patterns of activity observed during different epochs of the task can therefore be approximated with linear combinations of the mean and spatial components of the brain state patterns. For instance, maintenance related activity observed during trials in which the target appears in the right visual hemifield, can be approximated by adding the mean and spatial components of the maintenance brain state patterns, while maintenance activity observed during trials in which the target appears in the left visual hemifield can be approximated by subtracting the spatial component of the maintenance brain state from the mean component.

Taken together, the brain states characterize the patterns of engagement of canonical regions underlying the VME epoch (e.g., frontal eye fields), maintenance (e.g., prefrontal and frontal eye fields) and the non-visuomotor aspects of retrieval and response (e.g., preSMA) (Anderson et al., 1994; Fried et al., 1991).

We verified that the resulting brain state components captured the relevant whole-brain patterns of BOLD signal associated with specific task epochs, by projecting each whole-brain volume of a subject’s average trial time course, onto the complete set of brain state components using linear regression (see Materials and methods). We performed this operation separately for each task condition. By temporally ordering the regression weights associated with each brain state component from each TR of the average trial time course, we converted the whole-brain time courses of task activity into time courses of task-related brain state expression (Figure 3). This procedure is conceptually similar to principle component/independent component analyses in which whole-brain voxel-wise time series are converted into time courses of expression of whole-brain patterns. Here, however, the individual components, that is, brain states, have known behavioral and cognitive processes with which they have been empirically associated. These brain state patterns, although derived from only the long delay trials, also served as an effective basis for describing the whole-brain patterns of BOLD signal evoked during the short delay trials as well (Figure 3; lower panel).

Figure 3

Download asset Open asset

As expected, given our approach for constructing the brain states patterns, the time courses of brain state expression for all subjects exhibited similar temporal characteristics. However, it remained possible that there might be age-related differences in mean global gain, which would affect the average amplitude of brain state expression. To determine whether such differences were present, we performed omnibus F-tests on the brain state time courses for each task condition to assess the null hypotheses that all of the coefficients for Age*TR (Age*TR*TargetHemifield in the case of the spatial brain state components) were equal to zero. We observed minimal age-related differences in the average time courses of brain state expression. These differences were observed for the spatial, but not mean, component of the VME states across all trial types (all p<0.001). We did not detect any significant age-related differences in the expression of either mean or spatial components of the maintenance state (all p>0.14). Results for age-related differences in the time courses of expression for the retrieval states were mixed: We observed no omnibus age-related differences within either of the long delay conditions, but within the long presentation short delay trials we observed a small age-related difference in the expression of the spatial retrieval state (F(42,14784)=1.6; p=0.012). Post-hoc examination of the individual Age*TR and Age*TR*TargetHemifield coefficients at each time point in the trial revealed that this effect was driven by slightly greater expression of the state by adults, across right and left side targets, during the 5th, 6th, and 8th TRs. However, in our post-hoc analysis no single time point reached significance (minimum(p)=0.055). We also observed that the mean retrieval state was differentially expressed across age within the short presentation short delay trials (F(40,14112)=2.0; p<0.001). Post-hoc analyses revealed that this effect was driven by a slightly greater expression of the mean retrieval state by adults during this condition during the 9th-13th TRs, well after the occurrence of peak expression for this state.

From visual inspection it is clear that adults exhibit a slightly prolonged expression of the spatial component of the VME brain state during the different trials (seen most prominently in Figure 3; lower panel). We also wanted to know whether the peak amplitude of spatial VME expression differed with age. From each session we examined the amplitude of peak expression of the spatial VME state for each trial type. Because the sign of expression of the spatial VME state varies depending on target hemifield, we extracted the maximum value of positive expression for right side trials, and we extracted the minimum value of expression for left side trials. If adults expressed the spatial VME state to a greater extent than children and adolescents due to greater average gain, this would result in greater positive expression for right side trials and reduced (more negative) peak expression during left side trials. We therefore examined the Age*Target interaction term, which we found did not reach significance (t(2696)=1.59; p=0.111).

Combined, these results demonstrate that the set of brain state patterns provide a simplified low dimensional basis for describing BOLD signal changes evoked by the memory-guided saccade task. Importantly, age-related differences in the expression of the brain state patterns during task performance were minimal, and only the spatial component of the VME state exhibited consistent age-related differences in expression across trial types. Even here, however, the age-related differences were not ones of amplitude, but of duration, suggesting that mean global gain does not change during adolescent development.

We hypothesized that behavioral performance is affected by fluctuations in global gain signals when they occur around the time of a behavioral response. The signature of variability in global gain signaling would be ‘brain state variability’, or momentary fluctuations in the amplitude of expression of whole brain states of activity associated with ongoing neural processes.

However, due to the delayed and prolonged nature of the BOLD response, the activity measured near the time of the MGS consists of multiple superposed states of activity associated with visuomotor, WM maintenance, and retrieval processes. Each of these processes may be affected differently by global gain variability and variability affecting each process might differentially affect behavioral performance. We therefore examined the trial-wise relationship between behavioral performance and the expression of the mean and spatial components of the VME, maintenance and retrieval brain state patterns in an interval of time centered on the occurrence of each MGS.

After removing the mean trial responses from each voxel, we looked for remaining patterns of activity across the brain that matched each of the canonical brain states. To accomplish this, we projected the whole-brain pattern of BOLD signal residual values at each TR onto the set of canonical brain state patterns using linear regression. This approach converted the whole-brain residual time series into a time series of task-related brain state fluctuations, and allowed us to determine the extent to which each brain state component was over- or under-expressed at a particular point in time.

We divided the resulting fluctuation time series for each brain state component into snippets, intervals of time centered on each MGS and extending ±15 TRs (each TR=1.5 s) before and after (Figure 4a and b). Each snippet was associated with the RT and SE of a particular trial. After aligning snippets from each trial to the TR in which the MGS occurred (TR 0), we used regression models to measure the relationship between both trial-to-trial differences in RT (z-scored) and SE (z-scored and rectified) and trial-to-trial differences in the expression of each brain state component at different times relative to the MGS. In order to measure the contribution of brain state variability to behavioral variability, we first estimated the fraction of trial-to-trial behavioral variance accounted for by a null model that included terms for non-neural trial-to-trial factors (e.g. target eccentricity) and performance factors (e.g. RT when testing SE and vice versa). Then we measured the additional fraction of behavioral variability that could be explained by including the measures of brain state expression from each TR in the snippets. The difference in explained variance between the full and null models indicates the fraction of trial-to-trial behavioral variability uniquely associated with brain state variability occurring at a particular time relative to the execution of a MGS (Figure 5a and b [top panels]).

Figure 4

Download asset Open asset

Figure 5

Download asset Open asset

As hypothesized, the relationship between brain state expression and both RT and SE peaked around the time that the MGS was executed. For trial-wise RT, brain state/behavior associations were significant beginning with the TR when the MGS was executed, peaking 1 TR after the saccade, and lasting for a total of 6 TRs. The relationship between trial-wise SE and brain state expression was similar, but much less prominent and evident during only the third TR following the MGS.

Trials with faster RTs were associated with greater early expression of the mean VME and maintenance brain states (TRs 0–2), and reduced later expression of all mean states (TRs 3–5). Greater SE was associated with greater expression of the mean VME state (TR 3).

RT and SE also covaried with fluctuations in the amplitude of expression of the spatial components of the brain state patterns (bottom three rows of Figure 4a and b). The fastest RTs occurred when whole-brain patterns of task-related activity were biased in the direction of the target (VME: TRs 1–3; maintenance and retrieval: TRs 1–2). That is, the fastest right side trials were those in which brain states were expressed most ‘rightwardly’, and the fastest left side trials were those in which the spatial brain states exhibited the greatest ‘leftward’ expression. Interestingly, greater target hemifield appropriate expression of the spatial VME state (TRs 2–3) was associated with increased SE, indicating that increasing the amplitude of its expression does not improve all aspects of task performance. Increased hemifield appropriate expression of the spatial component of the maintenance states (TR 3), in contrast, was associated with reduced SE.

That greater expression of VME brain states was associated with faster RT and greater SE, prompted us to examine the behavioral data for signs of a speed-accuracy trade-off. We found a significant quadratic relationship between z-scored RT and SE at the trial level indicating that, within a session, excessively fast and slow responses were associated with reduced accuracy (t(16754)=4.64; p=3.52e-6). Collectively, these results demonstrate that, trial-to-trial, both reaction time and accuracy of subjects’ responses covary with the amplitude of expression of whole-brain patterns of task related activity.

Greater early expression and reduced later expression of the mean brain states for fast RT trials (represented by the transition from blue to red in some rows of the lower panel of Figure 4a), could possibly be explained by trivial correlation between the timing of a saccade and the latency of the expression of the brain state. That is, on trials with longer RTs, the BOLD activity would be shifted later, causing the brain state to appear under-expressed early and over-expressed late relative to the average time course. To determine whether this relationship could account for the observed correlation with behavior, we performed a set of simulations to compare the temporal patterns of BOLD signal residuals for fast and slow RT trials that would result from timing-, amplitude-, and timing and amplitude-based relationships (see Materials and methods).

As depicted by Figure 5, the contributions of timing and amplitude variability to brain state variability are distinguished by their effect on the temporal structure of the residual brain state time series. A purely timing based explanation of the mean brain state/RT relationship predicts that the integral of the mean brain state residual time series, for both fast and slow reaction time trials should converge to zero (Figure 5, column 1). A relationship between reaction time and brain state expression mediated by fluctuations in the amplitude of expression of the brain state patterns predicts that the same time integrals converge to non-zero values of opposite sign (Figure 5, column 2). A combination of timing- and amplitude-based relationships predicts an initial bifurcation of the time integrals of the fast and slow reaction time residuals that then a partial re-convergence (Figure 5, column 3). We found that the trial-wise relationship between RT and the expression of the mean VME brain state —the state associated with eye-movements— was inconsistent with both a purely amplitude-based mechanism, and purely timing-based mechanism and instead is likely to reflect a mixture of both the trivial time shifted BOLD response, and a true relationship between gain fluctuation and performance (Figure 5, column 4).

If the stabilization of behavior during adolescence is related to a reduction in brain state variability, then the proportion of BOLD signal variability associated with brain state variability will decrease with age. We examined the subset of TRs (0–5; the highlighted region Figure 4a) around each correct MGS that showed a significant relationship with trial-to-trial behavioral performance and computed SS_brain, the sum of whole-brain squared error associated with all brain state patterns across each TR as well as SS_error, the sum of the remaining squared error. For each session, we computed the ratio SS_brain/(SS_brain + SS_error), producing values we refer to here as total brain state variability, which corresponds to the fraction of residual whole-brain BOLD signal variability associated with brain state variability. We found that total brain state variability decreases with age (Figure 6a) after controlling for mean frame-wise displacement (FD)(Power et al., 2012; Siegel et al., 2014) and the number of correctly performed trials (t(333)=-3.35; p=9.0e-4).

Figure 6

Download asset Open asset

The measure of total brain state variability, SS_brain, can be simply decomposed into a sum of contributions from the mean and spatial components of VME, maintenance, and retrieval brain state patterns; SS_brain = SS_VME + S_SMaint + SS_Retrieval. We computed the ratio of each brain state component sum of squared error to (SS_brain + SS_error), and found that the stability of each of the three sets of brain state components exhibit distinct developmental trajectories (Figure 6b). The proportion of VME-related brain state variability did not significantly decrease with age (t(333)=-1.47; p=0.14). However, maintenance- and retrieval-related brain state variability both show significant age-related decreases (t(333)=-3.8; p=1.8e-4 and t(333)=-5.53; p=6.27e-8 respectively). Pairwise comparisons of slopes reveal that trajectories of VME and retrieval-related variability can not be significantly distinguished from one another (t(670)=0.60; p=0.54), while maintenance-related brain state variability decreases more slowly through adolescence compared to both VME and retrieval-related variability (t(670)= 2.44; p=0.014 and t(670)=-3.71; p=2.2e-4 respectively).

Next, we sought to determine whether a systematic relationship exists between in-scanner motion and measures of brain state variability. We reasoned that if brain state variability was unrelated to movement-related artifacts, then our finding that brain state variability was reduced in older subjects would still hold if we selectively sub-sampled our data so that we compared a group of adults who moved excessively to a group of children who moved relatively little. The biased subsampling routine that we employed is based on a mean-matching technique that has been described in detail elsewhere (Churchland et al., 2010). Briefly, we divided our data into two sets, split at the median age of our sample. We selectively drew samples from the two data sets such that, on average the older sample exhibited greater mean FD than the younger group. We found that reversing the relationship between motion and age does not significantly alter our finding that older subjects exhibited less brain state variability than younger subjects (Figure 7).

Figure 7

Download asset Open asset

Our analyses to this point have shown that trial-to-trial differences in behavioral performance are associated with brain state variability, which stabilizes during adolescent development. If the stabilization of behavior during adolescence is the result of reduced brain state variability, then subjects exhibiting the greatest longitudinal reduction in behavioral variability should be those that exhibit the greatest decline in brain state variability. To determine whether this is the case, we leveraged the longitudinal design of our dataset to examine how individual differences in the developmental trajectories of RT variability and saccade imprecision are related to individual developmental trajectories of brain state variability.

We selected a subset of 29 subjects for whom we had at least four complete sessions of data, and estimated individual developmental slopes of total brain state variability. We modeled developmental changes in total brain state variability as a linear effect of age after observing superior performance compared to a model in which total brain state variability was fit with a age^-1 term (simulated likelihood ratio test with 5000 iterations, age⁻¹ model: DF=9 AIC=−3116.2, Log-Likelihood=15567.1; linear age model: DF=9 AIC=−3120.2, Log-Likelihood=1569.1; p=0.015);. We also estimated individual regression weights for the developmental trajectories of reaction time variability and imprecision after controlling for task condition and mean reaction time. Here we used an age⁻¹ term to model individual trajectories.

Subjects exhibiting the greatest decreases in behavioral variability given their age have the greatest age⁻¹ regression coefficients. In contrast, subjects exhibiting the greatest decrease in total brain state variability have the smallest age regression coefficients. Evidence consistent with the hypothesis that the developmental stabilization of behavior is driven by a reduction in brain state variability would be the existence of a negative correlation between the brain state and behavioral regression coefficients. We observed just such a negative relationship for reaction time variability (r = −0.48; p=0.008) (Figure 8a). This result remained significant (p<0.05) when we limited our data set to subjects with 3–5 or more sessions as well. However, the within-subject relationship between brain state variability and saccade imprecision was not significant (r=0.28; p=0.14) (Figure 8b).

Figure 8

Download asset Open asset

Given the modest relationship between brain state variability and SE at the single trial level, we considered the possibility that we might be underpowered to detect the longitudinal relationship between total brain state variability and saccadic imprecision in our smaller sample size. We expanded our analyses to investigate whether individual differences in saccade imprecision were related to individual differences in brain state variability using our entire sample and including age⁻¹ terms as covariates. Here we did observe a significant positive relationship between total brain state variability and saccade imprecision (t(1344)=3.2; p=0.001) (Figure 8c).

We tested 152 subjects between the ages of 8 and 33. Subjects were initially recruited between the ages of 8 and 30 years and were scanned approximately annually for 1–10 years. Subjects were included based on two criteria: (1) Mean frame-wise displacement (FD) was less than 0.15 mm; and (2) at least 50% of the trials from each of the four trial types had to be measurably correct. Here, incorrect trials are those for which measurements of reaction time and endpoints for both visually- and memory-guided saccades were unavailable due to blink artifacts, noisy data, or transient loss of pupil- or corneal reflection-lock. After applying these exclusion parameters our dataset consisted of 126 subjects (60 female). We applied no further outlier control for our analyses. Participants and/or their legal guardians provided informed consent before participating in this study. Experimental procedures for this study complied with the Code of Ethics of the World Medical Association (1964; Declaration of Helsinki) and the Institutional Review Board at the University of Pittsburgh. Subjects were paid for their participation in the study.

Eye-movements were recorded in the scanner with an infrared camera system equipped with long-range optics and sampling at 60 Hz (Model R-LRO6, Applied Science Laboratories, Bedford MA). Subjects’ compliance with instructions was assessed and eye-movements were monitored via remote video during task performance. We used a nine-point calibration procedure to estimate the transformation from the eye-tracker's native encoding space to on-screen pixel location. Saccadic events were detected using an in-house suite of automated routines. Individual saccade candidate events were detected from local maxima in the eye-movement velocity trace. Saccade start and end times were determined by searching backward and forward in time in the velocity trace to find the sample where velocity dropped below 1/10th of the peak velocity (Gitelman, 2002).

Imaging data were acquired using a Siemens 3-Tesla MAGNETOM Allegra (Erlangen, Germany) system with a standard radio-frequency (RF) head coil at the Brain Imaging Research Center, University of Pittsburgh, Pittsburgh, PA. Structural images were acquired using a sagittal magnetization prepared rapid gradient echo (MPRAGE) T1-weighted pulse sequence with 224 slices with 0.7825 mm slice resolution. Functional images were acquired using a gradient echo echo-planar (EPI) sequence sensitive to blood-oxygen-dependent (BOLD) contrast (T2*) (TR=1.5 s, TE=25 ms, flip angle=70°, voxel size=3.125×3.125×4 mm slice resolution, 229 volumes). Twenty-nine slices per volume were collected with no gap and aligned to the anterior and posterior commissure (AC-PC) plane. Structural and functional fMRI data are available through the Dryad repository (Montez et al., 2017).

T1-weighted anatomical images were reconstructed from raw DICOM files and converted to NIFTI format. We estimated the bias field corrections using smoothed and highpass filtered anatomical data analyzed with FSLs fast algorithm. After bias field correction we constructed a skull stripped anatomical data set for the subject, which we used to estimate the 12 degree-of-freedom affine transformations that would align the subjects data with the MNI152 anatomical template. Finally, we computed the non-linear transformation that would bring the subject’s affine-aligned anatomical data set into registration with the MNI152 template. We saved final combined linear/non-linear transformation for later use in registering the subjects’ functional data to the standard space.

fMRI data were preprocessed using a combination of AFNI (Analysis of Functional NeuroImages, RRID:SCR_005927) and FSL software (FSL, RRID:SCR_002823). In our pre-processing pipeline, raw data was converted from DICOM format to NIFTI volumes and slice-timing correction was applied using AFNI tools. We performed motion estimation and correction in two phases. First we pre-aligned each frame of a subject’s functional data to a volume created by taking the temporal mean of the 4-D functional time series. Then, a second, ‘true’, average functional volume was computed from the pre-aligned functional data, producing a reference functional volume that was less affected by motion artifacts. We then aligned each frame of the original function time series to this second reference volume using sinc-function interpolation and estimating the time course of translational and rotational motion throughout the run. We used these estimated time series throughout our later analyses of the functional data.

Next, using FSL’s brain extraction tool, we stripped the skull and superfluous tissues from the subject’s motion corrected mean functional EPI images, afterward aligning the resulting mean EPI volumes to their anatomical MPRAGE volume using a six degree-of-freedom rigid-body transformation estimated using spline interpolation. To align each frame of the motion corrected EPI sequence to the subjects structural image, we applied the translation estimated in the previous step to each frame of the motion corrected functional time series and then removed the skull and extraneous tissues from each frame of the functional time series. Tissue remaining within the mean functional volume after the skull stripping procedure was removed by applying a dilated binary mask to the mean aligned functional volume that removed extreme voxels whose values did not reside in middle 98th percentile. We then removed voxel-wise temporal extrema using AFNI’s 3dDespike software.

To align a subjects functional data to a standard MNI152 (Montreal Neurological Institute; MNI) template in a single transformation step, we used FSL convertwarp, and applywarp functions to combine the estimated motion correction, functional-to-structural, and linear and non-linear subject-to-MNI152 transformations into a single operator, which we applied separately to each frame of the original slice time-corrected functional data.

We performed minimal spatial smoothing on the aligned functional data, using a SUSAN algorithm with a 5 mm FWHM kernel, followed by a conservative high-pass filtering of the voxel-wise time series, which removed or attenuated BOLD signal frequencies below 0.0083 Hz (corresponding to fewer than three cycles per task run). Finally, we rescaled all voxel values by a value defined to be 10,000 divided by the global median.

In our analyses, measuring brain state variability requires determining the spatial structure of canonical whole-brain patterns of BOLD signal associated with distinct task-related processes. The flowchart in Figure 9 depicts an outline of the processing steps used to transform individual preprocessed fMRI time series into the average time courses of brain state expression and brain state variability as well as the relationship between the major processing steps and certain key analyses.

Figure 9

Download asset Open asset

From each session's data we extracted eight voxel-wise average time courses of BOLD activity corresponding to each of the four task conditions when stimulus presentations occurred in wither the left or right visual hemifields. We estimated these time courses with a finite impulse response (FIR) regression model. FIR design matrices were constructed manually and applied to the voxel-wise time series using 3dDeconvolve (AFNI). All trials, including incorrect responses and blinks, for each stimulus type were modeled over an interval consisting of the duration (from initial stimulus presentation to the execution of the memory-guided saccade) plus an additional 22.5 s (15 TRs). The design matrix included nuisance regressors to account for the effects of signal drift, subject motion, and global signal changes as captured by white matter and cerebrospinal fluid (CSF) signals and their derivatives. Signal drift for each run was modeled as a third order Legendre polynomial time series.

Head motion was computed along six affine components corresponding to translation in the three cardinal directions and rotations about three orthogonal axes. In addition we computed a time course of total displacement for each session based on the Euclidean norm of the time derivative of the movement time series at each time point. To account for the prolonged effect of autocorrelated movement on the BOLD signal, we included temporally leading (−1TR) and lagging (+1–2 TRs) copies of each of the seven motion regressors (Friston et al., 1996). Each of the seven motion time courses therefore contributed four motion regressors to the deconvolution design matrix. After deconvolution, we scaled the resulting whole-brain average trial time courses at each voxel, normalizing them to the standard deviation of the regression residuals at the same voxel location.

Idealized voxel-wise trial time courses for the long delay conditions were estimated from the scaled average trial time course estimates. We modeled these separately for each condition and target hemifield using 3dLME (AFNI), a linear mixed-effects framework. Each time point was modeled as a separate categorical fixed effect and we did not include an intercept term in the model. To account for any bias due to the over representation of subjects who participated in more scans, we included subject identity as a random effect component in the regression model. For each trial type we computed the total displacement undergone by each subject's brain during the BOLD signal measurement intervals (trial durations plus 15 TRs) and included it as a fixed-effect component in the regression analysis. We calculated subjects' average age for all of their sessions and, after centering by the global mean, included it as a subject level fixed-effects regressor. We included a mean age by time interaction term to capture age-related differences in the voxel-wise time courses. We included the subjects' age at each session, after subject level mean-centering, as a second age-related random-effects regressor. Within a given voxel, a single whole trial time course may include independent contributions from visually- and memory-guided saccade events. To account for potential differences due to variability in the number of correct saccades, we included the proportion of unclassifiable and incorrect visually- and memory-guided saccades and their interactions with time as fixed-effect components of the model. We produced idealized trial time courses by generating the voxel-wise model estimates for a subject of mean age, mean in-scanner displacement, and perfect trial performance. This process generated four idealized whole-brain time series corresponding to both long-delay conditions in which targets were located in either the left or right visual hemifield. We used these idealized BOLD time series in our construction of the canonical brain states

Including lagged motion regressors during deconvolution accounts for the prolonged linear effects of in-scanner motion up to 2 TRs after a given movement (Friston et al., 1996). To control for any effects of motion that may continue beyond that time and bias measures of brain state variability, we developed a method to account for temporally prolonged linear motion artifacts from fMRI data based on motion template volumes that model the spatial pattern of artifacts associated with the linear effects of motion (Figure 10a). We normalized the regression coefficients associated with each lagged motion regressor at each voxel by the temporal standard deviation of the voxel’s post-deconvolution residuals and computed their resulting mean spatial patterns across all subjects. We used these whole-brain patterns of normalized regression coefficients to construct motion artifact templates. For each of the 28 templates (7 motion components with 4 temporal lead/lags), we subtracted the spatial mean of all voxel values and scaled the resulting volumes to a common vector magnitude. Then, using principle component decomposition, we found a set of 11 motion templates that captured >90% of the variability in the set, which were then converted back into 3D volumes.

Figure 10

Download asset Open asset

We verified that the set of motion templates captured components of linear motion related BOLD signal variability that remained after deconvolution by computing SS_motion, the sum of squared error associated with all motion templates across all TRs in the whole-brain residual BOLD time series. We then computed the ratio SS_motion/(SS_brain + SS_error) where SS_brain is the sum of squared error associated with brain state variability and SS_error is unexplained sum of squared error. We refer to this ratio as motion template variance. We found that estimates of mean FD and motion template variance are highly related (t(337)=10.3; p=1.06e-21), indicating that even with rigorous motion controls during deconvolution (see Materials and methods), significant linear motion artifacts remain (Figure 10b). We found that while brain state variability is not significantly related to mean FD (t(337)=1.39; p=0.17), it is significantly greater in subjects exhibiting greater motion template variability (t(337)=2.01; p=0.045), (c-d). The usage of motion templates as an additional control for motion-related artifacts arose from an abundance of caution and the need for high dimensional nuisance regressors that could be fit simultaneously to each volume alongside the set of brain state patterns. While these preliminary analyses suggest promise for the approach, we acknowledge that it has yet to be fully validated.

Although the relationship between brain state variability and movement was small and inconsistent across different estimates of motion, we sought a more rigorous control by comparing brain state variability in a group of high motion adults to low motion children and adolescents (see Results).

To remove trivial components of whole-brain BOLD signal variability associated with TR-to-TR shifts in the global mean across the entire imaged volume, we included a constant offset template, which consisted of a whole-brain binary mask. To account for TR-to-TR shifts in mean BOLD signal that are limited to the gray matter (the sort of variability which in some resting state studies is removed via global gray matter signal regression) we included a second binary template defined by the brain state mask (see Materials and methods Estimating canonical brain states).

In addition, we constructed a set of spatial gradient templates to account for other trivial modes of whole-brain BOLD signal variability consisting of simple linear spatial gradient patterns. We first created 3 spatial gradient volumes whose voxel values were equal to their x, y, and z coordinates relative to the volume's center of mass. We set each voxel that fell outside of a whole-brain MNI mask to zero. We then computed a set of 3 ‘interaction’ templates that corresponded to each pair-wise product of the spatial gradient templates. As a whole, this set of 19 templates constituted the set of spatial nuisance regressors that we used to capture and remove remaining unwanted spatial modes of whole-brain BOLD signal variability associated with motion and trivial variability in global signal.

We included the motion templates and spatial constant and gradient templates as additional nuisance regressors in all analyses that involved projecting whole brain volumes of BOLD signal onto the canonical set of brain state components.

We derived the canonical brain states from whole-brain patterns of BOLD activity extracted from different time points within the idealized trial time courses. To constrain the brain state patterns to swathes of gray matter that were reliably imaged across sessions, we constructed a mask defined by the intersection of a probabilistic MNI gray-matter mask (threshold≥0.5) and a mask constructed from the group average of the model R-squared maps from each session's initial deconvolution (threshold≥0.27, selected to removed unreliably imaged brain regions). Combined these masks served to eliminate white matter; CSF; the hind-brain; cerebellum, the inferior portion of which was not consistently imaged across sessions; artifact-prone basal forebrain; and infero-temporal lobes). The extent of the resulting brain state mask can be observed in (Figure 2b).

Each canonical brain state consists of a mean and spatial component. The mean component represents the average pattern of whole-brain BOLD activity evoked by a specific task epoch, regardless of target location. The spatial component represents the difference between the patterns of whole-brain BOLD activity evoked by specific task epochs for right side versus left side targets. To reduce the extent to which BOLD signal resulting from a trial's initial visually-guided saccade contaminated our estimate of maintenance- and retrieval-related brain states, we maximized the temporal distance between the task epochs by using time courses from long delay trials only.

The canonical mean and spatial VME brain states were derived from the whole-brain pattern of activity occurring within the brain state mask around the time of the visually guided saccade. To account for hemodynamic lag, we extracted four volumes of BOLD activity (corresponding to the left and right side targets from both long delay conditions) from the idealized long delay time courses six seconds (4 TRs) after the visually guided saccade was performed to ensure that most regions would have achieved their peak BOLD response (Handwerker et al., 2004). We constructed the mean VME brain state by computing the voxel-wise average of normalized BOLD activity across all four volumes. To construct the spatial VME state we separately computed the voxel-wise average of the VME volumes for right and left side targets. The spatial state was defined as the difference –right minus left– of the resulting volumes.

Canonical maintenance-related brain states were based on the patterns of normalized BOLD signal that occurred during the TR immediately prior to a subjects’ execution of the memory-guided saccade. Since brain states were estimated from the long delay trial types only, the maintenance brain state patterns were naturally separated from the BOLD response evoked by the initial stimulus presentation and the subjects' subsequent visually guided saccade. To completely remove any remaining stimulus and visuomotor contributions, as well as global average signal, we extracted all volumes from the idealized whole-brain trial time courses that occurred before or included a trial's initial visually guided saccade and regressed these patterns from the maintenance state patterns. States were then constructed for maintenance intervals as described for the VME states above.

Construction of the retrieval-related mean and spatial brain states followed a similar course as that used in the construction of the VME and maintenance-related brain states. Like the VME states, retrieval-related states were based on the normalized BOLD responses occurring during the 4^th TR after the execution of the memory-guided saccade. We removed all components of VME and maintenance activity by regressing out every pattern of activity within the brain state mask that occurred before the TR at which the subject began a memory-guided saccade.

The resulting mean and spatial brain states, as expected, exhibited a high degree of mirror symmetry across the midline plane. However, we observed noise and artifacts were introduced into the brain state patterns by accumulating error in our regression-based orthogonalization procedure. To correct for these artifacts, we leveraged symmetry of the brain state patterns by averaging them with corresponding brain state patterns derived from a left/right mirrored version of the idealized BOLD time series. To do this, we constructed left-right mirrored versions of the idealized trial time courses, performed an identical set of operations to define each brain state and combined them with the original set of brain state patterns. We first applied a mirror matrix to the idealized trial time course volumes. We aligned the mirrored volumes to the standard MNI template using a 12-degrees of freedom affine transformation. We then applied a non-linear transformation that warped the mirrored and aligned volumes to the standard MNI template. The final mean brain state components were constructed by computing the voxel-wise average of the mirrored and non-mirrored mean brain states. We inverted the sign of all voxel values of the mirrored spatial volumes before averaging them with the non-mirrored spatial volumes to produce the final spatial brain state components. In practice the mirroring procedure had only a minor effect on our results, visibly removing noise and improving the relationship that we observed between trial-to-trial behavioral performance and brain state expression (Figure 4).

We converted the average whole-brain trial time series and whole-brain residual time series into average time courses of brain state expression and variability respectively. For each TR, we extracted the whole-brain pattern of activity that we then vectorized and modeled using a linear regression. Our design matrix consisted of vectorized versions of the six brain states (the mean and spatial components of the VME, maintenance and retrieval states) as well as the 19 nuisance regressors templates described above. For each TR we extracted the regression weights for the six brain states, motion, and nuisance components and ordered them into a time series. When this procedure is performed on the whole-brain average trial time series, the result is a time course of expression of each of the brain states during a trial. When performed on the whole-brain residual time series, the result is a time course of brain state fluctuations, where positive values indicate that a particular brain state was present to a greater extent than average and negative values indicate that a state was expressed less than average. For each session, temporally z-scored the time series of variability for each brain state component.

For each session we separately transformed reaction time and saccadic error from each of the four main task conditions into z-scores. SE was rectified such that high SE values reflect greater error in memory-guided saccadic endpoints on a trial. We excluded all trials for which measurements of reaction time and endpoints for both visually- and memory-guided saccades were unavailable due to blink artifacts, noisy data, or transient loss of pupil- or corneal reflection-lock.

We related trial-to-trial variability in reaction time and accuracy to variability in the expression of each brain state across a range of times (±15 TRs) relative to the TRs that contained the memory-guided saccades from each trial (Figure 4). Using all correct trials across all sessions, we extracted our z-scored measurements of brain state fluctuation derived from the whole-brain residual time series. We then constructed a regression model that included terms for the measured values of each brain state at the relative TR. We also included terms for the spatial brain state interaction with target hemifield.

Each model contained terms that varied across trials but did not vary across relative TRs. These included terms for run number (coded as 1–3), target hemifield (coded as −1,1), target location (eccentricity, coded from least to most eccentric as 1–3), and the square of target location term. Because it is possible that RT and SE are correlated on a trial-to-trial basis, a true relationship between brain state variability and RT may result in a trivial relationship between brain state variability and SE, or vice versa. To account for this possibility, in the RT regression models, we included a term for trial-to-trial SE and its square. Similarly, for the SE regression models we included a term for RT and its square. This set of regressor terms served as a null model against which the full brain state model was compared. The trial-wise reaction time and accuracy models were fit using a linear mixed-effects framework (MATLAB) to account for the different numbers repeated measurements for many of the subjects. Subject identity was modeled as a random effect. We used the difference between the ordinary R² values for full and null models at each relative TR to assess the amount of unique behavioral variability accounted for by trial-to-trial fluctuations in the expression of different brain states. At each relative TR we compared the null and full modes using simulated maximum-likelihood estimation procedure with 5000 iterations (MATLAB).

To perform these simulations, we compiled a distribution of reaction times for all correct memory-guided saccades across our subject database. We then simulated 400 trials with reaction times selected to produce a distribution that was identical to the compiled empirical distribution. To simulate the simple timing-based effect of BOLD signal variability we generated an impulse function for each simulated trial, a vector where all but one element is equal to zero, where each consecutive element refers to 60 ms time bin after the extinction of the central fixation cross (the signal to perform a memory-guided saccade). For each draw from the reaction time distribution we generated new impulse function vector by inserting a 1 into the vector at the index corresponding to the reaction time on that trial. We convolved each of the 400 impulse functions with a canonical HRF modeled at the same 60 ms resolution. This produced a set of HRF time series whose time of peak amplitude varied with reaction time

Amplitude-based simulations were performed similarly but with two key differences: (1) for each trial we inserted 1 into all impulse function vectors at the same time index, corresponding to mean reaction time, for all trials; and (2) we added or subtracted from the 1 a linearly interpolated value between ±0.25 where+0.25 corresponded to the fastest reaction time and −0.25 corresponded to the slowest reaction time. Mixed amplitude and timing based simulations were a hybrid of the two described above. The index of the 1 for each trial's impulse function was selected to coincide with the reaction time on that trial. An additional amplitude modulation factor, as above, was added to the impulse index.

Separately for timing-, amplitude-, and timing and amplitude-based simulations, we computed the mean HRF time series across all trials and simulated residual time series by subtracting the mean HRF time series from the individual trial time series. Next we divided the simulated residuals in to fast and slow RT sets, defined by median split and calculated their average. Lastly, we computed the time integral of the mean residual time series for fast and slow trials.

To compare the simulated pattern of high temporal resolution residuals to the actual data we extracted equivalent snippets (beginning at the TR containing the MGS) of the mean VME brain state fluctuation time course for all trials and all subjects. We selected the mean VME state because of its close relationship with the visuomotor processes which makes it most likely to reflect a trivial timing based relationship between RT and brain state expression. As in the simulated data, for each subject we divided the snippets in to fast and slow RT trials based on a median split and then calculated the group average time series. Finally, we interpolated the resulting time series to a matched temporal resolution using shape preserving piece-wise cubic interpolation (MATLAB).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Human subjects: Participants and/or their legal guardians provided informed consent before participating in this study. Experimental procedures for this study complied with the Code of Ethics of the World Medical Association (1964; Declaration of Helsinki) and the Institutional Review Board at the University of Pittsburgh. Subjects were paid for their participation in the study.

1. Received: January 31, 2017
2. Accepted: August 3, 2017
3. Accepted Manuscript published: August 19, 2017 (version 1)
4. Version of Record published: August 31, 2017 (version 2)

© 2017, Montez et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.