Multiscale and multimodal network dynamics underpinning working memory

TL;DR

This study uncovers how multiscale neural dynamics involving frontoparietal and default mode systems, their subnetworks, and structural linkages influence working memory performance, integrating neuroimaging and computational modeling.

Contribution

It reveals the distinct roles of frontoparietal subnetworks and their gene expression, structural connectivity, and dynamic interactions in working memory, advancing mechanistic understanding.

Findings

Inverse correlation between FPS-DMS interaction strength and WM performance

Identification of two FPS subnetworks with different gene expression patterns

Computational model shows subnetwork activity modulates functional interactions

Abstract

Working memory (WM) allows information to be stored and manipulated over short time scales. Performance on WM tasks is thought to be supported by the frontoparietal system (FPS), the default mode system (DMS), and interactions between them. Yet little is known about how these systems and their interactions relate to individual differences in WM performance. We address this gap in knowledge using functional MRI data acquired during the performance of a 2-back WM task, as well as diffusion tensor imaging data collected in the same individuals. We show that the strength of functional interactions between the FPS and DMS during task engagement is inversely correlated with WM performance, and that this strength is modulated by the activation of FPS regions but not DMS regions. Next, we use a clustering algorithm to identify two distinct subnetworks of the FPS, and find that these subnetworks display distinguishable patterns of gene expression. Activity in one subnetwork is positively associated with the strength of FPS-DMS functional interactions, while activity in the second subnetwork is negatively associated. Further, the pattern of structural linkages of these subnetworks explains their differential capacity to influence the strength of FPS-DMS functional interactions. To determine whether these observations could provide a mechanistic account of large-scale neural underpinnings of WM, we build a computational model of the system composed of coupled oscillators. Modulating the amplitude of the subnetworks in the model causes the expected change in the strength of FPS-DMS functional interactions, thereby offering support for a mechanism in which subnetwork activity tunes functional interactions. Broadly, our study presents a holistic account of how regional activity, functional interactions, and structural linkages together support individual differences in WM in humans.

Figures

Equations

$$P\left(C_f|z,\mu ,{\sigma }^2\right)=\prod _{i=1}^N\prod _{j=1}^N\exp \left(C_{f,ij}\cdot \frac{{\mu }_{z_i{z}_j}}{{\sigma }_{z_i{z}_j}^2}-\frac{C_{f,ij}^2}{2{\sigma }_{z_i{z}_j}^2}-\frac{{\mu }_{z_i{z}_j}^2}{{\sigma }_{z_i{z}_j}^2}\right).$$

