Analyzing the resting state functional connectivity in the human language system using near infrared spectroscopy

We have evaluated the use of phase synchronization to identify resting state functional connectivity (RSFC) in the language system in infants using functional near infrared spectroscopy (fNIRS). We used joint probability distribution of phase between fNIRS channels with a seed channel in the language area to estimate phase relations and to identify the language system network. Our results indicate the feasibility of this method in identifying the language system. The connectivity maps are consistent with anatomical cortical connections and are also comparable to those obtained from functional magnetic resonance imaging (fMRI) functional connectivity studies. The results also indicate left hemisphere lateralization of the language network.


INTRODUCTION
Functional Near Infrared Spectroscopy (fNIRS) is an optical method for brain functional imaging. Using near infrared light with high penetration depth in human tissue, fNIRS can measure the neuronal activity in response to a task or an environmental stimulation through neurovascular coupling in superficial areas of the brain. fNIRS has been used for functional studies as a portable, less expensive and less restraining alternative to functional Magnetic Resonance Imaging (fMRI) in different task based brain functional studies. One of the more recent areas of interest in both fNIRS and fMRI is the study of interaction between different cortical areas through their intrinsic neuronal signaling (Biswal et al., 1995;White et al., 2009). This intrinsic signaling appears in the form of slow varying spontaneous fluctuations in the blood oxygen-level dependent (BOLD) signal in the absence of stimulation. These fluctuations are correlated between brain areas that are anatomically and functionally connected and have been used to map brain functional networks such as sensorimotor, visual, and auditory as well as higher networks such as language and attention (Beckmann et al., 2005;Zhang and Raichle, 2010). These networks are often mapped from data collected at rest and are therefore referred to as Resting State Functional Connectivity (RSFC) maps.
The RSFC analysis can also potentially identify changes in intrinsic neural activity as a result of disease in some neurological and psychiatric conditions. Changes in connectivity strength in different brain networks have been observed in conditions such as autism (Cherkassky et al., 2006), depression (Anand et al., 2005), Alzheimer disease (Li et al., 2002), and attention-deficit hyperactivity disorder (Tian et al., 2006).
Given the advantages of fNIRS, different brain networks have been investigated through RSFC using fNIRS. One of the most common methods for analyzing brain network connectivities using RSFC in fNIRS is cross correlation (Zhang and Raichle, 2010). In cross correlation, an fNIRS channel is selected as the seed channel and the correlations of the signal in all other channels with the seed channel are calculated. The objective is to find the cortical areas whose resting state fluctuations are similar to that of the seed channel. Cross correlation based functional connectivity has been used in conjunction with fNIRS to derive connectivity maps in different brain networks (White et al., 2009;Lu et al., 2010;Duan et al., 2012). Cross correlation and cross coherence of changes in fNIRS channels in particular have been used on infants to study functional connectivity. In a study by Homae et al., the functional connectivity in 3 month old infants with and without audio stimuli was investigated . The study showed that high correlation exists between neighboring channels and also their corresponding channels in the contralateral hemisphere during resting state. One drawback of correlation based connectivity is that it can be sensitive to detection of spurious connection as a result of presence of cross talk between channels, systemic interference or noise (Tass et al., 1998).
The phase of the oxygenation changes in fNIRS has been shown to contain information on functional connectivity. In a study by Taga et al., strong phase synchronization was observed between cortical regions during sensory processing in young infants . The phase information was also used in this study to determine the timing and flow of activity in the developing brain. The phase relationship between oxygenated and deoxygenated hemoglobin in resting state in infants was investigated in another study by Taga et al. (2000). Spatially synchronized oscillations were also observed in this study in the measurement channels over the occipital cortex (Taga et al., 2000). Phase based methods have also been used in other modalities such as EEG and MEG for functional analysis (Stam et al., 2007).
Previous studies have investigated the neural connectivity in infants using fNIRS. However, there have been limited studies on the neural systems involved in speech processing in the newborn with some level of inconsistencies in the data regarding the lateralization of the network. In this paper, we have investigated the application of a new method based on the phase relation between spontaneous oxygenated hemoglobin changes in fNIRS channels and have used it as a measure of functional connectivity to identify the speech processing neural systems in newborn infants. Compared to methods based on signal amplitude such as cross correlation, phase is much less sensitive to noise and interferences. It also does not require the assumption of stationarity for the signals. The phase synchronization is not equivalent to coherence or frequency synchronization and is an independent characteristic of the interrelationship between two processes (Tass et al., 1998). To evaluate the feasibility of this analysis method for detecting functional connectivity, we applied it to a study of processing of speech vs. non-speech in newborn human infants. This type of comparison is of particular value for the question being asked as there is a growing body of evidence on the brain areas involved in language processing in neonates, but less on the underlying connectivity.

