|
1.IntroductionFunctional near infrared spectroscopy (fNIRS)1–3 is a promising imaging technology for noninvasive, continuous monitoring of regional blood flow and tissue oxygen consumption. It utilizes low-intensity light at the near-infrared spectrum (between 600 and 900 nm), within which the optical absorbance of living tissues is small and hence, light can penetrate up to a few centimeters into the tissue. The two forms of the oxygen-carrying molecule of the blood, denoted as oxyhemoglobin () (when oxygen is bound) and deoxyhemoglobin (Hb) (when oxygen is delivered), have distinct spectroscopic characteristics within the near infrared (NIR) spectrum. Changes in optical absorption parameters are measured at two different wavelengths within the NIR spectrum and then, are converted to changes in Hb and concentrations using the Beer–Lambert law. The depth of NIR light penetration is proportional to the square root of the distance between the light source and the photodetector.4 With multiple source-detector separations (S-D), fNIRS enables concurrent monitoring of oxygenation changes at superficial and deep tissues across different body parts, such as hands, feet, and head. In a head fNIRS measurement with a typical S-D of 3 cm, NIR light penetrates sufficiently to reach the uttermost layer of the brain, the cortex. fNIRS is increasingly being used for detecting task-related oxygen consumption changes in the cortex. Mainly because of the portability and low cost of the measurements, fNIRS has become considerably popular for research in neuroscience, sports medicine, psychology, psychiatry, and rehabilitation. Over the past two decades, it has been applied to several clinical conditions too, such as schizophrenia,5–12 Alzheimer’s disease,13,14 epilepsy,15–18 and neonatal intensive care.19–27 The measured hemodynamic response by fNIRS is a realization of a naturally continuous physiological phenomenon. Typical fNIRS data consist of a set of discrete time points sampled every fraction of a second. In the presence of observational noise, the measured hemodynamic time series are often not smooth and fluctuating but one can assume that the true underlying trajectory is a smooth function. Moreover, due to the slow nature of the hemodynamic response and given the typical sampling rate of an fNIRS device, it is likely that the adjacent samples are correlated to some extent. It is, therefore, desirable to view an individual’s repeated measures over a time interval as a continuous function or curve. Traditionally, fNIRS studies involve linear regression and analysis of variance on features extracted from the recovered hemodynamic response. The time series of the evoked hemodynamic response can be directly estimated by linear deconvolution of the measured hemodynamic changes and the experimental timing function (a boxcar function that determines the timing of the experimental paradigm). An accurate estimate of the evoked response, however, requires a precise choice of the timing function.28 Another approach for describing the waveform of the evoked hemodynamic response, which is adapted from functional magnetic resonance imaging (fMRI) analysis, assumes a temporal shape for the so-called hemodynamic response function (HRF).29 Given this a priori assumption, it models the HRF using a set of basis functions (such as gamma or Gaussian functions).30 This technique has the advantage of reducing the number of unknowns by returning a smooth function with definite parameters. However, results may be biased by the choice of the canonical basis functions. The estimated evoked hemodynamic responses are averaged to enhance signal-to-noise ratio, provided the interstimulus intervals (ISIs) are sufficiently large to avoid overlap of HRFs.28 Analysis of fNIRS data has primarily been limited to discrete-time methods. Functional data analysis (fDA)31 is an alternative technique developed for handling a large number of data sampled over a continuum which is often time. In a functional domain, we study functional objects rather than sample points; therefore, the vector observations , where is the number of observations, are replaced by functions . The philosophy behind fDA is “to think of observed data functions as single entities, rather than merely as a sequence of individual observations.”31 Discrete observations can be converted to continuous functions using a linear combination of basis functions. The advantages of using such representation is twofold: (1) it provides a computational platform that allows storing and analyzing large amount of data with improved computational efficiency and flexibility as computational challenges arise due to huge number of measurements on each of small number of subjects, known as large small problem32 and (2) smooth functions allow study of the dynamics of the underlying processes through their derivatives. fDA has been applied to neuroimaging studies.33 Viviani et al. first proposed an fDA approach for exploratory analysis of fMRI images for a single subject case.34 They showed that compared to ordinary principal component analysis (PCA), the functional version of PCA could better visualize the variability of the data introduced by experimental alternations as it takes advantage of smooth functions. Later, Long et al. followed an fDA approach for dimension reduction of fMRI data for the estimation of noise covariance kernel.35 Our primary aim in this research was to introduce fDA methodology for “exploratory analysis” of fNIRS data to obtain a better insight into the physiological response. Given a set of Hb and functions collected from a number of individuals, we wish to characterize main directions of variability and investigate shared features by the two functions. We specifically would like to explore: (1) temporal variability within observations and (2) linear relationship and temporal association between Hb and time series. We employed two novel functional techniques: functional principal component analysis (fPCA) and functional canonical correlation analysis (fCCA).31 fPCA explores the variability between a set of observations. In a functional sense, it reveals when, in a time series, the maximum variation between several observations occurs. fPCA works without a priori knowledge of the experimental setting or any specific assumptions on the data. In other words, fPCA provides an initial assessment and helps understand variation and population structure in high-dimensional data.34 We hypothesized that fPCA can discriminate events in a set of fNIRS data during a functional task. fCCA examines modes of variation that two sets of functions share. With this method, we would like to investigate the interaction between Hb and time series and discover how these two functions covary over time. In other words, we want to know, during the course of an experiment, how variability in Hb data is related to the variability in data and what types of and how much variation they share. In this study, we use Hb and parameters measured during a painful task, namely cold pressor test (CPT). A CPT is performed by immersing a limb into a cold water (usually freezing water) container for a specific period of time. Three trials of a CPT were used to study the effect of repeated sustained noxious stimuli on the hemodynamics, and the results are published elsewhere.36 The rest of the paper is organized as follows: in Sec. 2, we will briefly describe the protocol and measurements and will present a description of the fDA framework; in Sec. 3, we will illustrate the application of fPCA and fCCA on fNIRS data collected during a CPT; lastly, in Sec. 4, we will review and discuss the results. 2.Materials and Methods2.1.ParticipantsTwenty healthy, right-handed individuals from the Drexel University community participated in this study after giving the informed consent form approved by the Institutional Review Board (IRB). We instructed subjects to avoid smoking and drinking any caffeinated or alcoholic beverages for at least 3 h prior to the experiment. 2.2.ProtocolSubjects performed three successive trials of a CPT. Only data from the first CPT were used for the present study. The experiment block consisted of four conditions: (1) a 30-s baseline when the subject relaxed; (2) a 2-min immersion of the right hand up to the wrist into a bucket of circulating tepid water () for adaptation; (3) a 45-s immersion of the same hand into a bucket of circulating ice water (); and (4) a 2-min poststimulus hand immersion in the tepid water for the hemodynamic recovery. 2.3.MeasurementsWe used a continuous wave fNIRS system designed and developed at Drexel University. The principle and instrumentation of fNIRS are described elsewhere.36 The fNIRS sensors consisted of one light source with two light-emitting diodes (LEDs) at 730- and 850-nm wavelengths and three photodetectors. Two detectors were placed at 2.8 cm from the LED making the “far” channels to investigate the hemodynamic response at intracranial layers and one detector was located at 1 cm from the LED making the “near” channel to measure the hemodynamic changes within the superficial extracranial tissues. Two fNIRS sensors with the same configuration were positioned symmetrically on the left and right sides of a subject’s forehead proximate to the anterior median line and were secured using a Velcro strap. The sampling rate of raw optical intensity measurements was 2 Hz. The optical density parameters for 730- and 850-nm wavelengths were calculated by taking the logarithm of the ratio of the detected light intensity during baseline to the detected light intensity during the task. The optical density time series were converted to changes in oxyhemoglobin () and deoxyhemoglobin (Hb) concentrations using the modified Beer–Lambert law.37 2.4.Functional Data Analysis2.4.1.Building continuous, smooth curves from discrete dataModern signal or image acquisition machines sample data points from a smooth, continuous process subject to observational noise. This can be expressed as where is the observed data, is the assumed noise-free latent data at time , and is the observational error that explains the roughness of the observed data.Prior to any functional analysis, the discrete data need to be converted to a continuous functional object. Then, a flexible method that gives control over the level of smoothing can be applied. Spline smoothing is the most popular and powerful approach with well-defined properties.38–40 Spline smoothing or regularization approach controls the degree of smoothness by a single parameter that penalizes the function’s roughness. This parameter, denoted as the smoothing parameter (), compromises between the goodness of fit and the smoothness. To formulate this method, we first need to define “goodness of fit” and “roughness” criteria. The goodness of fit is measured by the least squares criterion where is the fitted data and is the observed data.Roughness of a function can be quantified by the total curvature criterion where denotes the second derivative operator. A smaller implies a less variable function, whereas a larger indicates a rougher curve. The roughness criterion can be generalized to any arbitrary order of a function’s derivative.A penalized residual sum of squares is formed by putting the two opposing criteria in Eqs. (2) and (3) together where is the smoothing parameter that controls the weight of data fit and smoothness.As mentioned earlier, the smoothing parameter () controls the trade-off between the closeness to the observed values and the smoothness. The amount of imposed smoothness, however, cannot be arbitrary. If is close to zero, we obtain an estimate close to the data, and if is too large, we obtain an estimate equivalent to the linear regression estimate of the data and consequently, the details of the signal of interest may be diminished. Therefore, it is important to select a reasonably good smoothing parameter. An appropriate smoothing parameter may be chosen subjectively by visual judgment and prior knowledge of the process generating the data. An objective, data-driven method is also developed using the generalized cross validation (GCV) measure41 where is the number of observations, SSE is the mean square error, and is the trace of the “smoothing matrix.”42 The optimum minimizes function plotted against . GCV is the most popular procedure for selecting an optimal value for the smoothing parameter. However, like other cross-validation methods, GCV tends to under-smooth the data.2.4.2.Functional PCAPCA explores main modes of variation in high-dimensional data on a set of orthogonal, linearly uncorrelated directions. The functional counterpart of PCA (fPCA) visualizes variability among a set of functional observations, which theoretically have infinite number of data values. fPCA returns a set of orthogonal weight functions or principal components (PCs) spanning the same range as the functional data. Each PC accounts for a percentage of the total variability in data. Depending on how much variability one would like to explain, the first few PCs are retained to estimate the data. fPCA algorithmLike ordinary PCA, functions are first centered, i.e., the mean across all observations is removed from each observation, to avoid the first PC to mostly reflect the average response. Basically, fPCA of observations , finds the PC weight function for which the PC score defined as maximizes constrained toAfter computing the first PC weight function and the first PC score according to Eqs. (6) and (7), the second PC weight function and the second PC score are calculated with an additional constraint which requires the new PC weight function to be orthogonal to that computed on the previous step. This constraint requires the new PCs to reflect new directions of variability. After each step, the amount of explained variability decreases. A sequence of descending PC scores and the corresponding weight functions are computed stepwise. By definition, and are referred to as the “eigenvalues” and “eigenfunctions.”To compute the functional PCs, one can formulate the fPCA as the eigenanalysis of the bivariate covariance function: The functional eigenequation is then expressed as where is the ’th eigenvalue and is the ’th eigenfunction of the covariance function. Computer software is developed to solve the eigenequation [Eq. (10)] for pairs of and .Only the first few PCs are considered for further analysis and the remaining PCs are regarded as either irrelevant or are assumed to reflect observational error or noise. In this case, the problem is to figure out how many components are to be considered. A solution is to plot the base 10 logarithm of eigenvalues according to their size (against their indices ) and to identify the so-called elbow point at which the slope of the graph changes from steep to flat. Then, the researcher may decide to retain only the components before the elbow. This method is called the scree or elbow test and is yet considered subjective. A traditional approach includes the components with eigenvalues greater than average. Another method is to include as many as components that account for a specific cumulative percentage of total variability. 2.4.3.Canonical correlation analysisFunctional canonical correlation analysis (fCCA) is an exploratory technique that examines modes of variation that “pairs” of functions share.43 By fCCA, we would like to explore how variability in Hb and data is interrelated. Given pairs of observations , , fCCA estimates pairs of functions that most explain the interaction between and . As in fPCA, data are first centered to focus on the analysis of variation from the mean.43 and are termed “canonical weight functions” and are estimated such that the correlation between the following integrals is maximized are called “canonical variates” and are calculated by maximizing the squared canonical correlation defined as Similar to fPCA, one can obtain a nonincreasing series of squared canonical correlations by imposing the condition that successive canonical weight functions be orthogonal.42 However, imposing a strong roughness penalty on the computed weight functions and is a rule to obtain meaningful results. A cross-validation criterion for objective calculation of the roughness penalty is described by Leurgans et al.43 as follows: The roughness penalty for the smoothed weight functions is obtained by finding a value for the smoothing parameter () that maximizes the squared correlation of pairs of scores , where and are the smoothed canonical weight functions with the ’th curve omitted. 3.Results3.1.fNIR Data PreprocessingRaw optical intensity data for 730- and 850-nm wavelengths were first filtered using a finite impulse response lowpass filter with a cut-off frequency of 0.14 Hz to remove high-frequency noise, respiration, and heart pulsation artifacts. The filtered raw data were converted to changes in Hb and concentrations relative to the mean value of the optical intensity during the first 20 s of the baseline. Hb and data were registered in reference to the time when the hand was immersed in the ice water. All functional data analyses were conducted in MATLAB (R2011a, MathWorks, Natwick, Massachusetts) using the fDA package for MATLAB.42 Hb and were smoothed by imposing a penalty on the roughness of the second derivative of the data with a of . The smoothing parameter was chosen empirically since the GCV plot resembled a sigmoid function which was not informative. Figures 1 and 2 show the nonsmooth and smooth Hb and data for two far channels located on the right and left sides of forehead, respectively. 3.2.Results of fPCAThe overall mean function of Hb and across all subjects was first subtracted from each observation before performing the fPCA. Results of fPCA of Hb and data were obtained separately for the six channels. No significant variation between channels was observed in the calculated PCs. Thus, in favor of saving space, we show and discuss the figures for two left and right far channels shown in Figs. 1 and 2. However, overall results for other channels will be discussed in Sec. 4. We estimated the number of components based on the total explained variance. We desired to account for 95% of total variability; therefore, the first four PCs were considered. The four considered PC weight functions of are presented in Figs. 3(a) and 3(b) for the right and left far channels, respectively. The first PC always has the highest PC score and accounts for the highest variability; the remaining PCs are associated to descending PC scores and account for less amount of total variability. A visualization technique that may simplify interpretation of fPCA results is to plot the PCs as perturbations of the mean function.31 In doing so, a multiple of each PC is added and subtracted to the mean function. By definition, this factor is 0.2 times the root mean square of the difference between the PC and its overall average.31 Results for PC perturbations of are shown in Figs. 3(c) and 3(d) for the right and left far channels, respectively. PC1 is mostly associated with amplitude variation, particularly during the stimulus and poststimulus conditions. PC2 accounts for a considerably smaller percentage of total variability and corresponds to amplitude variation mostly during prestimulus condition. PC3 and PC4 account for a very small percent of total variability and reflect some small variations during the prestimulus and the last minute of poststimulus conditions. Figure 4 shows the same information as Fig. 3 for Hb data. One may identify some similarities between the PCs calculated for Hb and data. For both right and left channels, PC1 is associated with an overall amplitude variation after the baseline. PC2 and PC3 explain amplitude variation during the prestimulus and poststimulus and a phase shift in early poststimulus condition. PC4 corresponds to some already identified variations mostly during prestimulus condition. 3.2.1.Rotating PCsAs we saw in the previous section, each considered PCs did not give substantive “localized” information of the variability; that is, several PCs may explain variation during an experimental condition or a single PC may correspond to variation during several conditions. Thus, the interpretation of PCs is not straightforward. Although the original PC weight functions are computed to be optimal, it does not mean that there are no other alternatives that can do just as well. One solution is to rotate the original PC weight functions by a rotation matrix. Varimax is the most popular rotation method for both multivariate and functional analyses.31 The rotated PC weight functions by varimax remain orthogonal as the original PCs but they are not necessarily uncorrelated. Moreover, the rotated PC weight functions are not associated to descending PC scores, unlike the original PCs. However, the cumulative variability explained by the rotated PCs is the same as before the rotation. A varimax approach finds an orthogonal rotation matrix () to transform the PC weight functions: where is a matrix containing the considered PC weight functions. A varimax solution is obtained by maximizing the variance of the vector of ; where are the elements of . This is only achievable when the are relatively large or relatively small. As such, the rotated weight functions have maximized variance due to few large elements and many minuscule elements. Consequently, each PC is associated with one or a small number of modes of variation, which considerably enhances interpretability.Results of the varimax rotation of and Hb data are shown in Figs. 5 and 6, respectively. The weight functions rotated by the VARIMAX rotation still account for the same cumulative percentage of total variability but in different proportions. Each rotated PC weight function [Figs. 5(a), 5(b), 6(a), and 6(b)] shows variation during a brief time window compared with unrotated PC weight functions [Figs. 3(a), 3(b), 4(a), and 4(b)] that described variation over a long time span. This effect of rotation is more evident in the mean function perturbations [Figs. 5(c), 5(d), 6(c), and 6(d)]. For data (Fig. 5), the two most important PCs that account for the largest percentage of total variability identify the evoked response during the stimulus condition and the recovery during poststimulus condition. The PC weight function associated to the stimulus condition identifies slope and amplitude variation of the evoked response and a temporal shift in the peak of response across the 20 subjects. The PC weight function associated to the poststimulus condition corresponds to amplitude and timing variation of the hemodynamic recovery curve. The two least important PCs are associated to variations in the prestimulus and end of the poststimulus conditions. For Hb data (Fig. 6), the largest variability is observed during the poststimulus condition by PC2 and PC3 for both channels; one associated to variation in the hemodynamic settlement at the end of the poststimulus and the other one associated to variation of peak of the evoked response occurred with a delay after removal of the stimulus. The two least important PCs identified variation during the stimulus and prestimulus conditions. 3.3.Results of fCCAfCCA was applied to Hb and curves to explore their interaction and shared variation. We used the cross-validation criterion proposed by Leurgans et al. to find the optimum value of the smoothing parameter () for the canonical weight functions.43 Based on the cross-validation results and visual inspection of plotted canonical weight functions, a of for far channels and a of for near channels were selected. We present and discuss the first canonical weight functions, a.k.a. the leading weight functions, for two far channels and two near channels on the right and left sides of forehead (Fig. 7). Results for the other far channels were similar to the presented far channels. Figure 7 shows that the variation from the norm for Hb and is positively correlated throughout the experiment for all channels. For the right far channel, however, the dominant mode of covariation exists in the peak and timing of Hb and changes, with variation being smaller and much faster than Hb response. 4.DiscussionfDA characterizes the main features shared by a set of observations from a number of individuals. In this research, we proposed the fDA as a novel tool for the analysis of fNIRS data and showed that the fPCA is capable of decomposing the Hb and time series during a CPT task into the experimental conditions. Multivariate PCA has been previously used to identify the systemic interference and remove it from fNIRS measurements. Zhang et al.44 used spatial eigenvector-based analysis of baseline activity in diffuse optical imaging. They used dozens of optodes over a large area of the cortex and showed that baseline information can be used to estimate global physiological interference during a finger tapping task and a tactile stimulus. Virtanen et al. compared PCA and independent component analysis (ICA) performance in removing superficial physiological interference during hypocapnia and hypercapnia maneuvers.45 Instead of using a fine grid of source–detector pairs, they used two channels with short (1 cm) and far (3 cm) source-detector distances and showed that both PCA and ICA are capable of effectively removing the physiological activity from the fNIRS measurements. fNIRS data are naturally continuous time series that are observed in discrete time. The frequency response of the hemodynamics is relatively low compared to typical fNIRS devices sampling rate and the recorded sample points are likely to be correlated. Therefore, a more accurate and effective way to look at fNIRS data is to incorporate the information that is inherent in the time order and smoothness of fNIRS data. In traditional or standard PCA, the data are treated as vectors of discrete samples and permutation of components will not affect the analysis and thus, the time order of the physiological phenomena is irrelevant. Moreover, it is well-known that large number of measurements on each of a comparatively small number of subjects or experimental units, the so-called “large , small problem,” might lead to ill-posed problems32 and also, the estimation of maximum eigenvalue might be asymptotically biased.46 The fDA approach exploits smoothed curves and estimates a smoothed fixed covariance function to avoid the “curse of dimensionality.” The fixed covariance function solves the problem of large covariance matrices whose dimensions exceedingly increase with the sample size. Hence, fDA alternative solution for handling high-dimensional data is becoming very popular. Also, because the fDA methodology considers the entire curve as a single entity, there remains no issue regarding high correlation between repeated measurements. fPCA is an exploratory tool to capture the principal directions of variation in the population and dimension reduction. Our study provides a good understanding of the variation in the population structure of Hb and at different conditions of a CPT. Our observations are summarized as follows:
In a secondary exploration, we investigated the association between Hb and variability using fCCA. We were interested in the dominant modes of correlation between Hb and . Results demonstrated an overall positive correlation between Hb and variation through the CPT experiment. However, for a right far channel, there was an evident time shift in the Hb response relative to . In general, fCCA gives a broad view of the interaction between two variables. The identified patterns of covariation give important directions for further investigation of interaction between Hb and in a particular experimental setting. Feature selection for the preliminary study of CPT data published elsewhere was done empirically.36 By visual inspection of the Hb and data, we observed that the amplitude and slope of the evoked hemodynamic response were quite variable between subjects, and we chose two features based on these characteristics for statistical analyses.36 However, empirical search did not give any information regarding the power of these features for explaining variability without running comprehensive statistical analyses. However, fPCA aids a straight feature selection by characterizing the geometry of curves as a set of principal directions of variation. Exploratory functional analysis of CPT data gave us new insight into our dataset by: (1) explicit estimation of PCs explaining the highest variability in the dataset, (2) calculation of the percentage of total variability explained by each component, and (3) study the covariation of Hb and throughout a CPT experiment. We showed that fPCA could successfully discriminate different phases of the hemodynamic response during a CPT experiment and interestingly, it identified two components for the poststimulus condition: an early component associated to variation in the hemodynamic recovery and a late component associated to the hemodynamic settlement. More importantly, it described inter-subject variations in the dynamic of the response, such as slope of the evoked response and slope of the recovery curve. In addition to this qualitative information, fPCA provided quantitative measures of the importance and weight of each of these features based on the percentage of total variability explained by that feature; for instance, experimental manipulations that created the largest variation among the dataset were identified to be the stimulus and early poststimulus conditions. Extracting such detailed information from the raw data (Figs. 1 and 2) was impractical, if not impossible. For future studies, we will exploit these findings for further statistical analysis of CPT data as well as designing new protocols according to the extracted features (e.g., how much time would be necessary for the hemodynamic recovery/settlement after a cold stimulus). In the context of brain imaging studies, as Tian33 pointed out in her review article, fPCA can play a critical role in research involving a stimulus-response paradigm when the onset, duration, amplitude, and/or speed of the evoked brain activity is unknown or uncertain, such as in a decision-making process or an emotional reaction. As an exploratory technique, fPCA does not make any a priori assumption on the form and timing of the waves of neural activity and this key property makes fPCA application particularly relevant in neuroscience studies employing fNIRS as well as other neuroimaging tools such as electro-encephalography, fMRI, and magneto-encephalography.33 After exploratory steps, one can extract more pertinent information. For instance, fPCA is a good tool for computing proximities between curves in a reduced dimensional space. Using such a distance measure, one can investigate whether the curves are similar or can be split into several classes. Due to inherent interaction between the autonomic nervous system and central pain processing system, researchers try to assess the individuals’ perception of pain through physiological signs. Hemodynamic changes at the skin and major arteries have been investigated to study the sympathetic response to a noxious stimulus.47–49 In this exploratory study, we attempted to learn what feature of response to a cold noxious stimulus is most variable among control subjects. We showed that the PCs explaining the variability in the CPT data and their weight varied between Hb and parameters and channels. This information may be helpful to decide which parameter and/or channels are most associated with the hypothesized effect, for example, in brain mapping studies. With the CPT experiment, we previously showed that the amplitude and timing of the evoked response of was correlated with subjects’ report of the pain intensity.36 Through the results of this research, we also learned about importance of the dynamics of the evoked response to a CPT as revealed by fPCA. In another study, we used this information and considered the derivatives of curves during a pain tolerance test using CPT. Results confirmed that the amplitude and timing of the peak of the first and second derivative of the response right after hand immersion in cold water is a very robust feature. In future, we will use this knowledge to better understand the physiology of pain through explaining some of the observed variations by known factors such as age, gender, or information about whether the subjects are high/low tolerance to that specific stimulus. The observed modes of variation and covariation patterns can help improve design of fNIRS studies for pain research in particular and for other experiments in general. This study was preliminary but opens an interesting avenue for future research. Without a good understanding of population structure, it would not be feasible to formulate and study physiological phenomena and identifying variability is the first step in this direction. AcknowledgmentsAuthors would like to acknowledge Dr. Patricia A. Shewokis to suggest the functional data analysis for the cold pressor data. The first author would like to thank Dr. Shewokis for instructing her in the topic. ReferencesP. Rolfe,
“In vivo near-infrared spectroscopy,”
Annu. Rev. Biomed. Eng., 2 715
–754
(2000). http://dx.doi.org/10.1146/annurev.bioeng.2.1.715 ARBEF7 1523-9829 Google Scholar
A. VillringerB. Chance,
“Non-invasive optical spectroscopy and imaging of human brain function,”
Trends Neurosci., 20
(10), 435
–442
(1997). http://dx.doi.org/10.1016/S0166-2236(97)01132-6 TNSCDR 0166-2236 Google Scholar
M. FerrariV. Quaresima,
“A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,”
Neuroimage, 63
(2), 921
–935
(2012). http://dx.doi.org/10.1016/j.neuroimage.2012.03.049 NEIMEF 1053-8119 Google Scholar
M. S. Pattersonet al.,
“Absorption spectroscopy in tissue-simulating materials: a theoretical and experimental study of photon paths,”
Appl. Opt., 34
(1), 22
–30
(1995). http://dx.doi.org/10.1364/AO.34.000022 APOPAI 0003-6935 Google Scholar
K. Taniguchiet al.,
“Multi-channel near-infrared spectroscopy reveals reduced prefrontal activation in schizophrenia patients during performance of the kana Stroop task,”
J. Med. Invest., 59
(1–2), 45
–52
(2012). http://dx.doi.org/10.2152/jmi.59.45 JINVFI 1081-5589 Google Scholar
S. Shimoderaet al.,
“Mapping hypofrontality during letter fluency task in schizophrenia: a multi-channel near-infrared spectroscopy study,”
Schizophr. Res., 136
(1–3), 63
–69
(2012). http://dx.doi.org/10.1016/j.schres.2012.01.039 SCRSEH 0920-9964 Google Scholar
S. Koikeet al.,
“Association between severe dorsolateral prefrontal dysfunction during random number generation and earlier onset in schizophrenia,”
Clin. Neurophysiol., 122
(8), 1533
–1540
(2011). http://dx.doi.org/10.1016/j.clinph.2010.12.056 CNEUFU 1388-2457 Google Scholar
Y. Fujitaet al.,
“Asymmetric alternation of the hemodynamic response at the prefrontal cortex in patients with schizophrenia during electroconvulsive therapy: a near-infrared spectroscopy study,”
Brain Res., 1410 132
–140
(2011). http://dx.doi.org/10.1016/j.brainres.2011.06.052 BRREAP 1385-299X Google Scholar
M. Azechiet al.,
“Discriminant analysis in schizophrenia and healthy subjects using prefrontal activation during frontal lobe tasks: a near-infrared spectroscopy,”
Schizophr. Res., 117
(1), 52
–60
(2010). http://dx.doi.org/10.1016/j.schres.2009.10.003 SCRSEH 0920-9964 Google Scholar
K. Ikezawaet al.,
“Impaired regional hemodynamic response in schizophrenia during multiple prefrontal activation tasks: a two-channel near-infrared spectroscopy study,”
Schizophr. Res., 108
(1–3), 93
–103
(2009). http://dx.doi.org/10.1016/j.schres.2008.12.010 SCRSEH 0920-9964 Google Scholar
A. WatanabeT. Kato,
“Cerebrovascular response to cognitive tasks in patients with schizophrenia measured by near-infrared spectroscopy,”
Schizophr. Bull., 30
(2), 435
–444
(2004). http://dx.doi.org/10.1093/oxfordjournals.schbul.a007090 SCZBB3 0586-7614 Google Scholar
F. Okadaet al.,
“Impaired interhemispheric integration in brain oxygenation and hemodynamics in schizophrenia,”
Eur. Arch. Psychiatry Clin. Neurosci., 244
(1), 17
–25
(1994). http://dx.doi.org/10.1007/BF02279807 0940-1334 Google Scholar
A. J. Fallgatteret al.,
“Loss of functional hemispheric asymmetry in Alzheimer’s dementia assessed with near-infrared spectroscopy,”
Brain Res. Cogn. Brain Res., 6
(1), 67
–72
(1997). http://dx.doi.org/10.1016/S0926-6410(97)00016-5 CBRREZ 0926-6410 Google Scholar
C. Hocket al.,
“Near infrared spectroscopy in the diagnosis of Alzheimer’s disease,”
Ann. N. Y. Acad. Sci., 777 22
–29
(1996). http://dx.doi.org/10.1111/j.1749-6632.1996.tb34397.x ANYAA9 0077-8923 Google Scholar
E. WatanabeY. NagahoriY. Mayanagi,
“Focus diagnosis of epilepsy using near-infrared spectroscopy,”
Epilepsia, 43
(Suppl 9), 50
–55
(2002). http://dx.doi.org/10.1046/j.1528-1157.43.s.9.12.x EPILAK 0013-9580 Google Scholar
K. Haginoyaet al.,
“Ictal cerebral haemodynamics of childhood epilepsy measured with near-infrared spectrophotometry,”
Brain, 125 1960
–1971
(2002). http://dx.doi.org/10.1093/brain/awf213 BRAIAK 0006-8950 Google Scholar
E. Watanabeet al.,
“Noninvasive cerebral blood volume measurement during seizures using multichannel near infrared spectroscopic topography,”
J. Biomed. Opt., 5
(3), 287
–290
(2000). http://dx.doi.org/10.1117/1.429998 JBOPFO 1083-3668 Google Scholar
D. K. Sokolet al.,
“Near infrared spectroscopy (NIRS) distinguishes seizure types,”
Seizure, 9
(5), 323
–327
(2000). http://dx.doi.org/10.1053/seiz.2000.0406 SEIZE7 1059-1311 Google Scholar
R. Bokiniecet al.,
“Assessment of brain oxygenation in term and preterm neonates using near infrared spectroscopy,”
Adv. Med. Sci., 57
(2), 348
–355
(2012). http://dx.doi.org/10.2478/v10039-012-0050-6 AMSDCH 1896-1126 Google Scholar
M. Rangeret al.,
“Cerebral near-infrared spectroscopy as a measure of nociceptive evoked activity in critically ill infants,”
Pain Res. Manag., 16
(5), 331
–336
(2011). Google Scholar
P. Gilibertiet al.,
“Near infrared spectroscopy in newborns with surgical disease,”
J. Matern. Fetal Neonatal Med., 24
(Suppl 1), 56
–58
(2011). http://dx.doi.org/10.3109/14767058.2011.607673 JMNMAE 1476-4954 Google Scholar
G. E. Morenoet al.,
“Regional venous oxygen saturation versus mixed venous saturation after paediatric cardiac surgery,”
Acta Anaesthesiol. Scand., 57
(3), 373
–379
(2013). http://dx.doi.org/10.1111/aas.12016 AANEAB 0001-5172 Google Scholar
A. Milanet al.,
“Near-infrared spectroscopy measure of limb peripheral perfusion in neonatal arterial thromboembolic disease,”
Minerva Pediatr., 64
(6), 633
–639
(2012). MIPEA5 0026-4946 Google Scholar
R. M. Cerboet al.,
“Cerebral and somatic rSO2 in sick preterm infants,”
J. Matern. Fetal Neonatal Med., 25
(Suppl 4), 97
–100
(2012). JMNMAE 1476-4954 Google Scholar
G. GreisenT. LeungM. Wolf,
“Has the time come to use near-infrared spectroscopy as a routine clinical tool in preterm infants undergoing intensive care?,”
Philos. Trans. A Math. Phys. Eng. Sci., 369
(1955), 4440
–4451
(2011). http://dx.doi.org/10.1098/rsta.2011.0261 PTRMAD 1364-503X Google Scholar
W. BaertsP. M. LemmersF. van Bel,
“Cerebral oxygenation and oxygen extraction in the preterm infant during desaturation: effects of increasing FiO(2) to assist recovery,”
Neonatology, 99
(1), 65
–72
(2011). http://dx.doi.org/10.1159/000302717 NEONCC 1661-7800 Google Scholar
A. PetrovaR. Mehta,
“Near-infrared spectroscopy in the detection of regional tissue oxygenation during hypoxic events in preterm infants undergoing critical care,”
Pediatr. Crit. Care Med., 7
(5), 449
–454
(2006). http://dx.doi.org/10.1097/01.PCC.0000235248.70482.14 1529-7535 Google Scholar
T. J. Huppertet al.,
“HomER: a review of time-series analysis methods for near-infrared spectroscopy of the brain,”
Appl. Opt., 48
(10), D280
–298
(2009). http://dx.doi.org/10.1364/AO.48.00D280 APOPAI 0003-6935 Google Scholar
K. J. FristonP. JezzardR. Turner,
“Analysis of functional MRI time-series,”
Hum. Brain Mapp., 1 153
–171
(1994). http://dx.doi.org/10.1002/hbm.v1:2 HBRME7 1065-9471 Google Scholar
Y. ZhangD. H. BrooksD. A. Boas,
“A haemodynamic response function model in spatio-temporal diffuse optical tomography,”
Phys. Med. Biol., 50
(19), 4625
–4644
(2005). http://dx.doi.org/10.1088/0031-9155/50/19/014 PHMBA7 0031-9155 Google Scholar
J. O. RamsayB. W. Silverman, Functional Data Analysis, Springer, New York
(2005). Google Scholar
I. M. JohnstoneD. M. Titterington,
“Statistical challenges of high-dimensional data,”
Phil. Trans. R. Soc. A, 367
(1906), 4237
–4253
(2009). http://dx.doi.org/10.1098/rsta.2009.0159 PTRMAD 1364-503X Google Scholar
T. S. Tian,
“Functional data analysis in brain imaging studies,”
Front. Psychol., 1 1
–11
(2010). Google Scholar
R. VivianiG. GronM. Spitzer,
“Functional principal component analysis of fMRI data,”
Human Brain Mapp., 24
(2), 109
–129
(2005). http://dx.doi.org/10.1002/(ISSN)1097-0193 HBRME7 1065-9471 Google Scholar
C. J. Longet al.,
“Nonstationary noise estimation in functional MRI,”
Neuroimage, 28
(4), 890
–903
(2005). http://dx.doi.org/10.1016/j.neuroimage.2005.06.043 NEIMEF 1053-8119 Google Scholar
Z. Baratiet al.,
“Hemodynamic response to repeated noxious cold pressor tests measured by functional near infrared spectroscopy on forehead,”
Ann. Biomed. Eng., 41
(2), 223
–237
(2013). http://dx.doi.org/10.1007/s10439-012-0642-0 ABMECF 0090-6964 Google Scholar
B. Chanceet al.,
“A novel method for fast imaging of brain function, non-invasively, with light,”
Opt. Express, 2
(10), 411
–423
(1998). http://dx.doi.org/10.1364/OE.2.000411 OPEXFF 1094-4087 Google Scholar
G. WahbaS. Wold,
“A completely automatic french curve: fitting spline functions by cross validation,”
Commun. Stat., 4
(1), 1
–17
(1975). http://dx.doi.org/10.1080/03610917508548493 0361-0918 Google Scholar
B. W. Silverman,
“Some aspects of the spline smoothing approach to non-parametric regression curve fitting,”
J. R. Statist. Soc. B, 47
(1), 1
–52
(1985). Google Scholar
R. L. Eubank, Spline Smoothing and Nonparametric Regression, 1988). Google Scholar
P. CravenG. Wahba,
“Smoothing noisy data with spline functions,”
Numer. Math., 31
(4), 377
–403
(1979). NUMMA7 0029-599X Google Scholar
J. O. RamsayG. HookerS. Graves, Functional Data Analysis with R and MATLAB, Springer, New York
(2009). Google Scholar
S. E. LeurgansR. A. MoyeedB. W. Silverman,
“Canonical correlation analysis when the data are curves,”
J. R. Statist. Soc. B, 55
(3), 725
–740
(1993). Google Scholar
Y. Zhanget al.,
“Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,”
J. Biomed. Opt., 10
(1), 011014
(2005). http://dx.doi.org/10.1117/1.1852552 JBOPFO 1083-3668 Google Scholar
J. VirtanenT. NoponenP. Merilainen,
“Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals,”
J. Biomed. Opt., 14
(5), 054032
(2009). http://dx.doi.org/10.1117/1.3253323 JBOPFO 1083-3668 Google Scholar
P. HallH. G. MüllerJ. L. Wang,
“Properties of principal component methods for functional and longitudinal data analysis,”
Ann. Stat., 34
(3), 1493
–1517
(2006). http://dx.doi.org/10.1214/009053606000000272 ASTSC7 0090-5364 Google Scholar
M. MusicZ. FinderleK. Cankar,
“Cold perception and cutaneous microvascular response to local cooling at different cooling temperatures,”
Microvasc. Res., 81
(3), 319
–324
(2011). http://dx.doi.org/10.1016/j.mvr.2011.01.004 MIVRA6 0026-2862 Google Scholar
S. Roattaet al.,
“Effect of generalised sympathetic activation by cold pressor test on cerebral haemodynamics in healthy humans,”
J. Auton. Nerv. Syst., 71
(2–3), 159
–166
(1998). http://dx.doi.org/10.1016/S0165-1838(98)00075-7 JASYDS 0165-1838 Google Scholar
R. B. Paneraiet al.,
“Cerebral blood flow velocity response to induced and spontaneous sudden changes in arterial blood pressure,”
Am. J. Physiol. Heart Circ. Physiol., 280
(5), H2162
–H2174
(2001). 0363-6135 Google Scholar
|