Estimation of directional coupling between cortical areas using Near-Infrared Spectroscopy ( NIRS )

This study invesitigated the feasibility of measuring directional coupling between cortical areas with near-infrared spectroscopy (NIRS). Cerebral hemodynamic responses were recorded at the primary somatosensory cortex (S1), secondary somatosensory cortex (S2), and primary motor cortex (M1) regions of the rat barrel cortex during electrical stimulation of rat whiskers. Deoxyhemoglobin concentration changes were calculated from NIRS recordings and the Granger causality based on the multivariate autoregressive (MVAR) model was used to estimate the effective causal connectivity among S1, S2, and M1. The estimated causality patterns of seven rats showed consistent unidirectional coupling between the somatosensory areas and the motor areas (S1 and S2 → M1), which coincided well with our hypothesis because the rats’ motor function was completely anesthetized. Our preliminary results suggest that corticocortical directional coupling can be successfully investigated with NIRS. ©2010 Optical Society of America OCIS codes: (170.0170) Medical optics and biotechnology; (170.1610) Clinical applications; (170.2655) Functional monitoring and imaging References and links 1. K. J. Friston, L. Harrison, and W. Penny, “Dynamic causal modelling,” Neuroimage 19(4), 1273–1302 (2003). 2. L. Lee, K. Friston, and B. Horwitz, “Large-scale neural models and dynamic causal modelling,” Neuroimage 30(4), 1243–1254 (2006). 3. O. David, S. J. Kiebel, L. M. Harrison, J. Mattout, J. M. Kilner, and K. J. Friston, “Dynamic causal modeling of evoked responses in EEG and MEG,” Neuroimage 30(4), 1255–1272 (2006). 4. S. J. Kiebel, O. David, and K. J. Friston, “Dynamic causal modelling of evoked responses in EEG/MEG with lead field parameterization,” Neuroimage 30(4), 1273–1284 (2006). 5. L. Astolfi, F. Cincotti, D. Mattia, C. Babiloni, F. Carducci, A. Basilisco, P. M. Rossini, S. Salinari, L. Ding, Y. Ni, B. He, and F. Babiloni, “Assessing cortical functional connectivity by linear inverse estimation and directed transfer function: simulations and application to real data,” Clin. Neurophysiol. 116(4), 920–932 (2005). 6. B. Abler, A. Roebroeck, R. Goebel, A. Höse, C. Schönfeldt-Lecuona, G. Hole, and H. Walter, “Investigating directed influences between activated brain areas in a motor-response task using fMRI,” Magn. Reson. Imaging 24(2), 181–185 (2006). 7. J. R. Sato, E. A. Junior, D. Y. Takahashi, M. de Maria Felix, M. J. Brammer, and P. A. Morettin, “A method to produce evolving functional connectivity maps during the course of an fMRI experiment using wavelet-based time-varying Granger causality,” Neuroimage 31(1), 187–196 (2006). 8. W. Hesse, E. Möller, M. Arnold, and B. Schack, “The use of time-variant EEG Granger causality for inspecting directed interdependencies of neural assemblies,” J. Neurosci. Methods 124(1), 27–44 (2003). 9. D. W. Gow, Jr., J. A. Segawa, S. P. Ahlfors, and F. H. Lin, “Lexical influences on speech perception: a Granger causality analysis of MEG and EEG source estimates,” Neuroimage 43(3), 614–623 (2008). 10. A. Korzeniewska, C. M. Crainiceanu, R. Kuś, P. J. Franaszczuk, and N. E. Crone, “Dynamics of event-related causality in brain electrical activity,” Hum. Brain Mapp. 29(10), 1170–1192 (2008). 11. M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, “Phase synchronization of chaotic oscillators,” Phys. Rev. Lett. 76(11), 1804–1807 (1996). #122234 $15.00 USD Received 4 Jan 2010; revised 16 Feb 2010; accepted 18 Feb 2010; published 5 Mar 2010 (C) 2010 OSA 15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 5730 12. G. Pfurtscheller, and C. Andrew, “Event-Related changes of band power and coherence: methodology and interpretation,” J. Clin. Neurophysiol. 16(6), 512–519 (1999). 13. F. Mormann, K. Lehnertz, P. David, and C. E. Elger, “Mean phase coherence as a measure for phase synchronization and its application to the EEG of epilepsy patients,” Physica D 144(3-4), 358–369 (2000). 14. C. W. J. Granger, “Investigating Causal Relations by Econometric Models and Cross-Spectral Methods,” Econometrica 37, 414-& (1969). 15. L. Ding, G. A. Worrell, T. D. Lagerlund, and B. He, “Ictal source analysis: localization and imaging of causal interactions in humans,” Neuroimage 34(2), 575–586 (2007). 16. F. De Vico Fallani, L. Astolfi, F. Cincotti, D. Mattia, M. G. Marciani, S. Salinari, J. Kurths, S. Gao, A. Cichocki, A. Colosimo, and F. Babiloni, “Cortical functional connectivity networks in normal and spinal cord injured patients: Evaluation by graph analysis,” Hum. Brain Mapp. 28(12), 1334–1346 (2007). 17. M. A. Franceschini, D. K. Joseph, T. J. Huppert, S. G. Diamond, and D. A. Boas, “Diffuse optical imaging of the whole head,” J. Biomed. Opt. 11(5), 054007 (2006). 18. T. Austin, A. P. Gibson, G. Branco, R. M. Yusof, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three dimensional optical imaging of blood volume and oxygenation in the neonatal brain,” Neuroimage 31(4), 1426–1433 (2006). 19. S. Muehlschlegel, J. Selb, M. Patel, S. G. Diamond, M. A. Franceschini, A. G. Sorensen, D. A. Boas, and L. H. Schwamm, “Feasibility of NIRS in the neurointensive care unit: a pilot study in stroke using physiological oscillations,” Neurocrit. Care 11(2), 288–295 (2009). 20. T. Durduran, C. Zhou, B. L. Edlow, G. Yu, R. Choe, M. N. Kim, B. L. Cucchiara, M. E. Putt, Q. Shah, S. E. Kasner, J. H. Greenberg, A. G. Yodh, and J. A. Detre, “Transcranial optical monitoring of cerebrovascular hemodynamics in acute stroke patients,” Opt. Express 17(5), 3884–3902 (2009). 21. J. Sepulcre, J. C. Masdeu, M. A. Pastor, J. Goñi, C. Barbosa, B. Bejarano, and P. Villoslada, “Brain pathways of verbal working memory: a lesion-function correlation study,” Neuroimage 47(2), 773–778 (2009). 22. N. Shimoda, K. Takeda, I. Imai, J. Kaneko, and H. Kato, “Cerebral laterality differences in handedness: a mental rotation study with NIRS,” Neurosci. Lett. 430(1), 43–47 (2008). 23. F. Babiloni, F. Cincotti, C. Babiloni, F. Carducci, D. Mattia, L. Astolfi, A. Basilisco, P. M. Rossini, L. Ding, Y. Ni, J. Cheng, K. Christine, J. Sweeney, and B. He, “Estimation of the cortical functional connectivity with the multimodal integration of high-resolution EEG and fMRI data by directed transfer function,” Neuroimage 24(1), 118–131 (2005). 24. S. Aydore, and M. Mihcak, “On Temporal Connectivity of PFC Via Gauss-Markov Modeling of fNIRS Signals,” IEEE Trans. Biomed. Eng. (2009). 25. G. Paxinos, and C. Watson, The rat brain in stereotaxic coordinates (Elsevier Academic Press, Amsterdam; Boston, 2005). 26. C. I. Moore, S. B. Nelson, and M. Sur, “Dynamics of neuronal processing in rat somatosensory cortex,” Trends Neurosci. 22(11), 513–520 (1999). 27. F. Crespi, A. Bandera, M. Donini, C. Heidbreder, and L. Rovati, “Non-invasive in vivo infrared laser spectroscopy to analyse endogenous oxy-haemoglobin, deoxy-haemoglobin, and blood volume in the rat CNS,” J. Neurosci. Methods 145(1-2), 11–22 (2005). 28. S. Lee, D. Koh, K. Kwon, H. J. Lee, Y. Lang, H. C. Shin, and B. M. Kim, “Hemodynamic Responses of Rat Brain Measured by Near-infrared Spectroscopy During Various Whisker Stimulations,” J. Opt. Soc. Korea 13(1), 166–170 (2009). 29. S. Lee, M. Lee, D. Koh, B. M. Kim, and J. H. Choi, “Cerebral hemodynamic responses to seizure in the mouse brain: A simultaneous NIRS-EEG study,” J. Biomed. Opt. (to be published). 30. D. T. Delpy, M. Cope, P. Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33(12), 1433–1442 (1988). 31. J. C. Ye, S. Tak, K. E. Jang, J. Jung, and J. Jang, “NIRS-SPM: statistical parametric mapping for near-infrared spectroscopy,” Neuroimage 44(2), 428–447 (2009). 32. T. J. Huppert, R. D. Hoge, S. G. Diamond, M. A. Franceschini, and D. A. Boas, “A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans,” Neuroimage 29(2), 368–382 (2006). 33. G. Schwarz, “Estimating Dimension of a Model,” Ann. Stat. 6(2), 461–464 (1978). 34. J. Geweke, R. Meese, and W. Dent, “Comparing Alternative Tests of Causality in Temporal Systems Analytic Results and Experimental-Evidence,” J. Econom. 21(2), 161–194 (1983). 35. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Incomplete Beta Function, Student's Distribution, F-Distribution, Cumulative Binomial Distribution (Cambridge University Press, 1993). 36. J. Hamilton, Covariance-Stationary Vector Processes (Princeton University Press, 1994). 37. A. K. Seth, ““Causal connectivity of evolved neural networks during behavior,” Network-Comp Neural. 16(1), 35–54 (2005). 38. A. Schlogl, “A comparison of multivariate autoregressive estimators,” Signal Process. 86(9), 2426–2429 (2006). 39. C. C. Petersen, “The functional organization of the barrel cortex,” Neuron 56(2), 339–355 (2007). 40. S. Chakrabarti, and K. D. Alloway, “Differential origin of projections from SI barrel cortex to the whisker representations in SII and MI,” J. Comp. Neurol. 498(5), 624–636 (2006). 41. L. Petreanu, D. Huber, A. Sobczyk, and K. Svoboda, “Channelrhodopsin-2-assisted circuit mapping of longrange callosal projections,” Nat. Neurosci. 10(5), 663–668 (2007). 42. T. Schreiber, “Measuring information transfer,” Phys. Rev. Lett. 85(2), 461–464 (2000). 43. G. Nolte, A. Ziehe, V. V. Nikulin, A. Schlögl, N. Krämer, T. Brismar, and K. R. Müller, “Robustly estimating the flow direction of information in complex physical systems,” Phys. Rev. Lett. 100(23), 234101 (2008). #122234 $15.00 USD Received 4 Jan 2010; revised 16 Feb 2010; accepted 18 Feb 2010; published 5 Mar 2010 (C) 2010 OSA 15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 5731 44. A. Gallagher, M. Lassonde, D. Bastien, P. Vannasing, F. Lesage, C. Grova, A. Bouthillier, L. Carmant, F. Lepore, R. Beland, and D. K. Nguyen, ““Non-invasive pre-surgical investigation of a 10 year-old epileptic boy using simultaneous EEG-NIRS,” Seizure-Eur, J. Epilepsy 17, 576–582 (2008). 4