$$BC(i)=\{\begin{array}{cc}1-{\left(\frac{k_i(j)}{k_i}\right)}^2-{\left(\frac{k_i(l)}{k_i}\right)}^2& \text{for}{k}_i(j)+k_i(l)=k_i\\ & \\ {\left(\frac{k_i(j)}{k_i}\right)}^2+{\left(\frac{k_i(l)}{k_i}\right)}^2& \text{for}{k}_i(j)+k_i(l)<k_i\end{array}$$

$$\frac{d{u}_j}{dt}=u_j[a_j+i{\omega }_j-|u_j{|}^2]+\xi {\eta }_j(t),$$

$$\frac{d{x}_j}{dt}=\left[a_j-x_j^2-y_j^2\right]x_j-{\omega }_j{y}_j+G\sum _{i=1}^{n_o}{D}_{ij}\left(x_i-x_j\right)+\xi {\eta }_j(t),$$

$$\frac{d{y}_j}{dt}=\left[a_j-x_j^2-y_j^2\right]y_j+{\omega }_j{x}_j+G\sum _{i=1}^{n_o}{D}_{ij}\left(y_i-y_j\right)+\xi {\eta }_j(t).$$

$$R(t)=\frac{1}{n_o}\left|\sum _{j=1}^{n_o}{\text{e}}^{i{\varphi }_j(t)}\right|,$$

$${\mathrm{OutcomeVariable}}_{it}={\mathcal{B}}_{0i}+{\mathcal{B}}_{1i}{\mathrm{SubNetwork}}_{it}+e_{it},$$

$${\mathcal{B}}_{0i}={\gamma }_{00}+u_{0i},$$

$${\mathcal{B}}_{1i}={\gamma }_{10},$$

$${Measure1}_{it}={\overline{Measure2}}_i+{\mathrm{Subject}}_i+c({Measure2}_i)+e_{it},$$

Full text

**Multiscale and multimodal network dynamics underpinning working memory

Andrew C. Murphy1,2, Maxwell A. Bertolero1, Lia Papadopoulos3, David M. Lydon-Staley1, and Danielle S. Bassett1,3-6,*

1Department of Bioengineering, School of Engineering & Applied Sciences, University of Pennsylvania, Philadelphia, PA 19104, USA

2Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA

3Department of Physics & Astronomy, School of Arts & Sciences, University of Pennsylvania, Philadelphia, PA 19104, USA

4Department of Neurology, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA

5Department of Psychiatry, Perelman School of Medicine, University of Pennsylvania, Philadelphia, PA 19104, USA

6Department of Electrical & Systems Engineering, School of Engineering & Applied Sciences, University of Pennsylvania, Philadelphia, PA 19104, USA

To whom correspondence should be addressed: dsb@seas.upenn.edu*

Abstract

As a key component of executive function, working memory allows information to be stored and manipulated over short time scales. Performance on laboratory and ecological working memory tasks is thought to be supported by the frontoparietal system, the default mode system, and interactions between them. Yet, little is known about how these systems and their interactions might relate to, or perhaps even explain, individual differences in working memory performance. Here, we address this gap in knowledge using functional magnetic resonance imaging data acquired during the performance of a 2-back working memory task, as well as diffusion tensor imaging data collected in the same individuals as part of the Human Connectome Project. We show that the strength of functional interactions between the frontoparietal and default mode systems during task engagement is inversely correlated with working memory performance. Further, we demonstrate that the strength of functional interactions between these two systems is modulated by the activation of frontoparietal regions but not default mode regions. To gain a deeper understanding of these processes, we used an unsupervised clustering algorithm informed by an explicit model of network architecture to identify two distinct subnetworks of the frontoparietal system, and we subsequently demonstrate that these subnetworks also display distinguishable patterns of gene expression. We found that activity in one subnetwork was positively associated with the strength of the functional interaction between the frontoparietal and default mode systems, while activity in the second subnetwork was negatively associated. Using the diffusion imaging data, we then demonstrate that the pattern of structural linkages between these subnetworks explains their differential capacity to influence the strength of the functional interaction between the frontoparietal and default mode systems. To determine whether these correlative observations could provide a mechanistic account of large-scale neural underpinnings of working memory, we built a computational model of the system composed of simplified coupled oscillators. We demonstrate that modulating the relative amplitude of the subnetworks in the model causes the expected change in the strength of the functional interaction between the frontoparietal and default mode systems, thereby offering support for a candidate mechanism in which subnetwork activity tunes functional connectivity. Broadly, our study presents a holistic account of how regional activity, functional connections, and structural linkages together support individual differences in working memory in humans.

1 Introduction

Working memory supports the short-term storage of information, thereby facilitating its further manipulation and processing [7]. Individual differences in working memory performance have been associated with differences in the recruitment of distinct cognitive systems, and in the functional interactions between such systems [48, 52, 11]. Heavily implicated in working memory function is a group of regions in frontoparietal cortex involved in cognitive control, among other functions [62, 60], as well as a group of regions comprising the so-called default mode system, which is most notable for its strong activation in the human resting state. Interestingly, these two systems are thought to have opposite effects on working memory: frontoparietal activity is vital for directing attention to external stimuli [66], while default mode activity is important for internally-directed cognition [2]. Furthermore, these two systems tend to be in functional competition during tasks with high working memory load; in such tasks, the activity of the two systems is anti-correlated [21].

Despite extensive study, the functional role of competition between the frontoparietal and default mode systems remains unclear. Over the past decade, several hypotheses have been put forward. For example, one notable hypothesis suggests that competition between the frontoparietal and default mode systems during certain tasks enables maximally disjunctive levels of activity [21]. Put more simply, this inter-system competition might allow for a pattern of whole brain dynamics characterized by heightened activity in the frontoparietal system coupled with decreased activity in the default mode system. Explaining such a pattern of dynamics is particularly important in light of evidence that it favors improved working memory performance [6]. A second notable hypothesis suggests that competition may enable a toggling between (i) brain states eliciting default mode activity for introspective processes and (ii) brain states eliciting frontoparietal activity for externally-directed attentional processes [34]. A third hypothesis suggests that inter-system competition could prove useful for continuously modulating the brain’s response to task complexity; while activity in the frontoparietal system scales with the complexity of the task, activity in the default mode system decreases with the complexity of the task [12]. It is worth noting that – while of course conceptually interesting – none of these hypotheses offer particularly explicit mechanisms for competition.

Its functional role aside, the anatomical substrates and extent of the competition between the frontoparietal and default mode systems remain unclear. Structurally, it is well-known that the two systems are quite distinct in terms of their topological role within the connectome: the default mode system is part of the so-called rich club of the structural connectome, which is a set of densely interconnected high degree nodes, while the frontoparietal system is part of the so-called diverse club, which is a set of densely interconnected nodes with diverse connectivity across all putative cognitive systems [13]. It is intuitively plausible that such distinct placements within the larger whole-brain network could constrain or define the roles that each system can play in inducing certain types of dynamics [42, 25]. In addition to its anatomical substrates, the anatomical extent of the competition remains unclear. Is the entirety of the default mode system competitive with the entirety of the frontoparietal system? Are these systems in fact homogeneous or could they contain subnetworks with distinct dynamics? Initial evidence suggests that the frontoparietal system can be separated into two subnetworks: one more connected to the default mode system and one more connected to the dorsal attention system [29]. Similar evidence suggests that the default mode system can be separated into several subnetworks [3] supporting distinct processes in the internal generation of thought and the local processing of information [20]. A careful investigation into the nature of these system-specific subnetworks may inform more mechanistic models of competition and its relevance for working memory performance.

Here, we seek to perform a deeper study of individual differences in working memory function and its relation to multimodal neural phenotypes including regional activity, inter-regional connectivity, and structural linkages. We also seek to complement data-driven analysis of empirical measurements with biologically-motivated computational modeling to probe the validity of our explanations and posited mechanisms. Specifically, in this study we use functional magnetic resonance imaging data collected from 644 healthy adult human participants in the Human Connectome Project during the performance of a 2-back working memory task. To address potential structural drivers of our findings, we also use diffusion tensor imaging data acquired in the same participants. We begin with a broad investigation into the relations between behavior, activity within the frontoparietal and default mode systems, and competition between the frontoparietal and default mode systems. We then use a model-based machine learning algorithm to uncover functional groups or subnetworks within the entire frontoparietal system, which we subsequently find to display distinguishable patterns of gene coexpression. We show that subnetwork activity tracks the strength of the functional interaction between the frontoparietal and default mode systems, in a manner that is consistent with the underlying pattern of structural linkages between them. Lastly, – drawing on methods from complex systems physics and non-linear dynamical systems theory – we build a computational model of the system and its dynamics as a collection of coupled oscillators, with parameters and coupling architecture informed by the biological evidence we uncover. We use the model to probe which features of the biology are able to produce the observed phenomena, and to test our hypotheses about the relationships between activity, anatomical connectivity, and functional interactions.

2 Materials and Methods

2.1 Imaging data acquisition and preprocessing

2.1.1 Data acquisition

For each subject in the Human Connectome Project (HCP) S900 release [81], we extracted the task-based functional magnetic resonance imaging data acquired during the performance of the n-back working memory task, a resting state functional magnetic resonance imaging scan, a high-resolution anatomical scan, and a diffusion tensor imaging scan. In this release, 644 subjects contained all 4 data types, all four resting-state scans, and BedpostX diffusion data. Participants were 346 females, mean age (std) = 28.6 (3.68) years. No additional exclusion criteria were applied. All analyses were performed in accordance with the relevant ethical regulations of the WU-Minn HCP Consortium Open Access Data Use Terms. Informed consent was obtained in writing from all participants.

The acquisition parameters for each data type are as follows. The parameters for the acquisition of high-resolution structural scan were: TR = 2400 ms, TE = 2.14 ms, TI = 1000 ms, flip angle = 8 deg, FOV = 224 $\times$ 224 mm, voxel size = 0.7 mm isotropic, BW = 210 Hz/Px, acquisition time = 7:40 minutes. Functional magnetic resonance images were collected during both rest and task with the following parameters: TR = 720 ms, TE = 33.1 ms, flip angle = 52 deg, FOV = 208$\times$180 mm, matrix = 104$\times$90, slice thickness = 2.0 mm, number of slices = 72 (2.0 mm isotropic), multi factor band = 8, echo spacing = 0.58 ms. Diffusion tensor images were collected with the following parameters: TR = 5520 ms, TE = 89.5 ms, flip angle = 78 deg, refocusing flip angle = 160 deg, FOV = 210$\times$180, matrix = 168$\times$144, slice thickness = 1.25 mm, number of slices = 111 (1.25 mm isotropic), multiband factor = 3, echo spacing = 0.78 ms, $b$-vaues = 1000, 2000, and 3000 s/mm2.

2.1.2 Working memory task and associated measure of behavior

We focused our analyses on data acquired during the n-back task, due to its reliable recruitment of the executive system [69]. The working memory task was presented at two different levels of difficulty: 0-back and 2-back. For both levels, subjects were presented with a stream of images taken from the following four categories: faces, places, tools, or body parts. The latter images presented body parts that were whole (non-mutilated); no images contained nudity. In the 0-back condition, subjects were meant to respond positively during every image presentation. In the 2-back condition, subjects were meant to respond positively if the present image was identical to the image presented two previously. The task was divided into two runs, each run composed of 8 task blocks and 2 fixation blocks [9]. Fixation blocks lasted for 15 seconds each. Each task block consisted of 10 trials, where a stimulus was presented for 2 seconds, followed by a 500 ms ITI (2.5 seconds total per trial) [9]. Each block begins with a 2.5 second cue indicating the 0-back or 2-back condition. During the 0-back condition, working memory loads are minimal. Within each run, half of the task blocks used the 2-back paradigm, and half used the 0-back paradigm. Additionally, within a run, each stimulus type (images from a single category) was presented in a different block. To estimate behavioral performance, we calculated the accuracy of responses across all image categories separately for 0-back and 2-back conditions. We chose to focus on accuracy due to its interpretability [11]. However, we also considered d-prime [69, 47] and demonstrate that our main results hold when using this metric in place of accuracy (see Supplementary Materials).

2.1.3 Processing of functional magnetic resonance imaging data

For both resting-state and task functional connectivity, CompCor, with five principal components from the ventricles and white matter masks, was used to regress out nuisance signals from the time series. In addition, the 12 detrended motion estimates provided by the Human Connectome Project were regressed out from the time series. The mean global signal was removed and then time series were band-pass filtered from 0.009 to 0.08 Hz. Finally, frames with greater than 0.2 mm frame-wise displacement or a derivative root mean square (DVARS) above 75 were removed as outliers. Sessions composed of greater than 50 percent outlier frames were not further analyzed.

We chose to regress out the global signal from the time series because it has been shown to remove motion signal and global scanner noise. We also note that the mathematics of global signal regression does not necessitate a specific spatial distribution of negative correlations [58], and our claims regard the relative strengths of connectivity between networks rather than their sign. Moreover, the processing pipeline used here has been suggested to be ideal for removing false relations between connectivity and behavior [72]. Finally, we note that our main results hold when using wavelet coherence as a measure of functional connectivity (see Supplementary Materials); this measure is bounded between 0 and 1.

2.1.4 Preprocessing of diffusion tensor imaging data

For the diffusion imaging, the Human Connectome Project applied intensity normalization across runs, the TOPUP algorithm for EPI distortion correction, the EDDY algorithm for eddy current and motion correction, gradient nonlinearity correction, calculation of gradient $b$-value/$b$-vector deviation, and registration of mean $b$0 to native volume T1w with FLIRT. BBR+bbregister and transformation of diffusion data, gradient deviation, and gradient directions to 1.25 mm structural space were also applied. The brain mask is based on the FreeSurfer segmentation. The BedpostX (Bayesian Estimation of Diffusion Parameters Obtained using Sampling Techniques) output was then calculated, where the ‘X’ stands for modeling crossing fibers. Markov Chain Monte Carlo sampling was used to build probability distributions on diffusion parameters at each voxel. The process creates all of the files necessary for running probabilistic tractography. Using the Freesurfer recon-all data computed by the Human Connectome Project, the fsaverage5 space cortical parcellation was registered to the subject’s native cortical white matter surface and then transformed to the subject’s native diffusion volume space. From these data, we derived seeds and targets for probabilistic tractography, which we ran with the FSL probtrackx2 algorithm using 1000 streams initiated from each voxel in a given parcel.

2.1.5 Whole brain functional parcellation

We parcellated the brain into 400 discrete and non-overlapping regions of interest using the Schaefer atlas (fslr32k surface) [70]. Notably, the Schaefer atlas was originally developed in the same HCP data that we study here, and it yields a functional demarcation of both the default mode and the frontoparietal systems. Of course other functionally defined atlases exist, but they are less ideal for our purposes for several reasons; the Power atlas [64] does not provide full cortical coverage, and the Gordon [41] and Brainnetome [32] atlases are lower spatial resolution including 333 and 246 regions, respectively. The Schaefer atlas provides an assignment of each region to one of 17 putative cognitive systems: two visual, two somatomotor, two dorsal attention, two salience/ventral attention, one limbic, three frontoparietal, three default mode, and one temporo-parietal system. To ensure that the granularity of the data was consistent with the granularity of our hypotheses, we collapsed these 17 systems into 8 systems by combining individual systems that belonged to the same cognitive system; that is, we combined the two visual systems into a single system, the two somatomotor systems into a single system, the two dorsal attention systems into a single system, the two salience systems into a single system, the three frontoparietal systems into a single frontoparietal system, and the three default mode systems into a single system.

2.2 Analysis of functional magnetic resonance imaging data

2.2.1 Estimation of functional activation

To measure functional activation during task performance, we employed a general linear model (GLM). We collated time series from each brain region during either the 0-back condition or the 2-back condition. For each of the two conditions, we formulated a regressor that indicated the volumes in the time series during which the subject was engaged in the task (not in a fixation block). Treating the 0-back and 2-back time series separately, we convolved this regressor with the canonical hemodynamic response function using the SPM (Wellcome Centre for Human Neuroimaging, London, UK) function spm_hrf, to create a matrix $M_{conv}\in\mathcal{R}^{N\times T}$, where N is the number of regions and T is the number of time points. We collated observed regional time series to create a matrix $M_{TS}\in\mathcal{R}^{N\times T}$. We then used the MATLAB function mldivide to solve the element-wise multiplication equation $M_{TS}=\gamma\times M_{conv}$ for the maximum likelihood estimate (MLE) of the $\gamma\in\mathcal{R}^{1\times N}$ coefficients, where we took the $\gamma$ estimates to be regional activity. We solved this MLE separately for each of the two conditions, providing distinct estimates of regional functional activation for 0-back and 2-back.

2.2.2 Estimation of functional connectivity matrices

We used the preprocessed data to construct functional connectivity matrices reflecting functional interactions between regions and systems. Specifically, we extracted processed time series from each of the 400 regions in the Schaefer atlas (Fig. 1 A). Next, we calculated the Pearson correlation coefficient between each pair of regional time series (Fig. 1 B). We chose to use the Pearson correlation to represent functional connectivity due to its widespread use in the neuroimaging literature, as well as its ease of interpretability [83], but we also demonstrate robustness of our results to other measures of functional connectivity in the Supplementary Materials. We collated all inter-regional estimates of functional connectivity into a single $400\times 400$ connectivity matrix, $C_{f}$ (Fig. 1 C), which we then treated as the formal encoding of a network model of brain function [10]. To be explicit, in this network model regions are represented by network nodes, and functional connections are represented by weighted edges, where the weight of the edge between node $i$ and node $j$ is given by the Pearson correlation coefficient between the time series of region $i$ and the time series of region $j$. Finally, we averaged the estimates of functional connectivity within systems, and between pairs of systems, to construct a system-by-system connectivity matrix (Fig. 1 D).

2.2.3 Identification of subnetworks using a weighted stochastic block model

To partition the frontoparietal system into two functionally disjoint groups, we employed the weighted stochastic block model, which is a powerful community detection method. This method is complementary to the more widely used modularity maximization methods [59], but is noted to have increased flexibility and sensitivity to a more diverse set of network architectures [33, 16, 57, 14, 16]. Briefly, the weighted stochastic block model is a generative model that places each of the $N$ nodes of network $C_{f}$ into one of K communities. This placement is accomplished by finding a network partition $z\in\mathcal{Z}^{N\times 1}$ where $z_{i}\in\{1,2,\dots,K\}$ and $z_{i}$ denotes the membership of node i. Assuming that network edge weights are normally distributed, following [1] and [16], the generative model takes the following form:

[TABLE]

Here we have introduced model parameters $\mu\in\mathcal{R}^{K\times K}$ and $\sigma^{2}\in\mathcal{R}^{K\times K}$, where $\mu_{z_{i}z_{j}}$ and $\sigma^{2}_{z_{i}z_{j}}$ parameterize the weights of normally distributed connections between community $z_{i}$ and community $z_{j}$, and $C_{f,ij}$ denotes the $ij$th element of the network $C_{f}$. Furthermore, $P(C_{f}|z,\mu,\sigma^{2})$ is the probability of generating the observed network $C_{f}$ given the parameters. This model is fit to $C_{f}$ in order to estimate the parameters $z$, $\mu$, and $\sigma^{2}$. We fit the model using MATLAB code provided in [1] and freely available at http://tuvalu.santafe.edu/~aaronc/wsbm/.

For each subject, we fit the weighted stochastic block model to the $m\times m$ subgraph of the adjacency matrix representing functional connections between the $m=61$ regions of the frontoparietal system. We selected $K=2$ a priori due to prior evidence that the frontoparietal system can be separated into two distinct components [29]. The implementation generated a single maximum likelihood partition of regions into functional communities for each subject. Then, we pooled partitions across subjects and used a consensus similarity method [30] to identify a single partition that is most similar to all others, where similarity is quantified by the $z$-score of the Rand coefficient [80]. To assess the statistical significance of this partition, we performed a non-parametric permutation test. Specifically, for each subject, we calculated the log-likelihood of the weighted stochastic block model fitting the final consensus partition to their individual functional network. As a null, we calculated the log-likelihood of the weighted stochastic block model fitting a random permutation of the final consensus partition to their individual network. We assessed the difference between the true and null data using a multilevel model (see Section 3.2).

2.3 Analysis of diffusion tensor imaging data

2.3.1 Estimation of structural connectivity matrices

After performing probabilistic tractography, we applied the same 400-region Schaefer atlas. Next, we calculated the proportion of streams seeded in a voxel in one region that reached another region. We chose to use the proportion of streamlines to represent structural connectivity due to the inhomogeneity of the region sizes. We collated all inter-regional estimates of structural connectivity into a single $400\times 400$ connectivity matrix, $C_{s}$, which we then treated as the formal encoding of a network model of brain structure [10]. Similar to the model of brain function, in this structural network model, regions are represented by network nodes, and structural connections are represented by weighted edges, where the weight of the edge between node $i$ and node $j$ is given by the proportion of streams seeded at region $i$ that reach region $j$. Finally, we averaged the estimates of structural connectivity within systems, and between pairs of systems, to construct a system-by-system connectivity matrix, akin to the one that we constructed from the functional data.

2.3.2 Calculation of boundary controllability

We posited that structural connections between systems would play an important role in the functional coupling between the frontoparietal and default mode systems. Specifically, we hypothesized that the formal nature of that role was one of boundary controllability, more commonly studied in the field of control and dynamical systems theory [61]. Boundary control is a quantifiable metric describing the notion that the topological location of a region within a structural network partially governs that region’s influence on the function of modules or communities in the network [15, 42]. Intuitively, if region i has strong structural connections to regions j and l, the activity of region i influences the functional connection between regions j and l. Boundary control assesses the positioning of a region between two other regions, and can be calculated for region i with respect to its control over regions j and l:

[TABLE]

Here, $k_{i}$, the degree of region i, is the sum of all region i’s structural connections. The variables $k_{i}(j)$ and $k_{i}(l)$ are the strength of region i’s structural connections to regions j and l, respectively. A region with high boundary control is predicted to more effectively modulate the functional connection between region j and region l, although boundary control does not assess whether the modulation will increase or decrease the connection strength.

3 Gene coexpression analysis

To determine whether the two subnetworks that we identified in the frontoparietal system displayed distinct patterns of gene expression, we used gene expression data from six post-mortem brains available from the Allen Brain Institute. We focused our analyses on 16699 genes that had previously been identified as relevant for brain function [67]. Data for these specific genes were available in 338 of the 400 brain regions. We assigned the anatomical location of each probe to one of 338 a priori defined parcels. For each parcel and each gene, we calculated the mean expression of that gene across all probes, after subtracting the mean expression of each probe for that gene [50]. Collectively, these calculations generated a data matrix of size 338 (parcels) by 16699 (mean expression across probes in that parcel for a given gene).

Gene coexpression between parcel i and parcel j is measured by the Pearson correlation coefficient r between gene expression values of parcel i and gene expression values of parcel j. To assess subnetwork specificity of gene coexpression, we first randomly sampled 500 genes and constructed a $338\times 338$ gene coexpression matrix. Second, from this coexpression matrix, we calculated the mean gene coexpression both within and between subnetworks. Between subnetwork coexpression is the mean of the gene coexpression values for pairs of nodes for which one node in the pair is located in subnetwork (A) and the other node in the pair is located in subnetwork (B). Third, we calculated the ratio of within- to between-network gene coexpression, where a ratio $>1$ indicates greater within-subnetwork coexpression than between-subnetwork coexpression. Finally, we compute the difference between this ratio and that expected in a nonparametric null model in which we randomly permuted subnetwork membership 1000 times, and recomputed coexpression ratios for each permutation. Next, we repeat the above process 10000 times, selecting a different random sample of 500 genes each time, generating a distribution of differences indicating whether the true ratio is larger than the mean of the null model ratios. To complete our statistical analysis, from this distribution we compute the probability of the difference being less than 0. We chose the approach of taking 10000 samples of 500 genes rather than taking a single sample of all the genes in order to demonstrate that the result is robust to the particular selection of genes, and to protect against our true ratio being influenced by outliers. For additional analyses controlling for the distance between regions, see the Supplemental Materials.

3.1 Dynamical network model

As motivated more fully in the Results section, we propose that two subnetworks of the frontoparietal system work in functional opposition to either couple or decouple the frontoparietal system from the default mode system. Specifically, we suggest that the functional connection between the frontoparietal and default mode systems is governed by the relative activation amplitudes of the two frontoparietal subnetworks. To further probe these relationships and their mechanistic underpinnings, we built a coarse-grained network in which each unit represented a particular brain system, and we simulated system dynamics with a canonical coupled oscillator model. More specifically, network activity was modeled by the normal form of a supercritical Hopf bifurcation (also referred to as the Stuart-Landau model), which describes the transition between a state of low activity and a state of oscillatory dynamics [46, 51]. We chose this model because (i) it permits the independent manipulation of oscillator amplitudes (Fig. S1), allowing us to further investigate the empirically observed relationships between activity and connectivity, and (ii) it is often used to model large-scale brain activity [35, 36, 55, 26, 71]. Following [26, 71], the local dynamics of the $j^{th}$ unit are given by the following equation:

[TABLE]

where $u_{j}=\rho_{j}e^{i\theta_{j}}=x_{j}+iy_{j}$, $\eta_{j}$ is drawn from a normal distribution to add Gaussian noise to the system, and $\xi$ scales the noise. In Eq.(3), the term $a_{j}$ is commonly called the bifurcation parameter. When $a_{j}<0$, the system goes to a low-activity fixed point and when $a_{j}>0$, the system obeys a stable limit-cycle solution with angular frequency $\omega_{j}$ and signal amplitude governed by $a_{j}$.

Following [26, 71], we model a network of interacting components by separating the system into its real and imaginary parts, and we link different components via the following set of coupled differential equations:

[TABLE]

and

[TABLE]

Here, G is the coupling strength and $n_{o}$ is the total number of oscillators, which in our case is four (default mode network, dorsal attention network, frontoparietal subnetwork (A), and frontoparietal subnetwork (B)). As suggested by [71] we set $\xi=0.02$, and we took $x_{j}$ as the oscillatory signal of interest. To estimate the frequency parameters, we empirically calculated the peak frequency of each of the 4 systems during the resting state, fit a normal distribution to the peak frequency values, and drew from the normal distribution. Similarly, to establish $D_{ij}$, the coupling of the network nodes, we calculated the mean structural connectivity estimated from diffusion tractography between system $i$ and system $j$ across all subjects. We integrated the equations using a time step of 0.01 seconds for 6 minutes, which was the approximate length of the task scans.

To estimate reasonable values for the coupling and bifurcation parameters, we conducted a parameter sweep (Fig. S2), and computed both the root mean square of the time series, and the synchrony among all oscillators using the Kuramoto order parameter, defined at time point t as:

[TABLE]

where $\phi_{j}(t)$ is the instantaneous phase of oscillator j at time t. To get a summary statistic for the entire time series, we took the mean of R across time (Fig. S2 E). The instantaneous phase was computed by taking the Hilbert transform of the unfiltered time-series.

For our baseline working point, we selected $a=-0.075$ and $G=0.1$ for all units, where the Kuramoto order parameter has an intermediate value, signifying a realistic dynamical regime between a state of no synchrony and a state of complete synchrony among all oscillators. Furthermore, at this working point the root mean square of the time series is higher than the noise level (Fig. S2 F). For each point in our grid defined by $a$ and $G$, we averaged results from 25 initializations of the system. Importantly, the bifurcation parameter of a node in this four oscillator system is linearly related to the root mean square of the time series values of that node (Fig. S1 B).

3.2 Statistical analysis

At several points, we calculate the statistical difference between outcome variables of the two sub-networks. To that end, we initially take a parametric approach, and then we confirm all of our findings using a non-parametric permutation-based approach. In all visualizations of statistical relationships, subject effects have been regressed out from the dependent variable.

3.2.1 Subnetwork connection differences

In testing our hypotheses, we often asked questions of the following form: Does the strength of the connection between subnetwork (A) and the default mode system differ from the strength of the connection between subnetwork (B) and the default mode system? For questions of this form, we used a multilevel model where each outcome variable (e.g., a measurement of connection strength) has attributes encoding subnetwork membership, task run, and subject identity. The multilevel model framework [75] accounts for the nested nature of the data (multiple scans nested within subject). We specified the model as:

[TABLE]

where $\mathrm{OutcomeVariable}_{it}$ is the outcome variable (i.e. connection strength) for person $i$ on run $t$; $\mathcal{B}_{0i}$ indicates the level of the outcome in subnetwork (A); $\mathcal{B}_{1i}$ indicates differences in the level of outcome associated with subnetwork (B) versus subnetwork (A); and $e_{it}$ are residuals.

Person-specific intercepts (from Level 1) were specified (at Level 2) as:

[TABLE]

and

[TABLE]

where $\gamma$ denotes a sample-level parameter and $u$ denotes residual between-person differences that may be correlated, but are uncorrelated with $e_{it}$. The multilevel model was fit with lme in R using maximum likelihood estimation. In the case of many outliers, we treat our data with robust models, rather than standard linear models. Robust models down-weight points of data, where the most outlying points are down-weighted most severely. Specifically, we implement robust multilevel models using robustlmm in R [49]. We note in the text whenever a robust multilevel model is used.

3.2.2 Repeated mesures correlations

Unless otherwise noted, we use a repeated measures correlation when examining the association between two continuous variables [8]. The repeated measures correlation accounts for non-independence among observations (due to multiple runs per subject) by using a form of analysis of covariance (ANCOVA) to adjust for between-person variance. The model is specified as:

[TABLE]

where $\mathrm{Measure1}_{it}$ is the value of variable one for subject $i$ during measurement occasion $t$, $\overline{\mathrm{Measure2}}_{i}$ is the mean value of the second variable in the $i$-th participant, $\mathrm{Subject}_{i}$ is a unique identifier for each participant, and $c(\mathrm{Measure2}_{i})$ is the covariate for the $i$-th participant and is equal to $\mathcal{B}(\mathrm{Measure2}_{it}-\overline{\mathrm{Measure2}}_{i})$, where $\mathcal{B}$ is the slope coefficient of the covariate. Like a Pearson correlation coefficient ($r$), the repeated measures correlation ($rrm$) is bounded by -1 to 1, and represents the strength of the linear association between two variables. The repeated measures correlation was estimated using the rmcorr package in R [8].

3.2.3 Nonparametric permutation-based approach

In addition to the multilevel linear model, we employ a complementary permutation-based approach. We begin with vectors $Y_{1}\in\mathcal{R}^{1\times 4n}$, $Y_{2}\in\mathcal{R}^{1\times 4n}$, and $S\in\mathcal{R}^{1\times 4n}$, and we wish to test whether there is a difference in the means of $Y_{1}$ and $Y_{2}$ against a null model. To construct the null model, for subject $i$ we find the two entries $J=(j_{1},j_{2})$ for which $S=i$. We then randomly reassign elements $J$ between $Y_{1}$ and $Y_{2}$, and repeat this procedure for all subjects to construct a null $Y^{null}_{1}$ and null $Y^{null}_{2}$ where we would expect the means to be equal. We then calculate the mean difference $d_{null}=\textup{mean}(Y^{null}_{1}-Y^{null}_{2})$. We re-permute and recalculate the mean 10000 times to establish a null distribution of the difference, and we determine a $p$-value for the true effect by calculating the proportion of null differences that are greater than the observed difference.

4 Results

4.1 Connectivity between default mode and frontoparietal systems tracks working memory performance.

We first sought to determine whether the strength of the connection between the default mode and frontoparietal systems tracks performance on a 2-back working memory task. For each participant, we calculated the average functional connectivity between all relevant pairs of brain regions, where one region of the pair was located in the default mode system and the other region of the pair was located in the frontoparietal system. We then estimated the relationship between behavioral accuracy and between-system strength. We found that the strength of the between-system connection was negatively correlated with performance on the task across subjects (repeated measures correlation, $r=-0.4731$, $p<0.001$, $DF=558$; Fig. 2A).

To assess task-specificity, we repeated the analysis using the imaging and performance data from the 0-back condition. During this condition, minimal working memory load exists. We found no significant relationship between behavioral performance and the strength of the connection between the default mode and frontoparietal systems (repeated measures correlation; $r=-0.0293$, $p=0.487$, $DF=562$; Fig. 2B-C). To assess the robustness of the connectivity-behavior relationship, we used an alternate behavioral measure (see Supplementary note 1 and Fig. S4 ), an alternate network parcellation (see Supplementary note 2 and Figs. S5, S6), and an alternate measure of functional connectivity (see Supplementary note 3 and Fig. S7). The results of these additional analyses serve to confirm our main findings.

4.2 The activity of the frontoparietal system is correlated with the strength of the functional connection between the frontoparietal and default mode systems.

After observing a statistical relation between working memory performance and the functional connectivity between the default mode and frontoparietal systems, we sought to better understand its potential drivers. We began by testing the hypothesis that system activity drives inter-system connectivity. We defined regional activity as the $\beta$-weight from a GLM fit to the regional BOLD magnitude during the 2-back task. To assess the activity of a system, we averaged those GLM $\beta$-weights across all regions in the frontoparietal system or across all regions in the default mode system. We observed that the frontoparietal system was more active than the default mode system (frontoparietal: mean $=2.74$, 95% CI [2.51, 2.98]; default mode: mean $=-3.09$, 95% CI [-3.33, -2.85]). We confirmed this difference in activity between the two systems using a multilevel model ($\beta=-7.4501$, $p<0.0001$, $t(1769)=-35.7$, $SE=0.20847$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$); see Section 3.2.

In assessing the relationship between system activity and inter-system connectivity, we found that the activity of the default mode system is unrelated to the strength of the connection between the default mode and frontoparietal systems (repeated measures correlation; $r=-0.0441$, $p=0.296$, $DF=562$; Fig. 2D). In contrast, we found that the activity of the frontoparietal system was negatively correlated with the strength of this connection (repeated measures correlation; $r=-0.0859$, $p=0.042$, $DF=562$; Fig. 2E). Given that frontoparietal activity is correlated with the strength of the inter-system connectivity, and that the inter-system connectivity is correlated with behavioral performance, it stands to reason that the frontoparietal activity should also be correlated with behavioral performance. Indeed we found that the activity of the frontoparietal system is positively related to behavioral performance (repeated measures correlation; $r=0.1019$, $p=0.016$, $DF=558$; Fig. S3 A), and this relation was not observed when considering activity of the default mode system (repeated measures correlation; $r=0.0384$, $p=0.363$, $DF=558$; Fig. S3 B).

4.3 Subnetworks of the frontoparietal system differentially modulate the strength of the functional connection between the frontoparietal and default mode systems.

The results reported in the previous section are correlative, and to press further we must posit a testable causal model. We propose that the dynamics of the frontoparietal system directly modulate the strength of the functional connection between the frontoparietal and default mode systems. Our candidate mechanism assumes that the frontoparietal system is composed of two distinct, non-overlapping sub-systems. Subnetwork (A) is posited to display dynamics that are correlated with the dynamics of the default mode system, while subnetwork (B) is posited to display dynamics that are anticorrelated with the dynamics of the default mode system (Fig. 3A). Further, we propose that the relative activity magnitudes of these two systems interact to tune the strength of the inter-system (frontoparietal - default mode) connection. When subnetwork (A) is highly active relative to subnetwork (B), the strength of the inter-system connection will be positive; conversely, when subnetwork (B) is highly active relative to subnetwork (A), the strength of the inter-system connection will be negative.

To assess the validity of this conceptual model, we first test the assumption that the frontoparietal system is composed of two distinct, spatially non-overlapping sub-systems. For each subject, we applied a weighted stochastic block model with $K=2$ to the subgraph of the adjacency matrix representing functional connections between frontoparietal regions, and we then extracted a group-representative partition using a consensus similarity method [30] (see Methods). The two subnetworks are shown on the cortex in Fig. 3B. To assess the statistical significance of this partition, we performed a non-parametric test, permuting the association of regions to subnetworks (see Methods). Against the null model, we found that this final consensus partition had a significantly higher log-likelihood using a multilevel model ($\beta=-492.57$, $p<0.0001$, $t(1769)=-71.0$, $SE=6.9371$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$).

To further validate their biological distinctness, we sought to determine whether the two subnetworks showed distinguishable patterns of gene expression. To this end, we quantified the average magnitude of gene coexpression for pairs of regions for which both regions were located within a single subnetwork. We also quantified the average magnitude of gene coexpression for pairs of regions for which one region of the pair was located in one subnetwork and the other region of the pair was located in the other subnetwork (see Methods). When tested against a non-parametric null model, we found that gene coexpression within subnetwork (A) was significantly higher than gene coexpression between subnetwork (A) and subnetwork (B) ($p=0.0057$). Similarly we found that gene coexpression within subnetwork (B) was higher than gene coexpression between subnetwork (B) and subnetwork (A) ($p=0.0111$). For complementary analyses demonstrating that this result cannot be explained by distances between regions, see Supplementary note 4. These findings support the notion that subnetwork (A) and subnetwork (B) are biologically distinct sectors of the frontoparietal system.

Next, we sought to test our proposition that – of the two frontoparietal subnetworks – one displays dynamics that are correlated with the dynamics of the default mode system, while the other displays dynamics that are anticorrelated with those of the default mode system. For each subject and each subnetwork, we calculated the strength of the functional connection between all pairs of regions for which one region of the pair is located in the subnetwork and the other region of the pair is located in the default mode system. Consistent with our hypothesis, we found that the activity of subnetwork (A) is positively correlated with the activity of the default mode system (mean $r=0.042$, 95% CI: [0.038, 0.046]; Fig. 3B left), while the activity of subnetwork (B) is negatively correlated with the activity of the default mode system (mean $r=-0.082$, 95% CI: [-0.088, -0.076]; Fig. 3B right). Using a multilevel model, we found that these correlations are significantly different ($\beta=-0.12469$, $p<0.0001$, $t(1769)=-37.4$, $SE=0.0033272$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$). To further unpack these findings, we also considered the relation between these two subnetworks and the dorsal attention system, which has been described as antagonistic to the default mode system in working memory tasks [5, 76, 28]. Consistent with this account, we found that the activity of subnetwork (A) was negatively correlated with the activity of the dorsal attention system (mean $r=-0.14$, 95% CI: [-0.18, -0.11]; Fig. S8 A right), while the activity of subnetwork (B) was positively correlated with the activity of the dorsal attention system (mean $r=0.12$, 95% CI: [0.11, 0.12]; Fig. S8 A left). Using a multilevel model, we found that these correlations are significantly different ($\beta=0.13057$, $p<0.0001$, $t(1769)=46.8$, $SE=0.0027885$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$).

To address the final proposition in our model, we sought to determine whether increased subnetwork (A) activity would lead to a stronger positive functional connection between the frontoparietal and default mode systems, while increased subnetwork (B) activity would lead to a stronger negative functional connection between the two systems. Because the summation of the subnetwork timeseries should reflect the complete frontoparietal timeseries, we also reasoned that when subnetwork (A) is more active (higher amplitude) relative to subnetwork (B), the frontoparietal signal would be more similar to the dynamics of subnetwork (A) than to the dynamics of subnetwork (B). To test these expectations, we began by quantifying subnetwork activity using the root mean square (RMS) of the subnetwork timeseries. We then used a single robust linear model to explain the strength of the functional connection between the default mode and frontoparietal systems by a linear combination of the activity of subnetwork (A) and the activity of subnetwork (B). Consistent with our hypothesis, we found that an increase in subnetwork (A) activity corresponded to a stronger positive functional connection between the two systems (estimate of regression coefficient $\beta=0.006535$, 95% CI: (0.00519, 0.00788)); Fig. 3D), while an increase in subnetwork (B) activity corresponded to a stronger negative functional connection between the two systems (estimate of regression coefficient $\beta=-0.0112$, 95% CI: (-0.01273, -0.0097)); Fig. 3E). In a complementary analysis, we used a single robust linear model to explain the strength of the functional connection between the dorsal attention system and the frontoparietal system by a linear combination of the activity of subnetwork (A) and the activity of subnetwork (B). We found that an increase in subnetwork (A) activity corresponded to a stronger negative functional connection between the two systems ($\beta=-0.00553$, 95% CI: (-0.00677, -0.00486)), while an increase in subnetwork (B) activity corresponded to a stronger positive functional connection between the two systems ($\beta=0.00989$, 95% CI: (0.00848, 0.011316)).

