Brain state flexibility accompanies motor-skill acquisition.

Learning requires the traversal of inherently distinct cognitive states to produce behavioral adaptation. Yet, tools to explicitly measure these states with non-invasive imaging - and to assess their dynamics during learning - remain limited. Here, we describe an approach based on a distinct application of graph theory in which points in time are represented by network nodes, and similarities in brain states between two different time points are represented as network edges. We use a graph-based clustering technique to identify clusters of time points representing canonical brain states, and to assess the manner in which the brain moves from one state to another as learning progresses. We observe the presence of two primary states characterized by either high activation in sensorimotor cortex or high activation in a frontal-subcortical system. Flexible switching among these primary states and other less common states becomes more frequent as learning progresses, and is inversely correlated with individual differences in learning rate. These results are consistent with the notion that the development of automaticity is associated with a greater freedom to use cognitive resources for other processes. Taken together, our work offers new insights into the constrained, low dimensional nature of brain dynamics characteristic of early learning, which give way to less constrained, high-dimensional dynamics in later learning.

Introduction

The human brain is an inherently adaptive (Mattar et al., 2016), plastic (Della-Maggiore et al., 2015) organ. Its fundamental malleability supports changes to its architecture and function that are advantageous to human survival. Importantly, such changes can occur on multiple time scales: from the long time scales of evolution (Kirschner and Gerhart, 1998; Clune et al., 2013) to the shorter time scales of multi-year development (Gu et al., 2015b), or even short-term learning (Ellefsen et al., 2015; Hermundstad et al., 2011). Notably, even in the shortest time scales of learning, adaptation can occur over multiple spatial scales (Mattar and Bassett, 2016), from the level of single neurons (Richardson et al., 2012) to the level of large-scale systems (Tunik et al., 2007). Moreover, this adaptation can affect functional dynamics (Heitger et al., 2012; Krakauer et al., 2005, 2004; Grefkes et al., 2004) or can evoke a direct change in the structure of neuroanatomy, driving new dendritic spines (Xu et al., 2009), axon collaterals (Chklovskii et al., 2004), and myelination (Sampaio-Baptista et al., 2013).

Malleability, adaptability, and plasticity often manifest as a variability in quantitative statistics that describe the structure or function of a system. In the large-scale human brain, such statistics can include measures of neurophysiological noise (Garrett et al., 2014, 2013; Breakspear and McIntosh, 2011) or changes in patterns of resting state functional connectivity (Deco et al., 2009, 2011, 2013). More recently, dynamic reconfiguration of putative functional modules in the brain – groups of functionally connected areas identified using community detection algorithms (Porter et al., 2009; Fortunato, 2010) – has been used to define a notion of network flexibility (Bassett et al., 2011b), which differs across individuals and is correlated with individual differences in learning (Bassett et al., 2011b), cognitive flexibility (Braun et al., 2015b), and executive function (Braun et al., 2015b).

Indeed, in the context of motor skill learning, dynamic network techniques have proven to be particularly advantageous for longitudinal designs, where data is collected from the same participants at multiple time points interspersed throughout the learning process (Bassett et al., 2013b, 2015; Wymbs and Grafton, 2015). Using a 6-week longitudinal design where participants trained motor sequences while undergoing functional magnetic resonance imaging, motor sequence learning was found to be associated with both increasing and decreasing motor system activity, with sequence-specific representations varying across multiple distinct timescales (Wymbs and Grafton, 2015). With a network modeling approach based on coherent activity between brain regions, the same dataset revealed the existence of a core-periphery structure that changes over the course of training and predicts individual differences in learning success (Bassett et al., 2013b). More recently, these changes were shown to reflect a growing autonomy between sensory and motor cortices, and the release of cognitive control hubs in frontal and cingulate cortices (Bassett et al., 2015). Yet despite these promising advances, dynamic network reconfiguration metrics are fundamentally unable to assess changes in the patterns of activity that are characteristic of brain dynamics, as they require the computation of functional connectivity estimates over extended time windows (Telesford et al., 2016; Bassett et al., 2013a).

To overcome this weakness, we use an alternative technique inspired by network science to identify temporal activation patterns and to assess their flexibility (Medaglia et al., 2018; Chen et al., 2016). Leveraging the same longitudinal dataset from the above studies, we begin by defining a brain state as a pattern of regional activity – for instance, estimated from functional magnetic resonance imaging (fMRI) – at a single time point (Gu et al., 2015a; Liu et al., 2013; Liu and Duyn, 2013) (see Fig. 1). Time points with similar activity patterns are then algorithmically clustered using a graph-based clustering technique (Porter et al., 2009; Fortunato, 2010), producing sets of similar brain states. Finally, by focusing on the transitions from one state to another, we estimate the rate of switching between states. This approach is similar to techniques being concurrently developed in the graph signal processing literature (Huang et al., 2016; Goldsberry et al., 2016), and it allows us to ask how activation patterns in the brain change as a function of learning. We address this question in the context of the explicit acquisition of a novel motor-visual skill, which is a quintessential learning process studied in both humans and animal models. As participants practice the task, we hypothesize that the brain traverses canonical states differently, that characteristics of this traversal predict individual differences in learning, and that the canonical states themselves are inherently different in early versus late learning.

Fig. 1

Schematic Depicting Construction of Adjacency Matrices

(a) Blood-oxygen-level-dependent (BOLD) signal from functional magnetic resonance imaging (fMRI) data was acquired from healthy adult subjects. (b) We calculated the mean BOLD magnitudes in each of 112 cortical and subcortical regions as a function of time. (c) The regional time series is represented in matrix format, and (d) the correlation between matrix columns (TRs) is used to create a time-by-time adjacency matrix. The ij element of this matrix measures the similarity between the regional pattern of BOLD magnitude between TR i and TR j. Adjacency matrices representing time-by-time networks form the fundamental data structure on which community detection algorithms function. We maximize a modularity quality function informed by these matrices to extract network communities: groups of TRs that show similar regional patterns of BOLD magnitudes. (e) Due to the near-degeneracy of the modularity landscape, this procedure is repeated 100 times per matrix. (f) Across these 100 partitions of TRs (nodes) into groups (communities), we construct a representative or “consensus” partition (g) that summarizes the significant structure in the original matrix.

