We performed a grey-box linear systems analysis of the fNIRS-based CVR measure (change in tHb) to anodal HD-tDCS in healthy humans. The hemodynamic response is crucial for dosing tDCS due to its role in action53, which can lead to inter- and intra- subject variability in neuronal responses40,54. Furthermore, tDCS can be a promising method to evoke regional CBF5 to ameliorate hypoperfusion in cerebrovascular diseases, including facilitating cognitive rehabilitation. The hemodynamic response to tDCS current density in the brain can be captured based on the tHb changes using our parameterized grey-box linear model. In the hemodynamic response, low-frequency oscillations can have frequencies between 0.01 and 0.2 Hz, where the top of that range is related to the conventional hemodynamic response function of 5 seconds – TF1(s) in Fig. 5. Lower frequency oscillations < 0.1Hz are related to vasomotor and myogenic activity in the terminal arteriole and micro-vessels27. The current study investigated very low-frequency oscillations between 0.01 and 0.05 Hz that can originate from arterioles under neurogenic innervation27,26. Our physiologically detailed model captured the related non-linear calcium dynamics that is of significance since the effect of tDCS on coupling between the neuronal and hemodynamic low frequency oscillations has been shown in our prior work42 that may represent tDCS effects on the cerebral astrocytes. Guhathakurta and Dutta25 have postulated that tDCS electric field spread in the highly conductive CSF directly affects the pial arteries and arterioles that contain perivascular nerves within their adventitial layer. Since pial arteries start the pressure-driven blood pathway to the cortex (reviewed in Schmid et al.13) so direct electric field effects on the volume of pial arteries can lead to an initial dip in the blood volume (and tHb) in the capillary bed that can be seen in few healthy subjects in Fig. 1A. Here, tDCS effects on the vascular neural network6,7,11 under “neurogenic hypothesis” is in contrast to the “metabolic hypothesis” that causes initial dip in oxy-Hb where vasoactive signals are generated around the capillaries that may be transmitted upstream to the penetrating artery and pial artery (see Fig. 7C) if the vasoactive signal intensity is large enough.
In the current study, fNIRS data consisted of changes in tHb at the ipsilateral (to HD-tDCS) and the contralateral sensorimotor cortex, where the biological criteria26 were formulated in the frequency band of 0.01–0.05 Hz55. So, tHb time series data were down-sampled to 0.1Hz to study the slow transient response during HD-tDCS. We applied a system identification approach using a physiologically constrained linear model to capture the fNIRS-based CVR to anodal HD-tDCS in healthy humans where the pathway from perivascular potassium to vessel circumference (i.e., Pathway 3) presented the lowest MSE (median < 2.5%), as shown in Fig. 4. Pathway 3 also gave the lowest MSE (median < 0.5%) when fitted to the whole 10 min of tHb time series during HD-tDCS (see Supplementary Table S4 and Supplementary Figure S2 in the supplementary materials). Therefore, the primary mechanism of action for the HD-tDCS-induced CVR between 0.01 and 0.05 Hz is postulated to be perivascular potassium that can be a significant determinant of the local pial vessel diameter56. Figure 10 shows all the pathways contributing to HD-tDCS-induced CVR depending on the vessel circumferencesensitivity to those pathways. Also, pathway response time will be different, e.g., tDCS-induced changes in the synaptic potassium (i.e., Pathway 1) had the slowest effect on vessel circumference. Recent studies showed that tDCS-induced alterations in cerebral CBF could only partially be related to the cortical excitability changes8. Therefore, our mechanistic investigation of the vascular response to tDCS using a physiologically detailed model (differential equations in the supplementary materials, section A: Model Equations for Neurovascular Compartmental Dynamics) of the four pathways, namely, via synaptic space, astrocytic space, perivascular space, and smooth muscle cell space (depicted in Figs. 9 and 10) provided novel evidence supporting the postulated role of perivascular potassium in CVR to tDCS. Pathway 1 considered modulation of synaptic mediators by tDCS current density, which then mediates the vessel circumference's effect via perivascular and astrocyte compartments through nested neurovascular dynamics. This is based on conventional neurovascular coupling mechanism under the effect of tDCS current density on the neurons, where neuronal activity drives the hemodynamic response. Pathway 2 considered the change in the membrane potential of the astrocytic compartment due to tDCS current density, which then drives the vascular response (vessel circumference). Pathway 3 considered modulation of the potassium currents in the perivascular space by tDCS current density. Pathway 4 considered the direct influence of tDCS current density on the vasculature's smooth muscle cells leading to the vascular response (vessel circumference). Pathway 4 was based on prior evidence that the voltage-gated potassium channels on the SMC can respond to the electric field57,58, especially those in the pial vasculature. In our proposed model for Pathway 4, SMC's membrane potential was modulated by the tDCS current density that modulated the voltage-gated potassium channel currents. In this study, we found the potassium currents, either in the synaptic compartment (Pathway 1) or in the perivascular compartment (Pathway 3), as the primary vasoactive signal related to CVR between 0.01 and 0.05 Hz. The vasoactive signal modulation by tDCS current density was represented mathematically as a first-order filter as cascaded transfer function – see Fig. 10. Here, after the minimization of the cost to find the optimal parameters, the cost needs to be evaluated based on Chi-Square Goodness-of-Fit that determines the quality of fit59 for such nested models. Table 2 shows the reduced dimension transfer functions for the four model pathways. Here, reduced dimension grey-box linear model for the Pathway 3 with 9 poles and 2 zeros (all free parameters) provided the best Chi-Square Goodness-of-Fit of 0.0078.
In this study, our goal was to investigate the role of various neuronal and non-neuronal pathways leading to the CVR to tDCS, as demonstrated experimentally by our prior works42,60,61, through an objective physiologically constrained grey-box model that is summarized by a block diagram in Fig. 10. The physiologically detailed model was primarily derived from Witthoft and Karniadakis38. The first model pathway considered the effect of tDCS on the synaptic space's vasoactive agent (potassium ions) where the linear analysis at the model's initial condition was represented by seventeen states, i.e., differential equations (details in supplementary materials, section A: Model Equations for Neurovascular Compartmental Dynamics). The transfer function order denotes the number of model poles where the reduced dimension grey-box model for the first pathway was represented by a twelfth order system that included a first-order linear filter for the tDCS’s vasoactive effects – see Table 2. The reduced dimension grey-box model for the astrocytic-driven (Pathway 2) and the perivascular potassium-driven (Pathway 3) pathways were represented by eleventh and ninth order systems, respectively. For Pathway 4, the input path considered the influence of tDCS current density on SMC gated ion channel current, and the reduced dimension grey-box model for this pathway was represented as a seventh order system. The impulse response function of the Pathway 1 defining conventional hemodynamic response function peaked around 5 sec (see Fig. 5(A)) that was comparable to known hemodynamic responses62. This provided insights into the temporal dynamics where the vessel response through an astrocytic pathway or perivascular potassium pathway peaked around 2 sec found comparable to known capillary responses13. Our study attempted to identify the role of various neuronal and non-neuronal pathways using a reduced dimension lumped grey-box model of NVU based on the vascular response (using tHb) to anodal HD-tDCS in healthy humans. However, the linearized physiologically detailed model of the four pathways could not investigate the effects of the current density on low-frequency (0.05–0.2 Hz) oscillations driven primarily by the nonlinear calcium dynamics. In this study, the data was down-sampled to 0.1 Hz that removed the spontaneous rhythmic nonlinear calcium dynamics in the frequency range of 0.05–0.2 Hz63, including the ∼0.1 Hz hemodynamic oscillations in the fNIRS time series42. Here, ∼0.1 Hz hemodynamic oscillations can be related to the synchronization of the intermittent release of calcium within vascular SMC63 that was modeled in the physiologically detailed model, as shown by Fig. 2(B) (also see the equations 48–49 in the supplementary materials, section A: Model Equations for Neurovascular Compartmental Dynamics). Our prior works have found a significant cross-correlation between log (base-10) transformed electroencephalogram (EEG) band-power (0.5-11.25Hz) and fNIRS oxy-Hb signal in that low frequency (≤ 0.1Hz) range42. Investigation of ∼0.1 Hz oscillatory vessel response vis-à-vis neuronal response (EEG) to tDCS will require a parameterized coupled oscillator model. Here, the relation of the ∼0.1 Hz oscillatory vessel response vis-à-vis neuronal response may be related to the cortical excitability changes to anodal tDCS60 due to the involvement of calcium dynamics that needs future investigation using transcranial alternating current stimulation (tACS). Unlike tDCS, tACS can lead to physiological entrainment at frequency of stimulation which can provide physiological insights based on grey-box modeling.
In summary, our study presented a preliminary investigation of linear systems analysis using a physiologically-constrained grey-box model that was found useful to explore various pathways in a lumped NVU model related to CVR to HD-tDCS. Such grey-box linear systems analysis using fNIRS-tHb data from individuals with pathological conditions61 can elucidate dysfunction in various NVU pathways, e.g., due to the pathological infestation. Here, our proposed linear systems analysis using the grey-box model is amenable to pole-zeros analysis of the transfer function for various pathways of the neurovascular system in health, aging, and disease that can also be used to classify the dysfunction26. This is crucial since neurovascular coupling dynamics are complicated phenomena in humans, and it can be hard to uncouple neuronal and vascular effects without mechanistic model-based hypothesis testing. The association between neuronal activity and hemodynamic responses can also be indirectly examined through other functional neuroimaging techniques such as fMRI and PET. Here, the integration of tDCS with functional neuroimaging modalities holds immense promise for throwing light on the underlying neuromodulation processes of current density effects. Experimental studies have indicated that long-term tDCS tends to transform neuronal activity through an induced electric field modulation of the cortical neurotransmitters (like gamma-aminobutyric acid and glutamate) during tDCS24,64. Besides, there is evidence of studies on CBF modulation through tDCS5,6,8. Of the available functional neuroimaging technologies, fNIRS, a portable, noninvasive clinically available tool, allows monitoring the local cortical hemodynamic response with reasonable spatial resolution and better temporal resolution than the gold standard of fMRI. fNIRS provides a measurement of changes in the cerebral oxygenation (oxy-Hb and deoxy-Hb) and blood volume (sum of oxy-Hb and deoxy-Hb)65,66 that is a promising tool to evaluate CVR to tDCS as evident from several studies67,68. The modality has been extensively used in various brain diseases like epilepsy, stroke, Parkinson’s disease, and mild cognitive impairment over the last few decades69. It can provide an indirect measure of CBF during tDCS where the advantage over other neuroimaging modalities like fMRI and PET are: portability, better safety, higher temporal resolution, and cost effectiveness70. Also, fNIRS can be combined with tDCS with no electro-optic interference to measure hemodynamic response during electrical stimulation. Moreover, dosing of the tDCS-induced current density can be monitored for safety in a diseased state since the blood-brain barrier (BBB) dysfunction in severe cases71 can be worsened by an increased BBB permeability53. Therefore, our study provided fNIRS based rational approach to investigate the underlying mechanism of the hemodynamic response to tDCS that is crucial for the clinical translation of transcranial electrical stimulation approaches in cerebrovascular diseases due to its ease of use and low cost72,73. Also, fNIRS measurement of Cytochrome-C-Oxidase74 has been shown feasible that is important due to the relation of vascular density and the cytochrome oxidase activity13.
Indeed, various studies have evaluated the hemodynamic effects of tDCS using fNIRS in humans in published literature. Merzagora et al.75 assessed the changes in prefrontal cortical oxygenation related to tDCS using fNIRS in healthy participants at rest. A large increase in oxy-Hb was observed in the 10min period following anodal tDCS compared to baseline levels before tDCS. Muthalib et al44 showed that anodal HD-tDCS induced more significant increases in oxy-Hb significantly during 10min of stimulation (at 2mA) in the sensorimotor cortex region within the vicinity of the 4 x 1 HD-tDCS montage compared to the region outside this boundary. At the same time, there were minimal oxy-Hb changes in the contralateral non-stimulated sensorimotor cortex region. Yaqub et al.76 evaluated the prefrontal cortex resting-sate intra-hemispheric and inter-hemispheric connectivity changes induced by 10min (1mA) anodal 4x1 HD-tDCS in healthy participants. Also, network connectivity changes were observed by analyzing filtered oxy-Hb time-series signals using standard graph theory metrics. Compared to the pre-stimulation phase, Yaqub et al.76 observed that the oxy-Hb levels and the resting-state connectivity of the prefrontal cortex increased during and after the stimulation, and the connectivity changes were more significant in the vicinity of the 4x1 HD-tDCS electrodes. Sood et al.42 extended Dutta et al.61 study and presented an autoregressive model parameter estimator method using Kalman filter to evaluate the relationship between changes in the fNIRS oxy-Hb signal and the EEG bandpower signal during anodal HD-tDCS. The time-varying poles of the autoregressive model were found to be comparable in all the healthy participants. In contrast, the zeros of the model exhibited variations across the participants during HD-tDCS.
Limitations of this study included the methodical limitations of the fNIRS technique77. The fNIRS signal acquired with optodes placed on the scalp can represent different hemodynamic signal sources (cerebral versus extra-cerebral) and other physiological causes (neuronal versus systemic) that can be evoked by tDCS. Due to the lack of short-separation channels to perform short source-detector regression to remove extra-cerebral hemodynamics, we performed a data-driven principal component analysis to identify the extra-cerebral signal components that explained the most significant amount of covariance across all the 16 spatially symmetrically distributed fNIRS channels. We found that Pathway 3 presented the least MSE, as shown in Fig. 4, so the fNIRS-tHb signal is postulated to have a significant representation from the superficial pial vessels. Differential activation of the oxy-Hb and deoxy-Hb (see section D: Oxy and Deoxy-Hemoglobin Changes at the Stimulated Region of the supplementary materials) also indicated neuronal activation and “metabolic hypothesis” for the cerebral capillaries in the vicinity of the HD-tDCS electrodes, which was used as a marker of neuronal cause. During grey-box modeling, the physiologically detailed model was substantially simplified by model linearization that lost nonlinear system dynamics, where calcium dynamics may be necessary for understanding the effects of tDCS on the rhythmic oscillations in the frequency range of 0.05–0.2 Hz63. Such rhythmic nonlinear calcium dynamics are postulated to be affected by tDCS current density, which will be investigated in our future study by adding the nonlinear calcium dynamics to our linear model. In this study, the parameterized grey-box linear model was fitted to slow transients in 0.01 and 0.05 Hz of fNIRS tHb, since we performed downsampling to 0.1 Hz that removed oscillations in the frequency range of 0.05–0.2 Hz. So, our investigation considered the acute effects of anodal tDCS by considering potassium ions as the main vasoactive agent with lumped multi-compartmental model of NVU.