1.

Introduction

Functional near-infrared spectroscopy (fNIRS) has been proposed as a promising technique for hemodynamic imaging of brain activities.^{1, 2} The near-infrared light penetrates the skull and reaches the cortical layer of the brain. It is mainly absorbed by two principal chromophores in tissue, i.e., oxygenated hemoglobin $\left(\mathrm{Hb}{\mathrm{O}}_2\right)$ and deoxygenated hemoglobin (Hb) that have significantly different absorption coefficients in the near-infrared region (670 to 900 nm).³ Therefore, by using two or more wavelengths in this region, changes of cerebral $\mathrm{Hb}{\mathrm{O}}_2$ and Hb concentrations can be quantified. The fNIRS technique measures the hemodynamic signals of the brain similarly to the blood-oxygen-level-dependent (BOLD) signal in functional magnetic resonance imaging (fMRI).⁴ While the spatial resolution of fNIRS is limited^{5, 6, 7} compared to that of fMRI, the former has a remarkable temporal resolution down to milliseconds and has the advantages of being compact, low cost, and easy to use without complete body confinement. It has been widely accepted by the researchers that fNIRS and fMRI are complementary to one another in brain functional imaging.

Early studies in both fMRI⁸ and fNIRS⁹ pointed to a linear relationship between the stimulus parameters (i.e., duration and intensity) and the resultant cerebral, hemodynamic response (BOLD signal in fMRI, $\mathrm{Hb}{\mathrm{O}}_2$ and Hb concentrations in fNIRS). A linear system model,¹⁰ as defined by the principle of superposition with respect to the stimulus duration [ $(f(x1+x2)=f\left(x1\right)+f\left(x2\right))$ ] and of scaling with respect to the stimulus intensity [ $f\left(\alpha x\right)=\alpha f\left(x\right)$ , where α is a scalar constant], was applied to predict the cerebral response while the stimulus parameter was changed.^{11, 12, 13} However, recent studies in fMRI using various types of stimuli showed that considerable nonlinearity in BOLD signals did exist when utilizing a short period of stimulus (several seconds), and that such nonlinearity varied across cortices. For example, studies looking at the primary visual cortex reported that the BOLD signal was nonlinear for stimuli less than 4 sec,¹⁴ or nonlinear for stimuli less than 3 sec,^{15, 16} and linear above these durations. Two studies looking at the primary auditory cortex reported that the BOLD signal was nonlinear for stimuli less than 6 sec,¹⁷ or nonlinear for stimuli less than 10 sec,¹⁶ and linear above these durations. Soltysik ¹⁶ also reported that the BOLD signal in the primary motor cortex is nonlinear for stimuli less than 7 sec and linear above this duration. While the exact implication of this nonlinear behavior still needs to be further explored, such studies allow us to investigate the neuronal adaptation and other functional properties of neurons within a specific brain region.¹⁸ Also, knowing the ranges of stimulus duration and rate for the nonlinear response to occur would be useful in the experiment design, such as precluding the use of short intervals between experimental stimuli, etc.

The prefrontal cortex (PFC) is an executive area for cognitive activities that deal with memory, language, and problem solving.^{19, 20, 21} A number of studies on the roles of PFC during working memory tasks have been conducted on adults using fMRI,^{22, 23} and on adults and children using fNIRS.^{24, 25} Most of the higher cognitive functions have shown some use of working memory. Similar studies have also been conducted on infants by Baird ²⁶ to assess the cognitive maturation of infants associated with growth. However, less attention has been paid to the potential nonlinear behaviors associated with cognitive tasks. While two studies have been conducted to investigate the rate effects of verbal tasks and have reported more significant nonlinear behaviors in multiple cognitive cortices than other cortices,^{27, 28} knowledge on the stimulus duration for the nonlinear response to occur in cognitive cortices is still limited.

