Bayesian filtering of human brain hemodynamic activity elicited by visual short-term maintenance recorded through functional near-infrared spectroscopy ( fNIRS )

Functional near-infrared spectroscopy (fNIRS) is a neuroimaging technique that measures changes in oxy-hemoglobin (ΔHbO) and deoxyhemoglobin (ΔHbR) concentration associated with brain activity. The signal acquired with fNIRS is naturally affected by disturbances engendering from ongoing physiological activity (e.g., cardiac, respiratory, Mayer wave) and random measurement noise. Despite its several drawbacks, the so-called conventional averaging (CA) is still widely used to estimate the hemodynamic response function (HRF) from noisy signal. One such drawback is related to the number of trials necessary to derive stable HRF functions adopting the CA approach, which must be substantial (N >> 50). In this work, a pre-processing procedure to remove artifacts followed by the application of a non-parametric Bayesian approach is proposed that capitalizes on a priori available knowledge about HRF and noise. Results with the proposed Bayesian approach were compared with CA and with a straightforward band-pass filtering approach. On simulated data, a five times lower estimation error on HRF was obtained with respect to that obtained by CA, and 2.5 times lower than that obtained by band pass filtering. On real data, the improvement achieved by the present method was attested by an increase in the contrast to noise ratio (CNR) and by a reduced variability in single trial estimation. An application of the present Bayesian approach is illustrated that was optimized to monitor changes in hemodynamic activity reflecting variations in visual short-term memory load in humans, which are notoriously hard to detect using functional magnetic resonance imaging (fMRI). In particular, statistical analyses of HRFs recorded during a memory task established with high reliability the crucial role of the intraparietal sulcus and the intra-occipital sulcus in posterior areas of the human brain in visual short-term memory maintenance. ©2010 Optical Society of America OCIS codes: (300.6340) Spectroscopy, infrared; (200.4560) Optical data processing; (170.2655) Functional monitoring and imaging; (170.3890) Medical optics instrumentation. References and links 1. D. A. Boas, M. A. Franceschini, A. K. Dunn, and G. Strangman, “Noninvasive Imaging of Cerebral Activation with Diffuse Optical Tomography,”in Vivo Optical Imaging of Brain Function. E. D. (CRC Press), Chap 8, pp. 193–221, (2002) 2. S. C. Bunce, M. Izzetoglu, K. Izzetoglu, B. Onaral, and K. Pourrezaei, “Functional Near-Infrared Spectroscopy,” IEEE Eng. Med. Biol. Mag. 25(4), 54–62 (2006). #133907 $15.00 USD Received 24 Aug 2010; revised 12 Nov 2010; accepted 14 Nov 2010; published 3 Dec 2010 (C) 2010 OSA 6 December 2010 / Vol. 18, No. 25 / OPTICS EXPRESS 26550 3. H. Obrig, and A. Villringer, “Beyond the visible--imaging the human brain with light,” J. Cereb. Blood Flow Metab. 23(1), 1–18 (2003). 4. S. P. Koch, S. Koendgen, R. Bourayou, J. Steinbrink, and H. Obrig, “Individual alpha-frequency correlates with amplitude of visual evoked potential and hemodynamic response,” Neuroimage 41(2), 233–242 (2008). 5. R. J. Cooper, N. L. Everdell, L. C. Enfield, A. P. Gibson, A. Worley, and J. C. Hebden, “Design and evaluation of a probe for simultaneous EEG and near-infrared imaging of cortical activation,” Phys. Med. Biol. 54(7), 2093– 2102 (2009). 6. A. Gibson, and H. Dehghani, “Diffuse optical imaging,” Philos. Transact. A Math. Phys. Eng. Sci. 367(1900), 3055–3072 (2009). 7. G. Jasdzewski, G. Strangman, J. Wagner, K. K. Kwong, R. A. Poldrack, and D. A. Boas, “Differences in the hemodynamic response to event-related motor and visual paradigms as measured by near-infrared spectroscopy,” Neuroimage 20(1), 479–488 (2003). 8. S. G. Diamond, T. J. Huppert, V. Kolehmainen, M. A. Franceschini, J. P. Kaipio, S. R. Arridge, and D. A. Boas, “Dynamic physiological modeling for functional diffuse optical tomography,” Neuroimage 30(1), 88–101 (2006). 9. V. Kolehmainen, S. Prince, S. R. Arridge, and J. P. Kaipio, “State-estimation approach to the nonstationary optical tomography problem,” J. Opt. Soc. Am. A 20(5), 876–889 (2003). 10. S. Prince, V. Kolehmainen, J. P. Kaipio, M. A. Franceschini, D. A. Boas, and S. R. Arridge, “Time-series estimation of biological factors in optical diffusion tomography,” Phys. Med. Biol. 48(11), 1491–1504 (2003). 11. A. F. Abdelnour, and T. J. Huppert, “Real-time imaging of human brain function by near-infrared spectroscopy using an adaptive general linear model,” Neuroimage 46(1), 133–143 (2009). 12. Y. H. Zhang, D. H. Brooks, M. A. Franceschini, and D. A. Boas, “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt. 10(1), 011014 (2005). 13. M. L. Schroeter, M. M. Bücheler, K. Müller, K. Uludağ, H. Obrig, G. Lohmann, M. Tittgemeyer, A. Villringer, and D. Y. von Cramon, “Towards a standard analysis for functional near-infrared imaging,” Neuroimage 21(1), 283–290 (2004). 14. S. Tak, K. E. Jang, J. W. Jung, J. Jang, and J. C. Ye, “General Linear Model and Inference for Near Infrared Spectroscopy using Global Confidence Region Analysis,” in Proceedings of IEEE Conference on International Symposium on Biomedical Imaging (ISBI), pp. 476–479 (2008). 15. 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). 16. K. E. Jang, S. Tak, J. Jung, J. Jang, Y. Jeong, and J. C. Ye, “Wavelet minimum description length detrending for near-infrared spectroscopy,” J. Biomedical Opt. 14, 034004-(1–13) (2009). 17. G. Morren, U. Wolf, P. Lemmerling, M. Wolf, J. H. Choi, E. Gratton, L. De Lathauwer, and S. Van Huffel, “Detection of fast neuronal signals in the motor cortex from functional near infrared spectroscopy measurements using independent component analysis,” Med. Biol. Eng. Comput. 42(1), 92–99 (2004). 18. C. B. Akgül, A. Akin, and B. Sankur, “Extraction of cognitive activity-related waveforms from functional nearinfrared spectroscopy signals,” Med. Biol. Eng. Comput. 44(11), 945–958 (2006). 19. A. V. Medvedev, J. Kainerstorfer, S. V. Borisov, R. L. Barbour, and J. VanMeter, “Event-related fast optical signal in a rapid object recognition task: improving detection by the independent component analysis,” Brain Res. 1236, 145–158 (2008). 20. G. Taga, K. Asakawa, A. Maki, Y. Konishi, and H. Koizumi, “Brain imaging in awake infants by near-infrared optical topography,” Proc. Natl. Acad. Sci. U.S.A. 100(19), 10722–10727 (2003). 21. R. Sitaram, H. Zhang, C. Guan, M. Thulasidas, Y. Hoshi, A. Ishikawa, K. Shimizu, and N. Birbaumer, “Temporal classification of multichannel near-infrared spectroscopy signals of motor imagery for developing a brain-computer interface,” Neuroimage 34(4), 1416–1427 (2007). 22. S. Cutini, P. Scatturin, E. Menon, P. S. Bisiacchi, L. Gamberini, M. Zorzi, and R. Dell’Acqua, “Selective activation of the superior frontal gyrus in task-switching: an event-related fNIRS study,” Neuroimage 42(2), 945–955 (2008). 23. H. Kojima, and T. Suzuki, “Hemodynamic change in occipital lobe during visual search: visual attention allocation measured with NIRS,” Neuropsychologia 48(1), 349–352 (2010). 24. G. Sparacino, S. Milani, E. Arslan, and C. Cobelli, “A Bayesian approach to estimate evoked potentials,” Comput. Methods Programs Biomed. 68(3), 233–248 (2002). 25. R. Luria, P. Sessa, A. Gotler, P. Jolicoeur, and R. Dell’Acqua, “Visual short-term memory capacity for simple and complex objects,” J. Cogn. Neurosci. 22(3), 496–512 (2010). 26. J. J. Todd, and R. Marois, “Capacity limit of visual short-term memory in human posterior parietal cortex,” Nature 428(6984), 751–754 (2004). 27. Y. Xu, and M. M. Chun, “Dissociable neural mechanisms supporting visual short-term memory for objects,” Nature 440(7080), 91–95 (2006). 28. M. Cope, and D. T. Delpy, “System for long-term measurement of cerebral blood and tissue oxygenation on newborn infants by near infra-red transillumination,” Med. Biol. Eng. Comput. 26(3), 289–294 (1988). 29. M. L. Schroeter, S. Zysset, F. Kruggel, and D. Y. von Cramon, “Age dependency of the hemodynamic response as measured by functional near-infrared spectroscopy,” Neuroimage 19(3), 555–564 (2003). 30. M. L. Schroeter, S. Cutini, M. M. Wahl, R. Scheid, and D. Yves von Cramon, “Neurovascular coupling is impaired in cerebral microangiopathy--An event-related Stroop study,” Neuroimage 34(1), 26–34 (2007). #133907 $15.00 USD Received 24 Aug 2010; revised 12 Nov 2010; accepted 14 Nov 2010; published 3 Dec 2010 (C) 2010 OSA 6 December 2010 / Vol. 18, No. 25 / OPTICS EXPRESS 26551 31. A. Duncan, J. H. Meek, M. Clemence, C. E. Elwell, P. Fallon, L. Tyszczuk, M. Cope, and D. T. Delpy, “Measurement of cranial optical path length as a function of age using phase resolved near infrared spectroscopy,” Pediatr. Res. 39(5), 889–894 (1996). 32. M. A. Franceschini, V. Toronov, M. E. Filiaci, E. Gratton, and S. Fantini, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express 6(3), 49–57 (2000). 33. M. Okamoto, H. Dan, K. Sakamoto, K. Takeo, K. Shimizu, S. Kohno, I. Oda, S. Isobe, T. Suzuki, K. Kohyama, and I. Dan, “Three-dimensional probabilistic anatomical cranio-cerebral correlation via the international 10-20 system oriented for transcranial functional brain mapping,” Neuroimage 21(1), 99–111 (2004). 34. S. Cutini, P. Scatturin, and M. Zorzi, “A new method based on ICBM152 head surface for probe placement in multichannel fNIRS,” Neuroimage (to be published). 35. M. R. Nuwer, G. Comi, R. Emerson, A. Fuglsang-Frederiksen, J. M. Guérit, H. Hinrichs, A. Ikeda, F. J. Luccas, and P. Rappelsburger, The International Federation of Clinical Neurophysiology, “IFCN standards for digital recording of clinical EEG,” Electroencephalogr. Clin. Neurophysiol. 106(3), 259–261 (1998). 36. A. K. Singh, M. Okamoto, H. Dan, V. Jurcak, and I. Dan, “Spatial registration of multichannel multi-subject fNIRS