To test these hypotheses, 20 healthy adult human participants practiced a set of ten-element motor sequences. To assess the change in brain activity related to behavioral change, we acquired BOLD time series and defined a brain state to be a pattern of BOLD magnitudes across regions at each time point. We then quantified the similarities in brain states across time and used network-based clustering algorithms to find recurrent brain states independent of their temporal order (Fig. 1e–f). We observed three to five "brain states", with two anti-correlated "primary states" occurring more frequently than the rest. We also observed that "state flexibility" – the flexible switching among all brain states – increases with task practice, being largely driven by contributions from brain regions traditionally associated with task learning and memory. Moreover, individuals with higher state flexibility learned faster than individuals with less switching between brain states. These results demonstrate that the global pattern of brain activity offers important insights into neurophysiological dynamics supporting adaptive behavior, underscoring the utility of a state-based assessment of whole brain dynamics in understanding higher order cognitive functions such as learning.

Materials and methods

Experiment and data collection

Ethics statement

In accordance with the guidelines set out by the Institutional Review Board of the University of California, Santa Barbara, twenty-two right-handed participants (13 females and 9 males) volunteered to participate and provided informed consent in writing. Separate analyses of the data acquired in this study are reported elsewhere (Bassett et al., 2013b, 2015; Wymbs and Grafton, 2015).

Experimental Setup and Procedure

Head motion was calculated for each subject as mean relative volume-to-volume displacement. Two participants were excluded from the following analyses: one failed to complete the entirety of the experiment and the other had persistent head motion greater than 5 mm during MRI scanning. The 20 remaining participants all had normal or corrected vision and none had any history of neurological disease or psychiatric disorders. In total, each participant completed at least 30 behavioral training sessions over the course of 6 weeks, a pre-training fMRI session, and three test fMRI sessions. The training used a module that was installed on the participant's laptop by an experimenter. Participants were given instructions on how to use the module and were required to train at minimum ten days out of each of 3 fourteen day periods. Training began immediately after the pre-training fMRI session and test scans were conducted approximately fourteen days after each previous scan (during which training also took place). Thus a total of 4 scans were acquired over the approximately 6 weeks of training.

Training and trial procedure

Participants practiced a set of ten-element sequences in a discrete sequence-production (DSP) task, which required participants to generate these responses to visual stimuli by pressing a button on a laptop keyboard with their right hand (see Fig. 2). Sequences were represented by a horizontal array of five square stimuli, where the thumb corresponded to the leftmost stimuli and the pinky corresponded to the rightmost stimuli. The imperative stimulus was highlighted in red and the next square to be pressed in the sequence was highlighted immediately after a correct key press. The sequence only continued once the appropriate key was pressed. Participants had an unlimited amount of time to complete each trial, and were encouraged to remain accurate rather than swift.

Fig. 2

Trial Structure

(a) Each trial began with the presentation of a sequence-identity cue that remained on the screen for 2 s. Each of the six trained sequences was paired with a unique identity cue. A discrete sequence-production (DSP) event structure was used to guide sequence production. The onset of the initial DSP stimulus (thick square, colored red in the task) served as the imperative to produce the sequence. A correct key press led to the immediate presentation of the next DSP stimulus (and so on) until the tenelement sequence was correctly executed. Participants received a "+" as feedback to signal that a sequence was completed and to wait (approximately 0–6 s) for the start of the next trial. This waiting period was called the "intertrial interval" (ITI). At any point, if an incorrect key was hit, a participant would receive an error signal (not shown in the figure), and the DSP sequence would pause until the correct response was received. (b) There was direct S-R mapping between a conventional keyboard or an MRI-compatible button box (lower left) and a participant's right hand, so that the leftmost DSP stimulus cued the thumb and the rightmost stimulus cued the pinky finger. Note that the button location for the thumb was positioned to the lower left for maximum comfort and ease of motion.

Each participant trained on the same set of six different ten element sequences, with three different levels of exposure: extensively trained (EXT) sequences that were practiced for 64 trials each, moderately trained (MOD) sequences that were practiced for 10 trials each, and minimally trained (MIN) sequences that were practiced for 1 trial each. Sequences included neither repetitions ("11", for example) nor patterns such as trills ("121", for example) or runs ("123", for example). All trials began with a sequence-identity cue, which informed participants which sequence they would have to type. Each identity cue was associated with only a single sequence and was composed of a unique shape and color combination. EXT sequences, for example, were indicated by a cyan or magenta circle, MOD sequences by a red or green triangle, and MIN sequences by orange or white stars. Participants reported no difficulty viewing the identity cues. After every set of ten trials, participants were given feedback about the number of error-free sequences produced and the mean time to produce an error-free sequence.

Each test session in a laboratory environment was completed after approximately ten home training sessions (over the course of fourteen days) and each participant took part in three test sessions, not including the pre-training session, which was identical to the training sessions. To familiarize the participants with the task, we introduced the mapping between the fingers and DSP stimuli and explained each of the identity cues prior to the pre-training session.

As each participant's training environment at home was different than the testing environment, arrangements were made to ease the transition to the testing environment (see Fig. 2 for the key layout during testing). Padding was placed under the participants' knees for comfort and participants were given a fiber optic response box with a configuration of buttons resembling that of the typical laptop used in training. For example, the distance between the centers of buttons in the top row was 20 mm (similar to the 20 mm between the "G" and "H" keys on a Mac-Book Pro) and the distance between the top row and lower left button was 32 mm (similar to the 37 mm between the "G" and spacebar keys on a MacBook Pro). The position of the box itself was adjustable to accommodate participants' different reaches and hand sizes. In addition, padding was placed both under the right forearm to reduce strain during the task and also between the participant and head coil of the MRI scanner to minimize head motion.

Participants were tested on the same DSP task that they practiced at home, and, as in the training sessions, participants were given unlimited time to complete the trials with a focus on maintaining accuracy and responding quickly. Once a trial was completed, participants were notified with a "+" which remained on their screen until the next sequence-identity cue was presented. All sequences were presented with the same frequency to ensure a sufficient number of events for each type. Participants were given the same feedback after every ten trials as they were in training sessions. Each set of ten trials (referred to hereafter as trial blocks) belonged to a single exposure type (EXT, MOD, or MIN) and had five trials for each sequence, which were separated by an inter-trial interval that lasted between 0 s and 6 s. Each epoch was composed of six blocks (60 trials) with 20 trials for each exposure and each test session contained five epochs and thus 300 trials. Participants had a variable number of brain scans depending on how quickly they completed the tasks. However, the number of trials performed was the same for all participants, with the exception of two abbreviated sessions resulting from technical problems. In both cases, participants had only completed four out of five scan runs for that session when scanning was stopped. Data from these sessions are included in this study.