Directional coupling is sometimes referred to as causality, effective connectivity, or functional causality.Various algorithms to measure the directional coupling between different brain areas have been proposed and applied to both hemodynamics-based (fMRI) and electrophysiology-based (EEG, MEG, and ECoG) modalities.Among many algorithms, dynamic causal modeling (DCM) and multivariate autoregressive (MVAR) modeling have been the most widely used.DCM first assumes an underlying brain network model, regions of interests, and a priori connectivity probabilities, then estimates the effective connectivity values between the brain regions using stimulus functions describing experimental manipulations and measured time series as the inputs and outputs of the dynamic system, respectively.DCM has been successfully applied to fMRI studies [1,2] and EEG/MEG studies [3,4], mostly associated with cognitive tasks.However, DCM has a disadvantage in that it requires strong a priori assumptions on the number of directional connections and their directionality, which are difficult to validate in practical neuroimaging studies.On the contrary, an MVAR-based modeling technique called "Granger causality" does not require any a priori assumptions on the connectivity models.Granger causality, proposed by Granger in 1969 [14], infers directional influence between two cortical areas by quantifying improvement of prediction in one area by additional information of the other area.We can judge whether there exists the Granger causality from 'X' to 'Y' if the past information from both time series can significantly improve the prediction of the future of 'Y' rather than using the past information from 'Y' alone.The future time series have been usually estimated using MVAR modeling [5].Directed transfer function (DTF) and partial directed coherence (PDC) are frequency domain versions of Granger causality and have been generally applied to EEG/MEG time series analyses.During the last decade, MVAR-based modeling has been extensively applied to functional neuroimaging studies based on fMRI [6,7], EEG [8], MEG [9], and ECoG [10].The MVAR-modeling-based causality analysis is now being applied to clinical neuroscience as well as to cognitive neuroscience and several articles in these fields have already been published [15,16].
Near-infrared spectroscopy (NIRS) is a relatively new functional imaging modality in neuroscience.It uses NIR light to non-invasively estimate various functional parameters of a brain such as cerebral blood volume and tissue oxygenation by measuring absorption changes caused by hemoglobin level changes.Recently, NIRS has been actively applied to functional brain imaging [17,18], clinical neurophysiology [19,20] and cognitive neuroscience studies [21,22].NIRS has some advantages over other neuroimaging modalities in that it provides a relatively high temporal resolution, it is cheap to build, and is easily maintained.Despite these advantages, NIRS studies have focused only on classical functional neuroimaging and has not been widely used for investigating cortico-cortical directional coupling.Causality estimation with NIRS would have an advantage over the EEG or MEG-based analyses since it is not necessary to solve inverse problems, which are indispensable for EEG/MEG-based analyses because EEG or MEG topographies are not directly attributed to the underlying cortical regions [23].Aydöre et al. [24] first applied the MVAR-based causality estimation method to NIRS data.They demonstrated the feasibility of investigating functional connectivity with NIRS.In their experiment, the functional connectivity of the prefrontal cortex showed close relevance to the cognitive load, indirectly suggesting that the functional connectivity information estimated from NIRS data is physiologically meaningful.
In this study, we report another strong piece of evidence demonstrating that the corticocortical directional coupling can be investigated with NIRS.We applied the Granger causality method to frequency-domain NIRS data recorded from the rat brain during whisker stimulation and estimated the directional coupling among primary somatosensory cortex (S1), secondary somatosensory cortex (S2), and primary motor cortex (M1).

