Direct Evidence for the Feedforward Neurovascular Coupling Mechanism in Humans During Task Onset: An EEG-fNIRS-TCD Multimodal Imaging Study

TL;DR

This study shows that brain activity and blood flow in small vessels predict larger artery blood flow during tasks, supporting a one-way brain-blood flow mechanism in humans.

Contribution

First successful concurrent assessment of neuronal, microvascular, and arterial responses to quantify neurovascular coupling in humans.

Findings

Neuronal activity and microvascular oxygenation predicted upstream arterial blood flow during task onset.

The neurovascular coupling response was consistent across participants, with minimal influence from concussion history or mental health.

Sex, age, and cardiorespiratory fitness influenced baseline and peak hemodynamic values.

Abstract

What are the main findings?
Concurrent EEG, fNIRS, and TCD measurements showed that neuronal activity and microvascular oxygenation predicted upstream arterial blood flow during task onset, supporting a unidirectional feedforward neurovascular coupling response.The overall NVC response was consistent across participants, with minimal influence of concussion history, mental health, or learning disabilities, although sex, age, and cardiorespiratory fitness influenced baseline and peak hemodynamic values.

Concurrent EEG, fNIRS, and TCD measurements showed that neuronal activity and microvascular oxygenation predicted upstream arterial blood flow during task onset, supporting a unidirectional feedforward neurovascular coupling response.

The overall NVC response was consistent across participants, with minimal influence of concussion history, mental health, or learning disabilities, although sex, age, and cardiorespiratory fitness influenced baseline and peak hemodynamic values.

What is the implication of the main findings?
This study represents the first successful concurrent assessment of neuronal, microvascular, and intracranial arterial responses to quantify neurovascular coupling in humans; providing support for the feedforward neurovascular coupling model in which neuronal activation and microvascular changes drive upstream cerebrovascular responses during cognitive and motor tasks.The findings demonstrate the utility of multimodal imaging (EEG–fNIRS–TCD) for characterizing integrated neurovascular responses in humans, providing a framework for future studies of cerebrovascular and neurological health.

This study represents the first successful concurrent assessment of neuronal, microvascular, and intracranial arterial responses to quantify neurovascular coupling in humans; providing support for the feedforward neurovascular coupling model in which neuronal activation and microvascular changes drive upstream cerebrovascular responses during cognitive and motor tasks.

The findings demonstrate the utility of multimodal imaging (EEG–fNIRS–TCD) for characterizing integrated neurovascular responses in humans, providing a framework for future studies of cerebrovascular and neurological health.

This investigation assessed the neurovascular coupling response through integrated assessments of neuronal function [electroencephalography (EEG)], microvascular oxygenation concentrations [functional near-infrared spectroscopy (fNIRS)], and arterial responses [transcranial Doppler ultrasound (TCD)]. The NVC response was assessed in 113 participants (86 females, aged 19–40 years) during visual (“Where’s Waldo?”) and motor (finger tapping) tasks. Block-averaged, time–frequency power was computed from the EEG data, while hemodynamic response functions were obtained from the fNIRS and TCD metrics. Granger causality assessed the predictiveness between EEG-fNIRS-TCD waveforms for each participant and was converted into a percentage of individuals displaying a significant value. Linear models were computed to determine the influence of sex, concussion history, young adulthood age, cardiorespiratory fitness, and mental health/learning disabilities on NVC parameters. During the initial 10 s of task onset, unidirectional predictiveness was weak to very strong for EEG-TCD (range: 47–83%) and fNIRS-TCD (44–92%) relationships; however, very weak to weak predictiveness was seen for the E0EG-fNIRS (0–29%) relationship for both tasks. Aside from known sex-, age-, and fitness-based influences on baseline/peak hemodynamic values (p < 0.050), the addition of concussion history and mental health/learning disabilities had minimal influence on NVC responses (p > 0.050). The findings demonstrated a unidirectional feedforward mechanism from the neuronal and microvasculature to the upstream arteries during task onset.

Figures

Full text

1. Introduction

Neurovascular coupling (NVC), also known as functional hyperemia, describes the tightly regulated relationship between neuronal activity and cerebral perfusion, which ensures that active regions of the brain parenchyma receive sufficient oxygen and nutrients to meet metabolic demands during activation [1,2,3]. This response is critical for maintaining optimal brain function and relies on feedforward and feedback mechanisms [4]. The feedforward component is characterized by rapid adjustments in blood delivery to activated regions, while feedback mechanisms maintain cerebral homeostasis by modulating vascular tone in response to metabolic byproducts [4]. Previous work from the author group has demonstrated the feasibility of multimodal electroencephalography (EEG), functional near-infrared spectroscopy (fNIRS), and transcranial Doppler ultrasound (TCD) to assess the NVC response to motor- and visual-based tasks in a small pilot study (n = 15) [5]. However, a large-scale investigation characterizing NVC using concurrent assessments of neural and hemodynamic systems in humans has not been completed due to the difficulty of simultaneous integration.

Several investigations have noted that neuronal or hemodynamic variables are influenced by modifying variables such as biological sex, age, concussion history, and cardiorespiratory fitness [6,7,8,9]. For example, while females have smaller brains, they contain more gray matter density due to greater cortical folding of the brain parenchyma [6,7]. Functional differences have additionally been noted, with females displaying larger EEG delta, theta, alpha, and beta power amplitudes [10,11]. Moreover, long-term neuronal deficits and hemodynamic alterations have been noted in those with a history of concussion; however, it is unclear how these influence the NVC response [12,13,14]. Linear declines in cerebral perfusion and neuronal density are known to occur with increasing age starting in young adulthood, with some investigations denoting differences in NVC metrics in older populations [15,16]. However, regular physical activity and exercise have been shown to have protective effects, attenuating the rate of neuronal/perfusion decline [17,18,19,20]. Finally, it is currently unknown how mental health and learning disabilities influence the NVC response.

