Brain network dynamics during working memory are modulated by dopamine and diminished in schizophrenia

Dynamical brain state transitions are critical for flexible working memory but the network mechanisms are incompletely understood. Here, we show that working memory performance entails brain-wide switching between activity states using a combination of functional magnetic resonance imaging in healthy controls and individuals with schizophrenia, pharmacological fMRI, genetic analyses and network control theory. The stability of states relates to dopamine D1 receptor gene expression while state transitions are influenced by D2 receptor expression and pharmacological modulation. Individuals with schizophrenia show altered network control properties, including a more diverse energy landscape and decreased stability of working memory representations. Our results demonstrate the relevance of dopamine signaling for the steering of whole-brain network dynamics during working memory and link these processes to schizophrenia pathophysiology.


Participants and study design
All participants provided written informed consent for protocols approved by the Institutional Review Board of the medical faculty in Mannheim. For the first study including healthy controls and patients with schizophrenia, we included a total of 202 subjects (178 healthy controls, 24 individuals with schizophrenia, see Supplementary Table 1 For the second, pharmacological intervention study, 17 healthy individuals completed a subject-and observer-blind, placebo-controlled, randomized three-period cross-over study (see Table 2). Exclusion criteria included a regular consumption of drugs or history of drug or alcohol abuse; systolic blood pressure (SBP) greater than 140 or less than 90 mm Hg, and diastolic blood pressure (DBP) greater than 90 or less than 50 mm Hg; notable resting bradycardia (heart rate (HR) <40 bpm) or tachycardia (HR >90 bpm); use of any medication or herbal remedies taken within 14 days prior to randomization into the study or 5  One subject was excluded from the analysis due to an excessive body-mass index (BMI > 30).

N-back paradigm
The visual N-back paradigm is a well-established and reliable working memory task consisting of a high memory load (2-back) and an attention control condition requiring motor response (0-back) (5-7).
Specifically, a diamond-shaped stimulus containing a number from 1 to 4 was presented every 2 seconds (see Supplementary Figure 1). In the 0-back condition, subjects were required to press the button on the response box corresponding to the number currently displayed on the presentation screen. In the 2-back condition, subjects were required to press the button on the response box corresponding to the number presented two stimuli before the number currently displayed on the presentation screen. Stimuli were presented in alternating blocks of either 0-back or 2-back conditions. In each condition block, 14 stimuli were presented. Each condition block was repeated 4 times. Task performance was measured by accuracy (defined as the percent of correct answers) and reaction time (defined as the time span between stimulus onset and button press) for each condition separately.

Atlas construction
To combine structural and functional brain imaging data, we first constructed a brain atlas that equally well respects functional and anatomical features. We transformed a recently published multimodal atlas (8)  Combining the two atlases resulted in 374 regions that covered cortical and subcortical structures. A full list of regions included in the combined atlas can be found in Supplementary Table 3.

1.4.
Network Control Theory

Model assumptions
The framework of network control theory that we have employed here is based on linear system model that relies on several assumptions:  Linearity: Linearity implies that the system evolves linearly over time. This is probably not an accurate description of most brain dynamics, but non-linear dynamics can be locally approximated by linear dynamics for macroscale brain circuits (10,11).  Freedom from noise: Neural systems at all time and spatial scale are not noise free. Nevertheless, it seems reasonable to consider that the salient features of our model do no depend on noise. This is aided by our definition of brain states as meta states that are defined as statistical patterns of brain activation over repeated measurements.

Continuous versus discrete dynamical models
Linear dynamical system can be studies using different time-system, either continuous or discrete. A discrete-time system assumes that the system evolves in discrete time steps whereas a continuous-time system models continuously changing dynamics. The choice depends on which time assumptions best reflects the neural dynamics at study. Our choice was motivated by a) our definition of brain states as nondiscrete entities of dimensional brain activity summarized over extended brain region and b) the computational traceability. As we considered brain regions containing millions of neurons as network nodes, we assumed that each brain region`s state change is more heterogenous and therefore better represented as a continuous system. For reasons of computational feasibility, we used a discrete-system approach for the computation of the suboptimal trajectories.

Stabilization of the dynamical system
If not stabilized, a dynamical system could potentially increase infinitely over time. Such extreme brain states would be neurobiologically implausible due to the finite energy resources of the brain. We therefore chose to normalize the system by decreasing the average weight of the connectome such that it goes to zero over time: (1) norm = |λ( ) max |+ − Here, I denotes the identity matrix, A denotes the structural connectivity matrix and |λ(A)max| denotes the largest eigenvalue of the system. To normalize the system, we must specify the parameter c, which determines the rate of stabilization of the system. We choose c = 1, meaning that all modes decay and the system goes to zero over time. Within the range of brain states that converge to zero over time, we cannot make statements regarding whether any of these intermediate brain states are biologically plausible or are realized in human brains. Further work integrating experiment and theory is needed to more clearly define types of implausible states, and their respective mechanisms (e.g., metabolic, electrical, informational, or other physical constraints). A more in-depth and mathematical introduction and discussion can be found in Refs. (12)(13)(14).

Time horizon
The time horizon T specifies the time over which the control input is applied and the system can be pushed from one state to the other. It determines how quickly the system is required to converge and therefore small values might give the system insufficient time to reach the target state, making it hard to control. In theory, T is dimensionless if not coupled to external time domains. As we do not intend to model an evolving process in real time, we chose T = 1 to use a normalized time, in line with previous works (15), which allows the system to have adequate time to be controlled (12). For a systematic investigation of the influence of T on control processes in the context of brain networks, we refer the reader to (12), as well as to more general introductions of the control processes in linear dynamics (14,16).

On the relationship between BOLD signal and control theory measures
In the control theory framework, the control input (u) of a node and the state of that same node (x) are highly interrelated. For example, if we consider a simplified system consisting of only one node, then the control energy E necessary to change the state of that node from an initial state (x0) to a target state (xT) is basically a function of the squared difference in that node's state As our definition of brain states is based on estimates that depend on BOLD activity, in such a simplified system control energy would not give any additional information other than the usual contrast images 2back-0back. However, if we consider a more complex system with more than one node, and where all nodes are connected via either direct or indirect links as summarized in the connectivity matrix A, then the control energy of a single node is not a simple function of the squared differences in its state but additionally accounts for the influence of other connected neighbors.

On the use of control theory as a statistical framework
In our analysis, we apply control theory as a statistical and theoretical tool to answer questions based on the theoretical "dual-state" framework regarding neurobiological properties of brain function. Translating and transferring across these three levels (control theory as a statistical tool, dual-state theory as a nonlinear theoretical framework, brain imaging data defining meta-level brain states) is challenging and requires (reasonable) simplifications. The hypotheses that we aim to test are based on the dual-state theory framework, which also uses the terminology of brain states and energy. In this framework, states and transitions are based on non-linear dynamics, corresponding to attractor basins, which translate to stable reoccurring activation patterns in neuronal ensembles (17)(18)(19). Abstracting these concepts to largescale dynamics of brain macro-circuits provides the underlying basis for the idea that we aim to investigate here: relatively stable "meta"-level brain activation patterns as identified by neuroimaging (including all the caveats of the assumption of stationarity of brain activations measured by functional magnetic resonance imaging) populate a state-space for which we aim to identify the brain regions that are responsible for maintaining and shifting those activation patterns. To answer these cognitive neuroscience questions, we use network control theory as a toolkit that makes these questions computationally tractable in a linear dynamical system framework enabling us to quantify the associated "energy cost" of transitions on a brain region level. This effort requires certain (reasonable) assumptions, in particular to assume an equivalence between states defined by neuroimaging and states defined in the control theory framework, as well linear and continuous transitions between those states. Future work integrating biophysical models of taskinduced brain activity in combination with network control theory and tailored imaging paradigms is critically needed to provide further evidence for the assumed relationships (and distinctions) between actual data, network control tools, and the theoretical framework.

Polygenic co-expression index calculation
Previous publications have shown that gene sets defined using co-expression networks and selected for their association with the genes DRD1 and DRD2 provided replicable predictions of n-back-related brain activity and behavioral indices (20)(21)(22)(23). Weighted Gene Co-expression Network Analysis [WGCNA (24) (20). Both PCI-based predictions were significantly replicated in an independent post-mortem dataset, while controlling for ethnicity. The translational effect of these two scores on brain activity during nback has been assessed and replicated across multiple samples, which combined amount to approximately 600 participants (23,27,28).
It is important to note that these dopamine-related genetic effects are large in magnitude compared to those estimated by polygenic risk score approaches that focus on epidemiological data, rather than on molecular processes. The DRD2-PCI we developed (23)

Null models of structural brain networks
To study the impact of structural brain networks on control properties, we repeated the computation of control energy using a randomized null model of the individuals' structural brain networks that preserves the degree distribution and ensure fully connected networks. Null models were created using the group as between-subject factor: group by null_vs_data interaction, F(1,99) = 0.289, p = 0.592). These analyses suggest that human brain structural networks are in some form optimized to control brain state transitions independent of diagnostic status.

Null models of spatial activity patterns
To study the impact of the spatial distribution of activity patterns on control properties, we repeated the computation of control energy and spatially randomized individuals' brain activation patterns.
Randomization was done using the randperm function in matlab for the paired vectorized brain activation patterns (related to 0-and 2-back) followed by recomputation of control energy. This procedure was repeated 100 times and the averaged brain-wide control energy over all 100 iterations for each subjects was used in the subsequent analysis. In line with our expectation, control energy increased significantly for randomized networks in both groups (repeated measures ANOVA with model_vs_data and transition as within-subject factors, HC: main effect of model_vs_data, F(1,174) = 6.995, p = 0.009). Further analysis revealed no interaction with patient status (repeated measures ANOVA with null_model_vs_data and transition as within-subject factors and group as between-subject factor: group by null_vs_data interaction, F(1,99) = 3.904, p = 0.056). These analyses suggest that the spatial distribution of brain activity patterns is important for minimizing control effort, but individuals with schizophrenia have a differently, potentially less organized activity pattern than healthy controls.

Robustness to choice of parcellation scheme
To demonstrate the robustness of the results to our choice of parcellation scheme, we repeated our analysis using a recently published functionally defined atlas comprising a similar number of areas (29).
Specifically, we used the "Gordon" template (29) consisting of 333 regions that are functionally derived from resting-state connectivity analyses. Data were reprocessed using the same pipeline as for the main analysis and all parameters were kept identical in the subsequent analysis. Notably, we replicated all main results (see Supplementary Table 1), indicating that our reported findings are robust to the choice of parcellation scheme.

Robustness to choice of edge definition
To demonstrate the robustness of our results to our selection of connectivity measure, we repeated our analysis using the number of streamlines normalized by the respective size of the regions to construct structural connectivity matrices (15). All parameters were kept identical in the subsequent analysis. All main results could be replicated (see Supplementary Table 4), indicating that our findings are robust to the choice of edge definition.

Impact of medication and duration of illness on control properties
In patients, the potential relationship between control energy and stability, antipsychotic drug dose

Pharmacological validation using Risperidone
To demonstrate the robustness of our pharmacological intervention of dopaminergic signaling, we additionally analyzed the data of the Risperidone condition in the same subjects. Risperidone also preferentially targets D2 receptors, but also affects D1, adrenergic, serotoninergic and histaminergic pathways. Using the same models and covariates as in the main analysis, we detected a trend-wise

Null results for gene score and imaging associations
As mentioned in the main text, D1 receptor expression-related gene scores predicted stability of both states

Relation to previous gene score and imaging associations
To demonstrate the added value of our analysis, we extracted the BOLD parameters from the voxels D2 expression-related gene score predicted the control energy of both state transitions (0-to 2-back: b = -0.075, p = 0.042; 2-to 0-back: b = -0.139, p = 0.062, age, sex, brain stability, difference in brain activity, first 5 genetic PCA components and brain activity at [-29 53 24] as peak voxel reported in Fazio et al. as covariates of non-interest).
All our previous associations remained significant when controlling for the BOLD parameters, suggesting that our data explains different variance.

Comparison to more conventional SPM analyses
To demonstrate the added value of our control metrics, we performed the following two more conventional SPM analyses: 1) Association of control energy and stability with PCI scores for DRD1 and DRD2: We performed a conventional GLM analysis for Placebo versus Amisulpride, using a paired t-test of the 2-back>0back contrast in SPM12. We could not detect any significant clusters, either for the Placebo > Amisulpride nor the opposite contrast, even at a lenient threshold of P < 0.001 uncorrected. These results suggest that network control theory can detect biologically meaningful effects that cannot be detected by more conventional SPM analyses.

2) Stability and Control Energy in Schizophrenia:
We performed a conventional SPM analysis for individuals with schizophrenia versus healthy control, using the 2-back > 0back contrast. At a lenient threshold of P < 0.001 uncorrected, we detected two main significant cluster of voxels in the common region-of-interests for the N-back task

Suboptimal trajectories
As mentioned in the main text, the variability in suboptimal trajectories was greater in schizophrenia (rm-ANOVA: main effect of group: F(1,98) = 4.789, p = 0.031, controlling for age, sex, DTI tSNR, brain state activity difference). These results remained significant after additionally accounting for the stability of both states and for the control energy of both transitions (rm-ANOVA: main effect of group: F(1,95) = 11.2, p = 0.001).

Supplementary Figure 1: N-back task design
Design of the N-back task: Stimuli were presented in blocks of either 0-back or 2-back conditions. There was no additional control or resting condition. In the 0-back condition, subjects were instructed to press the button on the response box corresponding to the number currently displayed on the presentation screen.
Here, the red numbers below the screen images indicate correct responses. In the 2-back condition, subjects were instructed to press the button on the response box corresponding to the number presented two stimuli before the number currently displayed on the presentation screen. Here, the red numbers below the screen images indicate the correct responses. Each condition block lasted 30 seconds and was repeated four times in an interleaved manner.