Animal preparation
All animal preparation and experimental protocols followed the institutional guidelines of Korea University.Seven adult female rats (Sprague-Dawley, 210-280 g) were anesthetized with urethane (20%, 1.25 g/kg body weight, i.p.) and mounted on a stereotaxic frame.Prior to surgery and optical recording, the rat was tested to ensure the lack of motor response.Following a one-inch middle incision, the scalp was clearly removed and skull was exposed.For the NIRS recordings, we defined coordinates of the S1 (anterior-posterior: −2.0 mm, lateral: + 3.0 mm), S2 (anterior-posterior: −1.0 mm, lateral: + 6.0 mm) and M1 (anteriorposterior: + 2.0 mm, lateral: + 4.0 mm) area from bregma point on the right hemisphere [25].In order to stimulate the whole whisker, two biostimulator needle electrodes (Bcotron, Digitalstimulator BDS 301, South Korea) were inserted at B4 and E1 sites of the left whisker area [26].

NIRS and biostimulator setup
A commercial eight-channel frequency-domain NIRS system (Imagent, ISS, Il., U.S.A.) was used to measure the cerebral hemodynamic responses to whisker stimulation.This system uses two wavelengths (690 nm, 830 nm), and each channel contains two optical fibers with 400 µmcorediameter (FT-400EMT, Thorlabs, N.J., U.S.A) for 690 nm and 830 nm.The optical fibers are arranged as shown in Fig. 1(a) and the picture of the actual experiment is shown in Fig. 1(b).The source-detector separation for each channel was approximately 5 -6 mm and the corresponding penetration depth is estimated as 1.5 -2 mm, which is sufficient to cover the cortex area [27].After the optical fibers are placed at appropriate lications, dental cement (Ketac Cem, 3M, U.S.A) was applied to the rat skull to maintain a constant boundary condition between optical fibers and skull [28,29].The NIRS signals were recorded at 30 Hz sampling rate.The rat whisker was stimulated for 10 s with 10 Hz, 10 V electrical pulse trains, followed by 90 s of resting time for recovery, as shown in Fig. 2.These procedures were repeated approximately 40-60 times.