In the research area of fNIRS, little has been reported on the nonlinear relationship between the stimulus parameters and the resultant hemodynamic parameters derived from fNIRS. A possible reason is that the nonlinear behavior at the cortices exists only within a short stimulus duration, which has not been the focus in most of the previous fNIRS studies. Thus, we conducted a study in the PFC in response to the anagram tasks using a multichannel fNIRS imager. The duration (i.e., the executing time) of the anagram task was variable, ranging from several seconds to more than one minute. This study provided us with a broad view on the variability of the PFC in response to cognitive, problem-solving tasks. We hypothesized that the nonlinear behavior also existed between the short-time stimulus and the fNIRS-derived readings. This hypothesis is based on the fact that the BOLD signal in fMRI and hemoglobin concentrations derived from fNIRS share similar hemodynamic processes occurring in the cortices during brain activities. Specifically, we examined the PFC hemodynamic responses in three discrete subgroups (based on their task durations), and utilized the superposition principle of the linear system model to determine the range and extent of the potential nonlinearity. A better understanding of nonlinearity in fNIRS signals will allow us to investigate brain adaptation and to extrapolate neuronal activities from hemodynamic signals.²⁹

2.

Materials and Methods

2.2.

Principle of Near-Infrared Spectroscopy Imaging

Optical density (OD) is conventionally defined as $\mathrm{OD}=\log [I_b∕I\left(t\right)]$ , where $I_b$ is the baseline optical intensity and $I\left(t\right)$ is the optical intensity at time $t$ . The principle of continuous-wave NIRS imaging is based on the modified Beer-Lambert law,^{1, 30} in which the relationship between the change in optical density (ΔOD) due to light absorption perturbation and the relative changes of $\mathrm{Hb}{\mathrm{O}}_2$ and Hb concentrations from baseline can be expressed as:

Eq. 1

$\Delta \mathrm{OD}\left(\lambda \right)=r\cdot \mathrm{DPF}\left(\lambda \right)\cdot \{{\epsilon }_{\mathrm{Hb}}\left(\lambda \right)\cdot \Delta \left[\mathrm{Hb}\right]+{\epsilon }_{\mathrm{Hb}{\mathrm{O}}_2}\left(\lambda \right)\cdot \Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]\},$

where ε is the molar extinction coefficient $(\mu {\mathrm{M}}^{-1}\cdot {\mathrm{mm}}^{-1})$ of $\mathrm{Hb}{\mathrm{O}}_2$ or Hb,³ Δ[Hb] and $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ are the changes of Hb and $\mathrm{Hb}{\mathrm{O}}_2$ concentrations (μM), $r$ is the distance between the source and detector, and DPF is the differential path length factor indicating the lengthening of the average optical path length due to light scattering in tissue. Therefore, by measuring the changes of optical densities at two wavelengths, the changes of $\mathrm{Hb}{\mathrm{O}}_2$ and Hb concentrations can be calculated as:

2.

$\Delta \left[\mathrm{Hb}\right]=A\cdot \Delta \mathrm{OD}\left({\lambda }_1\right)-B\cdot \Delta \mathrm{OD}\left({\lambda }_2\right),$

$\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]=-C\cdot \Delta \mathrm{OD}\left({\lambda }_1\right)+D\cdot \Delta \mathrm{OD}\left({\lambda }_2\right),$

where $A$ , $B$ , $C$ , and $D$ are constant coefficients based on the molar extinction coefficients at the respective wavelengths.

The instrument used for prefrontal scanning was a LED-based imager (LEDi, NIM Incorporated, Philadelphia, Pennsylvania).^{31, 32} The system had a flexible probe to match the contour of the human forehead. The probe consisted of four LED diodes as light sources (at ${\lambda }_1=730\mathrm{nm}$ and ${\lambda }_2=850\mathrm{nm}$ ) and ten photodiodes as detectors that were symmetrically arranged in an area of $3.5\times 14{\mathrm{cm}}^2$ , conducing to 16 nearest source-detector pairs (i.e., channels) at 2.5-cm separation [see Fig. 1a ]. During the experiment, a Velcro band held the probe firmly on the forehead, and extended from ear to ear in axial direction and from hairline to eyebrow in a sagittal direction [see Fig. 1b]. The optically probed area should cover approximately the anterior PFC, a part of the dorsolateral PFC, and a part of the ventrolateral PFC. The four LEDs flashed in sequence; the reflected light from the brain was detected with the nearest photodiodes of each LED and converted into digital signals using an analog-digital converter (ADC) card in the control box. Then the digital data were sent to a laptop computer through a serial port. The sampling rate across all 16 channels was 3 Hz.

Fig. 1

(a) The LEDi system. The system utilized four light emitting diodes as light sources. The reflected light from the brain was detected with ten photodiodes, and then converted into digital signals using an ADC card in the control box. Then the digital data were sent to a laptop computer through the serial port. (b) The LEDi probe placed on a forehead. The probe covered an area from ear to ear in the axial direction and from hairline to eyebrow in the sagittal direction. The locations of 16 channels are indicated by numbers.