fNIRS DATA
The fNIRS data were collected from newborn infants at BC Children's Hospital, Vancouver, Canada, during a separate language perception study (May et al., 2012). Informed consent was obtained from parents when the experiment was being conducted. The study design was approved by the ethics committee of the University of British Columbia. The experiment design and setup are shown in Figure 1. A total of 19 subjects were used in the analysis (mean age 1.4 days, 7 females, age range 0-4 days) out of which two subjects were excluded due to severe artifacts in the signals and poor data quality resulting from optode displacement during data collection. During the experiment, audio stimuli were administered to subjects while the subjects were in state of quiet rest or sleep. The audio stimuli consisted of blocks of sentences in Spanish and Silbo-Gomero. Silbo-Gomero is a whistled language that is a surrogate language of Spanish (Carreiras et al., 2005;Rialland et al., 2005) It uses whistles rather than speech, and was developed by shepherds in the Canary Islands to communicate across long distances. Spanish and Silbo-Gomero were selected as both are unfamiliar to the infants, while one is a spoken language and the other is not. Each block was 15 s long followed by 25-35 s of silence. A total of 8 blocks for each stimulus were presented in which each block consisted of continuous speech. The total experiment time was 22-25 min. The subjects' brain hemodynamic response was monitored by a 24 channel fNIRS device (Hitachi ETG-4000 machine with 695 and 830 nm lasers at a power of 0.75 mW, interoptode distance of 3 cm and sampling rate of 10 Hz). Two chevron shaped optode holders were used for securing nine 1 mm fibers to the head. There were a total of 4 detector and 5 source fibers on each holder resulting in 12 recording channels per holder. Figure 2 shows the placement of optodes on the subject's head. Surface landmarks (ears or vertex) were used for the placement of the probe holder over the infants' perisylvian area of the scalp. Channels 11 and 12 in LH and 23 and 24 were ideally placed above the infant's ear. A stretchy cap was used to secure the holders on the infants' head.