Many investigations into these modifying variables have been conducted using unimodal neuroimaging techniques (isolated TCD, fNIRS, or magnetic resonance imaging), which only assess a singular aspect of cerebral functioning [21]. Therefore, the purpose of this study was to assess the directional relationship between neuronal function (bioelectric activity assessed via EEG), microvascular oxygenation changes (via fNIRS), and arterial changes (via TCD) during visual- and motor-based NVC in young adults. A secondary aim of this study is to assess the influence of sex, age, concussion history, and aerobic capacity on the NVC response. This comprehensive approach will enable an in-depth bidirectional assessment of the known feedforward and feedback loops in young adults. We hypothesize that multimodal assessment will provide objective concurrent evidence of feedforward NVC. Moreover, we hypothesize that NVC will be weaker in those with a history of concussion [8], older age [22], and lower aerobic capacity [23].

2. Materials and Methods

2.1. Ethics and Study Design

The University of Calgary Conjoint Health Research Ethics Review Board approved the study under protocol REB20-2112. Except for study registration within a public database, all procedures adhered to the principles outlined in the 2013 revision of the Declaration of Helsinki [24]. A convenience sample of healthy young adult volunteers was recruited using word of mouth and snowball sampling techniques. Data collection followed the guidelines and standards recommended in the Neurovascular Coupling White Paper and other recommendations [3,21,25]. A detailed description of all methods and the full study design flowchart can be found in the previously published protocol paper [26].

2.2. Participants

The participants completed a main multimodal neuroimaging visit (n = 113; 77% female; 23 ± 5 years, range: 19–40 years; body mass index: 24 ± 3 kg/m^2^) and a subsample additionally completed an optional maximal oxygen uptake test (VO_2max_) (n = 58; 77% female; 23 ± 6 years, range: 19–40 years; body mass index: 24 ± 3 kg/m^2^) (Table 1). The exercise test was completed within a median of 2 days (interquartile range: 1–3, range: 0–10) of the multimodal imaging visit. Testing for female participants occurred between days 3 and 10 of the early follicular phase, during which 11 used oral contraceptives, 33 used intrauterine devices, and 42 did not use hormonal contraceptive methods. All participants were free from chronic neurological, cerebrovascular, cardiovascular, and respiratory conditions.

Data collection occurred between 07:00 and 18:00 from March 2023 to May 2024 in a controlled environment: barometric pressure (666 ± 6 Torr), temperature (22 ± 1 °C), and humidity (22 ± 12 percent). For 12 h prior to data collection, the participants abstained from nicotine products, caffeine, cannabis, alcohol, and other illicit substances. The participants also avoided exercise for at least 6 h prior to each visit [27]. To control for ambient light impacting the fNIRS signals, testing was completed in a darkened room (lumens range: 30–100, equivalent to 2–8 candles; LUX LED Light Meter LM-50KL, LATNEX, Woodbridge, ON, Canada).

2.3. Instrumentation

An in-depth description of the multimodal imaging modality, along with a visualization of the montage, has been published elsewhere [26]. EEG electrodes and fNIRS optodes were held in place using a custom-designed headcap, with Cz positioned over the vertex of the head (i.e., the midpoint between preauricular points and the naison/inion) (Figure 1). EEG data were collected at 1000 Hz using a 16-channel Cyton Biosensing Board (OpenBCI, Brooklyn, NY, USA) and sent to an OpenBCI GUI. Ground and reference electrodes were positioned over the left and right mastoid processes, respectively [28]. A continuous-wave NIRS device (cwNIRS) was used to sample fNIRS data at 3.91 Hz from 16 sources and 16 detectors (NIRScout; NIRx Medical Technologies, Berlin, Germany), which created 40 channels (20 per hemisphere) and were collected using the NIRStar GUI (Figure 1) [29,30]. Two 2 MHz TCD ultrasound probes were held in place with an adjustable headframe (DWL Technologies, San Juan Capistrano, CA, USA) and positioned over the transtemporal windows (Figure 1) [29]. These probes insonated the middle cerebral and posterior cerebral arteries, which were confirmed using carotid compressions and simple visual tasks [31]. Beat-to-beat blood pressure was captured using a Finometer NOVA device (Finapres Medical Systems, Amsterdam, The Netherlands) and corrected using a height calibration unit (Finapres Medical Systems, Amsterdam, The Netherlands) [32,33]. Each PQRST waveform was captured using a 3-lead electrocardiography system, using lead II methodology (FE 201; AD Instruments, Colorado Springs, CO, USA). Breath-by-breath respiratory parameters were captured using an inline gas analyzer (ML206; AD Instruments, Colorado Springs, CO, USA). Physiological and TCD data were time-synchronized using the LabChart GUI (AD Instruments, Colorado Springs, CO, USA), which linearly interpolated the TCD (100 Hz) and Finapres NOVA (200 Hz) devices to 1000 Hz.

2.4. Experimental Protocols