2.3.

Data Screening and Analysis

Data screening and analysis were done using Matlab 6.5 (The Mathworks, Natick, Massachusetts). The data were processed for each individual test and analyzed channel by channel for 16 channels. The major factor to be considered in data processing was the physiological noises, which have received increasing attention in recent studies.^{33, 34, 35, 36} Besides the functional activation at the cortex, the fNIRS signals also comprise various oscillations related to the spontaneous physiological activities. The most pronounced spontaneous oscillation stems from the cardiac pulsation at about 1 to 2 Hz. It is a very fast oscillation compared to the functional activation and can be filtered out directly. In the lower frequency range, three standard frequency bands were defined³⁵: 1. the high frequency (HF) oscillation related to respiration, having a broad peak at the respiratory rate of around 0.3 Hz; 2. the low frequency (LF) oscillation at about 0.1 Hz resulting from vasomotion³⁷; and 3. the very low frequency (VLF) oscillation between 0.02 and 0.04 Hz. Little is known about the physiological interpretation of VLF oscillation. The magnitude of the physiological noises can be much bigger than that resulting from the functional activation, because these noises are systemic and embedded in both the cortex and overlying tissues.³⁶ The physiological noises are implicated in two different aspects: first, the spontaneous oscillations are distinctive characteristics of in-vivo optical reflectance. In data screening, we used a continuous wavelet transform (CWT) algorithm to decompose and examine these spontaneous oscillations to exclude the invalid data that usually resulted from poor probe contact on the scalp. Second, the HF, LF, and VLF oscillations could induce significant aliasing into data analysis, since they have very similar temporal characteristics to the functional activation. A practical strategy to reduce the physiological noises is to perform task-related averaging, since these noises are not temporally correlated. Two steps were followed during the data processing.

1. Data screening: a wavelet transform is a variable-scale representation of signal with a wavelet function that is called a mother wavelet.³⁸ The wavelet transform of a discrete signal $x\left(n\right)$ is the convolution of $x\left(n\right)$ with a scaled and translated mother wavelet $\psi \left(n\right)$ :

Eq. 3

$C_{\alpha ,k}=\frac{1}{\sqrt{\alpha }}\sum _{n=1}^{N}x\left(n\right)\psi \left[\frac{n-k}{\alpha }\right],$

where α is the scale factor, $k$ is the translational factor, and $C_{\alpha ,k}$ are wavelet coefficients. Briefly, by altering the scale factor α, the wavelet transform provides a decomposition of the signal at a certain frequency band. In this study, the fourth derivative of Gaussian function was used as the mother wavelet [see Fig. 2a ]. Normalized wavelet coefficients at integer scale α = 5, 15, and 50 were adopted to represent HF, LF, and VLF oscillations, respectively (see Table 1 ). If the LEDi probe contacted the scalp firmly and the reflected light from the brain was acquired properly, continual and systematic oscillations across all channels would be seen clearly after CWT decomposition. For instance, Figs. 2b, 2c, 2d, 2e show such a segment of undecomposed cerebral data (ΔOD, ${\lambda }_2=850\mathrm{nm}$ ) as well as its decomposed VLF, LF, and HF components, respectively. If some channels were in poor contact with the scalp during the test, the decomposed components in these channels were significantly asynchronous to the other channels, and frequently discontinued due to the motion artifacts. Consequently the data from this test would be excluded.

Fig. 2

(a) The fourth derivative of the Gaussian function as the mother wavelet in CWT; (b) the undecomposed, cerebral signal (ordinate: channel number); and (c) to (e) the decomposed CWT components at scale factor α = 50, 15, and 5, corresponding to the VLF, LF, and HF spontaneous oscillations. Continual and systematic spontaneous oscillations across all channels are observed with good probe contact on the scalp.

Table 1

The pseudofrequency and 3-dB passband of CWT at scale α = 5, 15, and 50.

α	5	15	50
$f_p$ (Hz)	0.30	0.10	0.03
Passband (Hz)	0.21 to 0.37	0.07 to 0.12	0.03 to 0.04