Calculation of hemoglobin concentration
A modified Beer-Lambert law was used to calculate the oxy-(HbO 2 ) and the deoxyhemoglobin (Hbr) concentrations as in a previous study [30]: where OD is the optical density, I Final is the measured light intensity, I Initial is the initial light intensity, ε represents the extinction coefficients of the oxy-and the deoxy-hemoglobin, B is the differential path-length factor, L is the source-detection separation, and λ is the wavelength of light (690 nm and 830 nm).Using the calculated ∆OD λ , we can obtain ∆[HbO 2 ] (oxyhemoglobin concentration) and ∆[Hbr] (deoxy-hemoglobin concentration) with the following equations: In this study, deoxyhemoglobin data was used to analyze the directional coupling of cortical areas because it is generally agreed that the concentration of deoxyhemoglobin is more tightly coupled to the fMRI BOLD signals [31,32], of which the causality analyses have already been studied in detail [1,2,6,7].

Granger causality
The directional coupling was estimated by computing the Granger causality [14] between each pair of time series acquired from different locations on the rat brain surface.The simplest way to exploit this technique is to use the linear prediction of the future values of data x t and y t with an autoregressive model, under the assumption that x t and y t are vectors of deterministic variables.The value of a time series x at a time slice t, x t , can be estimated using the following two different methods: where α and β are the AR-matrix coefficients, p is the model order, u t represents the uncorrelated Gaussian random variables with mean zero in a univariate signal, w t represents the uncorrelated Gaussian random variables with mean zero in the multivariate signal.M is the maximal model order, which was determined by using the Bayesian Information Criterion (BIC) [33].If y t causes x t , the prediction error obtained from univariate modeling, u t , should be larger than that obtained from multivariate modeling, w t .In this study, the Granger causality is evaluated using [34] 0 1 ) where and , where T is the number of observations.To assess the statistical significance of the estimated Granger causality, we adopted the F-test with null hypothesis, H 0 : β p = 0, that is, y t does not influence the generation of x t .If β p = 0 for all p = 1, 2, …, M, the Granger causality value F becomes zero as RSS 0 equals RSS 1 .On the contrary, if the null hypothesis is rejected, that is, F is sufficiently large; we can conclude that y t causes x t .The test for the hypothesis is performed using a traditional F-test resulting from an ordinary least squared regression for each equation.
To test the statistical significance of F, the cumulative F distribution is estimated first [35] and then the probability of the F value is calculated by p = 1 -CDF(F), where CDF represents the cumulative distribution function and is estimated using the 'cdf' function implemented in Matlab 7.6 (Mathworks, U.S.A), where p represents the probability of the Granger causality.For example, p = 1 means that there exists no causal interaction between the two time series; whereas p = 0 means that there is a strong directional influence from y t to x t .We rejected the null hypothesis when the probability p is less than 0.05.The above concepts can be readily expanded to multivariate cases by exploiting an MVAR model.MVAR coefficients were calculated using a multivariate linear regression method called ordinary least squares [36].The signals applicable for the Granger causality analysis must be covariance-stationary, which means that the input signals maintain consistent mean and variance [36].In the case of non-stationary signals, we applied first-order differentiation to the signals repeatedly until the signals become stationary [37].For more detailed information on the implementation of the Granger causality method, please refer to the literatures [7][8][9]14,34].