The results of the three tests described above serve to validate the formal structure of our model. Next we turned to an assessment of the relevance of this complex dynamical system for behavior. Specifically, we had observed previously that behavioral performance decreases as the correlation between the frontoparietal and default mode systems increases. Here we seek to explain that observation using the activity of the two subnetworks. Because subnetwork (A) is positively related to frontoparietal-default mode connectivity, which is itself negatively related to behavioral performance, we would expect subnetwork (A) activity to be negatively related to behavioral performance. Conversely, because subnetwork (B) is related negatively to frontoparietal-default mode connectivity, which is itself negatively related to behavioral performance, we would expect subnetwork (B) activity to be positively related to behavioral performance. To probe these relationships, we fit a single robust multilevel model with behavioral performance as the dependent variable, and subnetwork (A) activity and subnetwork (B) activity as independent variables. As expected, we found that subnetwork (A) activity is negatively related to behavioral performance ($\beta=-0.00271$, 95% CI: (-0.0041, -0.00132)), and subnetwork (B) activity is positively related to behavioral performance ($\beta=0.00273$, 95% CI: (0.00115, 0.00430)).

4.4 The structural role of the frontoparietal subnetworks

We next turn to an examination of what, if any, neuroanatomical support exists for the subnetwork-driven dynamics espoused in the previous section. Recent advances in network control theory have posited that changes in the activation of single brain regions can induce a propagation of activity along white matter tracts to affect distributed circuit behavior in a predictable fashion [42, 79]. Here we test this notion within the specific confines of our experiment, asking: Does the structural connectivity of frontoparietal subnetworks constrain how activity propagates to neighboring areas, thereby modulating the coupling between the frontoparietal and default mode systems? We hypothesize that subnetwork (A) is more structurally connected to the default mode system, while subnetwork (B) is more structurally connected to the dorsal attention system. This hypothesis is based on evidence that a higher number of white matter tracts between two regions can support stronger functional connectivity between them [73], allowing subnetwork (A) to strongly couple to the default mode system and allowing subnetwork (B) to strongly couple to the dorsal attention system.