This investigation was a subsection of a larger study additionally using the multimodal approach and the aforementioned equipment to assess dynamic cerebral autoregulation (squat-stand maneuvers), cerebrovascular reactivity (dynamic end-tidal forcing), and cardiovascular function (quiet sitting and standing) [26]. The current investigation included a motor finger-tapping task and a visual complex scene search (i.e., “Where’s Waldo/Wally?”). Pilot work has demonstrated the feasibility of the EEG-fNIRS-TCD approach to capture the NVC response using these tasks [5]. For the finger-tapping task, the participants completed 12 blocks consisting of 15 s of rest followed by 10 s of tapping their thumbs to their four fingers on both hands at a rate of 1 Hz (32–34). For the Waldo task, the participants completed 10 blocks of ~20 s with eyes closed and 40 s of active engagement with eyes open during the task [34,35]. Prior to completion, the participants were instructed to search for the five characters within the Waldo Universe (i.e., Waldo, Odlaw, Wenda, Wizard Whitebeard, and Woof’s tail), where a new puzzle was presented for each cycle (i.e., a total of 10) [36]. These tasks were employed to maximize engagement.

The fNIRS task instructions were automated to present the rest and finger tapping sections via PsychoPy [37]. With the Waldo task, the researchers manually input the triggers via PsychoPy, as the participants had their eyes closed and thus could not see visual cues. Additionally, this enabled the researchers to extend the resting portion in the presence of Mayer waves [38]. The PsychoPy platform also allowed for precise time synchronization across the three streaming platforms (OpenBCI GUI, NIRStar, and TCD) using a C-pod and a parallel port replicator. The methods paper for this project explicitly details the time synchronization process via PsychoPy, which allowed for the combination of datasets during offline analysis [26].

The participants completed a graded ramp exercise test on a cycle ergometer to quantify VO_2max_ [39,40]. The protocol included a 5-min warm-up at a fixed workload (40 W for females and 50 W for males), followed by a ramp protocol aimed to be 8–12 min in length [39,40]. For females, the workload increased by 15 W per minute; for males, it increased by 25 W per minute until volitional exhaustion. Standardized verbal encouragement was provided to ensure maximal effort.

2.5. Data Processing

The EEGLAB toolbox in MATLAB (version 2024.0) was used to filter, clean, and remove artifacts from the EEG data [41]. Specifically, a bandpass filter of 0.5 to 40 Hz was applied to reduce noise, particularly noise arising from nearby power lines (~60 Hz) and other neuroimaging and biomedical equipment. Artifact Subspace Reconstruction and Independent Component Analysis were subsequently used to identify, isolate, and correct artifacts arising from eye blinks, electrode movement, and muscle activity [41,42]. The Homer3 (version 1.80.2) software package was used to convert, clean, and remove motion artifacts from the fNIRS data [43]. The raw light voltages collected from the NIRScout were converted into optical densities and then transformed into chromophore concentrations (i.e., HbO, HbR, and THb) via the application of the modified Beer–Lambert law [43]. Signal quality was assessed using the Pollonini method, which evaluates the presence of cardiac pulsatility within each channel as an indicator of physiological signal integrity [44]. Compared to the scalp coupling index, which quantifies optode-scalp mechanical coupling based on wavelength correlation, the Pollonini approach verifies the presence of physiologically meaningful hemodynamic pulse oscillations [44]. Given the study’s focus on NVC and temporal dynamics, prioritizing physiological signal detection ensured that channels containing valid vascular information were retained for analysis. Using a cardiac pulsation cutoff of 0.6 Hz and a 2-s window, channels below this threshold were retained if a clear cardiac signal was present, as verified through frequency spectra and wavelet transform plots. This maximized the number of channels available for HRF analysis while maintaining data quality. Wavelet motion correction was then applied to identify and remove artifacts arising from head and/or optode movement [45]. Mean middle cerebral artery velocity (MCAv), posterior cerebral artery velocity (PCAv), and arterial pressure were determined by averaging data points from each individual pulsatile waveform [46]. From these traces, systolic values were determined from the peak of the pulsatile waveforms, while diastolic values were taken from the end-diastolic point of the waveform (i.e., the lowest value before the next upstroke). Partial pressure end-tidal carbon dioxisde (P_ET_CO_2_) was calculated as the peak expired value recorded during each breath. Heart rate was derived from the R-R intervals of the ECG, converted from milliseconds to beats per minute.

Oxygen uptake was measured breath by breath using a metabolic cart that analyzed the expired gas from the sampling chamber. The system continuously recorded minute ventilation, oxygen consumption, and carbon dioxide production. These variables were calculated based on the measured volumes of expired air and the concentrations of oxygen and carbon dioxide in the expired gas. The average respiratory exchange ratio (i.e., carbon dioxide expulsion to oxygen uptake) at the end of the exertion test was >1.10. Prior to each test, the metabolic cart was calibrated to ensure accuracy in flow rates and gas analysis. Oxygen uptake data were averaged over 15-s intervals to determine VO_2max_ throughout the exercise protocol [47,48].

2.6. Independent and Outcome Variables

The independent variables in the current investigation were sex, concussion history, relative VO_2_, age, and mental health and learning disability diagnoses. Diagnosed concussion history and the presence of a mental health or learning disability diagnosis were ascertained through reporting on the Athlete Background and Symptom Evaluation portions of the Sport Concussion Assessment Tool 6 [49,50,51]. Block averages were produced across all EEG, fNIRS, and TCD data, time-synchronized to task onset and averaged across all task trials. For finger tapping, neurodynamic and hemodynamic responses were produced 10 s before and after task onset, while for Waldo these responses were produced 10 s before and 30 s after task onset. The outcome variables were the NVC metrics obtained from the EEG-fNIRS-TCD multimodal approach, including baseline, peak, relative increases/decreases, absolute increases/decreases, area under the curve during task engagement (AUC).