Behavioral apparatus

The modules on participants’ laptop computers were used to control stimulus presentation. These laptops were running Octave 3.2.4 along with the Psychtoolbox version 3. Test sessions were controlled using a laptop running MATLAB version 7.1 (Mathworks, Natick, MA). Key-press responses and response times were measured using a custom fiber optic button box and transducer connected via a serial port (button box, HHSC-1×4-l; transducer, fORP932; Current Designs, Philadelphia, PA).

Behavioral estimates of learning

We defined movement time (MT) as the time between the first button press and the last button press for any single sequence. For the sequences of a single type, we fit a double exponential function to the MT (Schmidt and Lee, 1988; Rosenbaum, 2009) data in order to estimate learning rate. We used robust outlier correction in MATLAB (through the "fit.m" function in the Curve Fitting Toolbox with option "Robust" and type "Lar"): MT = D1e− + D2e−, where t is time, κ is the exponential drop-off parameter used to describe the fast rate of improvement (which we called the learning rate), λ is the exponential drop-off parameter used to describe the slow, sustained rate of improvement, and D1 and D2 are real and positive constants. The magnitude of κ determines the shape of the learning curve, where individuals with larger κ values have a steeper drop-off in MT and thus are thought to be quicker learners (see Fig. 3) (Dayan and Cohen, 2011; Snoddy, 1926). This decrease in MT has been an accepted indicator of learning for several decades (Heathcote et al., 2000) and various forms have been tried for the fit of MT (Newell and Rosenbloom, 1981; Heathcote et al., 2000), with variants of an exponential model being the most statistically robust choices. Importantly, this approach is also not dependent on an individual's initial performance or performance ceiling.

Fig. 3

Exemplar Learning Curves

To quantify learning rate, we fit a double exponential function to the movement time, defined as the time between the first button press and the last button press of a given sequence, as a function of trial number (see Methods). Fits are shown for the median learner, the fastest learner, and the slowest learner in the cohort.

fMRI imaging

Imaging procedures

Signals were acquired using a 3.0 T Siemens Trio with a 12-channel phased-array head coil. Each whole-brain scan epoch was created using a single-shot echo planar imaging sequence that was sensitive to BOLD contrast to acquire 37 slices per repetition time (repetition time (TR) of 2000 ms, 3 mm thickness, 0.5 mm gap) with an echo time of 30 ms, a flip angle of 90°, a field of view of 192 mm, and a 64×64 acquisition matrix. Before the first round of data collection, we acquired a high-resolution T1-weighted sagittal sequence image of the whole brain (TR of 15.0 ms, echo time of 4.2 ms, flip angle of 90°, 3D acquisition, field of view of 256 mm, slice thickness of 0.89 mm, and 256×256 acquisition matrix).

fMRI data preprocessing

Imaging data was processed and analyzed using Statistic Parametric Mapping (SPM8, Wellcome Trust Center for Neuroimaging and University College London, UK). We first realigned raw functional data, then coregistered it to the native T1 (normalized to the MNI-152 template with a resliced resolution of 3×3×3 mm), and then smoothed it with an isotropic Gaussian kernel of 8-mm full width at half-maximum. To control for fluctuations in signal intensity, we normalized the global intensity across all functional volumes. Using this pipeline of standard realignment, coregistration, normalization, and smoothing, we were able to correct for motion effects due to volume-to-volume fluctuations relative to the first volume in a scan run. The global signal was not regressed out of the voxel time series, given its controversial application to resting-state fMRI data (Murphy et al., 2009; Saad et al., 2012; Chai et al., 2012) and the lack of evidence of its utility in analysis of task-based fMRI data. Furthermore, the functional connectivity matrices that we produce showed no evidence of strong global functional correlations but instead showed discrete organization in motor, visual and non-motor, non-visual areas (Bassett et al., 2015).

General linear model

We performed a standard GLM analysis to quantify the degree to which brain regions showed a linear decrease in activity over the course of training. For each subject, the BOLD response was modeled using a single design matrix with parameters estimated using the GLM. We used an event-related design to model the expression of sequence-specific representations, with trial onset corresponding to the presentation of the sequence identity cue, 2 s before the presentation of the initial DSP target stimulus. This approach included both the preparation and the production of learned sequences. We constructed the design matrix for each subject using separate factors for each scan session (pretraining session and training sessions 13), exposure condition (MIN, MOD and EXT) and repetition (new or repeated trial). For a trial to be coded as a repeated event, the previous trial had to have been (i) of the exact same sequence and (ii) performed correctly. Repeated trials that followed error trials, as well as the error trials themselves, were modeled using a separate column in the design matrix. To account for nonspecific effects of sessions, blocking variables were included for each scan run. Potential differences in BOLD values due to MT-related kinematics were accounted for by using the MT from each trial as the trial duration for modeled events (Grinband et al., 2008). To control for the potential influence of the time between trials, we weighted each event by the time elapsed since the previous trial. Following center mean normalization, this column was added to the model as a covariate of non-interest. Events were convolved using the canonical hemodynamic response function and temporal derivative in SPM8. Using freely available software (Steffener et al., 2010), we then combined corresponding beta image pairs for each event type (hemodynamic response function and temporal derivative) at the voxel level to form a magnitude image (Calhoun et al., 2004). To generate linear contrast images at the individual-subject level, we multiplied magnitude images corresponding to conditions of interest by the appropriate contrast weight and then combined them through addition. We then averaged the beta weights of all voxels within a given region.

Network construction and analysis

Partitioning into regions of interest

We divided the brain into regions based on a standardized atlas (Bassett and Bullmore, 2006; Braun et al., 2015a; Bullmore and Sporns, 2009). There exist a number of atlases and the decision of which to use has been the topic of several recent studies on structural (Bassett et al., 2011a; Zalesky et al., 2010), resting-state (Wang et al., 2009), and task-based network architectures (Power et al., 2011). Consistent with prior graph-based studies of task-based fMRI (Bassett et al., 2011b, 2013b, 2014; Mantzaris et al., 2013), we divided the brain into 112 cortical and subcortical regions using the Harvard-Oxford (HO) atlas of the FMRIB (Oxford Centre for Functional Magnetic Resonance Imaging of the Brain) Software Library (Macke et al., 2011; Woolrich et al., 2009). For each participant and for each of the 112 regions, the regional mean BOLD was computed by separately averaging across all voxels in that area (see Fig. 1a and b).