We tested this hypothesis by calculating the strength of the structural connectivity between subnetworks and systems using diffusion imaging tractography (see Methods). Consistent with our hypothesis, we found that subnetwork (A) is more strongly connected to the default mode system (mean = 2.07, 95 %CI: (2.05, 2.08); Fig. 4A left) than is subnetwork (B) (mean = 0.676, 95%CI: (0.668, 0.685); Fig. 4A right). Using a multilevel model, we found that these correlations were significantly different ($\beta=-0.0013939$, $p<0.0001$, $t(1769)=-174.8$, $SE=7.9759\times 10^{-6}$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$). Similarly, we found that subnetwork (B) is more strongly connected to the dorsal attention system (mean = 1.54, 95% CI: (1.52, 1.56); Fig. 4B left) than is subnetwork (A) (mean = 1.03, 95% CI: (1.02, 1.04); Fig. 4B right). Using a multilevel model, we found that these correlations were significantly different ($\beta=0.00050413$, $p<0.0001$, $t(1769)=99.3$, $SE=5.0721\times 10^{-6}$, $n=2414$), and we found similar results using a non-parametric permutation test ($p<0.0001$). Further, across subjects we found that the less structurally connected the two subnetworks were, the better individuals tended to perform on the 2-back working memory task ($r=-0.0766$, $p=0.03$).