DATA ANALYSIS
In order to determine the phase relation between channels, we first extracted the phase of the signal in each channel using the Hilbert transform. The Hilbert transform converts a real valued signal to a complex one, known as analytic signal, whose real part and phase correspond to the original signal and its derived phase, respectively (Oppenheim et al., 1989). The Hilbert transform of signal x [n] in the frequency domain is defined as (Oppenheim et al., 1989): where X e jω is the Fourier transform of x[n] and sgn(ω) is the sign function having value of 1 for ω > 0 and -1 for ω > 0. The analytic signal can then be written as where y[n] is the inverse Fourier transform of Y(e jω ). We used the joint probability distribution of the phases across channels to describe their connectivity. A common model for probability distribution of phase which is the circular analog of the Gaussian distribution is the Von Mises distribution. The Von Mises probability density function (pdf) is defined as (Mardia and Jupp, 2000): where θ is an angle defined in the interval [−π, π) and I 0 (κ) is the modified Bessel function of order 0. The parameter κ is the equivalent of the covariance for the Gaussian distribution and μ is the expected value of the angle. The probability density function of the signal phase in channel m conditioned on that of channel n can therefore be written as: κ mn describes the intensity of the phase correlation between signals in channels m and n. In other words, it shows how much the prior information of θ n affects distribution of θ m . The first moment of the distribution given in Equation 3 can be calculated as (Mardia and Jupp, 2000) Using the first moment, one can estimate parameter κ by numerically solving the optimization problem arg min where N is the number of samples in the data segment and φ i is the measured phase of the signal at time point i. Parameter μ can then be estimated using A close relationship exists between parameters of the distribution and the phase locking value (PLV) which is a common measure used in EEG signal processing to detect functional connectivity through synchronization between channels. PLV is related to the distribution through the magnitude of the first circular moment of the phase distribution (Canolty et al., 2012) One advantage of using Von Mises distribution over phase locking for connectivity analysis is that once the parameters are estimated, one would have the distribution function and can resample from the distribution to determine the significance levels. Also, the preferred phase difference is not available in PLV.
To evaluate fNIRS functional connectivity, we first calculated the phase from all fNIRS channels using the Hilbert transform as described earlier. Ideally, we would be interested in phase relations between channels when subjects are not exposed to any type of stimulation to reveal intrinsic network activities. It was shown in an earlier study, however, that the infant brain produces no significant response in the language network to Silbo-Gomero stimuli (May et al., 2012). We therefore used the fNIRS data during Silbo-Gomero stimulation as an alternative to resting state. The BOLD spontaneous fluctuations are concentrated at frequencies of less than 0.1 Hz (Biswal et al., 1995). Therefore, the signals were first filtered with an infinite impulse response bandpass filter (IIR) between 0.02 and 0.08 Hz to extract spontaneous hemodynamic activities and reject other interferences. This frequency range is comparable to those used in other studies investigating RSFC using fMRI and fNIRS (Biswal et al., 1995;Duan et al., 2012).
fNIRS data in general can be contaminated with motion artifacts as the result of subjects' spontaneous movements. These artifacts create interference in the form of highly correlated phase changes in fNIRS channels, especially in spatially close channels. This interference results in a very high phase correlation and can obscure underlying phase connections between channels. Even though filtering of the motion artifacts is possible, in order to minimize possibility of introducing any inter-dependence between channels, no artifact removal procedure was applied. Instead, the channels for all subjects were inspected visually and artifact contaminated regions within the Silbo-Gomero stimulation window were marked. An artifact free segment of the data in each channel was then selected for the analysis and the phase of the selected signal segments were then derived using the Hilbert transform. Since the brain shows no response to this stimulation type, the stimulation onsets were ignored and the segments were selected independent of the stimulation onsets. The segments contained variable number of stimulation blocks and their length ranged from 50 to 220 s.
The channel with the highest activation in the grand average for the Spanish stimulation task during the original study in the left hemisphere was selected as the seed channel for the RSFC analysis. The joint phase distribution of the seed channel and all other channels was then estimated by calculating the phase difference between the seed channel and other channels and then estimating κ mn and μ mn using Equations 7 and 8. We used a simplex derivative-free method to solve Equation 7 and derive κ mn (Lagarias et al., 1998). The analysis was performed in MATLAB (Mathworks MA, USA) and the phase coupling estimation toolbox developed by Cadieu et al. was used for parts of the analysis (Cadieu and Koepsell, 2010) The analysis was performed on oxygenated hemoglobin (HbO 2 ) changes only. Previous studies on the application of fNIRS to detect language network activity and connectivity have shown that the oxygenated hemoglobin is more sensitive to regional cerebral blood flow changes than the deoxygenated hemoglobin with the equipment used here (Hoshi, 2007;Gervain et al., 2008;Lu et al., 2010).
To examine the validity and reliability of the connectivity information derived with this method, we divided the subjects randomly into two groups and evaluated the connections for each group, similar to the approach used in . We then compared the correlation of the connectivities between the groups.
Other studies have suggested that the language network is left lateralized. We verified this in the network derived using our method. The lateralization was quantified using (Binder et al., 1996;) where M is the total number of channels, κ i L is the value of κ is in which i is the channel number and s is the seed channel in the left hemisphere. κ i R is the value of the same parameter with the channel symmetric to i in the right hemisphere. The significance level of the calculated lateralization index is then evaluated. The lateralization index results in a number between −1 and 1 with more positive numbers indicating higher degrees of left lateralization.
As the final step, we defined 4 regions of interest (ROI), 2 inside the language network and 2 outside the network and evaluated the 1 http://redwood.berkeley.edu/klab/pce.html. connection strengths in these areas. In particular, channels 6 and 7 were selected inside the network, based on our prior knowledge that the physical area they cover is in the language network, and channels 1 and 12 outside the network with the optode configuration used in the current study. Channel 1 is over the frontal areas while channel 12 covers the temporal area. The choice of these channels as being outside the language area is justified by the fact that they showed no significant activation in response to native language or Spanish (May et al., 2012). Figure 3 shows the HbO 2 signal from channels 7 and 9, the seed, for a typical subject. Qualitatively, the histogram of the phase at different time points for both channels does not show a clear dominant phase range as shown in Figures 3B,C. However, the joint distribution histogram shown in Figure 3D has a sharp peak focused around the mean phase difference. This is also evident in the Von Mises distribution function plot derived from estimated parameters in each case (shown in red). The estimated distribution parameters are also indicated in the figure. In the case of the phase histogram for individual channels (Figures 3B,C), the values of estimated κ mn are much smaller than in the conditional distribution ( Figure 3D). This indicates a high phase relationship between the two channels and is interpreted as connectivity.

