2.1.

Optical Signal Conversion

The first preprocessing stage involves the conversion of the measured light attenuation into relevant hemoglobin concentration levels. While it is possible to perform the conversion at a later stage,^{23, 26, 27} fOSA requires the conversion of the absolute light intensities from the cw measurement at any two distinctive wavelengths to $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ concentrations before proceeding on to other signal processing. This enables the user to inspect the hemoglobin data as a “first-pass” check on signal size, optode positioning and contact stability, etc. The conversion algorithm used in fOSA follows the modified Beer-Lambert law,²⁹ which considers the light scattering between each light source-detector pair to be homogenous both spatially and temporally throughout the study period. While the detailed algorithms for the conversion can be found elsewhere,^{21, 27} in summary the derivation of $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ concentrations (in micromolar of hemoglobin per liter of tissue) for each time sample can be simplified in matrix form if a minimum of two-wavelength measurement is made:

2.2.

fOSA Data Processing Procedures

In principle, the three possible forms of noise that can cause interference on the optical signals include the instrumentation noise, experimental noise, and physiological noise. The characteristics of the system noise are often well-behaved and can be removed easily, and the experimental noise can often be minimized with a properly designed task paradigm. However, very often the physiological noise (e.g., spontaneous vasomotion) overlaps very closely in frequency with the expected activation-related brain response.^{32, 33} The major emphasis of the filtering and curve fitting procedures in the data processing is to reduce non-task-related effects and to correct for any observable noise from the physiological measurements. Since the software is designed to include basic data correction techniques that then allow a real-time assessment of the activation signals, fOSA incorporates three different types of digital filters with distinctive roll-off and passband characteristics. While the elliptic filter gives a steep roll-off, the Chebyshev filter produces a gradual roll-off with fewer ripples in the passband, and the Butterworth filter produces a “maximally flat” passband with optimum roll-off. fOSA enables the removal of noise from the known frequency spectrum using either the low-, high-, and/or band-pass digital filters, resampling the data points and/or smoothing of the data by means of moving average. To view the range of frequencies in the physiological signal, its power spectral density can be plotted in the frequency domain using a Fast Fourier Transform (FFT). The FFT function in fOSA is defined using Welch’s averaged periodogram method by dividing the whole time course into five Hanning windows (i.e., the window duration is set as one-fifth of the signal duration). With no overlap and no zero padding needed, the sampling frequency (of the instrument) is set by the user.

Inappropriate detrending can introduce artificial peaks in the lower power spectrum that could be misinterpreted as physiological oscillations.³⁴ The fitting procedure implemented in fOSA uses a predefined rest period to obtain a curve that best fits its time course and then subtracts it from the remainder of the data points, where the level of the residual depends on the selective fitting order. Another effective way to cancel out the low-frequency noise is to perform a block averaging over the specific trials. The direct effect of this is to remove any uncorrelated high-frequency noise. Block averaging has very little effect on the removal of slow variations and the overall mean. These effects are usually considered when the data is detrended from the baseline.

fOSA offers the possibility to view the processed signals after each procedure. Since all the operations are linear, the sequence of the processing does not affect the overall results. Depending on the optode configuration, the results can be displayed in the form of topographic images (with linear interpolation between the measurement channels) or spatially configured hemoglobin time courses ( ${\mathrm{HbO}}_2$ , HHb, or both). It is also possible to focus on a specific measurement channel for more detailed analysis. An option to visualize the topographic Hb maps either at a particular time slice or in the playback mode is available where the whole time series of topographic images is displayed in a movie format. There are a variety of display options to interpret the statistical results (p-values and t-values), which can be viewed in the form of a topographic map or in a numerical table. To allow the user to carry out further analysis using other software, fOSA provides an option that allows the results to be exported as spreadsheets. Alternatively, fOSA offers the flexibility to process a batch of optical signals with the same experimental protocols, by applying the same processes to the group data.