2. Data analysis: after data screening, the raw data from each valid test were filtered with a high-pass Butterworth filter (the cut-off frequency was at 0.02 Hz) to remove the long-term baseline drift. A low-pass Butterworth filter (the cut-off frequency was at 0.4 Hz) was also applied to remove the fast oscillation stemmed from cardiac pulsations. Then $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] were calculated based on Eqs. 1, 2. It is known that a shorter wavelength should have a little longer DPF because its scattering is relatively larger, and it was also noted³⁶ that inaccurate DPF values will cause a certain degree of cross talk between $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb]. Thus, we utilized two DPF values, DPF = 6.7 for ${\lambda }_1=730\mathrm{nm}$ and DPF = 5.7 for ${\lambda }_2=850\mathrm{nm}$ , according to the in-vivo measurements on the adult head.³⁹

As was mentioned before, the duration of each task, defined from the onset of the task to the following key-stroke entry, was variable. Hence the length of an entire anagram segment, which was the duration of the task plus the following 30-sec rest period, was also variable. To study the variability and the possible nonlinearity of the hemodynamic responses, the anagram segments were selected into three discrete subgroups according to their durations and expressed as (mean ± standard deviation): D1 = 2.0 ± 0.3 sec, D2 = 4.0 ± 0.7 sec, and D3 = 8.0 ± 1.0 sec. Each subject completed her anagram tasks with some short or long durations; therefore, each subject’s dataset consisted of D1, D2, and D3 subgroups. Data averaging was done within each individual subject to reduce the disturbance of the physiological oscillations: all of the data segments of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] belonging to the same subgroup were aligned at the time of task onset and averaged as the individual pattern for the respective subgroup (D1, D2 or D3). To ensure the physiological oscillations to be effectively reduced after averaging, only the subgroups consisting of seven or more anagram segments were adopted. Then at the group level, six- and four-letter anagrams were grand-averaged and analyzed separately. Notice that the activated channels (or brain regions) in the PFC during anagram might be different among subjects. Therefore, at the group level only the most active channels from each subject were picked and grand-averaged. The variation of the most active channels among subjects was investigated separately.

3.

Results

Among the six female subjects, three subjects (subjects 1 to 3) chose a six-letter anagram, and the other three (subjects 4 to 6) chose a four-letter anagram. After data screening, the number of valid tests for each subject is shown in Table 2 . The accuracy of their performance ranged from 61.4 to 72.0%, which did not show significant dependence on the difficulty level (six- or four-letter anagram) that they chose.

Table 2

The performance of six subjects having valid tests.

Two histograms of anagram duration from subjects 3 and 4 are shown in Figs. 3a and 3b, respectively, demonstrating that each individual’s histogram is commonly more concentrated on the left side, apparently having the exponential distribution. Some extremely long durations are observed from each subject, meaning that the corresponding tasks might be too difficult to solve and subsequently excluded from data analysis. The mean duration for each subject (as given in Table 2) showed significant dependence of duration on the difficulty level: subjects 1 to 3 having the six-letter anagram did take much more time to complete the tasks than subjects 4 to 6 having the four-letter anagram. For subjects 1 to 3, all of D1, D2, and D3 subgroups met the criterions of having seven or more tasks and were selected for data analysis. However, for subjects 4 to 6 in the four-letter anagram group, only subgroups of D1 and D2 were adopted for data analysis, since the D3 subgroup did not have enough anagram tasks (<7 tasks) as we predefined (see Table 3 ).

Fig. 3

The histograms of task duration from two subjects: (a) subject 3 choosing the six-letter anagram; and (b) subject 4 choosing the four-letter anagram. The bar width in each histogram is 3 sec. For each subject, task durations apparently obey the exponential distribution.

Table 3

The number of tasks completed within durations D1 = 2.0 ± 0.3 sec, D2 = 4.0 ± 0.7 sec, and D3 = 8.0 ± 1.0 sec and most active channels in each adopted subgroup for six subjects.

3.1.

Temporal Characterization