Block averages were produced across all EEG, fNIRS, and TCD data, time-synchronized to task onset and averaged across all task trials. The block averages for EEG data were computed within the alpha (8–12 Hz) and beta 1 (12–16 Hz) frequency bands [52,53], as pilot work demonstrated that these frequency bands had the largest changes associated with these tasks [5]. Moreover, given the regions of interest, these were averaged across central electrodes (C3, Cz, and C4) for the finger-tapping task and occipital electrodes (PO3, PO4, O1, Oz, and O2). Concentrating on two frequency bands and combining data across electrodes within neural regions reduced the number of comparisons and maximized statistical power.

For the fNIRS data, hemodynamic response functions were calculated for all chromophores (HbO, HbR, and THb) in the channels that met the signal quality criteria. Previous research has shown that averaging fNIRS metrics across a region can reduce the sensitivity of detecting subtle physiological differences between groups or populations [54]. To address this, a single channel within the five regions was chosen for further analysis, and identified as the channel with the largest HbO AUC for each task. This maximized the signal-to-noise ratio, preserving the ability to compare physiological responses across EEG, fNIRS, and TCD data. Similar to the neural responses, to maximize statistical power, only the motor fNIRS and MCAv data were included in the finger tapping analysis and the occipital fNIRS and PCAv were included in the Waldo analysis. Pilot work has shown these regions are closely associated with the respective neural and hemodynamic demands of these tasks.

To enable comparison across modalities with inherently different sampling frequencies, all signals were time-locked to task onset and analyzed within standardized temporal windows rather than at their native sampling rates. Block-averaged responses were computed for each modality, and temporal metrics (baseline, peak, and AUC) were derived from these synchronized time series. By analyzing task-evoked changes within predefined windows rather than millisecond-level fluctuations, cross-device comparisons were performed on physiologically meaningful response epochs that reflect the slower vascular dynamics relative to bioelectrical activity. For EEG, fNIRS, and TCD data, baseline values were calculated as the average signal over the 10 s preceding task onset. Peak values were determined as the maximum signal within 10 s after task onset for the finger-tapping task and 30 s for the Waldo task. For EEG and TCD, both relative (percent change from baseline to peak) and absolute (difference between baseline and peak) increases or decreases were calculated, as these metrics are in absolute units. For fNIRS, only absolute changes were used, as the recorded units are relative concentrations. The AUC was computed as the cumulative area of the signal above baseline during the task engagement period—10 s for finger tapping and 30 s for the Waldo task [54]. EEG signals were analyzed within physiologically relevant frequency bands, including delta (0.5–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), low beta (12–16 Hz), medium beta (16–20 Hz), high beta (20–30 Hz), and gamma (30–40 Hz). The EEG data were processed to align with the fNIRS timing by computing time–frequency power band analysis and averaging across 256 data points. The latter was also completed for the cerebral blood velocity (CBv) traces [46], ensuring that the temporal resolution of EEG and TCD data was comparable to fNIRS data while preserving critical timing information.

2.7. Sample Size Calculation

An a priori power analysis was conducted for a linear model with five predictor variables: sex, concussion history, age, relative VO_2max_, and mental health/learning disability composite [55]. Assuming a moderate effect size (f^2^ = 0.15), an alpha level of 0.05, and 95% power, the analysis indicated that a minimum of 89 participants would be required to detect significant effects. A second power analysis was conducted to determine that 98 participants would be required to detect a difference with a small effect size (f^2^ = 0.02), an alpha level of 0.05, and a desired statistical power of 80% [56]. Two power analyses were completed; for physiological findings, this ensures that both robust, easily detectable moderate effects and more subtle small effects are adequately captured. This is crucial when examining complex biological relationships, as even small physiological changes can have significant clinical implications or indicate early stages of adaptation or pathology [57,58].

2.8. Statistical Analysis

All statistical analyses were completed in R (version 4.4.1) [59]. Sex-based differences in baseline demographics and cardiorespiratory variables were assessed using chi-square tests for categorical data and unpaired t-tests with Cohen’s d effect sizes for continuous data. Cohen’s d thresholds included negligible (<0.2), small (0.2–0.5), moderate (0.5–0.8), and large (>0.8) [60,61].

For participants who did not complete the exercise test, missing VO_2max_ data were imputed using multiple imputation by chained equations with 10 permutations and 25 iterations per imputation, incorporating available data from EEG, fNIRS, and TCD, as well as sex, age, and concussion history [56]. Linear models were used to analyze baseline, peak, area under the curve (AUC), and relative and absolute increases. The independent variables included sex (nominal; ref: Female), concussion history (nominal; ref: Conc Hx-), relative VO_2_ (continuous), age (continuous), and the self-reported presence of mental health diagnoses and/or learning disabilities (nominal; ref: none). Data are displayed as means with 95% confidence intervals.

Temporal responses between imaging modalities were computed using Granger causality analysis [62,63]. Granger causality was employed to determine the directional interactions between neuronal signals (EEG time–frequency bands) and vascular responses (HbO, HbR, MCAv, and PCAv) [62,63]. This method determines whether one signal provides statistically significant predictive information about the future values of another signal, beyond what is explained by the latter’s own past values. This method is similar to, albeit less computationally intensive than, a transfer entropy approach, which quantifies directed information flow using information-theoretic measures and typically requires estimation of probability distributions from the data [64]. Vector autoregressive (VAR) models of order p were constructed for pairs of signals of interest (EEG-HbO, EEG-HbR, EEG-MCAv, and EEG-PCAv).