A classical statistical test approach used to compare the difference in signal means between periods of rest and activation is to use a t-statistic that derives its probability of likelihood based on the standard deviation of the mean differences. Two different t-statistical approaches have been incorporated into fOSA: a Student’s t-test method typically used in the analysis of functional OT data^{26, 28} and the implementation of the SPM-OT method that uses the SPM approach to make inferences about the spatial coherence of the functional OT data. The former approach provides a quick assessment of the significantly active brain area, as each measurement channel is treated independently from its neighbor and the epoch design means that the task and rest periods are specified. Given the sluggish temporal nature of the hemodynamic response,³⁵ fOSA considers the temporal variation of the response function and provides the user with the option to specify the task and rest periods for statistical comparison. While more complex analysis involving multiple levels of comparisons that are required in some studies can be conducted outside the fOSA environment, the current software provides a single-factorial comparison by allowing the choice of one of the following options:

The open-source software means that it is also possible to modify the existing scripts and customize a user-specified statistical test.

2.3.

Statistical Parametric Mapping of OT Data

During functional activation, the neurovascular responses observed from the multiple brain regions are often spatially and temporally correlated in some way.²⁰ In the conventional statistical test approaches used in most optical imaging studies, each measurement point is treated as independent of its nearest neighbors, assuming that the regional activated areas are localized and independent of systemic influences. These tests also require prior knowledge about the hemodynamic response in order to specify the most relevant periods during which to compare the changes in the ${\mathrm{HbO}}_2$ and HHb and cannot easily account for the intersubject variability in the time course of the hemodynamic response. While it is possible to determine the time course of the hemodynamic response, modeling of spatial corrections is also required given that OT measurement channels are not spatially independent. SPM offers the flexibility to analyze sequential functional neuroimaging data by modeling the hemodynamic response and correcting for the spatial smoothness of the functional data.³⁶ The methodology uses the mass univariate approach where a statistic value is calculated for every voxel and the resulting statistics are then assembled into a map, i.e., an SPM. A detailed description of SPM theory and its algorithms is given elsewhere.³⁷ This paper focuses on the specific SPM procedures necessary for the analysis of OT data. The statistical approach adapted by SPM involves two basic steps. First, statistics reflecting evidence against a null hypothesis of no effect at each voxel are computed, where the effect is specified by a so-called contrast vector.³⁸ An image resulting from these statistics is produced. This image can be three-dimensional (3-D—fMRI, PET), or two-dimensional (2-D—M/EEG, OT). Second, this image is assessed, using p-values, by performing voxel-wise tests, to determine whether there is activation. While there are various approaches to control the family-wise error rate³⁹ of this mass-testing of many voxels, the random field theory (RFT) has been shown to be appropriate for the analysis of functional neuroimaging data.⁴⁰ The Gaussian random field treats the statistical images as sampled versions of continuous, spatially resolved data. P-values, adjusted for multiple comparisons, are computed by estimating the probability that some peak in an image surpasses some user-specified threshold, given the null hypothesis of no activation. For typical, smooth images, it has been observed that the RFT provides much more sensitive tests than the Bonferroni correction. This is because the RFT is explicitly informed about image smoothness (i.e., spatial correlations), while the Bonferroni correction is not. It might be asked whether the RFT is applicable to OT data. To the best of our knowledge, OT data fulfill all the assumptions needed for the RFT. The results of the RFT are asymptotically valid (for higher and higher thresholds) for any type of correlation structure, as long as the local correlation is strong. (We arrive at this conclusion through personal communication with Tom Nichols and Stefan Kiebel; see also below.) This holds for OT data. Similar discussions about its applicability exist in the M/EEG field, where we observe that no definitive review has yet been written on the applicability of RFT to M/EEG data. (M/EEG is similar to OT data in that it is a spatial mixture of the expressions of underlying brain sources.)

2.4.

SPM Modeling and Inference

In the following, we will go through some basic techniques that SPM uses to model spatiotemporal data, by which we mean that both temporal and spatial correlations are modeled. The main point is that the assumptions that the SPM procedures involved hold not only for fMRI, PET, and MEG/EEG but also for OT data.

