Stress-induced aberrations in sensory processing predict worse cognitive outcomes in healthy aging adults

Neural responses to somatosensory stimulation

Statistical analysis of our MEG sensor-level oscillatory data revealed robust broadband synchronizations in sensors near the sensorimotor cortices from about 10 to 90 Hz during the first 50 ms following onset of the first stimulation, with responses to the second stimulation (500-550 ms) extending from roughly 10 to 75 Hz (p < .001, corrected; Figure 2). To evaluate the dynamics, we focused our beamformer image reconstruction analysis on the higher 30-75 Hz frequency range and two 50 ms time intervals in which the oscillatory response was the strongest following stimulation (i.e., 0-50 and 500-550 ms for stimulations 1 and 2, respectively). Of note, we limited our analysis to 30 Hz on the low end because this aligns with the traditional definition of the gamma band. In contrast, we restricted our analyses to 75 Hz on the high end as the relative power sharply decreased thereafter, especially in response to the second stimulus.

Beamformer images revealed robust increases in gamma activity in the contralateral hand region of the primary somatosensory cortex (Figure 2), with identical peak locations in response to the first and second stimulation. As described in the methods, these images were grand-averaged across all participants and both stimulations to extract virtual sensor time series from the peak voxel. We derived the relative power envelope for the 30-75 Hz band from the resulting baseline-corrected time series, and used these in subsequent multilevel models of life and biological stress on gamma activity. Importantly, paired-sample t-tests between stimulation response power revealed that the response to the second stimulation was significantly reduced compared to the first, indicative of significant gating of gamma oscillatory power in the primary somatosensory cortex across all participants (t(73) = 6.54, p < .001; Figure 2).

In contrast, evaluation of sensor-level time-domain responses revealed more temporally-extended clusters evolving shortly after each somatosensory stimulation (i.e., 25-125 ms and 525-625 ms for stimulations 1 and 2, respectively) in sensors near the contralateral sensorimotor strip (p < .001, corrected; Figure 3). sLORETA source images revealed robust increases in phase-locked, time-domain neural responses in the left primary somatosensory cortex contralateral to stimulation, with identical peak locations in response to the first and second stimulation. As described in the methods, these images were grand-averaged across all participants and both stimulations to extract time series from the peak voxel for subsequent multilevel models of life and biological stress on evoked data (Figure 3). Paired-sample t-tests of stimulation response power revealed significantly weaker time-domain responses to the second stimulus compared to the first, indicative of significant sensory gating in the time-domain across all participants (t(73) = 9.31, p < .001, Figure 3).

Oscillatory profiles of somatosensory gating are modulated by life and biological stress