Wavelet decomposition

Historically, wavelet decomposition has been applied to fMRI data (Bullmore et al., 2004, 2003) to detect small signal changes in nonstationary time series with noisy backgrounds (Brammer, 1998). Here, we use the maximum-overlap discrete wavelet transform, which has been used extensively (Achard et al., 2006, 2008; Bassett et al., 2006; Achard and Bullmore, 2007; Bassett et al., 2009; Lynall et al., 2010) to decompose regional time series into wavelet scales corresponding to specific frequency bands (Percival and Walden, 2006). Because our sampling frequency was 2 s (1 TR), wavelet scale 1 corresponded to 0.125–0.25 Hz, and scale 2 to 0.06–0.125 Hz. To enhance sensitivity to task-related changes in BOLD magnitudes (Sun et al., 2004), we examined wavelet scale 2, consistent with our previous work (Braun et al., 2015b; Bassett et al., 2013b, 2014). For a lengthier discussion of methodological considerations, see (Zhang et al., 2016).

Constructing Time-by-Time Networks

We were interested in studying the similarities between brain states as individuals learn. We defined a brain state as a pattern of BOLD activity across brain regions at a single instant in time (Gu et al., 2015a; Liu et al., 2013; Liu and Duyn, 2013). Using the wavelet decomposition as the input signal, we measured the similarities between these states in each trial block which was comprised of approximately 40–60 repetition times (TRs). We calculated the Spearman correlation of regional BOLD magnitudes between all possible pairs of time points (TRs). This procedure creates an undirected, weighted graph or network in which nodes represent time points and edges between nodes represent the correlation between brain states at different time points (see Fig. 1 c & d). Intuitively, this matrix – which we refer to as a “time-by-time” network provides the necessary information to uncover common brain states (Medaglia et al., 2018; Chen et al., 2016), and to study transitions between brain states, as a participant learns.

Isolating Brain States Using Community Detection

To uncover common brain states in the “time-by-time” network, we used a network-based clustering technique known as community detection (Fortunato, 2010; Porter et al., 2009). In particular, we chose a common community detection approach known as modularity maximization, where we optimize the following modularity quality function (Newman, 2006) using a Louvain-like (Blondel et al., 2008) locally greedy heuristic algorithm (Bassett et al., 2013a): where A is the time-by-time matrix, times i and j are assigned, respectively, to community g and g, the Kronecker delta δ(g) = 1 if g = g (and zero otherwise), γ is the structural resolution parameter, and P is the expected weight of the edge between regions i and j under some null model. Consistent with prior work (Mattar et al., 2015; Bassett et al., 2011b, 2013b, 2015; Braun et al., 2015b; Cole et al., 2014; Mantzaris et al., 2013), we used the Newman-Girvan null model (Girvan and Newman, 2002): where is the strength of region i and . Importantly, the algorithm we use is a heuristic that implements a non-deterministic optimization (Good et al., 2010). Consequently we repeated the optimization 100 times (Bassett et al., 2013a), and we report results summarized over those iterations by building what are known as consensus partitions (Bassett et al., 2015) (see Fig. 1). In order to do this, we construct a nodal association matrix A from a set of N partitions, where A is equal to the number of times in the N partitions that node i and node j are in the same community. Furthermore, we construct a null nodal association matrix A, constructed from random permutations of the N partitions. This null association matrix indicates the number of times any two nodes will be assigned to the same community by chance. We then create the thresholded matrix A by setting any element A that is less than the corresponding null element to 0. This procedure removes random noise from the nodal association matrix A. Subsequently, we use a Louvain-like method to obtain N new partitions of A into communities, where each of the N partitions is typically identical, and each of which is a consensus partition of the N original partitions.

Recurrent Brain States

Each community obtained in the aforementioned pipeline includes a set of TRs that show similar patterns of regional BOLD magnitudes, and could thus be interpreted as representing a single, repeated brain state in a single trial bock. We first sought to aggregate these brain states over trial blocks. To this end, we average the pattern of regional BOLD magnitudes across all TRs assigned to that community in that trial block. We then repeat community detection across all representative brain states found in the trial blocks to find sets of representative brain states for each subject at each scan. By averaging the pattern of BOLD magnitudes of the brain states in each set, we find a group of representative brain state for every subject at every scan.

Second, we sought to aggregate these subject-scan representative brain states over all scans to identify a group of representative brain states for each scan. We thus repeat community detection over the set of all subject-scan representative brain states, separated by scan, and we again average the pattern of regional BOLD magnitudes across all subject-scan representative brain states assigned to the same community. This final set of brain-states we consider to be scan-representative brain states for each scan of learning.

Finally, we sought to find analogous communities in each scan. Therefore, we repeated the community detection algorithm for these communities and interpreted two scan-representative brain states assigned to the same community as analogous. In summary, we repeatedly use this brain state isolation procedure hierarchically to first isolate representative brain states for each subject-scan combination, then for each scan, and finally to find brain states in each scan that are similar to one another.

State Flexibility Calculation

To examine how the pattern of traversals through brain states changed, we define a "state flexibility" metric (Medaglia et al., 2018). We specified state flexibility (F) to be the number of state transitions (T) observed relative to the number of states (S), or . Intuitively, state flexibility is a measure of the volatility versus rigidity in brain dynamics, directly representing the frequency of dynamic state changes.

Results

Time by time network analysis identifies frontal and motor states

Our first goal was to characterize the average anatomical distribution of BOLD magnitudes across all subjects and scans, to better understand the whole-brain activation patterns accompanying motor skill learning. To achieve this goal, we create a time-by-time network where nodes represent individual time points, and edges represent the Spearman correlation coefficient between the vector of regional BOLD magnitudes at time point i and at time point j (see Fig. 1). We represented the time-by-time network as a graph. From these graphs, we were able to find 3 recurrent brain states, of which, two were strongly anti-correlated (Pearson correlation coefficient r(446) = −0.4291, p = 1.6951 × 10−21). These anti-correlated states make up 95.67% of all time subjects spent learning, and are also the only states to be present in all scans. We therefore refer to these states as "primary states" and focus our analysis upon them. We refer to the first state as the "motor state," characterized by strong activation of the extended motor system and anterior cingulate, as well as simultaneous deactivation of the medial primary visual cortex (see Fig. 4 A, Table 1). We refer to the second state as the "frontal state," characterized by strong activation of a distributed set of regions in frontal and temporal cortices, as well as subcortical structures (see Fig. 4 B, Table 1).