Collectively, these data suggest that the frontoparietal subnetworks might be optimally positioned in the structural connectome to mediate coupling between the default mode and dorsal attention areas, a coupling that is negatively correlated with performance (Fig. S12). To more directly test this notion, we calculated the regional boundary control (Eq. 2) with respect to the default mode and dorsal attention systems (Fig. 4C). Of the 20 regions for which boundary control exceeded the 95-th percentile, 9 were located in the frontoparietal system. To evaluate statistical significance, we compared these results to those of a non-parametric null model, in which we randomly permute the association between boundary control values and brain regions. We found that the probability that 9 or more of the top 20 regions fell within the frontoparietal system was significant with respect to the null model ($p=0.0019$). In summary, the data suggest that the frontoparietal system is structurally positioned to effectively control the coupling between the default mode and dorsal attention systems.

4.5 A reduced dynamical model for probing the relation between system activity and connectivity.

In the previous sections, we found evidence consistent with (but not proving) the causal notion that increased activity of subnetwork (A) drives stronger coupling between the frontoparietal and default mode systems, while increased activity of subnetwork (B) drives anticorrelation between the two systems. While it is difficult to prove the validity of such a causal model in healthy human participants, we can gather additional supportive evidence from in silico experiments exercising a formal computational model of the dynamical system. Specifically, we modeled the dynamics of a reduced four node network, with nodes representing the frontoparietal subnetwork (A), the frontoparietal subnetwork (B), the default mode system, and the dorsal attention system (Fig. 5A). We consider each unit in the network to be an oscillator, with dynamics described by the normal form of a Hopf bifurcation and with frequencies randomly sampled from an empirically measured distribution (see Sec. 3.1 for details). We coupled the four systems according to the mean weight of the structural connections between them, as estimated from diffusion tensor imaging tractography, averaged across subjects (see Methods). The model has two free parameters: (i) the global coupling parameter, which tunes the general capacity for synchronization, and (ii) the bifurcation parameter of each oscillator, which tunes the amplitude of the oscillator time series (Fig. S1). Following a broad parameter sweep, we selected parameter values to ensure a realistic dynamical regime between a state of no synchrony and a state of complete synchrony among all oscillators (see Methods and Fig. 5B).