To test whether EEG activity Granger-caused vascular responses, the restricted model was defined as

$$[eqn]$$

and the unrestricted model as

$$[eqn]$$

where $[eqn]$ denotes the vascular signal at time $[eqn]$ and $[eqn]$ denotes the EEG band power, $[eqn]$ is the autoregressive coefficients, $[eqn]$ is the residual error term, and $[eqn]$ quantifies the influence of EEG activity on the vascular response.

Granger causality from EEG to the vascular signal was tested under the following condition:

$$[eqn]$$

Reciprocal models were constructed to test whether vascular signals Granger-caused EEG activity:

$$[eqn]$$ $$[eqn]$$

with the null hypothesis

$$[eqn]$$

Causality strength in both directions was assessed using an F-statistic comparing residual sum of squares between restricted and unrestricted models:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ denote the residual sum of squares of the restricted and unrestricted models, respectively, while $[eqn]$ represents the number of observations. A significantly lower prediction error in the unrestricted model indicated a directional predictive influence between modalities.

The Granger causality analysis was conducted for each individual’s block-averaged EEG, fNIRS, and TCD responses during the motor and visual tasks; however, these were divided into time bins to determine the predictive value of these signals during different phases of the tasks. For the Waldo task, time bins of −10–0, 0–10, 10–20, and 20–30 s were used based on stimulus onset with eyes open, while for finger tapping, time bins of −10–5, −5–0, 0–5, and 5–10 s were used based on the tapping onset. This was then converted into a percentage of participants who displayed a significant directional relationship within each time grouping. These were binned into quintiles of very weak (0–20%), weak (20–40%), moderate (40–60%), strong (60–80%), and very strong (80–100%) relationships. Alpha was set a priori at 0.05.

3. Results

3.1. Waldo Granger Causality Analysis

Figure 2 displays the EEG-fNIRS and fNIRS-EEG temporal relationships during Waldo. During the 10 s preceding task onset, alpha displayed a moderate relationship with HbO and ThB, while the remaining variables displayed very weak to weak relationships (Figure 2). During time bins 0–10 and 10–20 s, very weak to weak relationships were identified (Figure 2). During the 20–30-s interval, EEG-fNIRS had very weak relationships, while fNIRS-EEG had strong relationships (Figure 2).

Figure 3 displays the EEG-TCD and TCD-EEG temporal relationships during Waldo. During the 10 s preceding task onset and 20–30 s after task onset, both displayed strong to very strong bidirectional relationships (Figure 3). During the 0–10 s after task onset, EEG-TCD displayed strong to very strong relationships, while the TCD-EEG relationships were very weak (Figure 3). The bidirectional relationships were very weak to weak during the 10–20 s after task onset (Figure 3).

Figure 4 displays the fNIRS-TCD and TCD-fNIRS temporal relationships during Waldo. The bidirectional relationships were strong to very strong during the 10 s preceding task onset (Figure 4). During 0–10 and 20–30 s after task onset, fNIRS-TCD displayed strong to very strong relationships, while the reverse TCD-fNIRS relationships during these times were very weak to weak (Figure 4). During the 10–20 s after task onset, TCD-fNIRS relationships displayed moderate relationships for HbR, while HbO, THb TCD-fNIRS and all fNIRS-TCD relationships were weak (Figure 4).

3.2. Finger Tapping Granger Causality Analysis

Figure 5 displays the EEG-fNIRS and fNIRS-EEG temporal relationships during finger tapping. Strong to very strong relationships were observed for fNIRS-EEG for all chromophores during the 15–20 s after task onset (5–10 s after tapping offset) (Figure 5) and THb in the preceding 5 s before task onset (Figure 5). However, all other relationships were weak to very weak (Figure 5).

Figure 6 displays the EEG-TCD and TCD-EEG temporal relationships during finger tapping. During the 10–5 s preceding task onset, both bidirectional relationships were strong to very strong, while during the 5 s preceding task onset, both bidirectional relationships were very weak to weak (Figure 6). During finger tapping (0–10 s), the EEG-TCD relationships were weak to strong, while the TCD-EEG relationships were very weak (Figure 6).

Figure 7 displays the fNIRS-TCD and TCD-fNIRS temporal relationships during finger tapping. During the 10–5 s preceding task onset, the fNIRS-TCD relationships were very strong, while the TCD-fNIRS relationships were very weak (Figure 7). During the 5 s preceding task onset, fNIRS-TCD THb displayed very strong, HbO displayed strong, and HbR displayed moderate relationships with MCAv (Figure 7); however, the TCD-fNIRS relationships were very weak for THb, moderate to strong for HbO, and moderate for HbR (Figure 7). During the task, the fNIRS-TCD relationships were strong to very strong, while the TCD-fNIRS relationships were very weak (Figure 7).

3.3. Transcranial Doppler Ultrasound Hemodynamic Response