RESULTS
Using the method described earlier, the group level resting state functional connectivity map with channel 9 chosen as the seed channel was derived and is shown in Figure 4. The detected network includes the areas known to be associated with language network including the superior temporal gyrus and Broca's area. The maps are also in agreement with those obtained for the language network in adults using correlation based fNIRS connectivity studies .
The connectivity maps resulting from the 2 random subgroups are shown in Figure 5. The maps for the two subgroups cover similar areas in both hemispheres. The correlation between individual connections in the two subgroups is shown in Figure 6 (Pearson correlation r = 0.6). Figure 7 shows the results of ROI connectivity analysis where the connection strength between the seed channel and 2 channels in the language area (6 and 7) is compared with that with two channels outside the language system (channels 1 and 12). Analysis of variance indicates significant difference between the connections (ANOVA p < 0.01) inside and outside the language network. In particular, channels 6 and 7 connections are not different while they are both higher than that of channels 12 and 1 in the temporal and frontal areas, respectively. Results of multiple comparison test are shown in Table 1 (Tukey's test).
The results of the lateralization analysis are shown in Figure 8. The average lateralization index is 0.172 and is significantly different from zero (1 sample t-test, p < 0.001). Here, the lateralization index is also compared for all subjects between the language network and a control case. The control network is created by choosing channel 11 as the seed channel. There is a significant difference in lateralization index between the language and control network (paired t-test p < 0.01). These results suggest strong left lateralization in the detected language network.

DISCUSSION AND CONCLUSION
We evaluated use of phase synchronization to identify resting state functional connectivity in the language system in infants using fNIRS. We used joint probability distribution of phase between fNIRS channels with a seed channel in the language area to estimate phase relations and identify the language system network. Our results indicate the feasibility of this method in identifying the language system. The connectivity maps are consistent with anatomical cortical connections and are also comparable to those obtained from fMRI functional connectivity studies (Perani et al., 2011;Mahmoudzadeh et al., 2013). The results indicate left hemisphere lateralization of the language network. The brain networks connectivity reveals information about underlying anatomical areas involved in a particular task. In some disease conditions, changes in cortical connections occur (Zhang and Raichle, 2010). Application of connectivity estimating methods to fNIRS enables investigation of such changes in cases where use of fMRI is not possible, such as in infants and extends utilization of fNIRS in wider range of clinical applications.
Use of fNIRS for analysis of functional connectivity offers several advantages over more traditional fMRI based connectivity analysis. Collection of fMRI data from infants and young children under resting condition can be challenging. In contrary, fNIRS is easily applicable to even newborns. Also, in cases where subjects to be tested are immobilized and can not be transferred to an MRI scanner, portable fNIRS systems can replace fMRI for connectivity analysis. One limitation compared to fMRI is the limited penetration depth which means connectivity analysis will be limited to cerebral cortex. Use of phase information and in particular using the method described here has some advantages that can justify its use for functional analysis. In comparison with amplitude based methods such as cross correlation and cross coherence, Frontiers in Human Neuroscience www.frontiersin.org January 2014 | Volume 7 | Article 921 | 6 FIGURE 6 | Correlation between connection strengths in the two subgroups. channels 1 and 12). The bars indicate the standard error of the mean. phase-based methods do not suffer from spurious connections as a result of common source. The method described here can be generalized to estimate the phase distribution function in a multivariate approach (Canolty et al., 2012). In this way, the connectivity parameter between all channels is estimated at once. This eliminates the limitation of bivariate methods in detection of spurious connections. One limitation of fNIRS is the lack of information of anatomical structure of participants' brains. Due to the lack of an established standard system for processing infant brains, some procedure to localize the channels to the brain may be necessary. The probe set configuration we used in this study primarily allows measurement of activation from anterior temporal to posterior temporal areas, but given infants' small heads, the anterior temporal areas likely extend into frontal regions and the posterior temporal likely extend into parietal regions. Nonetheless, given the variability among infant head sizes and the lack of agreed upon stereotactic maps of infant heads, we chose the most conservative approach of using the general terms we have used, such as anterior and posterior temporal, to refer to the areas measured.