Fig. 4

Brain States Common Across Learning

We show the average activation magnitude of brain areas of the primary states. We refer to the first state shown in panel a as the "motor state" due to the strong activation of the extended motor system. In this state, we also observe high activity magnitudes in anterior cingulate, and low activity magnitudes in primary visual cortex along the medial wall. We refer to the second state shown in panel b as the "frontal state" due to the strong activation of frontal areas. In this state, we also observe high activation magnitudes in temporal cortex and subcortical structures, and low activity magnitudes in the extended sensorimotor system. In panel (c), we show the average BOLD signal across all regions in each scan and note that it does not significantly change (F(3, 669) = 1.17, p = 0.3221). In d we show that the primary state proportion does not significantly change across scans (F(3, 57) = 0.17, p = 0.9163). See the supplement for subcortical slices of both the motor and frontal state.

Table 1

Twenty regions with the greatest BOLD activity magnitude in the two primary brain states with corresponding MNI coordinates.

Motor State				Frontal State
Name	X	Y	Z	Name	X	Y	Z
R, supplementary motor area	7	−2	58	L, frontal medial cortex	−6	42	−18
L, supplementary motor area	−7	−2	56	L, caudate	−13	9	10
L, postcentral gyrus	−37	−28	53	R, frontal medial cortex	7	42	−18
L, superior parietal lobule	−29	−50	57	L, parahippocampal gyrus, posterior	−22	−34	−15
R, planum temporale	51	−26	12	L, hippocampus	−22	−5	−18
L, supramarginal gyrus, anterior	−56	−33	37	L, subcallosal cortex	−5	19	−14
L, precentral gyrus	34	−10	49	R, caudate	13	10	11
R, supramarginal gyrus, anterior	58	−27	37	R, parahippocampal gyrus, posterior	23	−31	−16
R, precentral gyrus	34	−10	49	R, subcallosal cortex	6	19	−14
R, superior parietal lobule	29	−48	59	R, hippocampus	21	−4	18
R, parietal operculum cortex	46	−28	23	R, parahippocampal gyrus, anterior	21	−7	−31
L, Heschl's gyrus	−45	−21	6	R, nucleus accumbens	9	12	−7
L, parietal operculum cortex	−45	−33	22	L, middle temporal gyrus, anterior	−59	−4	−22
L, globus pallidus	−19	−5	−1	L, parahippocampal gyrus, anterior	−22	−9	−30
R, central opercular cortex	48	−6	13	L, paracingulate gyrus	−8	36	22
R, supramarginal gyrus, posterior	53	−41	33	R, lingual gyrus	14	−63	−7
R, Heschl's gyrus	46	−18	6	L, inferior frontal gryus, pars triangularis	−48	28	9
L, central opercular cortex	−46	−8	14	L, parahippocampal gyrus (superior to 34,35)	−24	−23	−14
L, planum temporale	−49	−31	11	L, nucleus accumbens	−9	11	−7
R, postcentral gyrus	37	−27	53	R, middle temporal gyrus, anterior	59	−1	−26

While these two states were statistically present across the entire experiment, we did observe small fluctuations in the magnitudes of the regional activity of both states. Thus, natural questions to ask are (i) did either state became stronger or weaker with training? and (ii) did the frequency of primary states change with learning? To address the first question, we calculated the mean BOLD magnitude among all brain regions for each state. Using a repeated measures ANOVA, we found no significant differences among scans in either state (F(3, 669) = 1.17, p = 0.3221, see Fig. 4). This suggests that the activation of these two states did not significantly change – on average – with the level of training. To address the second question, we calculated the proportion of time spent in primary states in each scan. Using a repeated measures ANOVA, we found no significant differences among scans (F(3, 57) = 0.17, p = 0.9163). This suggests that the frequency of primary states remained the same during learning.

State flexibility increases with task practice

How does the brain traverse these states? Do individuals’ traversals change with learning? To examine how the pattern of traversals through brain states changes during learning, we defined a "state flexibility" metric. Following (Medaglia et al., 2018), we specified state flexibility (F) to be the number of state transitions (T) observed relative to the number of states (S), or . Intuitively, state flexibility is a measure of the volatility versus rigidity in brain dynamics, directly representing the frequency of dynamic state changes. We observed that state flexibility increased monotonically with the number of trials practiced (repeated measures ANOVA: F(9, 171) = 9.97, p = 3.04 × 10−12, see Fig. 5 right panels) and that differences between the sequences did not exist in the pre-training session (repeated measures ANOVA: F(2, 59) = 0.03, p = 0.97, see Fig. 5 left panels). This suggests that as subjects learned the sequences, regional patterns of BOLD magnitudes became more variable, indicating more frequent transitions between different brain states.

Fig. 5

Whole-Brain State Flexibility

(a) Right Panel: State flexibility shows a clear increasing trend as subjects complete more trials of the experiment (F(9, 171) = 9.97, p = 3.04 × 10−12 using a repeated measures ANOVA). The 50 trial box corresponds to the pre-training session for each of the three types while all other boxes correspond to a single scan from each of the three training sets. Left Panel: State flexibility in pre-training scans for sequences that would become extensively trained (EXT), sequences that would become moderately trained (MOD), and sequences that would become minimally trained (MIN). There was no effect of sequence type on state flexibility. (b) Right Panel: State flexibility shows a mild decreasing trend in the null model due to the reduced number of TRs (F(9, 171) = 2.6, p = 0.0078; null model was created by permuting the adjacency matrix while maintaining symmetry). The 50 trial box corresponds to the pre-training session for each of the three types while all other boxes correspond to a single scan from each of the three training sets. Left Panel: State flexibility in pre-training scans for sequences that would become extensively trained (EXT), sequences that would become moderately trained (MOD), and sequences that would become minimally trained (MIN). These results indicate that we find no effect of sequence type prior to training.

An important question to consider is whether this change in state flexibility is related to the length of time that participants take to complete the practice trials. Specifically, because the experiment is self-paced, the length of time to complete the sequences decreased as participants practiced; subjects became quicker with experience. To ensure that the length of time to complete a sequence was not driving the observed changes in state flexibility, we constructed a non-parametric permutation-based null model by permuting the adjacency matrix A uniformly at random while maintaining symmetry. Critically, this null model displayed a decrease in state flexibility with number of trials practiced (F(9, 171) = 2.6, p = 0.0078, Fig. 5), suggesting that neither the reduced length of time nor the correlation values themselves can explain the observed increase in state flexibility, but rather that the temporal structure of the data is required for the observed increase in state flexibility.