Introduction
Functional near-infrared spectroscopy (fNIRS) is an emerging neuroimaging technique which uses the near-infrared region of the electromagnetic spectrum to measure changes in blood oxygenation.This region (i.e., 650-950 nm) is generally poorly absorbed by biological tissues other than oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR).Among other important implementations, this particular property, in combination with different absorption spectra of HbO and HbR, allows cognitive neuroscientists to explore brain activity by monitoring online variations in HbO and HbR concentration in the cerebral (cortical) blood flow during the execution of a cognitive task.
Over the past decade, a growing number of researchers have used fNIRS as an alternative to older and more established neuroimaging modalities, such as functional magnetic resonance imaging (fMRI) and positron emission tomography (PET), mainly because of its potential to provide non-invasively physiological information, its high temporal resolution and the relative low cost of the equipment [1][2][3].Information provided by fNIRS is complementary to that provided by direct measures of electromagnetic effects of neurons' activation and, for this reason, fNIRS has recently been used for the simultaneous recording of hemodynamic and electromagnetic signals (i.e., electroencephalography-EEG and magnetoencephalography-MEG [4,5]; ).Despite its inherent high signal contrast, one limitation of fNIRS is the reduced spatial resolution with respect to fMRI, frequently solved via multiple fNIRS measurements in order to localize the signal source in the brain [6].
In functional brain imaging, the signal acquired with fNIRS is a mixture of evoked hemodynamic response function (HRF), several background physiological components (e.g., cardiac, respiratory, blood pressure Mayer wave) and measurement noise (artifacts).Several methods have been proposed in the literature to estimate HRF from fNIRS signal by reducing the effect of noisy components, like low-pass filtering [7], state space estimation using extended Kalman filter [8][9][10][11], principal component analysis [12], generalized linear model [13][14][15], wavelet-based algorithm [16] and independent component analysis [17][18][19].Although each of these methods is generally associated with increases in signal-to-noise ratio, the so-called conventional averaging (CA) technique is still probably the most widely used method to estimate HRF from the fNIRS signal [3,5,7,[20][21][22][23]. Succinctly, the HRF is determined by averaging the fNIRS recordings (trials) collected after N identical stimuli, with N being often in the order of several tenth.Estimation of the HRF is achieved by assuming both the independence of the background noise from the activity elicited by the to-beprocessed stimulus, and the difference in phase of the physiological components from stimulus to stimulus.Notably, CA is algorithmically blind to information about HRF and fNIRS signals that can be independently extracted from the optical signal.Furthermore, fNIRS signal is notoriously not stationary (i.e., some trials are more noisy and/or less reliable than others), which suggests that the above assumptions should be taken with great caution.
In this work, we describe an approach for the estimation of HRF from fNIRS data devised to overcome the problems mentioned in the foregoing paragraph.The proposed method includes a pre-processing phase devoted to both artifacts removal (e.g., anomalous drifts due to movements of the subject) and reduction of the impact of cardiac activity on fNIRS signal.This pre-processing step allows us to circumvent the necessity of developing a complex model of the data (e.g., as the combination of sinusoids to describe the physiological noise component in the fNIRS signal) which would require sophisticated, but complex to deal with, estimation tools, such as the extended Kalman filter used in [11].Then, a Bayesian nonparametric approach is used to improve the contrast to noise ratio (CNR).The strength of this latter algorithm, originally developed by Sparacino et al. [24] for the estimation of eventrelated potentials (ERP) in electroencephalography, is that, though relying on mild assumptions on the fNIRS signal, it improves substantially the CNR due to the generation of a suitable compromise between experimental data and a priori expectations (e.g., smoothness) available on the unknown HRF.The proposed method is a general-purpose and straightforward filtering technique, that requires only a small amount of a priori information and can be used in any fNIRS experiment.In the present work, we compared its performance with two general, simple and widely used methods, namely, CA and classical Butterworth band-pass filtering.
Results obtained with both simulated and real data indicated that the methodology proposed in this context improves the HRF estimate compared to both CA and band-pass filtering.On the empirical side, the present approach was used to investigate visual short-term memory in humans (VSTM).Specifically, the hemodynamic signal was acquired from subjects engaged in the cued variant of a change-detection task.Acquired data were filtered using the Bayesian filtering approach and the peak amplitudes of the obtained HRFs were analyzed to disclose the neural underpinnings of VSTM functions.In this experiment, the amplitude of the expected HRFs was comparable or even lower to that of noise and of physiological components, making a stable HRF estimation hard to achieve.Nonetheless, in line with findings of particularly recent EEG [25] and fMRI studies [26,27], it was shown that the application of the present filter allowed us to detect and localize activity in parietooccipital regions (i.e., in the intraparietal sulcus, IPS, and the intra-occipital sulcus, IOS) correlated with the maintenance of information in VSTM.Furthermore, results obtained following the adoption of the present Bayesian filtering approach shed light and substantiated empirically an expected lateralization effect (see section 3.4 Lateralization effects in visual short-term memory), that dovetails nicely with prior fMRI investigations of VSTM functions in humans.

Participants
Thirteen right-handed students at the University of Padova participated in the experiment after providing informed consent.Two participants were discarded from the analysis because of an inadequate behavioural performance (accuracy below 50%).Thus, the participants included in the analysis were eleven (6 females, mean age 23.36 years ± 2.73).All participants had normal or corrected to-normal vision, and normal color vision.No participant reported a prior history of neurological or psychiatric disorders, and none was under medication at the time of testing.This study was conducted in accordance with the guidelines approved by the Institutional Review Board (IRB).

Stimuli and procedure
Each subject was seated in a comfortable chair placed before a computer screen.A schematic illustration of the sequence of events on each trial of the present experiment is reported in Fig. 1.A plus sign was presented at the centre of the screen for 2 s.The plus sign was then replaced by a small circle for 600 ms, alerting participants about the incoming presentation of the first (memory) array of colored squares.The alerting circle, which remained on the screen throughout the entire trial, was then flanked by an arrow head pointing to the left or right side of the screen for 400 ms.The offset of the arrow head was followed by a 300 ms blank interval.The memory-array was exposed for 300 ms, and followed by a blank interval lasting 1400-1600 ms, jittered in steps of 20 ms.A second (test) array composed of colored squares placed in the same positions as the squares in the memory-array was then displayed for a maximum duration of 3000 ms.Following the onset of the test-array, participants had to press one of two appropriately labeled keys of the computer keyboard to indicate whether a change in color had occurred or not.On 50% of trials one of the squares in the cued half of the memory-array was replaced with a square of a different color in the test-array.In the other 50% of trials, the color of the squares did not change between memory-and test-arrays.In the example illustrated, the rightmost square, green in the memory-array, was changed to cyan in the test-array.The to-be-remembered colored squares could be two or four.The inter-trial interval varied from 15 to 20 s.Trials were organized in 6 consecutive blocks (with a short pause between successive blocks), and each block included 24 trials.

Instruments
The fNIRS signal was acquired with a multi-channel frequency-domain NIR spectrometer (ISS Imagent TM , Champaign, Illinois), equipped with 28 laser diodes (14 emitting light at 690 nm, and 14 at 830 nm) modulated at 110.0 MHz.The diode-emitted light was conveyed to the subject's head by multimode core glass optical fibers (heretofore, sources; OFS Furukawa LOWOH series fibers, 0.37 of numerical aperture) with a length of 250 cm and a core diameter of 400 µm.Light that scattered through the brain tissue was carried by detector optical fiber bundles (diameter 3 mm) to 4 photo-multiplier tubes (PMTs; R928 Hamamatsu Photonics).The PMTs were modulated at 110.005 MHz, generating a 5.0 KHz heterodyning (cross-correlation) frequency.To separate the light as a function of source location, the sources time-shared the 4 parallel PMTs via an electronic multiplexing device.Only two sources (one per hemisphere) were synchronously (t = 4 ms) active (i.e., emitting light) resulting in a final sampling period of 128 ms (f = 10 3 /128 = 7.8125 Hz), after a dual-period averaging.
Following detection and consequent amplification by the PMTs, the optical signal was converted into temporal variations (Δ) of cerebral oxy-hemoglobin (ΔHbO) and deoxyhemoglobin (ΔHbR) concentration, an index that is sensitive to age [29,30].These indices were therefore corrected for age based on the differential-pathlength factor (DPF [28]; for details see [22]) using the equations described in [31] (1) Horizontal EOG (HEOG) was recorded bipolarly from tin electrodes positioned on the outer canthi of both eyes, referenced to the left earlobe.HEOG activity was amplified, filtered using a highpass of 0.01 Hz, and digitized at the same sampling rate of the fNIRS signal.
Impedance at each electrode was maintained below 5 KΩ.HEOG activity was re-referenced offline to the average of the left and right earlobes.

Probe placement procedure
Sources and detectors were held in place on the scalp using a custom-made holder and velcro straps.Each source location comprised two source optical fibers, one for each wavelength.The distance between each source/detector pair (heretofore, channel) was L = 30 mm, so as to equate channels for optical penetration depth into the cortical tissue [32].As shown by [33], the average cortical surface depth varies across regions, measuring around 17-18 mm in the parieto-occipital cortex.Considering that the average cortical thickness is around 30 mm, we set the cortical depth at 20 mm.Thus, the position of each channel coincided with the 20 mm cerebral projection of the midpoint of the source-detector distance.This probe arrangement included 18 channels, providing 18 measurements for HbO and 18 for HbR.
The spatial arrangement of source/detector pairs on the scalp was determined using a new probe placement approach [34], based on the combined use of a physical model of the head surface of the ICBM152 template (ICBM152-PM) and a 3D digitizing software (Brain Sight, RogueResearch).
The sources on each hemisphere were numbered from 1 to 7. Left-hemisphere detectors were indicated with the letters A/B, and right-hemisphere detectors with the letters C/D.Channels A1/C1 and A2/C2 recorded activity from the superior IPS (sIPS), A3/C3 from the angular gyrus (ANG), A4/C4 from the IPS, A5/C5 from the posterior part of the superior parietal lobule (pSPL), B3/D3 from the ANG, B4/D4 from the region at the intersection between IPS and the intra-occipital sulcus (IOS), B6/D6 from regions adjacent to the lateral occipital cortex (LOC), and B7/D7 from the superior occipital cortex (SOC).The resulting spatial arrangement of sources and detectors on the head surface is illustrated in Fig. 2a.The channels overlaid onto the ICBM152 brain template are illustrated in Fig. 2b.Both Figs.2a and 2b were generated after remapping the stereotaxic points onto the ICBM152 template using MRIcron (http://www.sph.sc.edu/comd/rorden/mricron/; for details, see [24]).
Afterwards, in order to place precisely the probe on the scalp of participants, the biunivocal correspondence between 10 and 10 points [35] and the MNI coordinates of their cerebral projections [33] was used.The spatial arrangement of each set of landmarks generated a stringent spatial bind that allowed us to standardize the probe placement across participants.The degree of precision achieved with this procedure was comparable with that obtained using different approaches [33,36], yielding a worst-case average error compatible with the spatial resolution of the present fNIRS setup.

The pre-processing strategy
Concentration measurements were first band-pass filtered (pass band: from 0.01 Hz to 3 Hz) to further remove any slowly drifting signal components and other noise with frequencies far from the signal band.For each channel, the acquired signal was segmented into 15 s trials starting from the memory array onset.Trials associated with an incorrect response and/or with HEOG amplitude exceeding ± 30 µV during the interval from the onset of the arrow cue presentation and the offset of the memory array (8%) were discarded from analysis.Trials were divided into two conditions of the present experimental design depending on the position (left hemisphere or right hemisphere) of the considered channel: to-be-memorized colored squared displayed on the same side of the considered channel (ipsilateral condition) and displayed on the opposite side of the considered channel (contralateral condition).Each trial was zero-mean corrected by subtracting the mean intensity of the optical signal recorded during the 15 s period.
Two separate procedures were applied to remove artifacts present in the hemodynamic signal.First, a custom procedure based on Grubbs' test was separately applied on the aligned trials in each condition [37].Trials with one or more values exceeding the empirically established Grubbs' threshold (which was set to 0.05) were discarded from analysis (1%).Moreover, given that some artifacts could have not been detected by this procedure, the remaining trials were checked with a second method, which considered variations in concentration of the hemodynamic signal throughout the entire trial rather than considering single time-points as the Grubbs' test.The mean value and the difference between the maximum and minimum values (heretofore, range) were calculated considering all trials in a given condition.The mean value and the range were also calculated for each single trial.Single-trial mean and range values were compared with the mean values of all trials in that condition [38].Trials characterized by a range or mean value greater than the condition mean ± 2.5 SDs were discarded from analysis (2%).
Subsequently, a notch filter was applied to reduce the cardiac component of the signal, which was the cyclic physiological component that most affected the performance of the Bayesian filtering described below (specifically, since the cardiac component had a period of about 1 s, it influenced the estimation of the noise model).The centre frequency of the notch filter was set to the frequency corresponding to the maximum value of the power spectral density in the range 0.7 -1.5 Hz, and it was computed for each trial in order to take into account possible variations of heart rate during the experiment.

The Bayesian filtering approach
The complete details of the approach, at both the theoretical and implementative levels, are described in [20], where the method was originally proposed and applied to three different kind of auditory ERPs.Only a brief explanation is thus reported in the present manuscript.
Data of each of the N trials were described by an additive model like the following: where y, u and v were n-size vectors containing equally spaced samples of the measured fNIRS signal (after the application of the above described pre-processing procedure), the "true" HRF u and the noise v (time zero is the time in which the stimulus occurs).Then, each trial y of Eq. ( 2) was individually filtered within a Bayesian embedding by exploiting a priori 2 nd order statistical information on both u and noise v, different from one trial to another.As far as the information on v, a stationary AR model was employed.Hence, the a priori covariance matrix of v was: where A was the square n-dimensional Toeplitz matrix the first column of which was [1, a 1 , a 2 ,…, a p , 0,… 0] T , {a k } k = 1,…,p being the coefficients of the AR model, and σ 2 was the variance of the noise process which derives the AR model.This model was identified, for each of the N available trials, from data measured in an interval lasting 4 s and starting from 1.5 s before the stimulus, when HRF was not present.Model order was set to 4. The only physiological component that occurred entirely in this interval of 4 s was the cardiac component, with an amplitude comparable to that of the expected HRF, which was removed by the notch filter previously described in the subsection devoted to pre-processing.As far as the 2 nd order statistical information on u, containing samples of the unknown HRF, the strategy was to model a priori known smoothness as the realization of a stochastic process obtained by the cascade of d integrators driven by a zero-mean white noise process {ε k } with variance λ 2 .Therefore, the covariance matrix of u was: where F = Δ d , with Δ being the square n-dimensional lower-triangular Toeplitz matrix the first column of which was [1, 1, 0, …, 0] T .For instance, for d = 1 the unknown signal is described by a random-walk model: which, in a Gaussian setting, tell us that, given u k-1 , then u k will be with probability 99,7% in the range u k-1 ± 3λ, i.e., the lower λ, the smoother {u k }.The multiple integration of a white noise process is a model widely used in the literature to describe signals on which only qualitative information is available.In fact, the model is flexible (i.e., it can be employed for a large class of responses by suitably tuning d) and simple enough (only λ is unknown) to be easily identifiable from post-stimulus data.The optimally filtered trial was thus: where γ = σ 2 /λ 2 was estimated, independently trial-by-trial, by the so-called "discrepancy" smoothing criterion (Twomey, 1965).
Remark 1.The discrepancy criterion may occasionally fail, leading to oversmoothed or undersmoothed profiles.A "mean" value of γ (obtained from individual γ of trials with no over-or under-smoothing) has been used only in cases (i.e., 5% of total trials) in which unacceptable smoothing was obtained.Without this correction, estimates of the HRF were worse (not shown).
Remark 2. The outcome of Eq. ( 6) is equivalent to that of the standard noncausal Kalman smoother, where the a priori model of the state evolution is u(t + 1) = u(t) + w(t) and the measurement model is y(t) = u(t) + v(t) (relative to Eq. ( 2) respectively), with covariance matrix of the process noise w equal to [λ 2 ], and covariance matrix of measurement noise v equal to [σ 2 ].The estimated HRF (ū) was obtained from the average of the N trials filtered as in Eq. ( 6), belonging to the same condition: and it was finally baseline-corrected by subtracting from the overall hemodynamic response the mean intensity of the signal in the 0-500 ms interval from the onset.A prototype of the whole algorithm described in this section was implemented in Matlab © (version R2008b, The Mathworks, Natick, Massachusetts, USA) and run on a personal computer.

Results
The new methodology is first tested and compared with CA and with a raw band-pass filtering procedure, on synthetically generated data in order to assess the results in a situation in which the true HRF is known.The estimation error will be considered as merit criterion.Then, the application to real data will be illustrated.In this case, the pre-processing procedure of the signal was the same for all methods and the contrast-to-noise ratio will be used as indicator of the estimator performance.The profiles obtained in the empirical study were then analyzed with the specific purpose to detect hemispheric lateralization effects during the maintenance of lateralized visual information in visual short-term memory.

Synthetic data generation
Simulated data were generated to assess the performance of the developed algorithm.It is well known that in fNIRS measurements background signals from systemic physiology noise are additional signal to the functional hemodynamic response [10,11].In order to take into account these effects in the generation of simulated data, these fluctuations were expressed as a linear combination of three sinusoids at the specific physiological frequencies.The first sinusoid corresponded to the cardiac component, with a frequency (f 1 ) ranging from 0.85 to 1.35 Hz and an amplitude (a 1 ) of ± 200 nM.The second sinusoid represented the respiratory component, whose frequency range was [0.15 -0.35] Hz (f 2 ) and the amplitude (a 2 ) was ± 200 nM.The third sinusoid described the blood pressure Mayer wave, with a frequency ranging from 0.05 to 0.1 Hz (f 3 ) and an amplitude (a 3 ) set to ± 400 nM.Each sinusoid had a different phase (θ), different trial by trial, ranging from 0 to 2π.
The HRF, function of time t, was modeled by a linear combination of two gamma-variant functions Γ, time dependent, with a total of 6 variable parameters [39]: , where u true is the known HRF, α tuned the amplitude, τ i and υ i tuned the shape and scale, respectively, and β determined the ratio of the response to undershoot.The parameters were chosen in order to have a peak amplitude of 140 nM and a peak latency equal to 6.5 s, corresponding to the peak amplitude and latency of the HRF that was expected by experimental data.The measurement noise η was modeled as a white normal process with standard deviation tuned to bear the standard deviation of real trials.
Only physiological and measurement noise (i.e., heart beat, respiration and blood pressure Mayer wave) were added to simulated data; artifacts (e.g., due to movements of the subject or shifts of a source or a detector, causing short, noncyclic abrupt drifts) have not been added.
Thus, samples y(t) of each simulated trial were generated as in Eq. ( 2), where u contains the samples of u true in Eq. ( 8) and v is given by: These simulated trials were comparable to real trials relative to HbO concentration changes.15 simulated subjects were created.Each subject contained 50 simulated trials.

Assessment of the method
The Bayesian filtering approach was applied to the simulated subjects, while the preprocessing procedure was not performed because no artifacts were introduced in simulated data.
In order to give a quantitative measure of the improvement of the estimate obtained with the Bayesian filtering approach, the following quantity was defined:  (10) where ū was the estimate of the HRF of Eq. ( 7) and u true was the HRF used in Eq. ( 8) to generate the simulated data.The value of E was a sort of percentage estimation error.The index E was computed for both CA and the new method in simulated data.Furthermore, the proposed method was compared with a more conventional band-pass filtering (Butterworth, band-pass, from 0.01 to 0.3 Hz).The main disadvantage of band-pass filtering is that it can reduce hemodynamic responses as well as noise because these components overlap in terms of frequency spectra [40].The index E was computed for this method too.Two representative subjects are shown in Figs.3a and 3b respectively, in which true HRF and the HRFs estimated with CA, Bayesian filtering and band-pass filtering are reported.In both subjects the HRF obtained with the Bayesian filter is the one that best correspond to the true HRF.All the values of E are reported in Table 1.A remarkable improvement of the estimation error was obtained by Bayesian filtering (E = 7.09 ± 7.33, mean ± SD) with respect to both CA (E = 33.28 ± 17.19) and band-pass filter (E = 18.37 ± 13.87), denoting the good HRF estimates achieved by the proposed method.The improvement obtained by Bayesian filtering was 26.19 with respect to CA and 11.28 with respect to band-pass filter.Notably, an even more pronounced relative improvement (33.34 with respect to CA and 17.65 with respect to band-pass filter) could be seen if the number N of trials was reduced of 10% (N = 45), with an estimation error equal to 9.21, 42.55 and 26.86, obtained by Bayesian filtering, CA and bandpass filter respectively.This, and the fact that the noise components had an amplitude greater than that of the expected HRF, suggest that the use of the method could allow a reduction of the number N of trials needed to obtain an interpretable HRF estimate, with a consequent reduction of the duration of the experiment.The estimation error of the peak amplitude estimate (E peak ) was computed by the HRFs obtained with Bayesian filter and band-pass filter, while it was not computed by the HRFs obtained with CA because of the great amount of noise still present in the signal.The peak amplitude was the parameter used in the statistical analysis, and its correct estimation is crucial.In this case also, an improvement of the estimation error was obtained by Bayesian filtering (E peak = 4.24 ± 5.95, mean ± SD) with respect to band-pass (E peak = 11.90 ± 10.39).The lower values of the estimation error of the peak amplitudes with respect to the estimation error of the whole HRFs were due to the fact that the differences between true HRF and the estimated HRFs were greater in the last 5 s of the considered interval (lasting 15 s), where any peak was present.2), while the red dashed line depicts the filtered profile û obtained by Eq. ( 6).The yellow dashed line illustrates band-pass filtered data.High frequency physiological components and random measurement noise are notably reduced by the Bayesian filtering, while low frequency physiological components (i.e., blood pressure Mayer wave) is reduced by averaging all trials.Figure 4b shows a direct comparison between the average HRF ū (red dashed line) estimated by Eq. ( 7) in subject 9 from the available N = 37 trials (channel A2) and the HRF profile estimated via CA (blue line) and band-pass filtering (yellow dashed line).Similarly, Fig. 5a and Fig. 5b show representative filtered trial and estimated HRF relative to ΔHbR instead of ΔHbO.It is apparent that the proposed method is able to reduce noise and make possible a reliable estimation of peak amplitude.On the contrary, the estimate of the peak amplitude is not possible using the CA technique, because of the great amount of noise still present in the signal.

Application to real data
To assess the benefits of the proposed method in real data, the contrast to noise ratio (CNR) [41] was computed.The CNR was defined as the square root of the ratio of signal power to noise power.The signal power was calculated by integrating the power from the power spectral density over the "signal bands", the noise power was calculated by integrating the rest of the power spectrum.This definition of CNR is useful in real data, where is not possible to totally separate signal and noise, and it is preferable to time domain methods because of the periodic nature of our visual stimulation.The signal bands of our stimulation paradigm were 0.05-0.067Hz (fundamental band) and 0.1-0.133Hz (second harmonic).The noise bands were the rest of the spectrum.For ΔHbO, the mean CNR in single trial without Bayesian filtering was 0.86 after Bayesian filtering the CNR increases to 1.79, denoting a CNR improvement of 108%.Once estimated the HRFs, the mean CNR obtained with CA and with Bayesian filter were 1.42 and 2.30 respectively, denoting a CNR improvement of 62%.Worst CNR values were obtained for ΔHbR: 0.70 in single trial and 1.03 in estimated HRFs without Bayesian filtering, which increased to 1.48 in single trial and 1.78 in estimated HRFs using Bayesian filtering, denoting an improvement of 111% and 72% respectively.As regard the comparison with the band-pass filter, the CNRs obtained with the band-pass filter were similar to that obtained with the Bayesian filtering approach (ΔHbO, single trial: 1.78 vs 1.79, HRF: 2.30 vs 2.30).
CNR values obtained by CA, Bayesian and band-pass filtering are reported in Table 2. Furthermore, analysis on single-trial estimates of Eq. ( 6) showed a reduction of the standard deviation obtained by the Bayesian filtering approach with respect to the band-pass filtering.For each subject, the standard deviation of the trials of each condition was computed.The mean SD obtained by the proposed method was 248 nM in the contralateral condition and 240 nM in the ipsilateral condition.Using the band-pass filter, the SD values obtained were 318 and 306 nM in the contralateral and ipsilateral condition respectively.Obtained SD values for each condition and subject were submitted to a mixed ANOVA considering contralateral and ipsilateral as within-subject factors and the type of filtering (Bayesian or Band-pass) as between-subjects factors.The ANOVA revealed a significant effect of the type of filtering (F(1, 10) = 22.702, p = 0.001), reflecting reduced SD in the Bayesian filtering relative to the band-pass filtering.The reduced variability in single trials obtained by the proposed method suggests a more valuable reduction of noise.

Lateralization effects in visual short-term memory
The obtained HRFs were analyzed to investigate the lateralization effects in visual short-term memory.
Lateralization effects are expected based on the known architectural properties of the neural pathways subtended in the encoding and maintenance of visual information.Succinctly, due to the crossing of the neural pathways and the level of the geniculate nucleus [42], enhanced activation of parieto-occipital cortex contralateral to the side of visual stimuli presentation should be observed relative to ipsilateral activation (e.g., [43,44,25]).While such lateralization effects are normally found using electroencephalography, they can hardly be detected using fMRI, which reveals a strong concurrent ipsilateral activation that tends to conceal them (e.g., [45]).Lateralization effects were monitored at each recording channel.For each symmetrical channels' pair, ΔHbO and ΔHbR concentration values contralateral to the cued hemifield were calculated by averaging data recorded at left-sided channels when the to-be-memorized colored squares were displayed in the right visual hemifield and data recorded at right-sided channels when the to-be-memorized colored squares were displayed in the left visual hemifield.Ipsilateral ΔHbO and ΔHbR concentration indices were calculated with an analogous algorithm by averaging data at the complementary channels.The concentration of ΔHbT (ΔHbO + ΔHbR [46]; ) was calculated as an estimate of cerebral blood volume.For each condition, the mean value in the interval between 1 s before and 1 s after the maximum value of the hemodynamic response was considered, and a one tail t-test was performed to identify the channels showing a significant activation increase relative to the baseline.A second series of one-tail paired t-tests was conducted to compare contralateral and ipsilateral condition.The results of all statistical tests conducted on the concentration values were corrected for multiple comparisons using the false discovery rate method (FDR BH [47]; ).The q value specifying the maximum FDR was set to 0.1, such that no more than 10% false positive could be included, on average, in the set of significantly active channels submitted to statistical test.For ΔHbO and ΔHbT, all channels resulted activated in both conditions, showing that the whole parieto-occipital regions investigated were involved in VSTM.Only the channel at B4/D4 (corresponding to the IPS/IOS) showed a significant difference between the two conditions (ΔHbO: p = 0.0093; ΔHbT: p = 0.0084).The relative activation maps are graphically reproduced in Fig. 6, showing the lateralization effect for ΔHbO.A greater concentration values was found in the contralateral condition relative to values in the ipsilateral condition (mean ΔHbO value in nM ± SD of contralateral vs. ipsilateral trials: 160.87 ± 89.88 vs. 136.21± 86.47; mean ΔHbT value in nM ± SD of contralateral vs. ipsilateral trials: 119.74 ± 98.67 vs. 103.41± 100.1).For ΔHbR, not significant results were found by the statistical analysis, probably due to the poor CNR value.Figure 7 shows the mean response profile for ΔHbO-ΔHbR and ΔHbT in the IPS/IOS region in the contralateral and ipsilateral conditions.
These results suggest that the increase in BOLD response associated with increases in VSTM load (found in recent fMRI studies) was larger in contralateral cortex than the ipsilateral cortex.Such results are also in agreement with recent EEG findings.Equivalent analyses were carried out on band-pass filtered data.Such analyses revealed a less significant activation versus baseline for both contralateral and ipsilateral condition in all active channels.The p-value obtained in the contralateral condition are reported in Table 3, and analogous results were obtained in the ipsilateral condition.A less significant difference between the two experimental conditions was found analyzing low-pass filtered data recorded in the IPS/IOS region (the corresponding p-values are reported in Table 4).Data obtained with Bayesian filtering reveal a more significant activation versus baseline and a more significant difference between the two conditions.The analysis on ΔHbR (not reported) didn't reveal significant activation.
Results obtained with the band-pass filter suggests a reduction of both noise and HRF.The results obtained in simulated and real data reassert the capacity of our method to reduce noise preserving HRF.

Discussion
The fNIRS is an emerging neuroimaging method which can be usefully employed to provide, with limited invasivity and reasonable laboratory costs, crucial information for the study of cognitive processes.Here, the possibility of estimating the HRF from fNIRS signal was under investigation.Several methods, with different degree of sophistication and adaptability to general situations, were proposed in the literature, but CA is still a widely used method.In this work, after having developed a pre-processing procedure, we assessed the performance of a non-parametric Bayesian approach, originally developed and assessed in [24] for a broader class of auditory ERP.
A key feature of the method is that, under mild assumptions on the signals into play, it can significantly improve the accuracy of the estimates thanks to its ability to establish, for each individual trial, a suitable compromise between data and a priori expectations available on HRF smoothness.In particular, the proposed Bayesian filtering approach exploits models of the 2 nd order a priori statistical information on the background fNIRS noise and on the unknown HRF.While such a statistical description of the ongoing fNIRS noise is obtained, trial by trial, by fitting an auto-regressive model against pre-stimulus data, the a priori known smoothness of the unknown HRF is formalized by describing it as the multiple integration of a white noise process.This is a general and flexible way to give an a priori probabilistic description of a physiological signal, and it can be employ for a large class of fNIRS experiments, especially because the model has only one unknown parameter which, for each trial, can be estimated by a well established smoothing criterion.Good results are achieved in signal-to-noise ratio improvement and physiological noise reduction.A reliable estimation of the HRFs can be reached even if the number of trial is small (N50).Synthetic data clearly demonstrated the superiority of the approach over CA and other filtering methods of frequent use (such as the Butterworth band-pass filtering used in this paper).The trial-by-trial filtering strategy and, in particular, the estimation of the noise model from the first seconds of each trial is especially suited when the acquisition of the fNIRS signal requires a long time (tens of minutes), like in the proposed experiment.As a matter of fact, acquisition is structurally proned to be influenced by a number of factors, e.g., caused by movements of the subject, resulting in a general non-stationarity of the signal.Some of the more likely factors include changes in the optical coupling between the subject's head and the fNIRS sensors (the detectors are sensitive to the angle at which the remitted light is received and the reflectance of the skin surface depends on the angle of incidence) and changes in the physiological components (e.g., changes in heart frequency and in the amplitude of its component in the hemodynamic signal).Hence, a filtering approach that considers the non-stationarity of both the measurement noise and the physiological components of the fNIRS signal is fundamental.On the other hand, the method returns an average response of the HRF via Eq.( 7).Describing the HRF by an average profile is widely done in fNIRS signal processing [3,5,7,11,[20][21][22][23].As a matter of fact, fMRI investigations [e.g., 48] have shown that the hemodynamic response is quite consistent within subjects, although there might be some variability between sessions.Given that subjects performed the present experiment in a single uninterrupted session, providing an average HRF profile is acceptable.However, further development of the present work will include the design of a single-trial analysis technique, which will allow us to quantify HRF temporal variability.
The Bayesian filtering method is already published in [24].However, its application to fNIRS signal is original and not straightforward.The resulting overall method (including the pre-processing step) is simple and general.It is based on mild assumptions on the acquired signals, and it can be used in any fNIRS experiment.These features make it more easy-to-use and flexible than more sophisticated, but more demanding in terms of hypotheses, methods available in the literature.For instance, unlike PCA, the new method does not depend from the number of channels and from their location, and it doesn't assume the space-time separability of physiological components from HRF: these requirements are necessary for PCA, but they could make the analysis not robust [12,49].The only assumption the present method hinges on is the possibility to describe a-priori expected smoothness of the unknown HRF based on a well established a priori model (multiple integration of a white noise process), that makes the derivation of the single unknown parameter (i.e., λ 2 ) from the data almost immediate, by means of an automatic smoothing criterion.To highlight the opposite case, the method proposed in [8,11] is twined with a quite more complex a priori model of the data distribution, where, for instance, the unknown HRF is described by a waveform whose functional properties must all be specified a priori, with the exception of a subset of parameters which are to be estimated via the additional application of a nonlinear Kalman filter.This latter, per se, requires several parameters, e.g., initial states and covariance matrices to be empirically specified.To note, the HRF could not be modeled as in the aforementioned papers, for peak amplitude and latency were unknown.
The Bayesian approach was compared to two current standard methods (CA and classical band-pass filtering), demonstrating a sizable improvement on HRF estimation.On simulated data, a five times lower than that obtained by band-pass filtering.On real data, the improvement achieved by the present method is confirmed by an increase of 72% in the contrast to noise ratio (CNR) with respect to CA.It is worth noting that in real data we cannot totally separate signal and noise.Consequently, in the signal band used to compute the CNR may be present not only HRF, but also noise.This is probably the reason why the CNR obtained by the Bayesian approach is similar to that obtained by band-pass filtering.A more valuable reduction of noise obtained by the proposed method is suggested by the reduced standard deviation (SD) in single trial: a statistical test (ANOVA) revealed a significant effect of the type of filtering, reflecting reduced SD in the Bayesian filtering relative to the bandpass filtering.In addition, the more correct HRF estimate obtained by the proposed method with respect to that obtained by the band-pass filtering is suggested by the more significant difference between the HRFs of the two experimental conditions, that is confirmed by recent findings in fMRI studies.
The results obtained through the developed methodology were shown to be useful to provide new insights in the investigation of the VSTM.The peak amplitude of the HRFs of each channel and experimental condition can be reasonably estimated and analyzed.Results obtained with the Bayesian filtering approach denote a lateralization effect, that is hardly observable using fMRI but that is confirmed by EEG studies.Other experiments will be conducted using fNIRS and the developed filtering procedure to better understand VSTM mechanisms.
Further development of the present work will address the removal of the physiological components whose frequencies overlap the signal band, e.g., respiratory and, above all, blood pressure (Mayer) wave.

Fig. 1 .
Fig. 1.Sequence of visual events on each trial of the present experiment (see text for a full description).

Fig. 2 .
Fig. 2. Probe placement on the ICBM152 template (occipital view).(a) Sources (red circles) and detectors (black circles) overlaid on the head surface of the ICBM152-PM template; (b) Cerebral projections of sources (white circles) and detectors (black circles) overlaid on the ICBM152 brain template.

Figure
Figure 4adisplays an example of filtering for a representative trial of a representative subject (subject 1, channel A4, trial 5).The blue line denotes the raw signal of y in Eq. (2), while the red dashed line depicts the filtered profile û obtained by Eq. (6).The yellow dashed line illustrates band-pass filtered data.High frequency physiological components and random measurement noise are notably reduced by the Bayesian filtering, while low frequency physiological components (i.e., blood pressure Mayer wave) is reduced by averaging all trials.Figure4bshows a direct comparison between the average HRF ū (red dashed line) estimated by Eq. (7) in subject 9 from the available N = 37 trials (channel A2) and the HRF profile estimated via CA (blue line) and band-pass filtering (yellow dashed line).

Fig. 6 .
Fig.6.Occipital and top views of the statistical maps of contralateral vs. ipsilateral comparison, for ΔHbO.The maps have been overlaid onto the ICBM152 brain template.For illustrative purposes, the upper part of the figure shows examples of memory arrays with the hemifields including the to-be-memorized colored squares (contralateral to the hemisphere where the enhanced concentration of HbT and HbO was observed in the present study) cued by arrow heads.Note however that colored squares and arrow heads were never displayed synchronously during the experiment (see Fig.1).