SPM involves the construction of a statistical model and making inferences about these parameter estimates based on the spatiotemporal neuroimaging data and the predicted model. In order to capture the sluggish nature of the neurovascular response, SPM treats each scan as a dependent observation of the others and models the time-continuous data for each voxel. Here, each of the brain activation voxels is expressed as a combination of different basis functions that best represent the functional response using the general linear model (GLM). In simpler terms, the GLM expresses the observed response $Y$ as a combination of the explanatory variables $X$ plus a well-behaved error term $\epsilon$ (i.e., $Y=X\beta +\epsilon$ ).²² Translating the model into matrix form, each row in $Y$ is treated as a scan of time-series data. The design matrix $X$ is where the experimental knowledge about the expected signal is quantified. The matrix contains the regressors (e.g., designed effects or confounds), where each row represents each observation and each column corresponds to some experimental effects, including effects that are of no interest and pertain to confounding effects (e.g., physiological effects) to our analysis that are modeled explicitly. The associated parameter $\beta$ in each column of $X$ indicates the effect of interest (e.g., the effect of a particular cognitive condition or the regression coefficient of the hemodynamic response). The convolution model for the hemodynamic response takes a stimulus function encoding the supposed neuronal responses convolved with a defined canonical hemodynamic response function (HRF) to form a regressor that enters into a design matrix.²⁰ SPM also allows the modeling of latency and dispersion derivatives as additional regressors to its canonical HRF in order to treat the variability of hemodynamic responses due to different sorts of events. Other forms of explanatory data (e.g., blood pressure and/or heart rate measurements) can be defined as additional regressors in the design matrix. Unmodeled effects are then treated as residual errors. Often, the number of parameters (column) is less than the number of observations (row); hence, some method of estimating parameters to find the best fit for the observed response is needed. In SPM, the parameter estimation is achieved using a mixture of restricted maximum likelihood (ReML) and ordinary least squares (OLS) methods. In estimating the error covariance matrix, SPM assumes that the pattern of serial correlations is the same over all voxels (of interest) but that its amplitude is different at each voxel. To incorporate nonspherical error distributions, SPM uses the ReML procedure to estimate the model coefficients and error variance iteratively.⁴¹ In essence, the ReML method deals with linear combinations of the observed values whose expectations are zero. The relationship between experimental manipulations and observed data may consist of multiple effects, all of which are contained within the design matrix. To test for a specific effect, a contrast is defined to focus on a particular characteristic of the data. SPM allows the specification of different contrast vectors to the same design matrix to test for multiple effects without having the need to refit the model (i.e., univariate linear combinations of the parameter estimates).

Regional changes in error variance and spatial correlation in the data can induce profound nonsphericity in the error terms (as opposed to the assumption of identically and independently distributed error terms). This nonsphericity would require large numbers of parameters to be estimated for each voxel using conventional multivariate techniques. However, using SPM, the parameterization is minimized to just two parameters for each voxel—error variance and smoothness estimators. This is made possible because SPM uses the RFT to resolve the multiple comparison problem, which then implicitly imposes constraints on the nonsphericity implied by the spatially extended data. SPM improves the sensitivity of the statistical test by considering the fact that neighboring voxels are not independent by virtue of spatial correlations in the original data and employs the RFT approach to adjust the p-significance. RFT correction often expresses its search volume as a function of smoothness (or resolution elements that define the block of correlated pixels). It relies on the expected Euler characteristic (EC) that leads directly to the expected number of clusters above the threshold, and we want to derive this statistic height threshold once the smoothness has been estimated. A relatively straightforward equation on the expected EC to derive this threshold is given by Worsley that depends on the $\alpha$ -threshold and number of resels contained in the volume of voxels under analysis.⁴² In practice, the correction will need to take into consideration the search volume. Restricting the search region to a small volume (i.e., defining the region of interest) within the statistical map is likely to reduce the height threshold for given family-wise error rates. The smoothness is in turn calculated using the residual values from the statistical analysis. Kiebel have shown that it is possible to estimate the smoothness based on the residual components in the GLM with sufficient degrees of freedom $(>20)$ and the size of the spatial filter kernel.⁴³ Using the same analogy, the voxel-based inference can be extended to a larger framework involving cluster-level and set-level inference.

2.5.

SPM-OT Analysis