Regional contributions to state flexibility vary by function

Do regions contribute differentially to state flexibility? To answer this question we conducted a "lesioning" analysis, where we calculated state flexibility for each subject and scan while "lesioning out," or excluding, a single region. We then calculated the average difference between the true state flexibility for that subject and scan, and the lesioned state flexibility for each region across all trial groups. We normalized these values by subtracting off the mean effect of lesioning on flexibility.

To assess statistical significance, we created a matrix of the contributions to state flexibility for all regions and subjects. We found the contribution from seven regions to be significant (p < 0.05 for each region, df = 19) by computing a t-test between the ablated state flexibilities and the true state flexibility, while correcting for multiple comparisons across the 112 brain regions using a false discovery rate (α=0.05). By calculating the average of these contributions for each region, we identify (i) negative contributors, the removal of which increases state flexibility, and (ii) positive contributors, the removal of which decreases state flexibility. We find that the significant negative contributors to state flexibility are associated generally with motor and visual function (supplementary motor area, cuneus cortex, and the postcentral gyrus). In contrast, the significant positive contributors to state flexibility are associated with more integrative processing in hetermodal association areas (temporal occipital fusiform cortex, and planum polare on the temporoparietal junction) (see Fig. 6, Table 2).

Fig. 6

Schematic Depicting Ablasion Calculation and Regional Contribution to State Flexibility

(a) Blood-oxygen-level-dependent (BOLD) signal from functional magnetic resonance imaging (fMRI) data was acquired from healthy adult subjects. (b) Single regions were selected and "ablated" or ignored in the BOLD time-series. (c) State flexibility was recalculated for these new matrices and the difference between true state flexibility and the ablated state flexibility is defined to be the regional contribution to state flexibility for the ablated region. (d) We show the difference between the true state flexibility and ablated state flexibility for each region. We note that regions in the motor and visual cortex tend to have negative contributions to state flexibility while regions in the hetermodal association areas tend to have positive contributions to state flexibility.

Table 2

Significant contributors to state flexibility with corresponding MNI coordinates.

Most Negative Contributors				Most Positive Contributors
Name	X	Y	Z	Name	X	Y	Z
R, supplementary motor area	7	−2	58	L, planum polare	−45	−8	−9
L, superior parietal lobule	−29	−50	57	R, planum polare	46	−6	−8
L, postcentral gyrus	−37	−28	53	R, temporal occipital fusiform cortex	36	−50	−18
L, supramarginal gyrus anterior	−56	−33	37

State flexibility is correlated with learning rate

The results thus far indicate that state flexibility is an important global feature of brain dynamics that significantly changes as individuals learn a new motor-visual skill. Yet, they do not address the question of how such brain dynamics relate directly to changes in behavior. Therefore, we next asked the question: Are individual differences in state flexibility related to individual differences in learning rate? Here, we estimate the correlation of the learning rate (see Methods for definition) with the differences in state flexibility between sessions. All correlations are estimated using a linear mixed effects model that accounts for the effect of both subject and scan.

We observe a significant positive correlation between state flexibility and learning rate (ρ = 0.24, p = 0.041, Fig. 7), accounting for the effects of subject. That is, subjects will tend to be inherently better or worse at learning than other subjects. Therefore, we normalize for inter-subject differences in learning rate, and we find a significant correlation between state flexibility difference and learning rate. State flexibility difference is positively correlated with learning rate: intuitively, larger decreases in flexibility are associated with better learning. Furthermore, we observe a significant correlation between state flexibility difference and learning rate (p = 0.045) when accounting for the effects of scan. Learning rate tends to decrease as the number of scans increases. Therefore, we build this trend into our model, and find a significant correlation between state flexibility and learning rate. Extending our previous assertion, these results suggest that both the individual differences and the larger patterns of change are correlated with learning rate.

Fig. 7

Individual Differences in Learning Rate Correlated with State Flexibility

State flexibility difference refers to the difference in state flexibility between consecutive trial sets, ordered by the number of trials practiced). State flexibility differences for all regions were computed and were found to be significantly positively correlated with individual differences in learning rate (ρ = 0.24, p = 0.041), suggesting that the observed increase in flexibility is associated with the learning rate of subjects.

Comparison to traditional and other network analyses

Next, we were curious how our hierarchical clustering based approach compared to more traditional GLM univariate analysis. To answer this question, we calculate the mean BOLD magnitude in each of the primary states and compute the Pearson correlation with GLM β weights. We find that both the motor (ρ = −0.75, p = 7.83 × 10−22, see Fig. 8a) and frontal (ρ = 0.67, p = 6.29 × 10−16, see Fig. 8a) states are strongly correlated with the GLM β weights, but interestingly in the opposite direction. Motor regions that decrease in their task-related activity most had higher average BOLD signal, while frontal regions that decrease in their task-related activity most had lower average BOLD signal. These results demonstrate that the GLM β-weight is not directly related to the BOLD signal magnitudes driving the time-by-time network structure, as the relationship between these two variables can change drastically in magnitude and sign in different brain regions.

Fig. 8

Time-by-Time Network Analysis Compared to GLM

(a) GLM weights are strongly correlated with both the motor (ρ = −0.75, p = 7.83 × 10−22) and frontal (ρ = 0.67, p = 6.29 × 10−16) states as defiend through hierarchical clustering. Each data point represents a single brain region, with values of state-specific BOLD expression (x-axis) or GLM beta-weight (y-axis) averaged over subjects. (b) GLM weights are strongly correlated (ρ = 0.32, p = 4.86 × 10−4) with contributions to state flexibility as defined through state flexibility ablation analysis.

We were then curious how our ablation analysis, which identified individual regions’ contribution to state flexibility, compared to the same GLM beta weights. We again find a significant correlation between GLM beta weights and contribution to state flexibility (ρ = 0.32, p = 4.86 × 10−4, see Fig. 8b), suggesting that regions whose activity decreased most with training also showed greatest contributions to brain state flexibility.