For the subgroups with the same task durations, similar hemodynamic patterns due to solving the anagram tasks were observed across the subjects. We selected the most active channels (see the next section to know how to determine this selection) on the left and right sides of the dorsolateral PFC from each subject, as listed in Table 3. Then, we plotted Figs. 4a, 4b, 4c to show the grand-averaged time courses of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] (mean ± SD, $N=3$ ) with task durations of D1 to D3, respectively. The data were taken from subjects 1 to 3 who performed six-letter anagrams. A transient “calm” period about 2 to 3 sec after the anagram onset $(t=0\sec )$ was observed in all the subgroups, followed by a significant deactivation in the cortex as $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ decreased and Δ[Hb] increased. This deactivation maintained a little longer in the subgroups having longer task duration. Then both $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] gradually recovered, eventually reaching a distinct overshoot in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and undershoot in Δ[Hb] in all three subgroups. Moreover, the overall variation of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ was always bigger than that of Δ[Hb] in the same subgroup. Similarly, Figs. 4d and 4e show the grand-averaged time courses of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] (mean ± SD, $N=3$ ) with task durations of D1 and D2, respectively, which were from subjects 4 to 6 who performed four-letter anagrams. A reduced deactivation, as compared to the respective subgroups with six-letter anagrams, was also observed in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb].

Fig. 4

The grand-averaged temporal profiles of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ (red curves) and Δ[Hb] (blue curves) concentrations in response to anagram execution: (a) to (c) correspond to subgroups of task duration D1 = 2.0 ± 0.3 sec, D2 = 4.0 ± 0.7 sec, and D3 = 8.0 ± 1.0 sec in six-letter anagrams (mean ± SD, subjects 1 to 3); (d) to (e) correspond to subgroups of task duration D1 = 2.0 ± 0.3 sec and D2 = 4.0 ± 0.7 sec in four-letter anagrams (mean ± SD, subjects 4 to 6). The gray area in each graph indicates the task duration. The bottom and top ${}^{*}$ symbols in each graph indicate the period of significant changes from the baseline $(\mathrm{p}<0.05)$ in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb], respectively. The variation of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ was measured as the height from the maximum decrease (P1) to the subsequent maximum overshoot (P2), as indicated in graph (a). (Color online only.)

3.2.

Spatial Characterization

To investigate the spatial patterns of PFC response, pseudocolor topographies of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] in the probing area are created based on the hemodynamic data taken from 16 channels. Figure 5 shows the PFC topographic series of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ with a 3-sec interval, ranging from 0 to 18 sec associated with an anagram task (0 sec represents the onset of task). Data were taken from subject 2, with the task duration of D3 = 8.0 ± 0.7 sec. Each hexagon in the graphs corresponds to the reading from the respective channel at the same location, with its color scale proportional to the temporal value at the given times. As shown in this set of figures, the maximum deactivations (the second and third frames of Fig. 5) were seen at channels 5 and 11 that were located on the left and right dorsolateral PFC, respectively, followed by big overshooting increases approximately occurring in the neighboring regions.

Fig. 5

PFC topographic series of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ with a 3-sec interval (ranging from 3 to 18 sec) that are associated with the anagram task (0 sec represents the onset of task). Data were taken from subject 2, subgroup D3 = 8 sec. Each hexagon in the images represents $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ reading from one channel. Its color scale is proportional to its temporal value of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ in μM right at the given times, as labeled to the left of the figures. The 16 channels are indicated by numbers shown at the top image. (Color online only.)

For each individual subgroup with different task durations, the most active channels during execution of the anagrams were determined based on the peak-to-peak amplitude of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ variation, i.e., the height from P1 to P2 as indicated in Fig. 4a. For instance, Figs. 6a, 6b, 6c show the PFC topographies based on the peak-to-peak amplitudes taken from subject 2 across three duration subgroups (D1 to D3). Channels 5 and 11 in the dorsolateral PFC consistently are the most active channels. Similar and consistent results were also observed from other subjects, as listed in Table 3. Such results clearly demonstrate that each subject kept a consistent pattern to solve the anagram tasks through the tests, activating the same brain regions in the PFC regardless of different durations spent on each task. Interindividual comparison showed that dorsolateral PFC (channels 3 and 5 on the left, and channels 11 and 13 on the right side of PFC) was the region most frequently engaged in anagram tasks. Furthermore, the overall spatial trends obtained at the two difficulty levels (four- and six-letter anagrams) show a similar pattern, while the most active channels slightly differ in amplitude from subject to subject. Therefore, it appears that this functional segregation depends mainly on the characteristics of the individual subject.

Fig. 6