The fact that both the fMRI and OT methods measure the temporal neurovascular changes induced by the increased oxygen demand of the neuronal activities and hence share a similar variation of the hemodynamic response makes it possible to utilize some of the existing functions in SPM designed for the analysis of fMRI data and apply them to the treatment of the optical signals. In addition to implementing the signal-processing features in fOSA, a conversion program named SPM-OT has been incorporated into fOSA that allows the OT data processed in fOSA to be analyzed in SPM like other neuroimaging modalities.⁴⁴ The approach used by SPM-OT to analyze the OT data is simple:

The third option produces images where the two-dimensional (2-D) measurement patches are embedded into a 3-D head space. This is important for co-registration of the OT analysis with other modalities like fMRI or EEG. Similar SPM approaches to analyze 2-D data embedded into 3-D space have been described for fMRI.^{45, 46} However, this is not the approach we take in this paper. We simply analyze the OT data on a 2-D patch, which can then post hoc be co-registered to brain space.

2.5.1.

Spatial preprocessing of OT data

A series of 2-D brain maps are generated for each time bin, where each measurement point is represented by a pixel corresponding to the ${\mathrm{HbO}}_2$ and/or HHb parameter of the OT data. Figure 1 shows one possible statistical map based on the standard configuration of the 24-channel Hitachi OT system (Fig. 2 ). To improve on the sensitivity and facilitate intersubject comparisons, an interpolation process using nearest neighbors is used to produce a finer grid. As the interoptode spacing (IOS) is often larger than the spatial resolution afforded by the measurement system (the Hitachi OT system has a fixed IOS of $30\mathrm{mm}$ ), it is more appropriate to break down each channel into smaller pixels to better identify any focal change in the hemodynamics. It is important to consider the spatial information that could be revealed using the existing OT technique and a smoothed map of the brain activity to facilitate the RFT correction for multiple comparisons. The position of each optode can be measured using an electromagnetic tracking device that provides the 3-D coordinates of the optode pixels (in red and blue, as shown in Fig. 1). The same interpolation technique can be applied to derive the 3-D positions of the channels (in yellow), where each point is assumed to be at the center between each source-detector pair.⁴⁷ To mask out any irrelevant pixels (outside the two squares as shown in Fig. 1), a “not-a-number” (NaN) value is assigned on the configured grid. This allows SPM to identify regions that are not spatially correlated during the analysis. In considering the typical neurovascular response, SPM-OT provides an option to resample the optical data to a level similar to that of the repetition time (TR) used in fMRI. Each SPM brain map is tagged sequentially, and it is vital that the images are generated in a correct sequence, as each of these maps is being selected based on its header information during the specification of the design matrix.

Fig. 1

Schematic diagram of the $3\times 3$ configuration of optodes from the Hitachi ETG-100 with 10 sources and 8 detectors. We derive a total of 50 pixels, which includes 18 interpolated pixels (in red and blue) and $24+8$ measurement points (in yellow). The surrounding pixels are padded with NaN (not-a-number) values so that SPM can perform a masking operation before the estimation procedure. (Color online only.)

Fig. 2

Experimental arrangement of the optode positions over the C3 (left) and C4 (right) motor areas for the finger tapping task.

2.5.2.

Modeling of the spatiotemporal OT data