For the Waldo task, no differences were noted in the PCAv hemodynamic responses analysis for sex (all p > 0.078), concussion history (all p > 0.169), age (all p > 0.054), relative VO_2max_ (p > 0.320), or the mental health/learning disability composite (all p > 0.067) (Table 2). For finger tapping, compared to females, males had lower baseline and peak MCAv metrics across the cardiac cycle (all p < 0.001) (Table 2). No differences in NVC metrics were noted between those with and without a history of concussion (all p > 0.095) (Table 2). Age was negatively associated with baseline and peak MCAv across the cardiac cycle (all p < 0.001), in addition to a larger mean and systolic relative increase and systolic absolute increase in MCAv (all p < 0.004) (Table 2). A greater relative VO_2max_ was associated with higher baseline and peak MCAv across the cardiac cycle (all p < 0.049), excluding peak mean MCAv (p = 0.079) (Table 2). No associations were found between the mental health and learning disability composite and any MCAv metric (Table 2).

3.4. Functional Near-Infrared Spectroscopy Hemodynamic Response

For the Waldo task, no differences were noted across the three chromophore hemodynamic response analyses for sex (all p > 0.084), concussion history (all p > 0.240), age (all p > 0.364), relative VO_2max_ (p > 0.115), or the mental health/learning disability composite (all p > 0.321) (Table 3). For finger tapping, males had a greater baseline HbR (p = 0.048) and ThB (p = 0.021) concentrations; however, no other sex differences were noted (p > 0.070) (Table 3). No other associations were found for concussion history (all p > 0.220), age (all p > 0.196), relative VO_2max_ (p > 0.121), or the mental health/learning disability composite (all p > 0.283) (Table 3).

3.5. Electroencephalography Neurodynamic Response

For the Waldo task, no differences were noted across both alpha and low beta power time–frequency responses for sex (all p > 0.237), concussion history (all p > 0.287), age (all p > 0.191), relative VO_2max_ (p > 0.146), or the mental health/learning disability composite (all p > 0.494) (Table 4). For finger tapping, age was associated with a greater low beta AUC (p = 0.017); however, no associations were found for sex (p > 0.237), concussion history (all p > 0.273), age (all p > 0.072), relative VO_2max_ (p > 0.190), or the mental health/learning disability composite (all p > 0.384) (Table 4).

4. Discussion

The purpose of the present study was to assess the NVC response concurrently at the neural, microvascular, and intracranial cerebral artery levels in a large cohort of young adults. While the NVC response describes the coupling of blood flow to neurovascular metabolism, the three assessment techniques employed in this investigation have not historically been measured simultaneously. The study revealed many novel finCdings. Firstly, the Granger causality analysis demonstrated that neuronal and microvascular activities were unidirectionally predictive of arterial changes during the initial 10 s of both visual- and motor-based tasks. Moreover, when a steady state was achieved, a bidirectional predictive relationship was identified, suggesting that different regions of the brain communicated to maintain hemostasis around the new set-point. Furthermore, the EEG-fNIRS Granger causality analyses demonstrated minimal relationships between variables, potentially due to the insufficient sampling of the fNIRS (3.90 Hz; 256 ms) to detect the immediate neuronal–microvascular signaling. Finally, the NVC response metrics, characterized through the absolute increase, relative increase, and AUC, were minimally influenced by sex, concussion history, cardiorespiratory fitness levels, age during the young adulthood period, and mental health disorders/learning disabilities. Collectively, these results demonstrate a regional signaling propagation from the capillaries up the arterial tree to influence blood delivery; however, the introduction of higher temporal resolution fNIRS devices would enable a more intricate neuronal–microvascular–arterial quantification.

To date, the relative contributions of the known mechanisms that moderate the NVC response are debated [1,2,4]. The NVC response is typically described as an intracranial arterial response that lags behind an increase in neural metabolism, because of both an increase in energy substrate consumption and increased accumulation of metabolic waste [1,2,4]. The periods of transition from neurological rest to activation are of particular interest for investigating the mechanistic underpinnings of the known NVC response. The present analysis demonstrated that the active neurological response to the visual “Where’s Waldo?” and motor finger-tapping paradigms, quantified using EEG, was predictive of the vascular response assessed by TCD (Figure 2 and Figure 5). This finding provides objective evidence to support the proposed theory of neural activation initiating large-scale compensatory changes in regional CBF. While this association was expected, the unique capacity of the present study is the evaluation of the influence of the microvasculature, which represents the missing causal link in nutrient delivery from the microvasculature to the neuron, as well as between the microvasculature and the upstream arteries. The fNIRS-TCD relationship demonstrated similar predictive capabilities, albeit stronger during the visual task (Figure 3 and Figure 6).

There is evidence for a number of compensatory and overlapping mechanisms that influence the dilation of capillaries and arterioles. For example, a meta-analysis by Hosford and Gourine [65] on the mediators of the NVC response identified that neural nitric oxide (NO) synthase reduced the NVC response by 64%, while other neurochemicals explained some of the response. The authors concluded that one-third of the signaling was unaccounted for [65]. Furthermore, NO has been proposed to potentially cause dilation of the capillaries via pericyte activation and adenosine accumulation, arterioles via cyclic guanosine monophosphate (cGMP) and soluble guanylate cyclase pathways, and arteries via vascular smooth muscle cells [66,67,68]. In neurons, NO is produced via neural NO synthase, activated by neural depolarization, which then acts on pericytes to mediate capillary dilation [68]. NO is also produced by endothelial cells which then act on the vascular smooth muscle cells ‘underneath’ them via cGMP. Hence, NO is a potent agent for the NVC response at all levels of the vasculature, though via different mechanisms. Finally, the release of potassium during neuronal activity hyperpolarizes endothelial cells and triggers retrograde vasodilation through endothelial inward rectifier potassium channels [69,70]. This upstream propagation via ionic currents through gap junctions between endothelial cells and myoendothelial junctions with smooth muscle cells ensures coordinated dilation of upstream vessels [71,72]. Hence, the NVC response is driven by a complex interplay between both neuronal and vascular signaling molecules.