Comparison of PFC topographies in terms of peak-to-peak amplitude in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ (the height from P1 to P2 in Fig. 4). Data were taken from subject 2 across three duration subgroups: (a) D1 = 2.0 ± 0.3 sec, (b) D2 = 4.0 ± 0.7 sec, and (c) D3 = 8.0 ± 1.0 sec (mean ± standard deviation). Channels 5 and 11 consistently are the most active channels across different temporal subgroups. The unit of the color bars is μM. (Color online only.)

3.3.

Test of the Superposition Principle

If the brain activation in response to certain stimulation follows a linear system model, the hemodynamic response to a long-duration task $f(x1+x2)$ can be predicted by the summation of the responses to several shorter tasks $f\left(x1\right)$ and $f\left(x2\right)$ , based on the principle of superposition: namely, $f(x1+x2)=f\left(x1\right)+f\left(x2\right)$ , where $x1$ and $x2$ represent two time durations and $f\left(x\right)$ denotes the function of hemodynamic response. Assuming that the cognitive load is statistically constant, or has constant amplitude, for all the anagram tasks at the same difficulty level, the test of superposition was performed to examine if the brain activation in response to anagram tasks in the PFC is linear. The data taken from the six- and four-letter anagram groups were analyzed separately, and only the most active channels in dorsolateral PFC were used, as given in Table 3.

In the six-letter anagram group, the first test for superposition was to use the response with a task duration of D1 = 2.0 ± 0.3 sec to predict the response with a duration of D2 = 4.0 ± 0.7 sec. To do so, the hemodynamic signals of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] with a duration of D1 were shifted by 2 sec and then added to the original signals to predict the signals with a task duration of D2 [see Fig. 7a ]. Figures 7c and 7d are generated in this way and clearly show that the added results of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] overpredict the actual respective measured signals, meaning that the superposition principle is not valid for the measured signals of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] if the shorter stimulation signals of 2 sec are used for prediction. Close inspection on the predicted and measured signals, as expressed by the maximum $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ decrease and maximum Δ[Hb] increase after anagram task [ $A_p$ and $A_m$ in Figs. 7c and 7d], also reveals that the predicted signals are approximately 52 and 32% larger than the actual measured signals in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb], respectively.

Fig. 7

Test of superposition principle using the six-letter anagram. (a) Schematic diagram to predict the response function, ${\mathrm{y}}^{\prime }\left(\mathrm{D}2\right)$ , at D2 = 4 sec using the measured signal, y(D1), at D1 = 2 sec. (b) Schematic diagram to predict the response function, ${\mathrm{y}}^{\prime }\left(\mathrm{D}3\right)$ and ${\mathrm{y}}^{″}\left(\mathrm{D}3\right)$ at D3 = 8 sec using the measured signals, y(D1) or y(D2) at either D1 = 2 sec or D2 = 4 sec. (c) and (d) show comparison between the measured $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ , Δ[Hb] signals at D2 = 4 sec (solid) and the predicted signals from D1 = 2 sec (dash). (e) and (f) show comparison between the measured $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ , Δ[Hb] signals at D3 = 8 sec (solid) and the predicted signals from D1 = 2 sec (dash). (g) and (h) show comparison between the measured $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ , Δ[Hb] signals at D3 = 8 sec (solid) and the predicted signals from D2 = 4 sec (dash).

In the second test for linearity, the hemodynamic signals of $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] with the duration of D3 = 8.0 ± 1.0 sec were predicted using both duration subgroups of D1 and D2, respectively [see Fig. 7b]. As shown in Figs. 7e and 7f, the predictions using subgroup D1 are much larger than the measured signals, having approximately 114 and 88% overprediction in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] amplitude, respectively. The predictions using subgroup D2 approach much closer to the measured signals, as shown in Figs. 7g and 7h, having only about 52 and 18% overprediction in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb], respectively. This means that the hemodynamic signals taken at the PFC are highly nonlinear when the task duration is shorter than 4 sec, but tend dramatically to be linear if the task duration is longer than 4 sec.

Similar superposition tests were performed in the four-letter anagram group. In this case the hemodynamic signals with D1 duration were used to predict the signals with D2 duration [see Fig. 7a]. Figure 8a shows 25% overprediction in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ amplitude, while the predicted time course in Δ[Hb] matched well with the measured signal, as shown in Fig. 8b. The latter means that the nonlinear effect is not obvious in Δ[Hb]. It is not surprising since Δ[Hb] signals are often weaker than $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ signals, so Δ[Hb] could be more easily subject to and suffering from physiological noises. The nonlinear effects shown in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ signals are expected to be more reliable than those in Δ[Hb].