Previous fMRI-OT studies have shown some similarities between the BOLD and ${\mathrm{HbO}}_2$ signals,^{9, 10} and given the better signal-to-noise ratio typically seen in ${\mathrm{HbO}}_2$ signals during functional activation, the signal has often been used as an indicator of regional cerebral blood flow.⁴⁸ As the NIRS measurement is able to provide a higher sampling rate than fMRI, in principle, it allows a wider range of the oscillatory spectrum to be captured. Specifying these measurable physiological oscillations in the design matrix enables a better modeling of the temporal hemodynamic response in SPM. During functional activation, the hemodynamic response of the ${\mathrm{HbO}}_2$ signal resembles that of a gamma function with a time to peak typically between 5 to $8\mathrm{s}$ and a slight undershoot before and after the bulk response.^{49, 50} A similar response function has been observed during functional activation in the fMRI-BOLD signal that is used in SPM as a temporal basis function for the modeling of the brain response.²⁰ Figure 3 shows the representation of the canonical hemodynamic response function (HRF) in SPM. To accommodate the sluggish feature of the HRF and the variability between different voxels due to different events, SPM allows a combination of different basis functions to model the observed response from the OT data with the inclusion of temporal and dispersion derivatives to give additional weightings to the canonical HRF.⁵¹ An advantage of using the OT technique to map the neurovascular response is that it is able to monitor both the $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ signals simultaneously. Hence in the analysis of OT data, SPM-OT offers the possibility to model the two chromophores (and possibly $\mathrm{\Delta }\mathrm{HbT}$ ) in the GLM. By assuming that the two responses take the form of a gamma function but in different directions during functional activation, SPM-OT offers an alternative route to infer about the functional OT data. Using the same SPM procedures in the model specification, Fig. 4b shows the inverse canonical HRF prescribed in SPM-OT to model the $\mathrm{\Delta }\mathrm{HHb}$ signal in an experimental data set. A positive t-statistic from the SPM analysis will hence indicate that the particular pixel location is significantly deactivated. However, it should be noted that while this option offers an alternative to compare each of the chromophores between different conditions, an inference between the two chromophores will not be possible since the response functions in these signals are thought to be different, with each function having different time latencies and amplitude changes.³⁵ These differences would require explicit (nonlinear) models that are not available with the linear methods in SPM.

Fig. 3

A combination of the three gamma functions (canonical HRF and temporal and dispersion derivatives) characterize the sluggish nature of the hemodynamic response function used to model the BOLD signal in SPM-fMRI analysis.

Fig. 4

(a) Design matrix for the deoxy-hemoglobin model for the $l$ -ftap (left finger-tapping, first three columns) and $r$ -ftap (right finger tapping, next three columns) tasks, and a constant error term for the design matrix. Each row represents a time slice of the topographic brain map. (b) The inverted canonical HRF (blue) plus temporal (green) and dispersion (red) derivatives model for the deoxy-hemoglobin over the five repeated blocks for the $l$ -ftap task (corresponding to the first three columns in the design matrix). (Color online only.)

Unlike in the SPM analysis of fMRI data where a global normalization is recommended to threshold the BOLD response from the extracerebral volume, this procedure is not necessary for the analysis of OT data. SPM estimates the size of the experimental effects from the residuals with the combination of the maximum number of explanatory regressors using restricted maximum likelihood (ReML) for each voxel, in which one or more hypotheses can be specified. The associated parameter estimates are the coefficients that determine the mixture of basis functions that best model the HRF for the voxel in question. A classical univariate test statistic can then be formed to examine each voxel, and the resulting statistical parameters are assembled into an SPM image. Here a within-subject effect is specified in the first-level analysis, and the intersubject comparison is carried forward to the second-level analysis. A precheck on the variance between optode positions is often necessary in the event of a group analysis.

3.1.

Comparison Using Simulated Data and Results

Validating algorithms on experimental in vivo data is difficult since there is an uncertainty over the true chromophore concentration changes; therefore, a set of simulated data was used to validate the software.⁵⁴ The artificial data set consists of a time series of Hb values generated using known values for the tissue absorption and scattering coefficients together with a diffusion theory model to calculate the expected light attenuation.⁵⁵ To generate the time series data, we have taken measured in vitro spectra of the ${\mathrm{HbO}}_2$ and HHb and combined these in time varying concentration levels to simulate the different experimental conditions.⁵⁶ Figure 5a (top) shows the simulated data together with an imposed linear drift and a known range of simulated noise frequencies between 0.1 and $1\mathrm{Hz}$ . The duration for both the task and rest periods were set at $20\mathrm{s}$ . Figure 5a (bottom) shows the simulated HHb concentration varying between 0 and $-5\mu \mathrm{mol}$ and ${\mathrm{HbO}}_2$ concentration varying between 0 and $10\mu \mathrm{mol}$ .

Fig. 5

(a) Simulated ${\mathrm{HbO}}_2$ (magenta) and HHb (green) signals compared with the ${\mathrm{HbO}}_2$ (red) and HHb (blue) results from fOSA after conversion and detrending. (b) The ${\mathrm{HbO}}_2$ (red) and HHb (blue) results from fOSA after filtering compared to the original simulated signals. The major differences are due to the effects of filtering and disappear if the comparison is made to a filtered version of the original signal. (Color online only.)