The result of the Granger analyses indicated that, as previously suggested, the microvasculature responds to local changes in neural metabolism, driving upstream regulation of the intracranial arteries to maintain homeostatic delivery of nutrients to the neurons [1,2,4]. There are three potential explanations/hypotheses for the coordinated response in the microvasculature and intracranial arteries. The first suggests that the neural control of endothelial cells in the intracranial arteries is activated during periods of microvascular dilation, causing a compensatory increase in flow to areas of increased microvascular flow [73]. Conversely, microvascular dilation in response to neural activity may result in a regional decrease in microvascular pressure, increasing the perfusion pressure in the upstream intracranial arteries, which would increase CBv [1,2,4]. Lastly, direct innervation of the vascular endothelium may co-occur with direct innervation of the microvasculature, but with different timing, rendering the microvascular response predictive of the TCD response, despite the absence of a direct causal mechanism. For example, larger cerebral arteries receive sympathetic innervation from the superior cervical ganglion, while cortical microvessels are innervated by noradrenergic, serotonergic, cholinergic, and GABAergic neurons from subcortical and local cortical sources. This suggests that both vascular compartments are directly influenced by neuronal activity through distinct pathways and temporal patterns [25,74].

To the authors’ surprise, the Granger causality analysis observed very little predictive association between the EEG-derived frequency band response and the microvascular responses to both tasks (Figure 1 and Figure 4). Theoretically, one would expect this to be similar to the EEG-TCD comparison but temporally quicker. Hence, the lack of findings could be in part explained by the fNIRS device sampling every 256 ms, missing the near-instantaneous signaling of the EEG-fNIRS relationship. Evidence for this proposition is the similar pattern of predictiveness for the EEG-TCD and fNIRS-TCD (Figure 2, Figure 3, Figure 5 and Figure 6), with arterial changes occurring ~1–2 s post-stimulus onset [5]. A sampling rate of 3.90 Hz would be sufficient to capture the microvasculature-to-arterial delay (Figure 3 and Figure 6); however, the rapid coupling of the microvasculature-to-neurological activity may have been missed due to the lower temporal resolution of the microvascular assessment (fNIRS). Finally, the frequency band assessments may be too non-specific to capture the localized response to the NVC task by the neural areas driving the NVC response [5]. Interestingly, delta band activity demonstrated stronger Granger relationships with the vascular signals compared with other EEG frequency bands. This may reflect the closer alignment between slower cortical oscillations and hemodynamic fluctuations, as low-frequency neural activity has previously been linked to vascular and metabolic processes underlying NVC [75].

After ~15–25 s with the NVC response, a new homeostatic set point is achieved (75). The feedback loop then becomes more prominent, given that the regulation requires smaller adjustments to maintain this new steady state [4]. The Granger causality relationship between variables demonstrated a bidirectional predictive relationship showing the crosstalk between signals to maintain homeostasis (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6). This relationship was more apparent in the Waldo task compared to the finger-tapping task, potentially due to the visual task eliciting greater neuronal/hemodynamic differences between the rest and active sections of each task (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6). This could also be attributable to regional differences in the NVC response, which warrants further consideration. Finally, the fact that bidirectional relationships were present between EEG-fNIRS signals during the steady-state portions of the tasks provides further support for the notion that the sampling rate was insufficient during the initial task activation.

There has been a recent emphasis on improving the generalizability and understanding of sex differences in physiological research by including male and female participants [76,77,78]. There are known sex differences in intracranial arterial CBv [22,79], cerebral hemoglobin concentration [80], and neurological function at rest [10,11], as well as in evoked potentials [81]. The present study detected a sex difference in baseline MCAv measures and in the peak MCAv response to the finger-tapping stimulus. The difference in the MCA NVC response is likely a byproduct of the elevated resting CBv detected in the female participants compared to the males, rather than a purely mechanistic difference in how the NVC response occurs between biological sexes (Table 2). Baseline sex differences were also noted for fNIRS, with males having higher THb and HbR concentrations compared to females during the finger-tapping task (Table 3). While this may point to hemodynamic sex differences within the microvessels, it could also be the result of a type I error. Nonetheless, it is established that males have a greater hematocrit concentration, leading to greater absolute concentrations of HbO, HbR, and THb [82]. As this investigation utilized a continuous-wave device, only relative concentrations were discernible. Hence, the quantification of relative concentrations would dampen the known absolute concentration differences. Finally, the current investigation highlighted a lack of sex-based differences in the alpha and low beta power-band frequencies. This could potentially be attributable to the task states washing out potential resting-state differences. Moreover, the only investigations to date that have utilized the “Where’s Waldo?” paradigm with EEG have been conducted by the current research group.

Concussion is characterized in part by a traumatic force impacting brain function yet often shows a negative result on a standard neuroimaging evaluation, such as computed tomography or magnetic resonance scan [83,84]. Nonetheless, acute functional alterations have been noted in the immediate days following concussion, which has prompted mainstream concern regarding the long-term neurological consequences of concussion [9,85,86]. This has further been amplified given the perceived potential associations with Alzheimer’s disease and related dementias [87], mild cognitive decline with age [88], depression and suicidality [89], and chronic traumatic encephalopathy [90]; however, the direct associations have commonly been overstated [91,92,93,94,95,96,97]. Previous research focused on the consequences of concussion includes the use of EEG [81,98,99,100,101], fNIRS [102,103,104,105] TCD [9,14,106,107,108], and fMRI [109,110,111,112,113,114] to assess impaired functional differences following concussion, with mixed findings. The current investigation noted no differences in those with and without a history of concussion; however, it is important to note that the present investigation focused on the NVC response. Therefore, future work into resting-state, cerebral autoregulation, and cerebrovascular reactivity domains is warranted to fully elucidate the long-term ramifications of concussion on cerebrovascular regulation.