Fig. 8

Comparison between the measured $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ , Δ[Hb] signals and the predicted signals from the four-letter anagrams (subjects 4 to 6). Specifically, the measured (a) $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and (b) Δ[Hb] with D2 = 4 sec task duration (solid curves) are predicted by adding two temporal signals with D1 = 2 sec task duration (dashed curves), respectively.

4.

Discussion

In this study we have investigated the temporal and spatial variability of hemodynamic responses i.e., $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb], in the PFC of six human subjects during anagram tasks, with respect to different durations and difficulty levels. From temporal evolution, we observed early deactivation in the PFC as $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ decreased and Δ[Hb] increased after the task onset; both $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and Δ[Hb] maintained for some time before returning to the baseline, followed by an overshoot in $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ and an undershoot in Δ[Hb]. It appears that the early decrease in Δ[HbO] is affected by the task difficulty level: the PFC response caused by the six-letter anagram shows a bigger, longer decrease than that caused by the four-letter anagram under the same task duration.

From spatial segregation, the dorsolateral PFC near the locations of channels 3 and 5 on the left forehead and channels 11 and 13 on the right forehead show predominant signals during the anagram tasks. This is expected, since the dorsolateral PFC is believed to be specialized for active manipulation and monitoring of information in verbal tasks.²² The slight variation between channels 3 and 5 as well as between channels 11 and 13 across subjects is believed to result from the individual head size. For some of the subjects, channel 16 at the right ventrolateral PFC showed stronger deactivation than channels 11 and 13; this observation is also reasonable, since the brain area under channel 16 is frequently involved in verbal tasks.²² Moreover, within each individual subject, good consistency in spatial segregations across different duration subgroups (D1, D2, and D3) is observed.

In this study, the principle of superposition among three subgroups has been specifically tested. The three durations were used to apply the principle of superposition to measured data, i.e., we wish to examine whether the hemodynamic responses will be doubled or tripled if the stimulation periods are doubled or tripled. Our data analysis from both the four- and six-letter anagrams has clearly shown a nonlinear behavior in hemodynamic signals in the PFC with the task durations shorter than 4 sec. The nonlinear behavior in $\mathrm{Hb}{\mathrm{O}}_2$ is more robust than that in Hb. To our current knowledge, this is the first investigation on the nonlinearity of cognitive response based on the principle of superposition. Our finding is complementary to the previous fMRI studies on the nonlinear rate effects in cognitive cortices based on the principle of scaling.^{27, 28}

Although fNIRS and fMRI measure similar cortical hemodynamics, several factors need to be further considered while comparing the nonlinear phenomena between the two technologies. First, the spatial resolution of fNIRS is relatively crude and is mainly determined by the measurement density of the optode array,⁷ being in the range of centimeters in human brain studies.^{5, 6, 7} Since the sampling volume in fNIRS is much bigger than the voxel size in fMRI, the nonlinearity observed in fNIRS might be equivalent to that obtained through multi-voxel averaging in BOLD signals at the corresponding cortical areas. Second, fNIRS is more sensitive to the hemodynamic signals in small vessels (i.e., arterioles, venules, and capillaries).⁴⁰ Given the strength of the static magnetic field, in contrast, the BOLD signals in fMRI may be associated primarily with large vessels and particularly draining veins.^{41, 42, 43, 44, 45} Therefore, a simultaneous comparison of nonlinearity between fNIRS signals and BOLD signal in fMRI has potential to better understand the nonlinear effect of cerebral hemodynamics at both micro- and macrovascular levels.

Studying and understanding nonlinearity of hemodynamic responses to specific stimuli is significant for optimizing experimental designs and data analysis. In general, it is believed that a longer or stronger stimulation signal will lead to a bigger hemodynamic response. However, this statement may not be completely true if the individual trials are spaced too closely, with fast-rate activation, due to the existence of nonlinearity. Large nonlinearities could result in misleading conclusions when using very short intervals between experimental stimuli. Recognizing the appropriate temporal range when the nonlinearity presents will help design better experiments and improve the sensitivity of measurements, as demonstrated in the field of fMRI.¹⁸