Results
With the NIRS signals obtained from seven rat brains, we calculated the deoxyhemoglobin concentration using Eq.(3). Figure 3 shows an example of temporal changes of the deoxyhemoglobin concentration acquired from a rat (rat no. 5) when electrical pulse trains at 10 Hz were applied for 10 s.Data were acquired from three different locations, S1, S2, and M1 as demonstrated in Fig. 1, and they showed clear hemodynamic responses to the electrical stimuli.The tendency of signal changes look similar between three locations, but we cannot infer any causal interactions or time delays simply by visual inspection.
We evaluated the Granger causality based on the MVAR model for each of the trials.Window size for the causality evaluation was set to be 30 s starting from the stimulation onset time (0 s).According to Schlögl's study [38], the accuracy of the MVAR estimation is somewhat dependent on the number of time samples and the variance of the prediction error becomes smaller than 10% of the signal variance if the number of samples is larger than 70.It was also suggested that different MVAR estimators do not show significant differences for sufficiently large number of samples (N > 400).In the present study, we used 900 samples for the MVAR estimation; therefore, we believe that our MVAR estimation should be sufficiently accurate.Figure 4 shows the averaged probability values of the 7 rats, where the strength of the Granger causality is presented using the mean probability values, different arrow colors, and different arrow thicknesses.Although some individual differences in the p-value level are observed, the causal interaction patterns look similar and consistent.Such individual differences in p-values are believed to be originated from several factors such as different signal-to-noise ratio, slightly different sensor locations, different anesthetic depths, and so on.When p = 0.05 was used to assess statistical significance, all the results showed statistically significant Granger causality values for S1 → M1 connections and 6 out of 7 rats showed meaningful Granger causality for S2 → M1.Most rats did not show statistically significant causality values of the opposite directional connections (M1 → S1 and M1 → S2), except for one case (rat 5, M1 → S2), which will be discussed later.In all cases, directional coupling from the somatosensory cortex to the motor cortex dominated the opposite directional coupling.We also estimated the Granger causality values of the NIRS data acquired from the same rats during the resting period and could not find any consistent directional couplings.
Since the rats were fully anesthetized to suppress their neuromuscular function, the feedback somatosensory activation originating from motor action was reduced and rarely affected the NIRS signals recorded at the somatosensory areas.According to the Petersen's report [39] describing rat's anatomical structure, long-range, glutamatergic, cortico-cortical synaptic connections exist between S1 and S2 as well as between S1 and M1, but the direct anatomical connection between S2 and M1 had not been reported.However, since the connectivity or causality estimation relies on the recorded neuronal signals, the functional connectivity patterns estimated for EEG and fMRI time series using either Granger causality or DTF do not necessarily coincide with the known anatomical connectivity [5][6][7].Therefore, in the sense of the functional connectivity, it can be said that both S1 and S2 activities caused the activation of M1.In our experiment, the causal interactions between S1 and S2 did not show consistent directional connectivity patterns, preventing us from making any conclusions.We believe that the anatomical distance between the two areas (S1 and S2) is so close that the deoxyhemoglobin concentration changes on one area may have contaminated those on another area due to the limited spatial resolution of the NIRS recordings.Since the anatomical distance between somatosensory areas (S1 and S2) and M1 is relatively larger than that between S1 and S2, the causal interactions between the somatosensory areas and the motor area could be observed clearly.