The results showed good agreement between the simulated data set and the fOSA results. There was a small offset between the two data sets that could be due to the wavelength-dependent effects of scattering and that is not currently taken into account in the fOSA application of the modified Beer-Lambert law (although this could be added by the user if desired). By applying a fifth-order Butterworth low-pass filter at $0.08\mathrm{Hz}$ and a linear curve fitting, we were able to effectively eliminate the simulated noise in this example. Figure 5b illustrates the comparison between the simulated and processed data after the averaging process. The slight difference in waveform shape is due to the effect of the filtering. Although not shown here, both signals share the same wave forms if the comparison is made with a similarly filtered version of the simulated signal. This test validated the correctness of the algorithms used to convert the light intensity signals as well as some of the basic signal processing routines.

3.2.

Comparison Using Experimental Data and Results

The application of the fOSA software was tested using published experimental data.²⁶ In this study, left ( $l$ -ftap) and right ( $r$ -ftap) finger-tapping tasks were conducted in two separate trials. The $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ were measured using a 24-channel OT system (ETG-100, Hitachi Medical Corporation, Japan).⁵⁷ Figure 2 shows the $3\times 3$ optode configuration used in the study, where the optodes were placed close to the sensorimotor cortex (C3 and C4 areas). In summary, the subject was guided to perform the finger-tapping tasks at a regular pace of $3\mathrm{Hz}$ , over a task period of $30\mathrm{s}$ and a rest period of $30\mathrm{s}$ . Each trial was repeated five times.

Figures 6a, 6c, 6e show the time course results for the $l$ -ftap task over channel 16, and Figs. 6b, 6d, 6f show the time course results for the $r$ -ftap task over channel 9. By applying the algorithms employed in the fOSA software, the two-wavelength light intensities were converted to relative $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ concentrations. Figures 6a and 6b show the $\mathrm{\Delta }\mathrm{Hb}$ concentrations after smoothing the data with a moving average with a window span of 10 samples and the data were downsampled to $1\mathrm{Hz}$ . Figures 6c and 6d show the concentration changes after modeling out the slow drift using a linear fit. Last, the data were averaged over the five repeated blocks to eliminate any existing high frequency components, as shown in Figs. 6e and 6f.

Fig. 6

Data from the study reported (Ref. 26) (data from subject 3). Channels 9 and 16 have been selected for illustration for the $r$ -ftap and $l$ -ftap tasks, respectively. (a) and (b) show the $\mathrm{\Delta }{\mathrm{HbO}}_2$ (red) and $\mathrm{\Delta }\mathrm{HHb}$ (blue) time course for the $l$ -ftap and $r$ -ftap tasks, respectively, after smoothing the data using the “moving average” technique. (c) and (d) show the time course results after linear fit for the tasks, and (e) and (f) show the averaged signals over the five repeated blocks for each task.

The results of the analysis on the experimental time course data processed using fOSA can be compared to those obtained using the analysis methods described by Sato ²⁶ in the study from which the data set was taken. In general, the hemodynamic response was observed in both analyses, i.e., an increase in $\mathrm{\Delta }{\mathrm{HbO}}_2$ with a slight decrease in $\mathrm{\Delta }\mathrm{HHb}$ after the stimulus onset, with the increase in $\mathrm{\Delta }{\mathrm{HbO}}_2$ maintained throughout the activation period. However, there is a slight variation in the time course that could be due to the processing procedures used. First, intensity conversion and block averaging were the only preprocessing steps done in the original experimental study by Sato. fOSA performed additional preprocessing steps: linear detrending and moving average to remove the slow drift and high-frequency components. Second, Sato used a “peak search” technique to select a $25\text{-}\mathrm{s}$ maximum absolute mean change (i.e., maximum averaged area under the curve within a $35\text{-}\mathrm{s}$ window). In their analysis, a paired t-test $(p<0.1)$ between the predefined rest and selected activation periods was conducted over the five repeated blocks to determine the significance for each task.