Next, we implemented the model to probe the relation between subnetwork activity and connectivity, focusing initially on the connectivity between the frontoparietal and default mode systems. In agreement with our empirical results, we found that increasing the amplitude of subnetwork (A) activity increased the correlation between the frontoparietal and default mode unit time series (Pearson’s correlation coefficient $r^{2}=0.157$, $p=0.001$; Fig. 5C). Similarly, and again in agreement with our empirical results, we found that increasing the amplitude of subnetwork (B) activity decreased the correlation between the frontoparietal and default mode unit timeseries (Pearson’s correlation coefficient $r^{2}=0.75$, $p<0.0001$; Fig. 5D). To assess the reliability of these results, we performed the same numerical experiments for a range of coupling and bifurcation parameter values, across which the effects remained robust (Fig. S10). Next, we examined the complementary relation between subnetwork activity and the connectivity between the frontoparietal and dorsal attention systems. Again, consistent with our empirical results, we found that increasing subnetwork (A) activity decreased the correlation between the frontoparietal and dorsal attention unit timeseries (Pearson’s correlation coefficient $r^{2}=0.65$, $p<0.0001$; Fig. S11 A), and increasing subnetwork (B) activity increased the correlation between the frontoparietal and dorsal attention unit time series (Pearson’s correlation coefficient $r^{2}=0.246$, $p<0.0001$; Fig. S11 B). Collectively, this pattern of results offers support for a candidate mechanism in which subnetwork activity tunes functional connectivity.

