Transcranial direct current stimulation (tDCS) has been shown to evoke hemodynamics response; however, the mechanisms have not been investigated systematically using systems biology approaches. We postulate that such a mechanistic understanding of the hemodynamics response, called cerebrovascular reactivity (CVR) to tDCS, can facilitate adequate delivery of the tDCS current density to the neurovascular tissue in cerebrovascular diseases. Our study presents a systems biology approach to evaluate CVR during high-definition (HD) tDCS using a physiologically constrained grey-box linear model. Our grey-box linear model was developed for the multi-compartmental neurovascular unit consisting of the vascular smooth muscle, perivascular space, synaptic space, and astrocyte glial cell. The physiologically detailed model generated vessel oscillations in the frequency range of 0.05–0.2 Hz driven mainly by non-linear calcium dynamics. Then, model linearization was performed to develop a grey-box linear model for evaluating the acute effects during the first 150 seconds of anodal HD-tDCS in eleven healthy humans using functional near-infrared spectroscopy (fNIRS) based measure of blood volume. The grey-box linear model was fitted to the total hemoglobin concentration (tHb) changes with optodes in the vicinity of 4x1 HD-tDCS electrodes. We found that the grey-box model pathway from perivascular potassium to the vessel circumference (Pathway 3) presented the best fit to fNIRS tHb time series with the least mean square error (MSE, median < 2.5%). Then, minimal realization transfer function with reduced-order approximations of the grey-box model pathways was fitted to the average tHb time series. Pathway 3, with nine poles and two zeros (all free parameters), provided the best Chi-Square Goodness of Fit of 0.0078. Therefore, our study provided a sound systems biology approach to investigate the hemodynamic response to tDCS that needs further investigation in health, aging, and disease. Future studies can leverage transcranial alternating current stimulation for frequency-dependent physiological entrainment of relevant pathways, e.g., oscillations driven by nonlinear calcium dynamics, that can provide additional insights based on our grey-box modeling approach.