Since the experimental data was measured using the 24-channel Hitachi ETG-100 system, the optodes configuration as shown in Fig. 1 was used to generate the spatial map for the SPM analysis. As the two tasks were conducted in separate trials, SPM-OT first combined the data sequences as one continuous vector for $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ concentration data. Note that SPM provides a flexible approach to estimate the experimental data; it is able to accommodate different data length with varying stimulus onset periods in its model. A canonical HRF (presented in SPM) plus its latency and dispersion derivatives were used as the basis functions for the $\mathrm{\Delta }{\mathrm{HbO}}_2$ response model. To illustrate the flexibility of the software and validate the activation response, a second similar but inverted canonical HRF plus its two derivatives were used for the $\mathrm{\Delta }\mathrm{HHb}$ response model. Figure 4a shows the specification of the design matrix. The $x$ axis corresponds to the six regressors for the two tasks plus a constant error term, and the $y$ axis represents the 700 topographic images generated by SPM-OT. A t-statistic was computed for each of the $l$ -ftap [Figs. 7a and 7c ] and $r$ -ftap [Figs. 7b and 7d] tasks. The reported t-results were p-corrected using family-wise error correction $\left({\alpha }_{0.05}\right)$ and analyzed at the “pixel” level. Statistical results from the analysis showed a significant effect over the right cortex at $t_{oxyHb\text{-}ch16}=19.02$ $({\mathrm{p}}_{\mathrm{corr}}⩽0.001)$ and $t_{deoxyHb\text{-}ch16}=12.79$ $({\mathrm{p}}_{\mathrm{corr}}⩽0.001)$ when the subject performed the $l$ -ftap task and at $t_{oxyHb\text{-}ch8\text{-}9}=17.26$ $({\mathrm{p}}_{\mathrm{corr}}⩽0.001)$ and $t_{deoxyHb\text{-}ch8\text{-}9}=5.68$ $({\mathrm{p}}_{\mathrm{corr}}⩽0.001)$ on the left cortex when performing the $r$ -ftap task. These results showed good agreement in terms of laterality with the published results²⁶ and previous OT studies.²⁷ To illustrate the significant activated regions on the left sensorimotor area by estimating the optode placement (not measured in the study), Fig. 8 shows the possibility of extracting the 2-D t-results from SPM and superimposing them onto a SPM brain map.

Fig. 7

(a) and (b) show the SPM t-results for the $l$ -ftap task for the $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ signals, respectively. (c) and (d) show the SPM t-results for the $r$ -ftap task for the $\mathrm{\Delta }{\mathrm{HbO}}_2$ and $\mathrm{\Delta }\mathrm{HHb}$ signals, respectively. Darker pixels correspond to higher significant t-values. A positive t-score for the ${\mathrm{HbO}}_2$ $\mathrm{signal}=\mathrm{significant}$ increase in the ${\mathrm{HbO}}_2$ signal, and a positive t-score for the HHb $\mathrm{signal}=\mathrm{significant}$ decrease in the HHb signal.

Fig. 8

The SPM t-results superimposed onto a canonical brain map showing the effects on the left cortex from the $r$ -ftap task.

3.3.

Application with Other OT Systems

fOSA has been designed to process and analyze data from both fixed and flexible array systems and is not restricted to use with the ETG 100 Hitachi system. To date, the software has so far been used to process OT data from two other cw systems—the ETG 4000 (Hitachi Medical Systems; 24 channels, wavelengths $695\mathrm{nm}$ and $815\mathrm{nm}$ )²⁸ and an in-house cw OT system (84 channels, wavelengths $775\mathrm{nm}$ and $850\mathrm{nm}$ ).⁵⁸ Figure 9 shows the optode array design of the in-house system. Eight sources and eight detectors are positioned to provide five different optode spacings ( $14.3\mathrm{mm}$ , $17.8\mathrm{mm}$ , $22\mathrm{mm}$ , $29\mathrm{mm}$ , and $34\mathrm{mm}$ ) with 42 channels in each array. This distribution of optode spacing is designed to provide some depth discrimination in the measured signals.⁵⁹ Software extensions to read data from other OT formats, including the OXYMON 24-channel OT system (Artinis Ltd; http://www.artinis.com) and Hamamatsu OT system (Hamamatsu Photonics KK) are under constant development.

Fig. 9

Array design of the distributed sources and detectors for the UCL in-house optical imaging system.