Another important reason to study nonlinearity is our attempt to study neural activities to be derived from hemodynamic measurements. Since fNIRS measures changes in hemoglobin concentrations in the brain, this technique does not directly detect any neural or neuronal activities. We expect that nonlinearity study may allow us to extrapolate or explore neural parameters in the following way. When a nonlinear pattern is detected with multiple stimulation durations in such a way as we showed earlier, one can mathematically model a neural adaption function and fit the model to the experimental response functions. A set of fitted parameters can then be obtained to reflect neural activities. Such approaches have been made by an fMRI study⁴⁶ and adapted in our recent fNIRS investigation.⁴⁷ While our expectation needs to be further confirmed, we believe that a good understanding of nonlinearity in fNIRS signals will permit the exploration of neuronal activities based on hemodynamic signals.

One noticeable observation in our study is the significant, early decreases in anagram response. The deactivations have maintained about 10 sec after the anagram task started and seem very unlikely to result from the initial dip,^{18, 48} which results mainly from a temporary lack of oxygen in the activated brain area, being caused by a rapid increase in neural activities. To date, such a decrease in cerebral blood flow and oxygenation, i.e., the negative BOLD response (NBR), has been frequently observed and reported in the fMRI field,^{41, 49, 50, 51} while NBR is opposite to the conventional thought that the brain activation leads to an increase in regional blood flow and oxygenation. The physiological origin of NBR is still an open question and has not been fully understood.

Several explanations for NBR have been given in a variety of fMRI studies using different brain stimulation protocols: one interpretation was associated to the passive hemodynamic changes around the activated brain area, reflecting an effect of blood draining from the neighboring areas, a so-called “blood steal.”⁴⁹ Another interpretation against blood stealing was reported in a study of the visual cortex,⁵² suggesting that NBR resulted from an active, attention-related suppression of neural activity. Indeed, NBR has been frequently observed in dorsolateral PFC during versatile video-game-like tasks;^{53, 54, 55} similar to the latter protocols, several fNIRS studies have further confirmed sustained decreases in $\mathrm{Hb}{\mathrm{O}}_2$ .^{56, 57, 58} To compare all of these reported studies with NBR, there was little in common except that all the tasks required a response to certain visual stimuli. Such visual stimuli exist in our anagram tasks. We hypothesize that the dorsolateral PFC deactivation in our study results from attention-related suppression, because manipulating or reordering the letters in an anagram requires sustained visual attention. However, we cannot completely rule out the possibility of blood steal based on our fNIRS results, because the probe size in our study was not able to cover the entire dorsolateral PFC.

Given all possible interpretations and hypothesis as discussed, a more recent animal study brings another view⁵¹ that negative BOLD signals do not necessarily mean decreased neuronal activity, but can stem from increased neuronal activity. NBR highly depends on the interaction and balance between hemodynamics and metabolism. Such a conclusion indicates the complexity and uncertainty on how to interpret NBR, which is still under current research exploration in the fMRI field. Similarly, a better understanding on the decrease in $\mathrm{Hb}{\mathrm{O}}_2$ signals seen in our anagram study needs to be further explored, while such PFC hemodynamic deactivation is unambiguously observed and in good agreement with multiple fMRI^{51, 52, 53} and fNIRS^{54, 55, 56} studies.

5.

Conclusion

In this study, we wish to prove our hypothesis that nonlinearity between short-term anagram tasks and corresponding hemodynamic responses can be observed by fNIRS at the PFC. With duration-variable anagram tasks, we observe a significant hemodymanic deactivation, as $\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\right]$ decreases and Δ[Hb] increases, during and shortly after the task executions. By using the superposition principle of a linear system, we demonstrate clear nonlinearity in hemodynamic responses in the PFC with task durations shorter than 4 sec. This study basically proves that fNIRS can be utilized to observe a nonlinear relationship between the stimulus parameters and the resultant hemodynamic parameters in the PFC when a short-term anagram task is given. The significance of studying nonlinearity by fNIRS is to possibly investigate brain adaptation and to extrapolate neuronal activities from hemodynamic signals, and thus to comprehend further about neurovascular coupling during brain activities. Furthermore, our observation of significant deactivation in early hemodynamic responses in the PFC is in good agreement with previous reports measured with fNIRS, and is analogous to the NBR that has been frequently observed in the fMRI field. While several explanations for NBR exist, understanding the decrease in $\mathrm{Hb}{\mathrm{O}}_2$ signals seen in our study remains to be further explored.