FIGURE 7 | ROI connections comparison between language area (represented by channels 7 and 6) and outside language area (represented by
Our results are comparable to similar studies in the literature. In particular Zhang et al. analyzed RSFC in the language system in adults using fNIRS . Their results indicated significant RSFC between left inferior frontal and superior temporal cortices which are associated with language system . Using fMRI, Fransson et al. studied resting state networks in the infants brain (Fransson et al., 2007). They observed similar networks in the bilateral temporal/inferior parietal cortex which encompasses primary auditory cortex (Fransson et al., 2007). In a study using fMRI on newborn infants, Perani et al. showed strong interhemispheric functional and structural connections between left and right temporal regions and also between the left and right frontal regions (Perani et al., 2011). Homae et al used the cross correlation and spectral coherence to investigate the functional connectivity in infants' brain . This study reports high correlation between adjacent channels and the corresponding contralateral channels. This is in agreement with our findings of connectivity between neighboring channels in the language area in the left hemisphere and corresponding channels in the right hemisphere using the phase-based method described here.
Other phase-based methods for analysis of functional connectivity have been reported earlier. In a study by Taga et al, the phase information of the oxygenated and deoxygenated hemoglobin was investigated in the presence of audio stimuli . The relationship between the stimulation and the evoked responses were used to evaluate the temporal features of the response. One important difference of our work is the absence of the stimulation evoked changes for detection of functional connectivity.
The presence of motion artifacts has a significant effect on the connectivity strength results using phase synchrony method presented here. The motion artifact results in in-phase changes across affected channels which may result in stronger phase correlation compared to those resulting from spontaneous neuronal activity. Therefore, care must be taken to ensure segments being processed do not include motion artifacts. Various methods for removal of motion artifacts in fNIRS have been proposed (Cooper et al., 2012;Brigadoi et al., 2014). In this particular study, we used visual inspection to remove motion contaminated data segments instead Frontiers in Human Neuroscience www.frontiersin.org January 2014 | Volume 7 | Article 921 | 7 of using an automated removal technique. An artifact removal method may introduce some level of distortion in the artifact-free regions of the signal and also may not fully remove the artifacts. Therefore, to avoid detection of connectivity as a result of such distortions, an automated approach was not employed. Saturated channels or channels that have lost coupling to tissue due to displacement will also have similar effect. Our fNIRS data was not collected during strict "resting" state. We used continuous blocks of data in which subjects were listening to a non-speech audio stimulation. No significant activation compared to baseline was observed on this data and was therefore used as the baseline. Some fMRI studies have used a similar approach for mapping RSFC. The study by Greicius et al. on default mode network in Alzheimer's disease patients for example, was performed during a low demand cognitive task (Greicius et al., 2004;Zhang and Raichle, 2010). An alternative approach in the case of task-related data is to regress out the task evoked response from the data before performing RSFC analysis (Fair et al., 2007;He et al., 2007).
The fNIRS signal is known to contain systemic interferences. This includes interference from cardiac pulsation, respiration, cardiovascular autoregulation and heart rate variability. The frequency band for connectivity analysis must be chosen such that it includes the relevant variations caused by the neurovascular coupling while rejecting the frequency bands containing these interferences. The cardiac interference in our study is around 2 Hz and the very low frequency interference (heart rate variability, cardiovascular autoregulation) is around 0.01 Hz. The respiratory fluctuation is around 0.2 Hz. The frequency band we chose for analysis (0.02-0.08 Hz) does not include these interferences and therefore, connectivity detection as a result of these interferences in unlikely.
Bivariate methods in general can result in non-existing spurious connections when there is a propagation of information from one channel to others. A pairwise measure of connectivity will result in detection of connection between all possible connection pairs, ie. direct or indirect. However, since we are only looking to find channels which belong to the same network, use of a bivariate measure of connectivity can be justified. Most methods for functional connectivity mapping in the literature based on fMRI or fNIRS also use seed based methods which is relying on the bivariate concept of finding coherence/correlation between channels with the seed channel (Greicius et al., 2004;Zhang and Raichle, 2010).
In summary, the results of this paper suggest that the proposed method can be used to reveal underlying connectivity patterns of cognitive functions in the resting state through phase relations between hemodynamic changes in different brain regions. The results also indicate a left lateralization in the detected network which suggests the language system may be left lateralized already in newborns. Most fNIRS studies have shown left lateralization of the language network in older infants and in adults. However, data have revealed inconsistencies regarding lateralization of this network in newborn infants, when examining activation in response to language stimulation (Pena et al., 2003;May et al., 2011;Perani et al., 2011;Sato et al., 2011). The results of this study in terms of the lateralization of the network provide further information and evidence for characterization of the language network in newborn infants. Thus, these results significantly add to the theory on the neural foundations of language acquisition by revealing that language processing networks are left lateralized even before language acquisition proper has begun.