While cerebral perfusion declines with age, which is associated with declining cognition, high physical fitness slows the decline in cerebral perfusion by roughly 10 years [17] and attenuates the odds of age-related cognitive decline by 41% [115]. The findings in Table 2 highlight that age and cardiorespiratory fitness had the largest influence on the intracranial arteries, with baseline and peak MCAv declining at a greater rate with age across young adulthood. However, fitness was associated with a higher MCAv, illustrating and confirming prior research showing that exercise in young adulthood is beneficial for augmenting cerebral perfusion. Despite these findings, the NVC response quantified through the absolute/relative increase and AUC was only nominally influenced. While fitness improves overall CBv indexes, the NVC response appears to be intact across the young adult spectrum. Previous work by Koep et al. [22] highlighted an age-by-sex interaction with a positive increasing relationship between PCAv and age in females only. This previous investigation was completed across all adulthood stages (18–85 years), while the current investigation was restricted to young adulthood (18–40 years). Hence, it appears that cardiorespiratory fitness and age may not need to be included as confounding factors within statistical models of NVC in young adulthood; however, if researchers are examining absolute CBv metrics, fitness is an important variable to include in statistical models when examining MCAv/PCAv metrics.

Finally, the linear model analyses found minimal influence of the mental health disorders/learning disability composite on the NVC response across all three imaging domains. A limitation of the included composite was that it grouped all participants together with a physician-diagnosed disorder/disability across several domains (e.g., major depressive disorder, anxiety, attention deficit hyperactivity disorder) as the study was underpowered to investigate each domain independently. Nonetheless, the collective use of these variables has been under-examined to date, and further work using scoring on continuous-based questionnaires (e.g., Patient Health Questionnaire-9, Generalized Anxiety Disorder-7) may be more sensitive than nominal groupings.

The primary limitation of the present investigation is the inherent difference in temporal resolution across the multimodal neuroimaging techniques employed. EEG data were acquired at millisecond-scale temporal resolution, reflecting rapid neuronal dynamics, whereas TCD and fNIRS capture hemodynamic responses at 100 and 3.90 Hz, respectively. These device differences in temporal properties required temporal alignment procedures, including averaging, resampling, and interpolation, to enable cross-modal comparison and Granger causality analyses. While this limits direct comparability of absolute timing between neuronal and vascular responses and may potentially lead to biased timing estimates, there is currently no commercially available platform that integrates EEG, fNIRS, and TCD within a unified acquisition system for clinical or research use. Therefore, the present investigation relied on independent acquisition systems that required post hoc synchronization and temporal alignment across modalities. Although this introduces inherent methodological challenges related to sampling frequency and timing precision, this is the first large-scale investigation to simultaneously integrate EEG, fNIRS, and TCD to examine NVC dynamics in humans, and the methodological trade-offs associated with multimodal synchronization represent a necessary step toward comprehensive characterization of neuronal–vascular interactions in humans.

Another limitation was the number of participants who did not complete the VO_2max_ testing. However, to reduce the influence of bias, missing data were predicted using multiple imputation by chained equations models that used all of the NVC and demographic data. This approach has been demonstrated to be highly robust for accounting for missing data [56,116,117]. There is some evidence that hormonal contraceptives may influence cerebrovascular function, which was not considered in the present study. However, to minimize this influence, females were tested during days 3–10 of their menstrual cycle, when hormones are known to be at their lowest levels [118,119].

To combine the assessment modalities, the spatial resolution of fNIRS and EEG was limited by the locations on the multimodal cap (i.e., centralized to the midline of the head). Thus, regional differences in the NVC response were not compared in the present study; however, pilot work demonstrated the largest responses within the occipital and motor regions for the visual- and motor-based tasks, respectively. Furthermore, while the study aimed to quantify the coupling of neurological activity with cerebral blood flow, intracranial arterial blood flow could not be directly assessed. CBv was a surrogate measure for flow, given that the diameter of the intracranial vessels could not be assessed. Thus, the influences of the endothelial NVC response were not assessed.

5. Conclusions

This study sought to isolate the feedforward and feedback loops of the NVC response in humans using a concurrent multimodal imaging approach of EEG-fNIRS-TCD. The Granger causality findings demonstrated a clear unidirectional feedforward mechanism between the EEG-TCD and fNIRS-TCD responses; however, no associations were found during task onset for the EEG-fNIRS relationship. This could be due to the lower temporal resolution of fNIRS (3.90 Hz), which was not sufficient to detect the neuronal-to-microvascular signaling that could occur within 256 ms of task onset. Future work utilizing fNIRS devices with greater temporal resolution and/or different analytical methods may be able to tease apart the precise underpinnings between local neuronal-vascular regions. Finally, aside from known arterial influences of biological sex, age, and cardiorespiratory fitness on CBv, these factors, in addition to concussion history and mental health/learning disabilities, had minimal influence on the NVC response across all three imaging domains. This increases the generalizability of previous NVC research findings, suggesting that future investigations may not need to consider these as confounding influences to control within the study design and/or statistical analyses.