Finally, we sought to compare this new analysis with previous-defined network flexibility measures derived from multilayer community detection, which have been estimated on this dataset in prior studies (Bassett et al., 2013a,b; 2015). We find a significant correlation (ρ = 0.27, p = 2.69 × 10−19, see Fig. 9) between our defined state flexibility and this previously-defined network flexibility measure. This suggests, intuitively, that the global flexibility is related to (but not exactly the same as) the average regional flexibility. Our state flexibility measure explains 7.28% of the variance in previously-defined network flexibility, suggesting that state flexibility also captures some unique characteristic of network dynamics. In the context of this study, these results provide insight into changes in patterns of activity at the TR level (brain state flexibility) that can accompany changes in patterns of functional connectivity at the level of 2-min windows (network flexibility). It could be interesting in future to use the higher temporal specificity of the brain state flexibility measure to understand fine-scale temporal variations in other tasks in which network flexibility has been implicated in behavior, including attention (Shine et al., 2016), linguistic processing (Chai et al., 2016), and working memory (Braun et al., 2015b).

Fig. 9

Time-by-Time Network Analysis Compared to Network Flexibility From Multilayer Community Detection

State flexibility, as defined through our time-by-time network analysis, is shown to be correlated (ρ = 0.24, p = 2.69 × 10−19) with measures of network flexibility, accounting for 7.28% of its variance.

Discussion

In this work, we studied task-based fMRI data collected at 4 time points separated by about 2 weeks during which healthy adult participants learned a set of six 10-note finger sequences. During learning, we hypothesized that the brain would show a change in the manner in which it traversed brain states. We defined a state as a pattern of BOLD magnitude across 112 anatomically-defined brain regions. We identified two canonical states characteristic of the entire period of task performance, which showed high activation of motor cortex and frontal cortex, respectively. Interestingly, we observed that the flexibility with which participants switched among these canonical states and other less common states was lowest early in training and highest late in training, indicating the emergence of state flexibility. We find that the positive contributors to state flexibility are associated with integrative processing while the negative contributors are associated with motor and visual function. Finally, we observe that changes in state flexibility were correlated with learning rate: increasing state flexibility was correlated with higher learning rates.

Extensions of graph theoretical tools to the temporal domain

Over the past decade, tools from graph theory have offered important insights into the structure and function of the human and animal brain, both at rest and during cognitively demanding tasks (Sporns, 2010). In these applications, the nodes of the graph are traditionally thought of as neurons or brain areas, and the edges of the graph are defined by either anatomical tracts (Hagmann et al., 2008; Bassett et al., 2011a) or by functional connections (Achard et al., 2006; Bassett et al., 2006). Yet, the tools of graph theory are in fact much more general than these initial applications (Bollobas, 1979, 1985). Indeed, recent extensions have brought these tools to other domains – from genetics (Fulcher and Fornito, 2016; Conaco et al., 2012; Arcila et al., 2014) to orthopedics (Murphy et al., 2018) – by carefully defining alternative graph representations of relational data. As a concrete example, a graph can be used to encode the relationships between movements or behaviors, by treating a movement as a node, and by linking subsequent movements (or actions) by the inter-movement interval (Wymbs et al., 2012). Similarly, a graph can be used to code the temporal dependencies between stimuli, by treating a stimulus as a node, and by linking pairs of stimuli by their temporal transition probabilities (Schapiro et al., 2013; Karuza et al., 2016).

While these applications may initially seem vastly different, they in fact all share a common property: that entities are related to one another by some facet of time. Here, by contrast, we construct the edge-vertex dual of this more common form. We ask: How are times related to one another by some other entity? Specifically, we study how the brain state in one time point is related to the brain state in another time point, and we define a brain state as the vector of activation magnitudes across all regions of interest (Medaglia et al., 2018; Chen et al., 2016). The notion that a pattern of activation reflects a brain state is certainly not a new one (Gu et al., 2015a). In the context of fMRI data, a common approach is to study the multi-voxel pattern of activation in a region of interest to better understand the representation of a stimulus (Kubilius et al., 2015; Chadwick et al., 2012). In the context of EEG and MEG data, the pattern of power or amplitude in a set of sensors or a set of reconstructed sources is frequently referred to as a microstate (von Wegner et al., 2016). The composition and dynamics of these microstates have shown interesting cognitive and clinical utility, predicting working memory (Muthukrishnan et al., 2016) and disease (Geschwind et al., 2016). Yet, while patterns of activation are acknowledged as an important representation of a brain or cognitive state, little is known about how these states evolve into one another. Recent advances have made this possible by coding the relationships between brain states in a graph (Medaglia et al., 2018). Here we capitalize on these advances to extract the community structure in such a graph, to identify canonical states, and to quantify the transitions between them during task performance. It will be interesting in future to broaden the analytical framework applied here to study other properties of the graph – including local clustering and global efficiency – to better understand how the brain traverses states over time.

Brain states characteristic of discrete sequence production

Using this unusual graph theory approach in which network nodes represent time points and network edges represent similarities in brain states across two time points, we were able to identify two canonical brain states that characterized the task-evoked activity dynamics across the entire experiment, extending across 6 weeks of intensive training. The most common state, perhaps unsurprisingly, was characterized by high BOLD magnitudes in regions of the extended motor cortex, including the bilateral precentral gyrus, left postcentral gyrus, bilateral superior parietal lobule, bilateral supramarginal gyrus, bilateral supplementary motor area, bilateral parietal operculum cortex, and bilateral Heschl's gyrus (Dayan and Cohen, 2011). This map is consistent with the fact that this is an intensive motor-learning paradigm (Bassett et al., 2015; Wymbs and Grafton, 2015) in which participants acquire the skill necessary to perform a sequence of 10 finger movements over a short period of time. The second most common state was composed of a frontal-temporal-subcortical system, containing the anterior middle temporal gyrus, medial frontal cortex, parahippocampal gyrus, caudate, nucleus accumbens, planum temporale, and hippocamus. These areas are thought to play critical roles in sequence learning (Exner et al., 2002; Vakil et al., 2000; Bischoff-Grethe et al., 2000) facilitated by higher-order cognitive processes including reward learning (Hikosaka et al., 2014; Kim and Hikosaka, 2015), cognitive control and executive function (Stuss, 2011; Alvarez and Emory, 2006), predicting nature and timing of action outcomes (Brown, 2011; Grinband et al., 2011; Rushworth et al., 2005), and subcortical storage of motor sequence information (Lehericy et al., 2005). This system is particularly interesting because it displayed a competitive relationship with the motor state, with a strongly anti-correlated activation profile, suggesting that frontal-subcortical circuitry affects control by transient, desynchronized interactions.

State flexibility, task practice, and learning rate