Discussion
We have investigated if the cortico-cortical directional coupling in the rat brain could be estimated with NIRS.MVAR-based Granger causality estimation was applied to frequencydomain NIRS data recorded from seven rat brains during electrical stimulation of whiskers, and estimated the directional connectivity among S1, S2, and M1 locations on the right hemisphere [Fig.1(a)].The estimated causality patterns showed consistent directional coupling between the somatosensory cortex and the motor cortex, suggesting that the corticocortical causal interactions can be investigated with NIRS recordings.In our experiment using seven rats, the directional connectivity from S1 or S2 to M1 was observed for the majority of the cases (13 out of 14 connections).The opposing directional connection (rat 5, M1 → S2) was observed only for one case.
Recent neuroanatomy studies related to the functional mechanism of the rat barrel cortex [40,41] have shown that the primary somatosensory barrel cortex is tightly interconnected with the ipsilateral secondary somatosensory cortex and ipsilateral primary motor cortex.Since callosal projections have been reported to be less prominent than other connections [39], inter-hemispheric connections were not considered in this study.
Motor-to-somatosensory connections may be observed when the rats' motor functions were not fully anesthetized or if the anesthetic depth was shallow.Since we could not use a device to monitor the anesthetic depth, we simply applied the dose in proportion to their body weights.We believe that the opposite directional connectivity from M1 to S2 may be related to the depth of anesthesia.More sophisticated experiments using indices measuring the depth of anesthesia such as arterial oxygen saturation and arterial blood pressure will be performed in our future studies.
Various methods have been proposed to estimate the directional coupling among simultaneously recorded neural signals, including hemodynamic responses.DTF and PDC have been extensively applied to EEG and MEG time series analyses [5], but such frequency domain analysis methods are not adequate to be applied to hemodynamics-based modalities such as fMRI and NIRS, of which specific frequency components or frequency bands are not generally meaningful.The Granger causality used in this study has been the most widely used measure to estimate the directional coupling among fMRI BOLD signals acquired at different regions of interest.Recently, entropy transfer (or directed information transfer: DIT) [42] and phase slope index (PSI) [43] have been proposed to quantify the direction coupling.Although the new indices have not been sufficiently verified, they can be nice alternatives to solve the 'stationary' problem of the MVAR-based methods.We plan to apply these new indices to our NIRS data sets in future experiments and compare the results with this study.
An important clinical application of the cortico-cortical directional coupling estimation has been the localization of epileptogenic sources.Ding et al. [15] applied an EEG source localization algorithm called first principal vectors (FINEs) algorithm to ictal EEG data and investigated causal interactions among the source time series to estimate ictal onset zones and propagation patterns, which can be used as important information for pre-surgical evaluations.However, to investigate ictal propagation with EEG, inverse problems should be solved because the EEG topographies do not directly reflect the underlying brain activities.To solve the EEG inverse problems, complex pre-processing procedures such as artifact rejection, manual or semi-automatic head model generation, etc. are required.Moreover, ictal EEG data are generally highly contaminated by muscle artifacts, making it difficult to estimate accurate neuronal source locations.Difficulty in estimating the number of EEG sources is also one of the crucial problems of EEG dipolar source localization.
On the other hand, NIRS data directly reflect the underlying cortical activity changes and are relatively more resilient against muscle artifacts since light propagation does not interfere with electrical signals caused by muscle excitation.NIRS provides hemodynamic information rather than the electrical information, allowing for an extended recording time and can be recorded simultaneously with EEG [44].Therefore, the NIRS-based causality analysis can be a promising tool for noninvasively estimating ictal propagation patterns if simultaneously recorded EEG data are used to determine ictal onset [45].