5 Discussion

5.1 Functional coupling between default mode and frontoparietal systems tracks individual differences in working memory.

In a remarkably diverse range of tasks, the frontoparietal system tends to increase in activity while the default mode system tends to decrease in activity [21, 27]. A common and intuitive interpretation of these findings is that the two systems exist in competition, consistently displaying anti-correlated dynamics with one another. From a cognitive perspective, it has been proposed that such competition is indicative of opposing goals of the two systems [17, 43], and may enable toggling between the cognitive states supported by each system [28]. Importantly, the degree of competition between the frontoparietal and default mode systems has been linked to individual differences in performance on a range of tasks including the categorization of emotion in faces [82], the flanker task [21], and several tasks taxing human reasoning abilities [44]. Interestingly, attenuation of the competition between these two systems is related to decreases in executive function occurring in older age [77].

The nature of the functional connection between the frontoparietal and default mode systems appears to be flexible and task-dependent [45, 31]. Externally-directed tasks tend to elicit competition while internally-directed tasks tend to elicit cooperation [19, 78, 63]. For example, the two systems show competitive dynamics during visuo-spatial planning, and cooperative dynamics during autobiographical planning [76]. Likewise, the two systems show cooperative dynamics during a finger-tapping task, and competitive dynamics during a movie watching task [38]. Similarly, the two systems have been found to show cooperative dynamics during a goal-directed simulation task requiring imagination [39], as well as during the production of internal trains of thought [74]. It has been proposed that the default mode system is critically involved in self-referential tasks [4], and that coupling to the frontoparietal networks allows attention to be directed internally to these tasks, rather than to external stimuli [74, 4]. Our study extends this prior literature by demonstrating the existence of competitive dynamics between the two systems during the performance of an n-back working memory task, where the degree of competition is related to individual differences in behavioral performance.

5.2 Frontoparietal activity is related to its coupling with the default mode system.

While functional connectivity is estimated from activity time series, the generic relationship between activation and coupling remains far from understood. In the rather non-generic context of working memory tasks, prior studies have shown that the activity of the default mode decreases [53, 6], consistent with a suppression of mind-wandering and other self-referential processes that can hamper external goal-directed behavior [5]. In line with these findings, our results demonstrate that during n-back performance, the default mode system is less active on average than the frontoparietal system. Interestingly, we do not observe a relationship between the suppression of default mode activity and individual differences in performance (Fig. S3 B), consistent with a previous study of group-averaged n-back performance in patients with chronic pain and healthy matched controls [18]. Indeed, it has been suggested that a threshold may exist, below which working memory performance is unaffected by any further decreases in default mode activity [18]. Although our data do not directly support the threshold interpretation, as we observed no significant nonlinearities in the relationship between default mode activity and task accuracy, it would be interesting in future work to further test this hypothesis in groups with markedly lower performance, for example young children or individuals with executive dysfunction.

In contrast to default mode activity, frontoparietal activity directly tracked individual differences in behavior, with higher frontoparietal activation being associated with greater task accuracy (Fig. S3 A). The frontoparietal system is known to support a range of executive functions including trial-by-trial control, enabling the allocation of attention to trial-specific information [65, 22]. Such processes are heavily recruited during the n-back task, where successful performance is dependent upon flexibly redirecting attention to similarities and differences between the current stimulus and prior stimuli. More generally, the frontoparietal system is thought to act as a flexible hub during task performance [24], altering its connectivity with the default mode system in a task-dependent fashion [34, 76]. Indeed, it has been proposed that the frontoparietal system is functionally (and structurally) interposed between the default mode and dorsal attention systems [22], modulating their coupling in a task specific way [28]. Hellyer and colleagues proposed that increased frontoparietal activity promotes persistent stable states, whereas increased default mode activity may allow for transitions between cognitive states, which in turn would be undesirable during tasks requiring directed and fixed attention [45]. Thus, an active frontoparietal system may also help to maintain a persistent competitive relationship with the default mode, aiding in improved behavioral performance.

Lastly, we investigated the relationship between activity in the frontoparietal and default mode systems and inter-system coupling. It has been hypothesized that competition between the two systems may enable maximally disjunctive computational activities [21]. While we found that more inter-system competition was related to increased frontoparietal activity, we did not find that functional connectivity between the two systems was related to default mode activity. This pattern of results suggests that competition may not always support maximally disjunctive activities, and instead provides further evidence that – in the activity range observed here – default mode activity does not relate to either functional coupling with the frontoparietal system or behavioral performance on the n-back working memory task.

5.3 Subnetworks of the frontoparietal system modulate coupling to the default mode system.

The frontoparietal system flexibly alters its functional connections dynamically according to current task demands [24], perhaps controlling the strength of connectivity between cognitive systems [11]. In order to support such a broad range of cognitive states, it has been suggested that the frontoparietal system may be composed of subnetworks, where each subnetwork subserves a specific cognitive state [34, 29]. Yet the mechanics of this subservience is far from understood. We proposed that the frontoparietal system is composed of two disjoint subnetworks with anti-correlated dynamics, where subnetwork (A) displays activity that is correlated with the activity of the default mode system, and where subnetwork (B) displays activity that is anti-correlated with the activity of the default mode system. Using an unsupervised clustering algorithm informed by an explicit model of network architecture, we demonstrated that the frontoparietal system is decomposable into two subnetworks with distinct patterns of functional connections. Prior work by Dixon and colleagues has found a similar division of the frontoparietal system [29]. The authors employed a hierarchical clustering method to establish a data-driven partition of the frontoparietal system into two components, the first component being more functionally connected to the default mode, similar to our subnetwork (A), and the second component being more functionally connected to the dorsal attention system, similar to our subnetwork (B). Our findings critically extend these prior observations by providing evidence that by supporting competition between the two systems, subnetwork (B) may be critical to working memory performance, while subnetwork (A) may be less involved in working memory, and more closely linked with introspective processes.

In addition to their functional distinguishability, we also demonstrated that these two subnetworks displayed distinct patterns of gene coexpression. In line with this observation, it is interesting to note that prior work has suggested that cortical regions responsible for different cognitive functions can express different genes [40], and that gene coexpression provides a partial explanation for patterns of functional connectivity [67]. In agreement with these findings, our results demonstrate more similar patterns of gene expression within subnetworks than between subnetworks, an effect that cannot be explained by inter-regional distance. Notable prior work has also suggested a link between structural connectivity and gene expression [54], also supported by our finding that the two genetically dissimilar subnetworks display differing patterns of white matter connectivity. It is interesting to speculate that these genetic dissimilarities are partially responsible for the subnetworks’ differential capacity to tune the coupling between the frontoparietal and default mode systems via relative changes in activity amplitudes.

Observing a similar functional division of the frontoparietal system [34], Fornito and colleagues found that increased connectivity between frontoparietal areas and the default mode was correlated with improved task performance on an introspective recollection task. Here in a non-introspective task characterized by competitive coupling between the frontoparietal and default mode systems, we observe that increased connectivity between the frontoparietal subnetwork (B) and the default mode system was associated with poorer performance. This difference in results can be explained by the observation that the frontoparietal system is a flexible hub [24], mediating the connection in particular between the dorsal attention and default mode systems [22]. During introspective tasks it may be advantageous for the frontoparietal system to be more connected to the default mode and less connected to the dorsal attention system, buffering the internal trains of thought from external stimuli [74]. Conversely, during externally-directed tasks (like the n-back), it may be advantageous for the frontoparietal system to be more connected to the dorsal attention system and less connected to the default mode system, buffering externally directed attention from internal thoughts and mind-wandering [5]. In line with this argument, we found that during performance of the 2-back task, subnetwork (B) is positively connected to the dorsal attention system and negatively connected to the default mode system, while the opposite is true for subnetwork (A).