Beyond the anatomy of the states that characterize extended training on a discrete sequence production task, it is also useful to study the degree to which those states are expressed, and the manner in which one state moves into another state. The two primary states that we observed characterized 95.67% of all time points, indicating their canonical nature. Temporally, the brain frequently switched back and forth between these two states, with less frequent traversal of other non-primary states. We quantified this switching using a brain state flexibility measure (Medaglia et al., 2018), and observed that flexibility increased significantly over the course of the 6 weeks of training. Moreover, brain state flexibility was negatively correlated with learning rate, being lowest early in training when behavioral adaptivity was greatest. These results suggest that consistent activation patterns characterize early training, when participants must learn the mapping of visual cues to motor responses, the use of the button box, and the patterns of finger movements. Later in learning, when the skill has become relatively automatic, participants display more varied progressions of activation patterns (higher brain state flexibility), potentially mirroring the greater freedom of their cognitive resources for other processes (Shamloo and Helie, 2016).

Importantly, these results offer a complement to prior efforts to quantify network flexibility based on estimates of functional connectivity (Bassett et al., 2011b) where the nodes are brain regions and the edges are temporally defined correlations between those regions. Network flexibility appears to peak early in finger sequence training (Bassett et al., 2011b), followed by a growing autonomy of motor and visual systems (Bassett et al., 2015). In combination with our results, these prior data suggest that there may be distinct time scales associated with brain variability at the level of activity (where variability may peak late) in comparison to the level of connectivity (where variability may peak early). Such a hypothesis could be directly validated in additional studies that reproduce the results we present here. The apparent separation in time scales of these processes over learning also supports the growing notion that the information housed in patterns of activity can be quite independent from information housed in patterns of connectivity (Bassett et al., 2015; Cao et al., 2016; Siebenhuhner et al., 2013). For example, earlier studies have demonstrated that patterns of beta weights from a GLM do not necessarily map onto patterns of strong or weak functional connectivity (Bassett et al., 2015), the temporal dynamics of an activity time trace do not necessarily map onto patterns of functional connectivity (Siebenhuhner et al., 2013), and phenotypes indicative of psychiatric disease can be identified in functional connectivity while being invisible to methods focused on activity (Esslinger et al., 2009). Together, these studies indicate that activity and connectivity can provide distinct information regarding the neurophysiological processes relevant for cognition and disease. They also in principle support the possibility of differential time scales of flexibility in activity and connectivity as a function of learning.

Comparison to traditional analyses

Our unique time-based clustering approach to identify brain states also offers additional information beyond traditional analysis methods while remaining consistent with previous results on the same dataset. We found that across-subject mean BOLD signal during the motor and frontal states was strongly correlated with the across-subject mean beta-weights from a more traditional GLM univariate analysis. That is, regions with high across-subject mean BOLD signal in the frontal state tended to show high across-subject mean beta-weights, while the opposite was true for the motor state. These observations suggest a strong relationship between the two measures of brain function, where the directionality of that relationship is significantly dependent on which state (frontal or motor) the brain is in. Due to the state dependence of this relationship, we cannot claim that beta weights drive state flexibility. The differences between the two measures become intuitive when one considers the fact that the two measures capture information on different temporal and spatial scales of the experiment. Specifically, regions switch relatively rapidly between the motor and frontal states, which remain fairly conserved throughout learning, whereas the GLM beta weights measure decreases in task-related activity over the entire time course. Additionally, we find that our state flexibility results are consistent with previous network flexibility results from more traditional multilayer community detection, but that there still exists a large amount of variance that remains unexplained by either statistic, an observation that is consistent with the notion that there exists some unique information in both. State flexibility and network flexibility intuitively operationalize quite different concepts; while network flexibility measures regional allegiance to functional modules (measuring functional stability), state flexibility quantifies brain state allegiance to canonical brain states (measuring temporal stability). Leveraging both measures simultaneously in the future may allow for a more holistic understanding of brain dynamics. Such an approach would draw information from regional activity and interregional connectivity; both well established metrics of brain dynamics that are rarely used in concert.

Methodological considerations

There are several important methodological and conceptual considerations pertinent to this work. The first consideration we would like to discuss is a relatively philosophical one. It pertains to our use of the term "brain state". It is important to disambiguate the use of brain state as a quantifiable and quantified object, defined as the pattern of activation magnitudes over all brain areas (strung out in a vector (Gu et al., 2015a)), and other more conceptual notions of mental state or cognitive state. These latter notions can be difficult to quantify directly from imaging data, even if they may have relatively specific definitions from both psychological and clinical perspectives (Richman and Unoka, 2015; Martin and Santos, 2016; Eddy et al., 2016). It will be important in future uses of our brain state detection and characterization technique to maintain clarity in the use of these terms.

The second important consideration relevant to this work is that the data that we study here was collected with a traditional 2 s TR. It would be very interesting to test for similar phenomena in the high-resolution BOLD imaging techniques available now, for example using multiband acquisitions. Such higher sampling could offer heightened sensitivity to changes in brain state flexibility related to individual differences in learning. Moreover, they could provide enhanced sensitivity to variations in brain state flexibility across different frequency bands, particularly higher frequency bands that have been shown to be sensitive to shared genetic variance (Fornito et al., 2011).

Finally, on a computational note, it is important to emphasize that the results described here are obtained via the application of a clustering technique (Porter et al., 2009; Fortunato, 2010) to identify brain states from the temporal graph. Importantly, the technique that we use – based on modularity maximization (Newman, 2006) – is a hard partitioning algorithm that seeks to solve an NP-hard problem using a clever heuristic (Blondel et al., 2008). Although modularity maximization can accurately recover planted network modules in synthetic tests (Lancichinetti and Fortunato, 2009; Bassett et al., 2013a), it does have important limitations (Lancichinetti and Fortunato, 2011; Good et al., 2010). Therefore it would be interesting in future to examine the sensitivity of results to other clustering techniques available in the literature.

Conclusion

In summary, in this study we seek to better understand the changes in brain state that accompany the acquisition of a new motor skill over the course of extended practice. We treat the brain as a dynamical system whose states are characterized by a recognizable pattern of activation across anatomicaly defined cortical and subcortical regions. We apply tools from graph theory to study the temporal transitions (network edges) between brain states (network nodes). Our data suggest that the emergence of automaticity is accompanied by an increase in brain state flexibility, or the frequency with which the brain switches between activity states. Broadly, our work offers a unique perspective on brain variability, noise, and dynamics (Deco et al., 2009; Breakspear and McIntosh, 2011; Garrett et al., 2013, 2014), and its role in human learning.