Fig. 1 .
Fig. 1.(a) Arrangement of optical fibers (L1, L2, L3 are the locations of the source fibers and D represents a detector).M1, S1 and S2 represent primary motor cortex, primary somatosensory cortex and secondary somatosensory cortex, respectively.(b) Picture of the actual experiment.

Fig. 2 .
Fig. 2. The paradigm of electrical stimulation.Pulses at 10 Hz repetition rate lasted for 10 s and the resting time was set for 90 s.This process was repeated approximately 40-60 times.

Fig. 3 .
Fig. 3. Representative plots of changes in deoxyhemoglobin concentration at the (a) S1, (b) S2, and (c) M1 locations on the right hemisphere of a rat brain.10 Hz pulse trains lasted for 10 s.The negative offset of all the signals is part of the physiological noise, which is not meaningful.The probability values of the Granger causality, p, were then evaluated for 6 connections (S1 → M1, M1 → S1, S1 → S2, S2 → S1, S2 → M1, M1 → S2) and averaged for all trials.Figure4shows the averaged probability values of the 7 rats, where the strength of the Granger causality is presented using the mean probability values, different arrow colors, and different arrow thicknesses.Although some individual differences in the p-value level are observed, the causal interaction patterns look similar and consistent.Such individual differences in p-values are believed to be originated from several factors such as different signal-to-noise ratio, slightly different sensor locations, different anesthetic depths, and so on.When p = 0.05 was

Fig. 4 .
Fig. 4. The averaged Granger causality (probability) values of seven rats.The strength of the Granger causality is presented using the probability values, arrow colors, and arrow thicknesses.Smaller probability values, thicker arrows, and darker colors represent stronger causality along the direction.Each p-value was averaged from 40 -60 trials.