5.4 Structural connectivity of the frontoparietal subnetworks.

Motivated by prior work demonstrating that structural and functional connectivity share topographic similarities [73, 23], we extend our study to multimodal data to better understand the potential structural drivers constraining the manner in which activity in frontoparietal subnetworks impinges on functional coupling in other systems. We found that subnetwork (A) had fewer structural connections to the dorsal attention system than it had to the default mode system, and also than subnetwork (B) had to the dorsal attention system. Similarly, we found that subnetwork (B) had fewer structural connections to the default mode system than it had to the dorsal attention system, and also than subnetwork (A) had to the default mode system. Collectively, these data suggest that structural connectivity may play a role in constraining the functional dynamics of the two subnetworks. To better understand this role, we draw on tools from network control theory, which provides a mechanistic framework to link network structure to functional network dynamics [42, 15, 56]. Regions with high levels of boundary control are theoretically posited to have the capacity to steer the brain to different states by coupling and decoupling cognitive systems [42]. Our results demonstrate that the frontoparietal subnetworks are situated within the white matter architecture in a manner that can drive the system’s coupling with the dorsal attention and default mode systems. The effective activity of this hub of network control may be cognitively advantageous by, during externally-directed tasks, buffering the attentional systems from the internally-directed processes of the default mode, and likewise during internally-directed tasks buffering the default mode system from external stimuli.

5.5 A reduced dynamical model for probing the relation between system activity and connectivity.

We proposed a model wherein the coupling between the frontoparietal and default mode systems is modulated by activity levels in the oscillatory dynamics of two frontoparietal subnetworks. To test the validity of this putative mechanism, we employed a simplified dynamical model of brain system activity using oscillators coupled by empirically-determined structural connectivity [26, 71, 55]. This model allowed us to directly test our hypotheses by allowing us to alter subnetwork amplitude, and subsequently to observe the resulting synchronization (functional connectivity) between the frontoparietal and default mode systems. The results from these simulated experiments agree with our hypotheses and demonstrate that, when reducing the system to a network of 4 coupled oscillators, the modulation of subnetwork amplitude governs intersystem coupling in a way that matches the empirical discoveries. This finding not only provides a mechanistic explanation for our results, but also more generally illuminates how the frontoparietal system may mediate the effective functional coupling between the dorsal attention and default mode systems.

It has been proposed that the frontoparietal system is functionally interposed between the dorsal attention and default mode systems, altering their functional coupling in a task-specific way. In particular, during externally directed tasks the frontoparietal system may engage with the dorsal attention system and disengage with the default mode system [74]. Conversely, during internally directed tasks the frontoparietal system disengages with the dorsal attention system and engages with the default mode system [5]. Together, these complementary processes are thought to effectively segregate external stimuli from internal trains of thought during tasks that require more focus on one of the two. Our results demonstrate that this complementarity may be accomplished by the differential activation and deactivation of two frontoparietal subnetworks. Indeed, our results suggest that increasing subnetwork (A) activity engages the frontoparietal system with the default mode and increasing subnetwork (B) activity disengages the frontoparietal system from the default mode. Importantly, we found that the opposite relationship is also true when considering the functional coupling between the frontoparietal and dorsal attention systems: increasing subnetwork (A) activity disengages the frontoparietal system with the dorsal attention system and increasing subnetwork (B) activity engages the frontoparietal system with the dorsal attention system. The frontoparietal system is thought to flexibly reconfigure functional connections with the dorsal attention and default mode systems in a task-dependent way [28], and our results offer evidence for a mechanism by which this flexible reconfiguration may be achieved.

5.6 Methodological Limitations

Several methodological limitations are pertinent to this work. In the text, the frontoparietal system was divided into two subnetworks, creating group-level subnetworks. First, it should be noted that the subnetworks could also have been studied at the level of single individuals, although such granularity could hamper the ability to draw group-level conclusions. Second, these subnetworks were defined using a community detection algorithm based on a specific generative model. There exist several methods to find communities within networks, each with its own set of underlying assumptions. As a result, careful consideration must be taken when selecting the appropriate method of community detection. The flexibility of the WSBM makes it the most reasonable choice given the data used here.

While the employed dynamical model indeed recovered empirical results and allowed us to postulate mechanistic explanations, it is important to point out some of its methodological limitations. First, as the most straightforward way to test our main hypotheses, we considered only four major brain systems: the dorsal attention system, default mode system, and the two subnetworks of the frontoparietal system. However, it is crucial to be aware of the fact that these systems are embedded into a larger network, and hence their behavior is determined not only by intrinsic parameters and coupling to one another, but also by interactions with other brain systems. While a comprehensive computational study of large-scale brain dynamics is beyond the scope of this work, it would indeed be interesting in future endeavors to attempt to systematically understand – via modeling – the role that other brain systems might have in modulating the activity and functional connectivity of systems related to working memory.

Another limitation of our model revolves around the scale at which it operates. In particular, because our main empirical findings concern system-level dynamics – rather than dynamics at the scale of individual neurons or parcels – we assumed, for simplicity, that the different units in our computational model represented different brain systems. This coarse-grained approach is beneficial for a number of reasons. For example, it simplifies our analysis, and allowed us to focus explicitly on macroscopic, system-level drivers of various results, and therefore directly compare output from the model to corresponding empirical findings. However, although informed by experimental data, it is critical to acknowledge that such a setup is a great simplification of the true system, and allows little room for understanding how observations at the level of brain systems arise from interactions between more microscopic structural components.

Along this same vein, the dynamics of each brain system were described inherently phenomenologically, via a Hopf bifurcation model [51], which has been utilized in studies concerned with the interplay between network structure and dynamics in general [37, 68] and also in computational studies on brain network dynamics more specifically [26, 71, 55]. In particular, such dynamics indeed capture the oscillatory nature of observed brain system activity, but do not embody a biophysically precise description of neuronal activity. Therefore, the model cannot attempt to describe the emergence of brain system dynamics from dynamical processes on smaller scales. Although in this investigation we have chosen, as a first step, to employ a canonical model with few parameters that favors simplicity and interpretability, building and analyzing more realistic and detailed, multi-scale models is indeed an important area of ongoing research.

6 Conclusion

Our study provides evidence for a mechanism by which the dynamics of the frontoparietal system may drive working memory performance. Two distinct subnetworks within the frontoparietal system play a role in modulating the functional coupling between the frontoparietal and default mode systems during the performance of an n-back working memory task, which may help buffer the externally-directed attentional system from internal trains of thought, and lead to improved behavioral performance. We found that the position of the two subnetworks within the white matter scaffolding constrains the distinct function of each: one is structurally tied to the dorsal attention system (and drives the frontoparietal system into competition with the default mode), and the other is structurally tied to the default mode (and drives the frontoparietal system into cooperation with the default mode). We extend these descriptive observations by building a computational model instantiating and validating the putative mechanism, and we bolster our findings with corroborating differences in gene expression in the two subnetworks. Together, our findings contribute to a holistic view of working memory by linking activity, functional connectivity, structural connectivity, and gene expression, and present one way of understanding how these four modes work in concert to support cognitive processes necessary for working memory.

Acknowledgments

We thank Karolina Finc, Lorenzo Caciagli, and Ursula Tooley for helpful comments on an earlier version of this manuscript. Our team acknowledges support from the John D. and Catherine T. MacArthur Foundation, the Alfred P. Sloan Foundation, the ISI Foundation, the Paul G. Allen Foundation, the Army Research Laboratory (W911NF-10-2-0022), the Army Research Office (Bassett-W911NF-14-1-0679, Grafton-W911NF-16-1-0474, DCIST- W911NF-17-2-0181), the Office of Naval Research, the National Institute of MentalHealth (2-R01-DC009209-11, R0-MH112847, R01-MH107235, R21-MH106799, R01-MH113550), the National Institute of Child Health and Human Development (1R01-HD086888-01), National Institute of Neurological Disorders and Stroke (R01-NS099348), and the National Science Foundation (BCS-1441502,BCS-1430087, NSF PHY-1554488 and BCS-1631550). The content is solely the responsibility of the authors and does not necessarily represent the official views of any of the funding agencies.