Abstract
Recent theoretical advances have motivated the hypothesis that the periaqueductal gray (PAG) participates in behaviors that involve changes in the autonomic control of visceromotor activity, including during cognitively demanding tasks. We used ultra-high-field (7 tesla) fMRI to measure human brain activity at 1.1 mm resolution while participants completed a working memory task. Consistent with prior work, participants were less accurate and responded more slowly with increasing memory load—signs of increasing task difficulty. Whole-brain fMRI analysis revealed increased activity in multiple cortical areas with increasing working memory load, including frontal and parietal cortex, dorsal cingulate, supplementary motor area, and anterior insula. Several dopamine-rich midbrain nuclei, such as the substantia nigra and ventral tegmental area, also exhibited load-dependent increases in activation. To investigate PAG involvement during cognitive engagement, we developed an automated method for segmenting and spatially normalizing the PAG. Analyses using cross-validated linear support vector machines showed that the PAG discriminated high versus low working memory load conditions with 95% accuracy in individual subjects based on activity increases in lateral and ventrolateral PAG. Effect sizes in the PAG were comparable in magnitude to those in many of the cortical areas. These findings suggest that cognitive control is not only associated with cortical activity in the frontal and parietal lobes, but also with increased activity in the subcortical PAG and other midbrain regions involved in the regulation of autonomic nervous system function.
SIGNIFICANCE STATEMENT Functional neuroimaging in humans has shown that cognitive control engages multiple corticostriatal networks and brainstem nuclei, but theoretical advances suggest that the periaqueductal gray (PAG) should also be engaged during cognitively demanding tasks. Recent advances in ultra-high-field fMRI provided an opportunity to obtain the first evidence that increased activation of intermediate and rostral portions of lateral and ventrolateral PAG columns in humans is modulated by cognitive load. These findings suggest that cognitive control is not solely mediated by activity in the cortex, but that midbrain structures important for autonomic regulation also play a crucial role in higher-order cognition.
Introduction
Working memory and cognitive control are supported by cortical areas that are part of the frontoparietal and the dorsal attention networks (Fox et al., 2006; Spreng et al., 2013) and interconnected subcortical areas in the thalamus and dorsal striatum (Parent and Hazrati, 1995). Although these networks are well characterized during working memory (Wager and Smith, 2003; Owen et al., 2005), cognitive tasks also engage autonomic responses linked to energy mobilization (Critchley et al., 2003; Thayer et al., 2012). Recent research has focused on the energetic demands associated with cognitive control measured both metabolically (Gailliot and Baumeister, 2007) and in terms of subjective effort and computational costs (Shenhav et al., 2017). Integrating findings from these viewpoints, researchers have proposed that dorsal cingulate cortex, anterior insula, and ventromedial prefrontal cortex are involved in weighing the costs and benefits of engaging in cognitive control (Shenhav et al., 2013; Westbrook and Braver, 2015). Consistent with this view, these regions are also involved in maintaining energy balance in the face of changing environmental demands, referred to as allostasis (McEwen and Wingfield, 2003; Sterling, 2012; Barrett and Simmons, 2015; Satpute et al., 2019). Many of these regions modulate multiple brainstem nuclei that are thought to be involved in cognitive control because of their neuromodulatory influence on cortex, including the dopamine-rich substantia nigra (SN) and the ventral tegmental area (VTA) (Westbrook and Braver, 2016), which may play a role in several computations involved in behavioral control (Montague et al., 2004), including effort allocation (Niv et al., 2007) and updating and gating access to working memory (O'Reilly and Frank, 2006; D'Ardenne et al., 2012).
Another brainstem structure that is well suited to regulate the demands of engaging cognitive control is the periaqueductal gray (PAG). The PAG's profile of connectivity positions it to be important in behaviors and mental states that involve motivated performance. The mammalian PAG receives inputs from portions of the anterior cingulate cortex (Ongür and Price, 2000; Johansen-Berg et al., 2008) that are also involved in visceromotor regulation. Tract-tracing studies demonstrate that different PAG columns project to different autonomic nuclei in the brainstem (Cameron et al., 1995) that are critical for implementing allostasis (Sterling and Laughlin, 2015). Furthermore, different PAG columns receive topographically organized inputs from different portions of medial prefrontal cortex, including monkey homologs of aMCC (i.e., dorsal cingulate cortex; An et al., 1998), which are thought to be important for cognitive control (MacDonald et al., 2000; Shenhav et al., 2016). One human fMRI study of working memory at 3 tesla (van Ast et al., 2016) has shown increases in midbrain activation during cognitive control, with effects in multiple nuclei including the PAG. Human studies of PAG connectivity suggest strong connectivity with the aMCC in particular (Kong et al., 2010) and that different PAG subregions have different patterns of subcortical connectivity (Coulombe et al., 2016). Recent monkey studies have also confirmed that projections from multiple medial frontal areas project to sympathetic effectors in peripheral organs (Dum et al., 2016). Thus, although there are differences in the anatomy and connectivity of prefrontal–PAG pathways across species, and different species may perform cognitive tasks differently (Carruthers, 2013), this ensemble of connections suggests that a fundamental role of the PAG may be to maintain allostasis in the face of demanding cognitive tasks.
The goal of the present study was to use 7 tesla fMRI to test whether the PAG plays a role in cognitive demand. Prior work has yet to test this question due to methodological limitations. Conventional neuroimaging methods blur neural signals of interest, which can both lead to misidentification of midbrain signals, and an inability to localize individual PAG columns. Advances in human neuroimaging methods at higher field strengths partially overcome these issues, revealing PAG activity with improved spatial precision (Satpute et al., 2013; Faull and Pattinson, 2017; Sclocco et al., 2018). We scanned participants and contrasted brain activity during three-back and one-back working memory tasks. We hypothesized that the PAG would exhibit greater activation during three-back trials than during one-back trials. Our imaging protocol localized effects to different PAG columns, which are associated with distinct patterns of behavioral and autonomic responses (Carrive and Morgan, 2012). In particular, the ventrolateral PAG is associated with sympathetic inhibition (Johnson et al., 2004), whereas the dorsal PAG is associated with sympathetic activation (Dean et al., 2016). Thus, we investigated whether activity during cognitive demand engaged particular subregions of the PAG. To compare PAG engagement with that of structures with established roles in working memory (O'Reilly and Frank, 2006), we examined other brainstem and subcortical regions including the SN and VTA (anatomically defined in Pauli et al., 2018), in addition to frontoparietal and dorsal attention networks, which have been reliably linked to working memory (Smith et al., 2009).
Materials and Methods
Participants.
This study included 24 participants (Mage = 24.15 years, SD = 6.33 years, eight female). All recruited participants were between the ages of 18 and 40 years, were right-handed, had normal or corrected to normal vision, were not pregnant, were fluent English speakers, had no known neurological or psychiatric illnesses, and were recruited from the greater Boston area. Participants were excluded from the study if they were claustrophobic or had any metal implants that could cause harm during scanning. All participants provided written informed consent and study procedures were completed as approved by the Partners' Healthcare Institutional Review Board.
Experimental paradigm.
Participants completed a visual N-back working memory task during fMRI scanning (based on the design in Gray et al., 2003 and van Ast et al., 2016). The task was administered during a single scanning session and pseudorandomly alternated between blocks of either one-back or three-back conditions with 10 trials in each block. The session included a total of 12 blocks that were presented in ABBA or BAAB order (counterbalanced across subjects) and were each followed by a 25 s rest period. Each block began with a cue indicating the current task (one-back or three-back). The task was designed with fixed proportions of 20% targets and 80% nontargets (12.5% of which were lures) in each block. Lure trials were defined as two-back matches in the one-back blocks and lures in three-back blocks were either two- or four-back matches. The proportion of lure trials was the same for both one-back and three-back blocks to equate the requirement for resolving interference.
On every trial, a letter (Q, W, R, S, or T) was presented at the center of the visual field for 2 s followed by a fixation cross for 2 s. Participants were instructed to respond with a button press when the letter on the screen matched the one presented n trials ago. The task was administered in MATLAB (RRID:SCR_001622, The MathWorks), using the Psychophysics Toolbox extensions (RRID:SCR_002881, Kleiner et al., 2007). Visual stimuli were projected so participants could view stimuli on a mirror fixed to the head coil used for data acquisition. Responses were recorded using an MR-compatible button box. Response times and hit rates for target trials were computed using this signal and served as the primary behavioral outcomes. Before and after the N-back task, participants completed a series of self-report items (using a five point bipolar visual analog scale) to indicate the extent to which they felt awake (vs sleepy), energetic (vs fatigued), engaged (vs bored), pleasant (vs unpleasant), calm (vs tense), and relaxed (vs restless). They were also asked to indicate the extent to which the task would be (was) demanding, whether they would have (had) enough resources to complete the task and whether they would (did) perform well. Changes in self-report measures were computed by subtracting ratings made before the task from those made after its completion. Before scanning, participants completed a practice session on a laptop PC. These practice sessions mirrored the task used during scanning, including alternating blocks of one-back and three-back trials. A total of 48 trials were completed during practice. The goal of the practice session was to ensure that participants were adequately familiar with the N-back task.
fMRI acquisition.
Gradient-echo echo-planar imaging BOLD-fMRI was performed on a 7 tesla Siemens MRI scanner. Functional images were acquired using GRAPPA-EPI sequence: echo time = 28 ms, repetition time = 2.34 s, flip angle = 75°, number of slices = 123, slice orientation = transversal (axial), phase encoding = A > P, voxel size = 1.1 mm isotropic, gap between slices = 0 mm, field of view = 205 × 205 mm2, GRAPPA acceleration factor = 3; echo spacing = 0.82 ms, bandwidth = 1414 Hz per pixel, partial Fourier in the phase encode direction: 7/8. A custom-built 32-channel radiofrequency coil head array was used for reception. Radiofrequency transmission was provided by a detunable band-pass birdcage coil.
Structural images were acquired using a T1-weighted EPI sequence: echo time = 22 ms, repetition time = 8.52 s, flip angle = 90°, number of slices = 126, slice orientation = transversal (axial), voxel size = 1.1 mm isotropic, gap between slices = 0 mm, field of view = 205 × 205 mm2, GRAPPA acceleration factor = 3; echo spacing = 0.82 ms, bandwidth = 1414 Hz per pixel, partial Fourier in the phase encode direction: 6/8. This sequence was selected so that functional and structural data would have similar spatial distortions (Renvall et al., 2016) to facilitate coregistration and subsequent normalization of data.
fMRI preprocessing.
Preprocessing was performed using SPM12 software (RRID:SCR_007037, Wellcome Department of Cognitive Neurology) and included spatial realignment (Friston et al., 1996), affine coregistration of the mean functional and structural data with six degrees of freedom (Friston et al., 1995), and a first-pass spatial segmentation and normalization of anatomical data using the Computational Anatomy Toolbox (CAT12 Toolbox; Dahnke et al., 2013; Gaser and Dahnke, 2016). The deformations computed on anatomical data were then applied to functional data. First-level models were estimated for each subject using these normalized functional data. After model estimation, two different additional normalization procedures were applied: (1) a whole-brain normalization to MNI152 space (Fonov et al., 2011) using a group-specific DARTEL template created from the T1-weighted structural EPI data to provide a greater degree of anatomical overlap and more precise tissue segmentations and (2) a PAG-specific normalization that used model residuals and tissue-specific information from segmentation.
To localize activity within the PAG, we developed an automated procedure to segment and normalize functional data. The procedure was based on manual methods used in a previous 7 tesla imaging study examining PAG responses to emotional pictures (Satpute et al., 2013). This procedure involved identifying the cerebral aqueduct by: (1) finding voxels that have large model residuals for each participant and have a similar spatial profile across participants; (2) creating a mask of the PAG by dilating the aqueduct using a sphere with a 2-voxel (2.2 mm) radius and restricting its spatial extent to gray matter voxels identified during segmentation of brain tissue into gray matter, white matter, and CSF (Ashburner and Friston, 2005); (3) creating a custom group template of the PAG using DARTEL (which was approximately in MNI space due to the first pass normalization described above); and (4) warping functional data within the PAG to MNI152 space for evaluation and visualization. The consistency of this segmentation and normalization procedure was evaluated by computing the average Dice coefficient (percentage overlap) between all pairs of subjects. The group-level segmentation of the PAG was further subdivided into five subregions using k-means clustering using on voxel locations in MNI space as input. The PAG template and software for implementing this procedure are available on GitHub at https://github.com/canlab.
fMRI analysis.
To estimate brain activation during the N-back task, the preprocessed functional time series were modeled using general linear models as implemented in SPM12. Separate models were estimated for each subject; they included separate regressors for the one-back and three-back blocks, which were specified as boxcar functions with a variable duration based on the reaction time on every trial. These regressors were convolved with SPM's canonical hemodynamic response function. Six nuisance regressors modeled the effect of participant movement based on motion parameters estimated during realignment (translation in the x, y, and z directions in addition to roll, pitch, and yaw). A simple contrast between the three-back and one-back conditions served as the primary outcome of interest. After model estimation, contrasts of parameter estimates were smoothed using a 4 mm FWHM Gaussian kernel. Group-level effects were implemented using one-sample t tests. Thresholds for parametric maps were selected using multiple-comparisons correction based on the false discovery rate (FDR q < 0.05; Benjamini and Hochberg, 1995; Genovese et al., 2002) with a minimum extent of 5 voxels.
Complementing mass univariate assessment, we additionally conducted multiple confirmatory tests based on extant functional parcellations. Differences in brain activity for the three-back versus one-back contrast were computed within functional parcellations that have previously been linked to working memory and executive function. These include large-scale resting-state networks from Yeo et al. (2011), a parcellation of the striatum based on functional coactivation during a wide variety of tasks (Pauli et al., 2016), and anatomically defined subcortical and midbrain structures involved in reinforcement learning and decision making (Pauli et al., 2018). After voxelwise analyses, one-sample t tests were conducted for each region of interest (ROI) (i.e., seven resting-state networks and five striatal zones).
To investigate how PAG activation changed during cognitive processing, we developed a probabilistic atlas of the PAG using subject-specific segmentation of the aqueduct and surrounding gray matter and normalization to group (and MNI) space with DARTEL. Creating a group-specific PAG template through iterative normalization via DARTEL produced individual segmentations that were more similar to one another (mean pairwise Dice coefficient = 0.709, 95% CI = 0.667–0.740) than segmentations based on conventional nonlinear deformations to MNI space (mean pairwise Dice coefficient = 0.381, 95% CI = 0.347–0.413). Differences in the segmentation across subjects were primarily located near the inner and outer boundaries of the PAG, as opposed to differences in overall shape or rigid alignment.
To localize brain responses within individual PAG columns, we used an automated segmentation and normalization procedure to create a population-level template of the PAG. We parcellated the group-specific PAG template (described above) into multiple distinct columns (dorsal PAG, left and right ventrolateral PAG, and left and right lateral PAG) by performing unsupervised clustering of voxels based on their rostrocaudal position and angular position relative to the cerebral aqueduct (estimated using principal components analysis). To exclude superior portions of the dorsal raphe nucleus, which is just anterior to caudal portions of the aqueduct, we implemented k-means with a six-cluster solution and excluded the cluster in this area (see Fig. 1). Unlike parcellations based on functional activation or connectivity, which change over time and are context dependent (Salehi et al., 2018), this procedure provides an unbiased, spatially informed method of subdividing the PAG that provides parcels that approximately correspond to the functional columns identified in animal research.
To quantify the extent to which patterns of PAG activity differed between the three-back and one-back conditions, we developed linear support vector machine (SVM) classifiers using a leave-one-subject-out cross-validation scheme (implemented in MATLAB). The goal was to differentiate patterns of PAG activation during the one-back and three-back conditions versus rest. In this framework, classification models were trained on data from all but one subject, and data from the remaining subject (i.e., one activation map for the three-back condition and another for the one-back condition) were used for testing and estimating out-of-sample performance. We used a simple classification function (i.e., linear kernels with a default hyperparameter value of C = 1) to ease interpretation of models and to minimize overfitting (Norman et al., 2006). This approach is comparable to other linear methods commonly used in the neuroimaging literature (Misaki et al., 2010). This procedure was repeated until each subject had been used for testing (21 folds total). Signal detection metrics (i.e., sensitivity, specificity, positive predictive value, and area under the curve or AUC) were computed using cross-validated distances from the SVM hyperplane. Bootstrap resampling (b = 1000 samples) was used to identify which voxels reliably contributed to classification, using a normal approximation for inference.
We conducted repeated random subsampling to estimate observed power for detecting differences between the three-back and one-back conditions. The primary goal of this analysis was to differentiate between well established effects in cortex and brainstem effects that may not be as robust (due to challenges in registration or increased physiological noise). To minimize bias due to circularity (Kriegeskorte et al., 2009), spatially defined PAG subregions and regions predictive of the term “working memory” in the Neurosynth database (RRID:SCR_006798; Yarkoni et al., 2011) were selected as ROIs. For the Neurosynth meta-analysis, significant regions in the “working memory” reverse inference map were thresholded at an extent >200 voxels. For each ROI, distributions of the mean signal, t statistic, and associated p-values were estimated using repeated random subsampling (without replacement). Power was estimated as the proportion of random samples that exhibited a significant effect at α < 0.05.
Peripheral physiological recording.
Peripheral autonomic nervous system activity was measured using an AD Instruments PowerLab data acquisition system with MR-compatible sensors. Data were acquired during the entire N-back task (i.e., both during task blocks and intermediate rest periods). A pulse transducer (AD Instruments) was placed on the index finger to measure heart rate and recorded at 1 kHz with a 0.9 to 5 Hz band-pass filter. A respiratory belt with a piezo-electric transducer (UFI) was placed around the torso at the level of the sternum to measure changes in respiratory rate; this signal was recorded at 1 kHz with a 0.5 Hz low-pass filter. Skin conductance was measured using wired Ag/AgCl finger electrodes, with sensors containing isotonic paste with signals amplified via an FE116 GSR amplifier and recorded at 1 kHz.
Peripheral physiological analysis.
Physiological time series data were analyzed using custom scripts in MATLAB to calculate respiration rate, heart rate, and skin conductance separately for the one-back and three-back conditions. Peaks in the finger pulse and respiratory signals were identified using a sliding window analysis to identify heart and respiratory rate. Starting at the onset of scanning, peaks in the finger pulse signal were identified in the 10 s of data preceding every point of data acquisition. Thirty-second windows were used to calculate respiration rate and mean skin conductance. Separate general linear models were estimated for each autonomic measure, each with stimulus onsets and durations identical to those used for fMRI analysis. These models accounted for the effect of task (three-back vs one-back) in addition to linear changes over time (to account for habituation).
fMRI model comparisons.
To evaluate the functional significance of PAG activation, a series of regression models were developed that evaluated relationships between individual differences in PAG activation (indicated by distance from the SVM hyperplane; see fMRI analysis section) and variation in behavioral performance and autonomic activity. All models included predictors for the main effect of condition (one regressor), subject-specific means (19 regressors), and a constant term. Additional models were created by adding a single regressor characterizing individual differences in hit rate, response time, mean skin conductance during each condition, mean heart period, and mean respiration rate. A total of six models were constructed: a model that added all possible behavioral and autonomic predictors (the full model, 26 unique regressors), and models adding hit rate (22 unique regressors), response time (22 unique regressors), skin conductance (22 unique regressors), heart period (22 unique regressors), and respiration rate (22 unique regressors). Evidence for each model was computed by computing the Bayesian information criterion (BIC) for each model, computing differences between models, and estimating BIC weights (wBIC), which estimate the strength of evidence for each model in the set (Wagenmakers and Farrell, 2004). Parameter estimates for highly probable models were inspected to determine which behavioral and autonomic factors best explained differences in PAG activity.
Results
Behavioral and peripheral physiological findings
Consistent with a large literature using the N-back task, participants performed worse on the three-back compared with the one-back condition, with hit rates lower in the three-back condition [median = 0.750, interquartile range (IQR) = 0.271] compared with the one-back condition (median = 1, IQR = 0) (z = −3.90, p < 0.0001, Wilcoxon signed-rank test). Also, response times were slower for the three-back condition (median = 0.906 s, IQR = 0.284 s) compared with the one-back condition (median = 0.620 s, IQR = 0.147 s) (z = −3.92, p < 0.0001). Self-report measures indicated that participants felt more engaged (mean difference = 0.682, SD = 0.780, 95% CI = 0.336–1.038, p = 0.0005) and energetic (mean difference = 0.364, SD = 0.848, 95% CI = 0.167–0.786, p = 0.004) after the task relative to a prescanning baseline. Self-report measures of expected versus actual performance indicated that they performed worse than anticipated (mean difference = −0.429, SD = 0.811, 95% CI = −0.798 to −0.060, p = 0.025).
Autonomic nervous system measures indicated that participants had similar respiration rates during the three-back condition (median = 17.412 breaths/min, IQR = 5.465 bpm) and the one-back condition (median = 17.090 breaths/min, IQR = 4.452 bpm) (z = 0.486, p = 0.627). Mean skin conductance was lower during the three-back condition (median = 2.787 μS, IQR = 3.794 μS) compared with the one-back condition (median = 2.994 μS, IQR = 3.830 μS) (z = −2.314, p = 0.021). No difference in heart period (median = 0.783 s, IQR = 0.121 s) was observed between conditions (z = 0.714, p = 0.475). These behavioral and peripheral physiological findings showed increasing cognitive demand in the three-back compared with the one-back condition without substantial impact on peripheral autonomic activity.
Functional neuroimaging findings
A focused analysis contrasting activity within the PAG during the three-back and one-back conditions (Fig. 1a) revealed a cluster spanning intermediate and rostral lateral/ventrolateral PAG in the right hemisphere that exhibited cognitive load-dependent increases in activity (peak t-statistic = 4.57, MNIx,y,z = [4, −28, −6], extent = 67.5 mm3, voxelwise FDR q < 0.05; Fig. 1b). This localization was corroborated using a data-driven spatial clustering of PAG voxels based on their rostrocaudal and angular positions (see Materials and Methods; Fig. 1c): Voxels that exhibited the greatest load-dependent effects were located in the right hemisphere within the lateral (mean difference = 0.371, SD = 0.354, t = 4.814, p = 0.0001, 95% CI = 0.210–0.532) and ventrolateral columns (mean difference = 0.366, SD = 0.417, t = 4.021, p = 0.0007, 95% CI = 0.176–0.555).
Load-dependent brain activity during the N-back task. a, Anatomical segmentation and normalization of periaqueductal tissue, including the superior portion of the dorsal raphe and PAG. Rendering depicts a group-based probabilistic atlas of human PAG overlaid on the ICBM152 template in MNI space. The color map indicates the proportion of subjects with overlapping PAG segmentations at each voxel. b, Load-dependent modulation of PAG activity during the N-back task. The color map reflects t-statistics from a one-sample t test conducted on contrasts of three-back vesus one-back conditions, showing greater activation during the three-back task in lateral/ventrolateral PAG. c, Parametrization and unsupervised clustering of PAG voxels within the group-based PAG atlas (thresholded at 20% overlap) based on their spatial location (with rostrocaudal and angular dimensions identified using principal components analysis and conversion to polar coordinates) produce a five-cluster columnar organization. The PAG can be unrolled and displayed along these two principal dimensions (right). d, Classification of PAG activity along high and low levels of memory load. Inset plots depict effect size maps (Cohen's d) for the three-back and one-back conditions. Receiver operating characteristic curve shows high discriminability of high and low working memory load when testing on data from an independent subject (AUC = 0.972, two-alternative forced choice, p < 0.0001, permutation test). e, Mass-univariate whole-brain analysis reveals multiple cortical and subcortical regions that are more active as working memory load increases in the N-back task. Color map reflects t-statistics from a one-sample t test conducted on contrasts of three-back versus one-back conditions. Yellow (and orange) colors indicate greater activation during the three-back task, whereas blue colors indicate higher levels of activity during the one-back task.
In addition to testing for voxelwise differences in activity, we quantified the separability of response profiles in the whole PAG as a function of working memory load by training linear SVM classifiers to differentiate the one-back and three-back conditions. The classifier performed very well, separating one-back from three-back blocks with high sensitivity and specificity in new individual participants in cross-validated analyses (sensitivity = 95% (95% CI = 84%–100%), specificity = 95% (95% CI = 83%–100%), positive predictive value = 95% (95% CI = 83%–100%), AUC = 0.972 (Fig. 1e). Voxel weights differentiating between high and low memory load were moderately similar to those identified using voxelwise subtraction analysis (spatial Pearson correlation = 0.3196, SD = 0.0208, 95% CI = 0.3207–0.3421). The weights that most reliably contributed to classification (peak z-score = 3.70, MNIx,y,z = [3, −28, −8], extent = 8 mm3) were located in the same portion of the lateral/ventrolateral PAG as the peak univariate effect. This anatomical convergence across both analysis methods indicate that increasing activity in lateral and ventrolateral PAG reliably discriminates between different levels of memory load.
To characterize the functional importance of this particular profile of PAG activation, we performed post hoc model comparisons (using the BIC; see Materials and Methods) that evaluated which behavioral (i.e., hit rate and response time) and autonomic (i.e., skin conductance, heart period, and respiration rate) measures explained individual differences in PAG pattern responses above and beyond the average effect of memory load. This analysis revealed the strongest evidence for a model using all autonomic and behavioral measures to predict the degree of PAG activation (the full model; wBIC = 0.7268) followed by a model including response time as the main variable of interest (the response time model; wBIC = 0.2714). These models were far better accounts of PAG activation than other models (all wBIC < 0.01). Regression coefficients for response time were positive in both the full model (β̂ = 3.7865, SE = 1.2856; t14 = 2.9453; p = 0.0106) and the response time model (β̂ = 4.479, SE = 1.3122, t14 = 3.4132, p = 0.0031), indicating that greater l/vlPAG activation was associated with slower reaction times, independent of individual variation in autonomic reactivity, or the average effect of memory load. In addition to being linked to slower reaction times, the full model revealed a negative association between skin conductance and PAG activation (β̂ = −6.8308, SE = 2.903, t14 = −2.3530, p = 0.0338), with no other measures being individually significant.
A whole-brain analysis contrasting brain activity acquired during the three-back and one-back conditions demonstrated that multiple cortical and subcortical regions also exhibited load-dependent increases in activation (Fig. 1e). These regions include areas that are linked to allostasis and effort, including dorsal cingulate cortex extending into supplementary motor area and dorsal anterior insula, as well as other regions related to cognitive control including dorsolateral prefrontal cortex, superior frontal sulcus, caudate body, and inferior parietal lobule, all of which have been strongly implicated in working memory (Wager and Smith, 2003; Owen et al., 2005).
Consistent with previous imaging studies of working memory and cognitive control, brain activation at the level of large-scale networks (Yeo et al., 2011) showed load-dependent increases in the average activity of the frontoparietal (mean difference = 0.715, SD = 0.929, t = 3.528, p = 0.0021, 95% CI = 0.292–1.139) and dorsal attention networks (mean difference = 0.600, SD = 0.754, t = 3.646, p = 0.0016, 95% CI = 0.257–0.943), but no other networks. Examining activation differences within different striatal subdivisions (Pauli et al., 2016) revealed load-dependent differences specifically in the posterior caudate (mean difference = 0.319, SD = 0.510, t = 2.862, p = 0.0096, 95% CI = 0.086–0.551), consistent with evidence that this portion of the striatum is engaged during manipulations of executive function, cognitive control, and working memory (Pauli et al., 2016).
Given higher levels of noise in fMRI measures of brainstem activity (physiological, motion, and scanner-related) compared with the cortex, we also estimated effect sizes within PAG and other cortical regions commonly implicated in working memory. To estimate effect sizes for cortical, subcortical, and brainstem areas in an unbiased manner, we created functional ROIs of expected activation from an automated meta-analysis (conducted with Neurosynth; Yarkoni et al., 2011) using the term “working memory” (for a similar application, see van Ast et al., 2016) and the spatially defined subregions estimated from our PAG segmentation procedure (Fig. 2). This analysis revealed that the observed effects in vlPAG and lPAG in the right hemisphere were detected with similar levels of power as those in the cortex, exceeding 80% power with sample sizes as small as 15 subjects.
Power analysis using functionally independent ROIs. a, ROIs identified based on spatial clustering of PAG voxels (shown in the sagittal slices on the left) and automated meta-analysis of the term “working memory” using Neurosynth (Yarkoni et al., 2011). b, Power curves estimated using repeated random subsampling (up to the full sample of 21 subjects). The vlPAG and lPAG (in the right hemisphere) exhibit >80% power with <15 subjects (left). All cortical regions exhibit 80% or greater power with <15 subjects (right).
As a final exploratory analysis, we examined the involvement of midbrain and subcortical structures with high numbers of dopamine receptors (Pauli et al., 2018) because of their established involvement in cognitive control and working memory. Many of these subcortical structures, including the hypothalamus (Keay and Bandler, 2015), VTA (Omelchenko and Sesack, 2010), and SN (Carrive and Morgan, 2012), have reciprocal connections with the PAG. This analysis revealed load-dependent modulation of activity in a number of regions (Fig. 3): SN pars compacta (mean difference = 0.3477, SD = 0.4153, t = 3.8373, p = 0.0010, 95% CI = 0.1587–0.5368) and reticulata (mean difference = 0.2749, SD = 0.3193, t = 3.9460, p = 0.0008, 95% CI = 0.1296–0.4203), parabrachial pigmented nucleus (mean difference = 0.2496, SD = 0.4450, t = 2.5707, p = 0.0182, 95% CI = 0.0471–0.4522), VTA (mean difference = 0.4619, SD = 0.7228, t = 2.9288, p = 0.0083, 95% CI = 0.1329–0.7909), hypothalamus (mean difference = 0.1951, SD = 0.4138, t = 2.1611, p = 0.0430, 95% CI = 0.0068–0.3835), and mammillary nucleus (mean difference = 0.4824, SD = 0.6424, t = 3.4411, p = 0.0026, 95% CI = 0.1900–0.7748).
Load-dependent activity in dopamine-rich areas during a working memory task. a, ROI analysis reveals multiple subcortical and brainstem regions that are more active as working memory load increases in the N-back task. Black crosses reflect mean and SE of contrasts between the three-back and one-back conditions. Colored distributions reflect smoothed kernel density estimates of the data distribution. Individual circles indicate the value of the contrast for each subject. Asterisks denote regions that are significant at p < 0.05 uncorrected. b, Three-dimensional rendering of dopamine-rich areas (and neighboring landmark structures, i.e., extended amygdala, mammillary nucleus, and red nucleus) included in the analysis (Pauli et al., 2018). Areas that exhibited significant load-dependent effects are labeled. The parabrachial pigmented nucleus is not visible, as it is obscured by the SN (c). Mass-univariate analysis within the dopaminergic ROIs reveals multiple regions that show load-dependent effects. Color map reflects t-statistics from a one-sample t test conducted on contrasts of three-back versus one-back conditions. Orange colors indicate greater activation during the three-back task. Effects are thresholded at a voxelwise cutoff of qFDR < 0.05.
Discussion
The PAG is commonly thought to be involved in regulating systems of the body's internal milieu, particularly when dealing with moments of motivated behavior when energy expenditures are higher, such as during episodes of pain, sexual behavior, distress, and feeling threatened (Fanselow, 1991; LeDoux, 2012; Keay and Bandler, 2015). Our observations in this study broaden our understanding of the PAG's functional repertoire by demonstrating that changes in PAG activity are reliably associated with different levels of cognitive demand. These effects were large and robust: Differences in cognitive load were highly discriminable in terms of their associated profiles of PAG fMRI response at 7 tesla. Increased activation in intermediate and rostral portions of lateral/ventrolateral PAG predicted whether a participant was performing a one-back or three-back working memory task with 95% accuracy. This level of prediction shows that activation of subregions of the PAG, when measured in the context of a cognitively demanding task, is an effective indicator of cognitive load; indeed, it is as effective as frontostriatal circuits and dopaminergic pathways (Curtis and D'Esposito, 2003; Cools and D'Esposito, 2011; Chatham and Badre, 2015). Our findings support the view that the brain regions important for cognitive demand include the PAG and specific midbrain and striatal regions rich in dopamine, which are important in other (e.g., reinforcement learning) contexts.
Our brainstem findings provide a first look at how PAG subregions relate to increasing cognitive demand. We found that activity in the lateral and ventrolateral PAG in particular was modulated by increasing working memory load. Classically, lPAG/dlPAG has been implicated in “active” coping responses and vlPAG with “passive” coping in nonhuman animals when engaging in defensive behavior (Bandler et al., 2000). For example, stimulation of ventrolateral PAG produces freezing/immobility (Assareh et al., 2016) and reduced sympathetic output (Bago and Dean, 2001). Attention-capturing stimuli have also been shown to involve reduced sympathetic/increased parasympathetic activity when orienting; for example, viewing negative graphic images is associated with reduced heart rate and increased PAG activity (Hermans et al., 2013). In our prior work, too, we found preliminary evidence that activation in the right ventrolateral PAG during negative image viewing was associated with reduced subjective affective arousal (Satpute et al., 2013). The present observations of activation in the ventrolateral PAG and slight decreases in sympathetic nervous system activation (indexed by mean skin conductance) during the three-back task are partially consistent with the possibility of sympathetic disengagement during the more difficult three-back task. Nevertheless, given the variability in autonomic responses and the activation across both lateral and ventrolateral PAG regions that we observed, it could also be the case that three-back working memory performance here involved a mix of “challenge” (sympathetic engagement) and “threat” (sympathetic disengagement) responses (Bandler et al., 2000). Indeed, the precise pattern of PAG subregions involved may depend on an individual's unique response profile to increasing cognitive demand, a possibility that could be explored in a future study that is optimized to investigate individual differences (see Schaefer et al., 2006 for related work examining variation in amygdala responses).
The similarities in PAG activation observed across human imaging studies (i.e., the present results and vlPAG activation observed in Satpute et al., 2013) could reflect a shared mechanism involved in coordinating autonomic outcomes associated with passive coping. Ventrolateral PAG activity could be involved in inhibiting sympathetic outflow (Bago and Dean, 2001; Johnson et al., 2004; Hermans et al., 2013), which could lead to subjective changes in arousal, with complex effects on cognitive performance depending on the need for attentional orienting (proparasympathetic) versus metabolic activation (prosympathetic). A related functional role of vlPAG during working memory is to marshal cognitive resources, that is, to suppress competing motivations (Gear et al., 1999; Sprenger et al., 2012; Geuter et al., 2016) to engender “task focus.” This function may not require major changes in autonomic activity directed to the periphery (to the extent that these changes serve to prepare the body for skeletomotor action), but would be associated with diverting resources to support neural systems underlying cognitive performance. Consistent with this view, it has recently been shown that passive coping or “freezing” responses, which are considered to be mediated by lateral/ventrolateral PAG, also facilitate attention and perception (Roelofs et al., 2010; Lojowska et al., 2015). Alternatively, these overlapping effects could be the product of neighboring but distinct neural populations that are linked to different behavioral, autonomic, or experiential outcomes (e.g., subvocalization, autonomic nervous system changes during the three-back task, or the subjective feeling of mental fatigue). Future studies that independently manipulate these different constructs are necessary to determine whether PAG activity observed here is representative of a task-specific or domain-general process such as reducing action to facilitate cognitive or perceptual processing. Moreover, it remains to be determined whether incorporating these putative functions of the PAG (e.g., upregulating cognitive resources or decreasing motor activity) will improve the performance of computational models of working memory (O'Reilly and Frank, 2006).
Our observation of PAG engagement during cognitive control raises multiple questions regarding its involvement with frontoparietal networks, frontostriatal circuits, and dopaminergic pathways. Although we have demonstrated that the PAG is an important neural substrate for cognitive control, linking it to task behavior and autonomic activity, it is not clear whether it is better characterized as an integrated component of these large-scale brain systems or if it should be considered as a regulatory output concerned primarily with regulating visceromotor, autonomic, and behavioral endpoints. Although it is uncertain whether neurons project directly from the PAG to cortical areas involved in cognitive control, the PAG does project directly to the amygdala (Rizvi et al., 1991) and to the thalamus (Krout and Loewy, 2000), which in turn project to dorsal cingulate cortex and neighboring medial PFC (McDonald, 1991; Hoover and Vertes, 2007), raising the possibility that the PAG is a critical component of different cortical systems. Additional ultra-high-field work exploring the connectivity structure and dynamics of the PAG, for example, using resting-state (Coulombe et al., 2016) or task-based (Faull and Pattinson, 2017) measures of functional connectivity, is necessary to determine whether it should be considered a component of large-scale brain systems.
In addition to our findings in the PAG, we identified several dopamine-rich midbrain nuclei that have reciprocal connections to the PAG (Kirouac et al., 2004; Geisler et al., 2007), which exhibited activation that scaled with working memory load. A large body of work has implicated prefrontal dopamine in the maintenance and control of working memory (Sawaguchi and Goldman-Rakic, 1994; Romanides et al., 1999). More recent high-resolution imaging in humans (D'Ardenne et al., 2012) demonstrates that transient activity in the VTA and SN acts to gate information to the prefrontal cortex, consistent with computational models of working memory (Gruber et al., 2006; O'Reilly and Frank, 2006). Our findings of sustained activation differences in dopamine-rich midbrain areas during the N-back task are consistent with accounts that consider the tonic influence of dopamine in biasing action and modulating activity in frontal cortex (Westbrook and Braver, 2016). Although midbrain dopamine is associated with the effort required for securing rewards (Salamone et al., 2007), there is also evidence that the VTA is activated in response to stressors and other aversive stimuli (Brischoux et al., 2009), consistent with its role in regulating effort more generally (Niv et al., 2007). It is possible that the aversive nature of engaging cognitive control during the N-back task explains the involvement of these midbrain nuclei, rather than cognitive load per se (see Westbrook and Braver, 2015 for a related discussion). Regardless of the underlying cause, which is a topic of future research, our findings corroborate that VTA and SN activity is modulated by working memory load, further validating their involvement in cognitive control.
Considering the energetic costs required during cognitive tasks broadens the scope of behaviors that accounts of PAG function must explain. Our findings suggest the PAG may be involved in conveying information about the energetic costs of implementing cognitive control. In particular, the vlPAG, which contains dopaminergic neurons (Lu et al., 2006) and is innervated by the VTA (Suckow et al., 2013), may be an important component of a large-scale system that weighs the costs and benefits of engaging effortful control to maintain allostasis. PAG inputs from the brainstem about autonomic and energetic demands (e.g., via the nucleus of the solitary tract and the parabrachial nucleus) could be compared against neuromodulatory outputs from nearby midbrain structures (e.g., locus ceruleus, SN, and VTA) to ready the brain and body for anticipated demands. This reconceptualization of the PAG considerably broadens the possible behaviors in which it may be implicated and suggests that the neural basis for cognitive control is distributed across interconnected cortical, subcortical, and brainstem structures, rather than being predominantly localized to the cortex.
Footnotes
This work was supported by the National Institutes of Health (National Cancer Institute Grant U01 CA193632, National Institute of Biomedical Imaging and Bioengineering Grants K01 EB019474 and P41 EB015896, and National Institute on Drug Abuse Grant T32 DA017637-14). This work also used instrumentation supported by the NIH's National Center for Research Resources Shared Instrumentation Grant Program (Grant S10 RR023401).
The authors declare no competing financial interests.
- Correspondence should be addressed to Philip A. Kragel at philip.kragel{at}colorado.edu or Ajay B. Satpute at a.satpute{at}northeastern.edu