To investigate how indices of life and biological stress impacted sensory gating across the lifespan, a multiple regression of empirically-derived and sample-specific definitions of allostatic load, relative age acceleration/deceleration, and their interaction was computed on the somatosensory gating ratio derived from the gamma band relative time series. Of note, higher gating ratios are indicative of worse suppression of redundant information. For descriptive statistics regarding our empirically-derived, sample-specific measures of life and biological stress, see Figure 4. Interestingly, allostatic load, relative age acceleration/deceleration, and their interaction (i.e., cumulative stress) were significantly predictive of the integrity of sensory gating, such that increases in life, biological, and cumulative stress were predictive of higher gating ratios (i.e., worse suppression; Figure 5; allostatic load: b = 0.51, p = .012, 95% CI [0.18, 0.85], relative age acceleration: b = 0.43, p = .027, 95% CI [0.11, 0.75], cumulative stress: b = -0.02, p = .030, 95% CI [-0.04, -0.01].

Finally, to evaluate the mechanism of action by which life and biological stressors impact filtering in the primary somatosensory cortex, we performed structural equation modelling with allostatic load and relative age acceleration predicting gamma response power in response to the first and second stimulation, sequentially. As expected, response power to the first stimulation was significantly predictive of gamma activity in response to the second stimulation, such that increases in gamma power during the first stimulation were predictive of increased gamma power during the second stimulus (b = .60, p < .001, 95% CI [0.55, 0.65]). In regard to our measures of life and biological stress, allostatic load was significantly predictive of gamma power during the first stimulation, but not the second, such that greater levels of life stress were predictive of reduced neural response to the first stimulation in the pair (b = -.56, p = .043, 95% CI [-1.02, -0.11]). In contrast, relative age acceleration was significantly predictive of the neural response to the second stimulus in the pair, but not the first, such that greater age acceleration indices were predictive of greater gamma power during the second stimulation (i.e., less gating-related attenuation; b = 0.11, p = .014, 95% CI [0.01, 0.05]; Figure 5).

Evoked neural response is unaffected by life and biological stress

For correspondence with the previous literature, which has focused on evoked somatosensory processing, we computed the time-domain sLORETA source images for each participant. Time series data were then extracted from the grand-averaged peak voxel of the phase-locked response estimates to derive indices of somatosensory processing (i.e., gating ratio, source power). Next, we conducted multiple regressions to evaluate whether life, biological, and cumulative stress (i.e, their interaction) predicted changes in the gating of phase-locked neural responses. Interestingly, allostatic load (b = -0.001, p = .993, 95% CI [-0.15, 0.15]), relative age acceleration (b = -0.001, p = .987, 95% CI [-0.14, 0.14]), and their interaction (b = 0.000, p = .913, 95% CI [-0.01, 0.01]) did not significantly predict the gating ratio derived from our time-domain analysis, suggesting that the evoked neural responses to somatosensory stimulation were not sensitive to changes in life and biological stress, as defined in the current study (Figure 6).

For completeness, we conducted the follow-up structural equation model interrogating the role of life and biological stress on the time-domain response to the first and second stimulation response power, sequentially. Similar to our oscillatory analysis, increased evoked response power to the first stimulation was significantly predictive of increased neural response power to the second stimulation (b = 0.62, p < .001, 95% CI [0.54, 0.69]). Interestingly, while allostatic load was not significantly predictive of phase-locked responses to the first or second stimulation (ps > .158), relative age acceleration was significantly predictive of neural response to the second stimulus, but not the first, such that greater age acceleration (i.e., greater biological age relative to chronological age) led to stronger neural responses to redundant stimuli (i.e., less attenuation; b = 12.91, p = .011, 95% CI [4.58, 21.24]; Figure 6).

Stress-induced changes in sensory gating predict poorer cognitive outcomes

Finally, given the previous literature suggesting that bottom-up, pre-attentive gating is closely linked with top-down higher order cognitive functions (e.g., attention; (15)), we evaluated whether stress-induced changes in somatosensory gating across the lifespan significantly predicted cognitive function. To assess cognitive function, we used an extensive neuropsychological battery that probed six domains and global cognitive function, including learning, memory, processing speed, executive function, attention, motor function and global cognition. As described in the methods, domain-specific composite scores were computed by averaging the Z-scores of assessments comprising each respective domain. Of note, we used unadjusted Z-scores rather than demographically-normed ones due to the fact that age and years of education were important variables comprising our measures of relative age acceleration and allostatic load, respectively. Next, we calculated the predicted gating ratios accounting for levels of allostatic load, relative age acceleration and their interaction using the regression equation described above. This yielded a predicted value of stress-induced gating, with higher values indicative of worse sensory gating, accounting for levels of life and biological stress. Interestingly, linear regressions of predicted oscillatory gating ratios on each cognitive domain and global cognitive function scores revealed that higher gating ratios (induced by increased levels of stress) were significantly predictive of poorer global cognitive function (i.e., decreased Z-scores averaged across all domains (b = -2.08, p = .005, 95% CI [-3.50, -0.66])). In addition, stress-induced declines in sensory gating were also predictive of poorer cognitive function on learning (b = -3.78, p = .001, 95% CI [-5.94, -1.62]), memory (b = -3.03, p = .003, 95% CI [-5.00, -1.07]), attention (b = -2.17, p = .028, 95% CI [-4.10, -0.24]), and processing speed-specific composites (b = -2.07, p = .025, 95% CI [-3.86, -0.27]), suggesting that less filtering of redundant sensory information (accounting for elevated levels of stress) led to significantly worse performance on these cognitive domains (Figure 7).

Participant demographics and neuropsychological assessment

Seventy-four healthy adults (M_age = 43.6 years old, range: 22-72 years old, 35 females) were enrolled in this study. Exclusionary criteria included any medical illness affecting CNS function, any neurological or psychiatric disorder, history of head trauma, current pregnancy, current substance use, implanted ferromagnetic objects or extensive dental work, and cognitive impairment. Cognitive impairment was based on a thorough neuropsychological battery that assessed functionality across six domains (i.e., learning, memory, executive function, attention, processing speed, motor function). Participants who scored less than one standard deviation from the mean using demographically-normed scores per test were deemed cognitively impaired and excluded from the study. The battery included the following tests for each domain: learning (Hopkins Verbal Learning Test – Revised (HVLT-R) Learning Trials 1-3 [74]), memory (HVLT-R Delayed Recall and Recognition Discriminability Index [74]), executive function (Comalli Stroop Test Interference Trial [75], semantic verbal fluency [76], phonemic verbal fluency [76], and Trail Making Test Part B [76]), processing speed (Comalli Stroop Test Color Trial [75], Wechsler Adult Intelligence Scale (WAIS-III) Digit Symbol Coding [77], and Trail Making Part A [76]), attention (WAIS-III Symbol Search [77], and Comalli Stroop Word Trial [75]), and motor function (Grooved Pegboard, Dominant and Non-Dominant Hands [76, 78]). Z-scores were computed using raw scores and composite domain-specific scores were calculated by averaging the Z-scores of assessments that comprised each domain respectively (see above). Global cognitive function scores were computed by averaging all domain-specific composite Z-scores. Of note, for a subsample of the participants included in the current study (N = 29), the HVLT-R Learning, Recall and Recognition tests were administered faster than the recommended speed (i.e., approximately 1 word per 1 second rather than 1 word per 2 seconds). However, all statistical analyses including these data accounted for this administration difference, and the strength of these relationships did not differ when accounting for participant exclusions using Fisher Z comparisons (learning: Z = 1.29, p = .197; memory: Z = 1.52, p = .128). The University of Nebraska Medical Center Institutional Review Board approved the study and all participants provided written informed consent.

Experimental paradigm

Participants were seated in a nonmagnetic chair with their head positioned within the MEG helmet-shaped sensor array. Electrical stimulation was applied to the right median nerve using external cutaneous stimulators connected to a Digitimer DS7A constant-current stimulator system (Digitimer Limited, Letchworth Garden City, UK). For each participant, we collected at least 80 paired-pulse trials with an inter-stimulation interval of 500 ms and an inter-pair interval that randomly varied between 4500 and 4800 ms. Note that our use of a jittered inter-pair interval of 4500-4800 ms is commonly used to avoid habituation of the neural responses via anticipation of the upcoming paired-pulses [11, 13, 36, 38]. Additionally, our use of 500 ms as the inter-stimulation interval was chosen based on prior studies that have shown that 200-500 ms is the optimal inter-stimulation interval for maximizing gating effects in this paradigm [39]. Each pulse generated a 0.2 ms constant-current square wave that was set to a limit of 10% above the motor threshold required to elicit a subtle twitch of the thumb.

MEG data acquisition and coregistration with structural MRI

All recordings were performed in a one-layer magnetically-shielded room with active shielding engaged for environmental noise compensation. With an acquisition bandwidth of 0.1-330 Hz, neuromagnetic responses were sampled continuously at 1 kHz using an MEGIN/Elekta MEG system (MEGIN, Helsinki, Finland) with 306 magnetic sensors, including 204 planar gradiometers and 102 magnetometers. Throughout data acquisition, participants were monitored using a real-time audio-video feed from inside the magnetically-shielded room. MEG data from each participant were individually corrected for head motion and subjected to noise reduction using the signal space separation method with a temporal extension [79]. Each participant’s MEG data were coregistered with their structural T1-weighted MRI data prior to imaging analyses using BESA MRI (Version 2.0). Structural MRI data were aligned parallel to the anterior and posterior commissures and transformed into standardized space. After beamformer analysis (see below), each subject’s functional images were transformed into standardized space using the transform that was previously applied to the structural MRI volume and spatially resampled.

MEG preprocessing and sensor-level statistics

Cardiac and ocular artifacts were removed from the data using signal-space projection (SSP) and the projection operator was accounted for during source reconstruction [80]. Epochs were of 3700 ms duration, with 0 ms defined as the onset of the first stimulation and the baseline being the -700 to -300 ms window. Of note, we shifted our baseline away from the period immediately preceding stimulus onset to avoid potential contamination by any anticipatory responses, although there was no evidence of such anticipatory responses in our final analyses. Epochs containing artifacts were rejected based on a fixed threshold method, supplemented with visual inspection. On average, 73 trials per participant were used for further analysis.

Artifact-free epochs were further processed using two parallel pipelines. For the oscillatory analysis, epochs were transformed into the time-frequency domain using complex demodulation [81], and the resulting spectral power estimations per sensor were averaged over trials to generate time-frequency plots of mean spectral density. The sensor-level data per time-frequency bin were normalized using the mean power per frequency during the -700 to -300 ms baseline period. The specific time-frequency windows used for imaging were determined through a two-stage, data-driven approach involving statistical analysis of the sensor-level spectrograms across all participants and trials. First, paired-sample t-tests against baseline were conducted on each data point, with the output spectrogram of t-values initially thresholded at p < .05 to define time-frequency bins containing potentially significant oscillatory deviations. To reduce the risk of false positive results due to multiple comparisons, the time-frequency bins that survived that initial threshold were temporally and/or spectrally clustered with neighboring bins that were also significant, and a cluster value was derived by summing all of the t-values of all data points in the cluster. Nonparametric permutation testing (10,000 permutations) was then used to derive a distribution of cluster values and the significance level of the observed clusters were tested directly using this distribution. Based on this analysis, the time-frequency periods that contained significant oscillatory events across all participants were subjected to beamforming analyses. Of note, in the case of the broadband gamma oscillations, we focused on a window surrounding the peak of the response (i.e., greatest amplitude change from baseline) in order to optimize the signal to noise ratio. Note that the significant time-frequency extent of the gamma response extended beyond this 75 Hz band, especially following the first stimulation.

For the time domain (i.e., evoked) analyses, all artifact-free epochs per participant were averaged with respect to the onset of the first stimulation for each sensor in the array, and the resulting mean time series per sensor and participant were examined statistically to determine the specific time windows used for subsequent source analyses. Like the oscillatory analyses, we used a two-stage approach that included paired-sample t-test against baseline, followed up with cluster-based permutation testing to control for multiple comparisons (initial threshold: p < .05, permutations: 10,000). The phase-locked, time-domain period that significantly differed from baseline were used to guide subsequent time-domain source level analysis. Further details of this method and our processing pipeline can be found in recent papers [6, 11, 36, 38, 39].

MEG source imaging

Cortical oscillatory networks were imaged through the dynamic imaging of coherent sources (DICS) beamformer [82], which uses the cross-spectral density matrices to calculate source power for the entire brain volume. These images are typically referred to as pseudo-t maps, with units (pseudo-t) that reflect noise-normalized power differences (i.e., active vs. passive) per voxel. Following convention, we computed noise-normalized, source power per voxel in each participant using baseline periods of equal duration and bandwidth [83]. MEG preprocessing and imaging used the Brain Electrical Source Analysis (Version 6.1; BESA) software. Further details of our analysis pipeline can be found in Spooner et al., (2020) [39].

Normalized source power was computed over the entire brain volume per participant at 4.0 × 4.0 × 4.0 mm resolution for the time-frequency periods identified through the sensor level analyses. Prior to statistical analysis, each participant’s MEG data, which were coregistered to native space structural MRI prior to beamforming, were transformed into standardized space using the transform previously applied to the structural MRI volume and spatially resampled. The resulting 3D maps of brain activity were averaged across all participants and both stimulations to assess the neuroanatomical basis of the significant oscillatory responses identified through the sensor-level analysis, and to allow identification of the peak voxels per oscillatory response.

Voxel time series data (i.e., “virtual sensors”) were extracted from each participant’s data individually using the peak voxel from the grand-averaged beamformer images. To compute the virtual sensors, we applied the sensor weighting matrix derived through the forward computation to the preprocessed signal vector, which yielded a time series for the specific coordinate in source space. Note that virtual sensor extraction was done per participant, once the coordinates of interest were known. Once the virtual sensor time series were extracted, we computed the envelope of the spectral power within the frequency range used in the beamforming analysis. From this time series, we computed the relative (i.e., baseline-corrected) response time series of each participant to quantify indices of somatosensory processing, including the gating ratio (response power to stim 2/response power to stim1) and source power in response to both electrical stimulations.

To enhance comparability with previous work, source images of the time-domain averaged responses were computed using standardized low-resolution brain electromagnetic tomography (sLORETA; regularization: Tikhonov 0.01%) [84]. The resulting whole-brain maps were 4-dimensional estimates of current density per voxel, per time sample across the experimental epoch. These data were normalized to the sum of the noise covariance and theoretical signal covariance, and thus the units are arbitrary. Using the temporal clusters identified in the time-domain sensor-level analysis, these maps were averaged over time following each somatosensory stimulation (e.g., 25-125 ms and 525-625 ms). The resulting maps were then grand-averaged across the two stimulations to determine the peak voxel of the time-domain neural response to the stimuli across participants. From this peak, voxel time series were extracted in sLORETA units to quantify measures of somatosensory processing (i.e., gating ratio, source power).

Blood-based markers of age acceleration

To evaluate levels of biological stress on age-related changes in sensory gating, whole blood samples were collected from each participant using BD Vacutainer EDTA tubes to evaluate methylation metrics based on the Hannum, Horvath and consensus models of predicted biological age. DNA methylation analyses were conducted on the entire data set from a large study of aging adults (N > 180) reported in previous publications [11, 13, 38, 85–87] and closely aligned with epigenetic age estimations established in previous work [35].

Specifically, DNA was purified from whole blood using DNeasy blood tissue extraction kits (QIAGEN). Methylation analysis was performed using Infinium HumanMethylation450 BeadChip Kits (Illumina). Following hybridization, BeadChips were scanned using the Illumina HiScan System. All data were processed through the Minfi R processing pipeline [88]. Methylome data were downloaded from Hannum [26] and EPIC [89] (GEO: GSE40279 and GSE51032), and these data were processed together along with methylation data generated from the larger study mentioned above. Beta values were extracted and quantile normalized using Minfi; cell counts were estimated using estimateCellComposition and resulting normalized beta values were adjusted for cell types [35, 90]. All data was then normalized using a modified BMIQ procedure provided by Horvath [25]. The gold standard was set to the median beta observed in the Hannum study [26].

For the current study, the “consensus model” of predicted biological age (i.e., both Hannum and Horvath predictions) was used, as this has been shown to outperform either prediction model in isolation [35]. Importantly, our measure of relative age acceleration (i.e., acceleration or deceleration in biological age relative to chronological age) was computed using the residuals from a regression of the consensus model of predicted biological age on chronological age for the 74 participants eligible for the current study.

Definition of allostatic load

To index relevant markers of allostatic load (i.e., life stress) in the current sample, we conducted an exploratory factor analysis (EFA) to define a latent variable of allostatic load using a compilation of metrics that are known to contribute to overall life stress and health. We began with a set list of seven measurements and progressively removed individual variables based on poor loadings (λ > .20), and overall model fit. Criteria for good model fit included a non-statistically significant chi square, a root mean squared error of approximation (RMSEA) < .06, a comparative fit index (CFI) > .95, and a standardized root mean squared residual (SRMR) < .08 based on standards in the literature [91]. The best-fitting model was used to define a latent variable for which an allostatic load score was extracted per participant. Modeling was completed using Mplus (Version 8.1).

We began by fitting a latent variable defined by alcohol use (Alcohol Use Disorders Identification Test—Composition Score (AUDIT-C)), relative age acceleration, depression (Beck Depression Index total score), number of self-reported declines in activities of daily living (ADL), body mass index (BMI), household income, and years of education. Briefly, alcohol use was defined using the AUDIT-C [92], which contains items that reflect the amount and/or frequency of current and past alcohol use during the past year. Thus, participants with higher scores are defined as heavy alcohol users with scores of four or more for men and three or more for women being considered clinically-relevant heavy alcohol use. In regard to depression symptom severity, we used the Beck Depression Inventory-II (BDI-II [93]), which requires participants to answer questions regarding the presence of depressive symptoms within a two-week period, with higher scores indicative of greater symptom severity. Regarding ADL, these scores were determined based on self-report declines in the participant’s ability to carry out daily tasks (e.g., dressing, eating, personal hygiene, etc.). Similar to our other measures, higher ADL scores reflected greater perceived decline in performance of daily tasks. Income and education were reverse-coded such that increasing values reflected lower education and income, and thus were indicative of greater stress risk. Importantly, all measures included in the EFA were treated as continuous variables of potential life stress. The resultant latent variable yielded allostatic load scores in which higher values indicated more allostatic load (i.e., greater life stress). The initial EFA based on all seven life stress factors indicated a single-factor solution with excellent fit (χ²(14) = 14.92, p = .38; RMSEA = .03, 90% CI [.00, .12]; CFI = .96; SRMR = .06). However, AUDIT-C, relative age acceleration, and income all loaded poorly onto the factor (λ’s = -.009 to .032). Excluding these three variables yielded a superior EFA solution, still yielding a single factor (λ’s = .22 to .80) with excellent overall model fit (χ²(2) = 2.03, p = .36; RMSEA = .01, 90% CI [.00, .23]; CFI = .99; SRMR = .04). Thus, our empirically-defined sample-specific quantification of allostatic load was comprised of depression symptom severity, total years of education, perceived declines in ADL, and BMI, with higher values indicative of greater current life stress (Figure 1). An allostatic load score was extracted per participant using the structure of this final EFA solution.

Statistical analyses

To evaluate the predictive capacity of life and biological markers of stress on age-related declines in somatosensory gating, we conducted a series of regressions with allostatic load (i.e., life stress), relative age acceleration (i.e., biological stress) and their interaction (i.e., cumulative stress; allostatic load score x relative age advancement) as predictors of the gating ratio derived from oscillatory and time-domain source imaging techniques, separately (Figure 1). In addition, we conducted follow-up analyses to examine the mechanism of action by which life and biological stress modulated the age-related decline in sensory gating. Specifically, we examined whether life and biological stress differentially predicted the neural response power to the first and second stimulation in our paired-pulse paradigm using structural equation modeling in Mplus (Version 8.1). Finally, we aimed to examine whether stress-induced changes to bottom-up sensory gating were predictive of higher order cognitive function assessed outside of the scanner using an extensive neuropsychological battery. Essentially, we computed a predicted gating ratio per participant accounting for levels of allostatic load, relative age acceleration and their interaction (i.e., using the regression equation described above), with higher values indicative of worse suppression of redundant sensory input. Next, we conducted a series of linear regressions of predicted gating ratios on six cognitive domains and global cognitive function (Z-scores). Of note, we used unadjusted Z-scores to compute domain-specific composites rather than demographically-normed ones due to the fact that age and years of education were important variables comprising our measures of relative age acceleration and allostatic load, respectively.