1. Introduction

Diffuse Optical Imaging (DOI) has seen increased interest in recent years [1–5]. The emergence in the mainstream of this new modality has opened the way for a non-invasive method to image hemodynamics in vivo. The 600–1000 nm spectral range has low absorption thus allowing detection of photons diffusing and traveling through several centimeters of tissue.

As with other modalities, the first challenge faced is to separate the noise and physiological signal from the signature of interest. This is the focus of the present paper. Once this separation is done, tomography can also being pursued by using models of light propagation, but even with these models the unknown underlying absorption and diffusion heterogeneities lead to diminished spatial resolution and quantification. When applied to brain imaging, significant differences emerge when compared to functional Magnetic Resonance Imaging (fMRI): the acquisitions are done with much better temporal accuracy (10–100 Hz) and oxygenated hemoglobin can also be measured. As such physiology appears very prominently in diffuse optical data and can be better identified since the acquisition is typically faster diminishing aliasing artifacts. This separation of physiology from the hemodynamics originating from a given task is a challenge that has been solved in mainstream applications, by two main techniques. The first is block averaging whereby high frequency physiology as well as low frequency drifts are filtered out of the data by carefully keeping frequencies on which the task, convolved with a prototype response, has support. An average over each block is then done on this filtered data thereby synchronizing the activation and potentially eliminating other slow physiology (such as Mayer waves). The second technique uses either block or an event related experiment and a General Linear Model (GLM) where drifts are introduced in the model as regressors and modeled. Again, depending on the number of drifts used, the data can then be filtered by all regressors not correlating with the protocol [6].

More recently, these analysis methods have been refined in order to include the effects of non-stationary physiology on the estimators. For example, blind Principal Component Analysis (PCA) based filtering [12] has been proposed to reduce physiological signals. Separately, Kolehmainen et al. [9] have developed a method to perform dynamical state space estimation. A quasi-stationary technique is described in Prince et al. [11] by fitting sinusoids amplitudes and phases to model cardiac and respiratory nuisance signal. This work is extended in Diamond et al. [7] where physiological regressors, measured during the experiment, are included in the analysis. This recent work confirms that modeling physiology in the estimation process can be beneficial.

In all these cases, the goal is to provide an estimation of the response to a given task independently of the presence of physiology. A question remains however: is working in a time based referential (i.e. looking at signals in time) the best method to differentiate physiology? The finer techniques developed above may benefit from the usage of a basis that better separates physiology and evoked responses prior to performing estimation. As such it may be interesting to have a technique and associated basis that would not only estimate the responses, but also open the door for an understanding of the physiology. For example, it is expected that with the occurrence of the task, physiological sources alter their behavior and obscure the response estimation. This is the basis for dynamical modeling. Here a time-frequency analysis seems appropriate given that the different phenomena appear at different time and frequency scales. In this work we want to lay the bases for a time-frequency approach to the model proposed in [7]. In order to do so, we first need to establish the technique using more standard static analysis.

The first goal of this paper is to study optical signals within a continuous multi-resolution framework, and compare it with physiological signals measured in parallel [7]. The second goal is to explore the information originating from these wavelet representations, complex in our case, and see if it can help in characterizing the evoked response better. In doing so, we develop a new method for performing the estimation of the Hemodynamic Response Function (HRF) based on synchrony. Some issues concerning the oxy-hemoglobin (HbO) signal and physiology are addressed and we show how synchrony may help in filtering out physiology from the DOI data. Actually, the synchrony alone can be also used for defining a contrast between two conditions in a paradigm.

The paper is organized as follows. The next section introduces the notion of complex wavelets and applies it in Section 3 to physiological data. We then promote the wavelet based estimator introduced and relate this new approach with a GLM applied directly in the wavelet space. In section 4, we develop new approaches based on wavelet representation of signals and synchrony to extract activation information from diffuse optical data. First we discuss and evaluate performance with simulations. We then apply the method in two experiments with visual and motor paradigms.

2. Analytical wavelets

A continuous wavelet analysis consists of expanding a signal in terms of localized oscillations of different scales. The set of those oscillating modes is defined from a generative mother wavelet for which an infinite number of choices exist. The wavelet transform amounts to computing the projections of the signal s(t) over the whole set of wavelet functions defined by

where a and b are scale (dilation) and time (translation) parameters respectively. In the present work, we consider the Lusin wavelet [13] defined by

where n is any positive integer. The reason for choosing this wavelet will become clear in the following and is related to its most important property, its analyticity, shown by observing that the Fourier spectrum has a support strictly on the positive frequencies:

The parameter n in ψ_n is related to the number of vanishing moments of the wavelet. The wavelet is blind to polynomials of degree less or equal than n. When analyzing data, this property is useful since it removes polynomial drifts, leaving only the relevant fluctuations in the time-frequency analysis. A drawback is that a larger n leads to a wavelet having larger frequency spread as can be seen from the previous expression. This last expression is also used in the numerical implementation of the wavelet transform. For completeness, the Fig. 1 displays the wavelet for n=7 in the time domain as well as in the frequency domain. We observe the phase quadrature between the imaginary and real parts reflecting the analyticity of the wavelet ψ_n(t), the imaginary part of the wavelet is the Hilbert transform of the real part. The analytic wavelet allows the unambiguous definition of a phase of the complex wavelet coefficients which can be interpreted as an instantaneous phase.

 

Fig. 1. Lusin wavelet ψ₇ in the time domain. (A): real part (solid line) and imaginary part (dashed line). (B): ψ₇ in frequency. w₀ denotes the central frequency of the wavelet.

Download Full Size | PDF

The complex-valued wavelet coefficients of a real signal s(t) are given by the projections over the whole set of wavelets,

defined in the time-scale plane (a,b). Here and in the following, the bar will be used to denote the complex-conjugate values. Scales and frequency bands are in one-to-one correspondence through the definition of the wavelet. We use ω≈ω ₀/a, where ω ₀ is the central frequency of ψ _n (see Fig. 1) instead of the scaling factor a. The complex-valued wavelet coefficients can be displayed in a time-frequency plane, either in terms of their quadratic amplitudes (|w|²) or normalized amplitudes (a ^-1|w|²). This is the so-called scalogram of the signal. Besides this usual picture of the wavelet transform, which takes only into account the amplitude of the wavelet coefficients, we can also consider the sign of the coefficients. Let us define the sign of the complex wavelet coefficients as a unit vector in the complex plane:

In the following we will use the term sign for the previous expression ${\gamma }_{a,b}^{\left(s\right)}=e^{{i\varphi }_{\mathrm{ab}}}$ and also refer to the phase, ϕ_ab, to describe the phase synchrony. Essentially, the signs are the complex unit vectors spanned by the local phases of the signal. At each scale a and time b, this quantity represents the instantaneous phase of the signal. This complex vector will play an important role in the forthcoming time-frequency analysis of the DOI signals.

In the example of Fig. 2 (A,B) we show the scalogram of a signal composed of two time-shifted gaussian functions modulated by two harmonics at 0.12 Hz and 1 Hz. The vertical scale axis is conveniently converted to frequencies. The spectral range explored with the wavelet can be determined in terms of the sampling rate (upper bound frequency) and duration (lower bound frequency) of the signal. This example illustrates the relative accuracy of the wavelet localization of pure harmonics. The second example in Fig. 2 (C,D) displays the time-frequency representation of a theoretical hemodynamic response to a stimuli. This signal (hereafter denoted by HRF) is commonly modeled as a convolution between the stimuli paradigm B(t) and the canonical impulse hemodynamic response h(t) (SPM package www.fil.ion.ucl.ac.uk/spm/),

In both examples, the patterns displayed in the time-frequency plane localize the pieces of wave in both temporal and spectral dimensions. Notice that because of the fundamental uncertainty principle (see [19] for instance), simultaneous accuracy in both time and frequency cannot be achieved and may depend on the choice of the wavelet.

The wavelet transform is a redundant representation of the signal. Indeed, the time-frequency plane contains much more coefficients than originally provided in the sampled signal. Nevertheless, the reconstruction of the signal is possible [13,19,21]. In the present work, we use the Morlet reconstruction formula [20]:

where k_ψ is a normalization constant that will be estimated with a control signal.

To summarize this introduction, we emphasize that the wavelet transform provides a time-frequency representation of a signal containing information about local instantaneous phases, provided that the wavelets are analytical. It is further possible to synthesize a signal from this complex representation. However, because of the vanishing moments of the wavelet, this synthesis will not exactly reproduce the initial signal since polynomial content up to degree n are filtered out. In the application presented in this work, this implicit polynomial filtering is in fact beneficial since it removes undesirable drifts.

 

Fig. 2. Left column: Two modulated gaussians (A) at 0.12 Hz and 1.0 Hz and the corresponding scalogram (B) with ψ ₇. Note that the wavelet transform provides both temporal and frequency localizations of the two constituents of the signal. The figure displays amplitudes of the coefficients and the cone of influence. The wavelet coefficients outside the cone (close to the boundaries of the signal) are tainted by artifacts due to edge discontinuities of the signal. Coefficients inside the cone are sufficiently far from the edge to be free of these artifacts. Right column: A wave train (C) given by equation (6) for some block paradigm B and its time-frequency representation (D).

Download Full Size | PDF

3. Wavelet transforms of DOI measurements and physiological data

In order to assess whether the wavelet transform defined above is relevant to optical imaging, we first look at physiological data and correlate it with the measured DOI data in the complex wavelet domain. Then, moving to a motor task within the same framework, these preliminary observations will lead us to define a new technique to estimate the HRF.

Let us first recall the preprocessing stage of the DOI data that was acquired here using a CW Near Infrared (NIR) imager commercially available (Techen CW5). The optical probe is positioned on the scalp of the subjects and specific source-detector pairs (hereafter called optode pair) are kept for analysis. In general, we consider closest source-detector pairs to make sure there is significant SNR (see for instance Fig. 4 (A)). The raw data generated by the instrument was demodulated and downsampled to 10 Hz. In each of the cases presented below, no further filtering was done and the photon fluence at wavelength λ was then converted to variation in optical density through the Beer-Lambert law

Up to path length factors, a concentration measure is then obtained by using the extinction coefficients that characterize the two chromophores:

In the present work, we used λ ₁=690 nm and λ ₂=830 nm [15]. Although not applicable to tomographic imaging, the Beer-Lambert law is extensively used in the literature. Moreover since we are only interested in changes during a given task, one can argue that changes in the signal are added linearly (e.g. Born or Rytov approximation). Each change in concentration is thus computed for selected optode pairs.

3.1. Etiology of observed physiology

Before looking at evoked DOI signals, we consider data recorded while the subjects are at rest (10 subjects, data from [7,8]). This set of data not only contains optical signals originating from 3 pairs, but also physiological information measured in parallel. These consist of the heart beat, the instantaneous heart rate, the respiratory signal and the blood pressure.

By using the analytical wavelet defined above with n=7 (found to be a good compromise for exploring the spectral content of the signals), we can compare the time-frequency planes of physiological signals and variations of concentration in oxy-deoxy-hemoglobin, averaged over all subjects and optode pairs. This is done in Fig. 3. The temporal window is chosen so that boundary effects at the beginning and end of the acquisitions are absent. The frequency range is computed from the data and depends both on the acquisition rate (reduced to 10 Hz so we explored up to 2 Hz), and the length of the measures (297 seconds, i.e. the lowest frequency is 0.01 Hz). The results reveal some correlations that are expected. The observation of the heart beat (~1 Hz) in the HbO signal, a potential correlation between blood pressure and HbR. The most interesting thought is the heart rate which, when represented in the time-frequency plane, correlates with low-frequency oscillations, the Mayer waves [22], observed in the HbO signal.

3.2. Example: motor task

We have shown that various physiological compartments can be distinguished in the time-frequency plane. Let us now consider how neuronal activity is mapped in the same plane. Given that the hemodynamic response to stimuli is slow, we expect that only a few ranges of frequencies will be sensitive to the task and this is the underlying reason for filtering data. The interest here is to observe the interplay between the representations of the physiology and the task in measures of HbR and HbO. As we will move further in our analysis, we will find those observations helpful.

The task protocol used for this example is as follows. A single subject was asked to rest for the first 4 min, then to perform a 7s left hand finger movement every 14 s for a duration of 4 min, followed with a second rest period of 4 min. The instruction was presented 2 s before the first stimulus. Stimuli were indicated by displaying a red square on a screen. The activation in this case was expected to occur on the right motor cortex over which optode pairs were placed. The optode configuration is shown on Fig. 4(A). The results for HbO and HbR of some pairs, in particular the ones having the strongest response located over the motor cortex, are also shown in Fig. 4. In this figure we observe that the evoked response is represented in the scalogram. It is seen that without averaging or filtering, the response is present on the HbR channel but much less defined in the HbO channel. In the latter, we see in the same frequency band the strong presence of Mayer waves. Although there is higher intensity in HbO when the task is present, it is also seen that other contributions appear in the same frequency band for most of the pairs. This raises questions as to how HbO can be used as a signature for activation without removing those physiological contributions.

 

Fig. 3. Time-Frequency plane for HbO, HbR and four physiology signals at rest, group averaged wavelet power over 10 subjects. HbO and HbR signal power spectrum are averaged over 3 optode pairs. HbO channel exhibits theMayer waves around 0.1 Hz. A similar effect is seen in HbR channel but is in fact much smaller: the graphs are normalized individually and the HbR signal around 0.1 Hz is 20 times smaller than for HbO. We also observe a similar pattern in the time-frequency plane of the heart rate. The cardiac beats are visible in the HbO signal as well when compared to the heart beat. The respiratory signal does not show strongly in the optical data but its effect is expected to be smaller.

Download Full Size | PDF

The previous observations support the assertions that wavelets may help in analyzing optical data and distinguishing physiology from activation. As such, the observations are compelling but not really surprising, we know that the frequency support of different physiological phenomenon will enable their distinction. As described below, this approach becomes more quantitative when we exploit further properties of the analyticity of the wavelets.

3.3. Synchrony and phase-locking maps

Phase synchrony between two signals can be quantified with the circular average of the phase lag between them. It is formally defined as the amplitude of the time average of a complex random value (〈·〉 denotes the time average)

where ϕ_i and ϕ_k are the instantaneous phase of the signals s_i(t) and s_k(t) respectively. This quantity is between 0 and 1 and can be evaluated in various ways, but analytic wavelet analysis of the signals gives a natural framework to compute it at a particular frequency band and time. Indeed, the previous expression can be written in terms of the sign of the complex wavelet coefficients:

Here, the temporal average is done with N_T samples in the window of size T around t=b. The closer this value is to one, the more in-phase the signals are around time t=b. As an example, phase synchrony can be used to better quantify the correlations between DOI and physiological signals. This is done on Fig. 5 where we consider couplings between HbO (HbR) and other physiological measurements for the same set of data as in Fig. 3. We observe that, as expected, the HbO channel is tightly coupled with the cardiac dynamics and more specifically with the heart rate related to the Mayer waves [23]. Such coupling is less apparent in the HbR channel. This supports the observation that the stimulated hemodynamic response in the previous motor task experiment is more preeminent in the HbR signal than in the HbO channel, the latter being mostly dominated by cardiac and Mayer waves.

 

Fig. 4. Motor task (finger-tapping) paradigm, single subject. (A): optodes configuration composed with 4 sources and 8 detectors. (B): scalogram of the stimuli convolved with the SPM2 canonical hemodynamic response (expected response in the data). The time axis has been truncated to eliminate the cone of influence. The six other panels show the normalized scalograms (the single colorbar applies to all graphs) of the wavelet transforms for HbO (left) and HbR (right) channels at 3 selected optode pairs (4,9 and 13). Notice the clear activation in the HbR channels of optode pair 4, in correspondence with the stimuli and the presence of Mayer waves in the HbO channels somewhat obscuring the presence of evoked activation. Optode pair 13 has no activation but still sees the Mayer waves.

Download Full Size | PDF

The question is whether this information, which is independent of the signal amplitude, can be used to get better estimations. When looking at the case of evoked responses to a given stimulus, we can quantify the synchrony measures on these evoked signals. Consider for example the synchrony between the theoretical hemodynamic response (HRF) and the stimulus in one optical channel s_i(t),

where the brackets denote here a time average taken over epochs of repeated stimulations. Assuming the stimulus is a periodic sequence of N_e blocks, we can assume, if the response is linear, that γ ^HRF _a,b=γ ^HRF _a,b+mΔ where Δ is the period of stimulation and m indexes the successive epochs. Thus, the synchrony computed by summing over the epochs leads to the definition of a phase-lock index ρ_si which characterizes, at each point (a, b) of the time-frequency plane, the presence of the response evoked, independently from the amplitude of this response:

 

Fig. 5. Computation of the synchrony maps, i.e. ρ_i,k(a,b), for s_i=HbR (or HbO) and s_k one of the physiological signals. The top row shows the HbO correlations, the bottom row HbR (group averaged wavelet synchrony over 10 subjects, 3 optodes each).

Download Full Size | PDF

By summing over the blocks, this quantity is defined as the average of the wavelet coefficients signs over epochs of the stimuli. It is close to 1 for similar complex signs at each epoch. Conversely, a value close to 0 reveals the absence of a coherent repetition of a common and detectable constituent of the signal.

Let us now consider the following shrinkage procedure to recover the hemodynamic response:

g(ρ) is a smooth function between 0 and 1 and weights each contribution of the wavelet coefficient accordingly to the synchrony index ρ_s(a, b) (the translation parameter can be interpreted as time, i.e. b=t). The interest in this reconstruction procedure is that the index is measuring the instantaneous synchrony with the protocol independently of the strength of the underlying signal. Thus small signals may be detected more easily. Finally, spurious synchrony may emerge from the cardiac dynamics in a frequency range outside the spectral support of the expected response. We can thus consider a partial summation in equation (14) ignoring wavelet coefficients at high frequencies (as in down to 0.5 Hz):

Although the proper definition of g is not particularly fundamental, the choice of this function must keep the coherency of the coefficients in the definition of the signal to be reconstructed. The function g may be defined in different ways. Here we introduce a smooth thresholding of wavelet coefficients,

With this function illustrated in Fig. 8 (A), highest phase-locked coefficients are kept in the reconstruction. This reconstruction is original in two ways. First, it is based on the local phase content of the signal estimated by using analytical wavelets. Second, it is used with the continuous wavelet transformof the signal. This last pointmay be disputable because of the redundancy of the continuous wavelet representation. However, since our procedure is based on the fact that the expected signal must possess this coherent phase information, this procedure is valid. The forthcoming evaluations, either with simulation or real data, will confirm such an assumption.

3.4. Wavelet analysis of random stimulus

The previous description emphasized the use of phase information in a wavelet based analysis of DOI. On the other hand, we cannot apply the previous methodology when this HRF is randomly generated since block averaging and phase-locking is not possible. Nevertheless, the present approach can be adapted to this context.

Let us model the evoked hemodynamic response H(t) from some canonical model and the stimuli paradigm. The wavelet transform of one optode pair signal s is then a combination of the hemodynamic response, the physiology (denoted by Θ) and noise. Performing the wavelet transform (denoted W), we write:

At this point, let us assume we can identify and select wavelet coefficients that represent the evoked response with a very small contribution from physiology and other background noises. The selected coefficients are denoted by (a*_i,b*_i), i=1,…, S in the scalogram. With this, an estimate of the amplitude of the response, β, can be found by a least-square regression over this subset of coefficients

The solution is the well known normalized scalar product between the two S-dimensional vectors W_s=[w_s(a*₁,b*₁), …, w_s(a*_S,b*_S)]^t and W_H=[w_H(a*₁,b*₁), …, w_H(a*_S)]^t:

In this expression, W_H and W_s are vectors of selected wavelet coefficients for the HRF and the DOI signal respectively. The next step is to find a method to perform this selection. Previous sections have argued that phase synchrony is a discriminating measure in the characterization of the wavelet coefficients with respect to the paradigm. Following this path, a proposal is to select coefficients with a high degree of synchrony (greater than 0.9 for instance) provided they also have power on the evoked response. The response from the canonical HRF identifies the coefficients satisfying this last condition. For a given optode pair, the selection is thus the intersection of the set of coefficients having high synchrony with the stimuli and the set of dominant wavelet modes of the stimuli. In doing this procedure, we maximize the chance of keeping coefficients with a high signal-to-noise ratio such that the estimator of the form (19) is justified.

One point deserves attention: in the generic situation the HRF signal overlaps the physiological signal. In other words, our estimator may give β*≠=0 even for rest signals. This bias can be statistically estimated with a set baseline signals acquired before or after the experience. Denoting by $\overline{{\beta }_0^{*}}$ the average estimation on the baseline for some random stimuli, the corrected estimation for the response to the true stimulation will be given by ${\beta }^{*}-\overline{{\beta }_0^{*}}$ . This procedure is cumbersome since it requires a set of baseline data.

Another technique using a single baseline is to restrict the estimation by choosing only a subset of the entire set of selected wavelet coefficients. This subset can be chosen at random many times allowing the computation of a statistic for the estimator. Finally, a statistical test (t-score) for the presence of the evoked response in the signal:

where σ* and σ*₀ are the variance of the estimators, for the signal and the baseline respectively. Here we consider a t-statistics since we compare two conditions: the rest state versus the stimuli. High values of this score indicate a higher probability for brain activation since we discriminate against the baseline condition.

This section has provided examples and definitions surrounding the analytical wavelets and diffuse optical data. In the next two sections we will show how these measures can be used to reflect neuronal activity in both simulations and real data.

4. Hemodynamic evoked optical response : simulations

Before moving to real data, it is useful to discuss and evaluate the methodology on simulated data. The main idea here is to take optical data taken at rest and add simulated block activations with some strength β_c, c standing for HbR or HbO channel:

where s _rest(t) is a typical background signal for either HbO or HbR channels. In order to identify the protocol in the data, we can compute the pattern of the convolved protocol in the time-frequency plane and see if it emerges naturally in the wavelet coefficients. The results presented previously for the finger tapping experiment indicate that, if the stimuli is strong enough, it will emerge in the time-frequency plane, though more preeminently for HbR than HbO due to physiology. Here we simulated the same protocol with small values β _HbO=5×10^-6 and β _HbR=-1×10^-6 as illustrated in Fig. 6. With such (and weak compared to the previous finger tapping experiment) simulated activation, the time-frequency representation of DOI channels do not display clear response to the paradigm. This may show the limit of amplitude based methods. Therefore, we propose to exploit synchrony to extract this response from the complex wavelet representation.

4.1. Block average and phase-locking in the time-frequency plane

This simulation reveals many of the features that will be used later for estimation: first when looking at Fig. 6, we observe that the protocol alone emerges at very low frequency in the time-frequency plane. When we add the signal to the rest HbO and HbR traces (Fig. 6) we observe that the overall map (left columns of Fig. 6) does not prominently show the emergence of the protocol itself in the absolute values of the wavelet coefficients. This is expected since the added signal is not that strong when compared to the physiology.

 

Fig. 6. Simulation of an optode pair signal on a single run. (A): block design stimulation (black) and the hemodynamic response (in red). The HbR signal (B) and HbO signal (C) with the respective simulated activation (in red) added. (D): Wavelet coefficients of the HRF (red line in (A)). (E) and (F): Norm of the wavelet coefficients of the HbR channel and HbO channel respectively. The HbO channel is dominated by physiology. In figure (E), there are high intensity norms at high frequencies which have support on very small areas which is why the color map is not as well distributed. Wavelet transforms are normalized.

Download Full Size | PDF

Since averaging the epoch time series increases the SNR of the response, the scalogram can be equivalently averaged over the epochs. This is what is shown in Fig. 7(A) and Fig. 7(C) where the averaged power for both channels (HbO and HbR) has been computed. These averages also amount to performing the average in time and subsequently doing the wavelet transform.

The following observation emerges when looking at the data, even in the averaged case, the absolute value of the average coefficients is carrying strong components from physiology which can be seen in both HbR and HbO. In fact, the high power values observed in the HbO panel (Fig. 7(A)) are the sum of both Mayer waves and the protocol. Thus even filtering out the fast physiology will lead to a biased estimation. This situation is typical of optical experiments, Mayer waves often confound the responses that we are trying to observe because they have support over similar frequencies and have strong amplitudes. On the other hand the response we are looking for should be phase-locked to the stimulation while Mayer waves are not necessarily locked to the protocol. Thus if we move to the phase-lock measure (Fig. 7(B,D)), we see a different picture: the frequency region where the protocol is present has strong phase-lock index and this map may provide a way to gate the signal (or filter) in a coherent way to recover the response. Note that here this is true of both HbR and HbO channels where we observe an increase in synchrony in the region where the paradigm has frequency support. High phase-lock values also appear in the averaged map at higher frequencies shown by the presence of hot spots around 1 Hz. These are physiological artifacts that exhibit spurious synchronies. One way of removing those amounts to introducing a filter in frequency by selecting the phase-locked wavelet coefficients down to some frequency threshold above which we know that the hemodynamic response is absent. This is the main observation of the numerical experiment done here. The use of a phase-lock indice is expected to help differentiate better the physiology from the signal seeked. In order to better quantify this statement, we perform a set of simulations in the next section.

 

Fig. 7. Normalized wavelet power coefficients averaged over the blocks for HbO (A) and HbR (C) corresponding to the simulation of Fig. 6 on a single run. Corresponding phaselock maps: HbO (B) and HbR (D). In this case, the averaged power exhibits mostly the activation in the HbO channel whereas HbR channel only shows the evoked response through the phase-lock map.

Download Full Size | PDF

4.2. Simulation of wavelet coefficients selection and HRF estimation

To validate the reconstruction procedure, we performed a series of simulations by adding a simulated response to optical signals taken with subjects at rest. We followed the procedure described in Fig. 6, but extended it to the 10 subjects and 3 optode pairs in the data set (from [7, 8]) providing a comprehensive dictionary of signals for simulations. The same protocol as in Fig. 6 was simulated with the values β _HbO=5×10^-6 and β _HbR=-1×10^-6. For each case (N=30), we estimated the HRF by selecting coefficients according to (15) with cutoff function (16). As defined, this function vanishes for coefficients with phase-lock less than some threshold λ.

In order to set in place the methodology, two last points must be addressed. First we need to evaluate the normalization constant k_ψ, second, the value of λ needs to be computed. The constant k_ψ can be numerically computed by using the theoretical response evoked by the paradigm, H(t). Assuming that the wavelet transform of such a function is exactly invertible (i.e., absence of polynomial drifts), the reconstruction from the time-frequency plane should give it back:

Consequently, we have $k_{\psi }\approx \frac{max{H}^{*}\left(t\right)}{maxH\left(t\right)}$ .

To evaluate the best λ, simulations were performed for a set of values λ ∈[0,1]. In each case, we evaluated the estimation of the response s*(t) using the procedure described above. Then a regression fit of the wavelet-based estimations to the known response H(t) was performed to recover an estimated value of β. This is justified by the fact that with our reconstruction based on synchrony, the remaining physiology is now small compared to the evoked response, if any, and we can write

where δ is a small error. This is what is done to get Fig. 8(B) where we show estimated β*’s and related errors with respect to the threshold λ. It should be noted that the case λ=0 corresponds to block averaging (in the wavelet representation) without any phase-locking criteria nor selection. The best threshold corresponds to λ=0.3. Support for the assertion that phase-locked coefficients are informative with respect to the paradigm is also provided by the fact that the estimator improves when the selection criteria is introduced. This is observed in the HbO channel for λ between 0 and 0.3: wavelet coefficient shrinkage eliminates remaining physiological noises.

 

Fig. 8. (A): Cut-off function for keeping wavelet coefficients in terms of the synchrony index. (B): Estimation (and variance) of β in the HbO (red) and HbR (blue) channels in terms of the threshold λ in equation (16). λ=0 corresponds to the average in the time-frequency plane.

Download Full Size | PDF

Having chosen a threshold, the synchrony results can now be compared with other methods. In Fig.9 we present comparative estimation results for the General Linear Model technique described in [6] as well as a more standard filtering ([0.016,0.4]Hz) and bloc averaging technique. In the latter, we further fitted the value of the response to recover a value of β to provide a basis for comparison. To estimate the quality of our estimator, we also computed a Z-score on the the mean value with the measured variance (we used $\frac{{\beta }^{*}-\beta }{\frac{{\sigma }^{*}}{\sqrt{N}}}$ where β*, σ* are the computed estimate mean and standard deviation on the data set).

We observe that for HbO the complex wavelet approach performs better than the averaging and the GLMhaving a smaller error bar as well as bias. This is well represented in the associated Z score. For HbR, the averaging technique provides a smaller variance in the estimation, but it systematically introduces a bias in the estimator which leads to a slightly higher Z score than the complex wavelet approach. The good performance of the GLM in that case is due to the large variance and an estimator that has less bias. This is supported by the fact that the GLM model does not perform well in a bloc design with a small number of events. Thus the estimator performs better to recover HbO signals. For HbR, it is on par with the averaging technique. This is not surprising given that the synchrony measure was built to remove physiology which is not as present in HbR data.

 

Fig. 9. Estimated values of the responses for HbO and HbR respectively for N=30. In each case, three methods are compared: the General linear model (GLM), the averaging technique (AVG) and the complex wavelet technique (CPLX). In parenthesis: Absolute value of the Z scores for the above estimations given that we measure the mean of the distribution.

Download Full Size | PDF

The reconstructed signal can also be computed as described above, we provide here an example on a single estimation. The result is shown in the dotted-lines displayed in Fig. 10.

 

Fig. 10. Estimation of the hemodynamic response by wavelet shrinkage for a single case. HbR channel (A) and HbO channel (B). Block average (blue lines), synchrony based estimator (red lines) and true response (dot lines). The estimated β* are found equal to -0.9×10^-6 and 3.1×10^-6 for HbR and HbO respectively.

Download Full Size | PDF

4.3. Random Stimulus

The evaluation of the technique with random stimulus according to (19) has also been simulated. DOI data taking the form of equation (21) with different physiological backgrounds (N=30), and various response intensities, β_HbR ∈ [-4,0]×10^-6 and β_HbO ∈ [0,2]×10^-6 have been used. For a given random stimuli, our evaluation compares the GLM approach with the complex wavelet techniques described in section 3.4. Mean and variance of the error both for HbR and HbO channels are summarized in Fig. 11.

In comparison with the GLM method, the performance of the complex wavelet approach in the random case follows the same conclusion as in the bloc paradigm case. The regression on the selected wavelet coefficients improves the accuracy of detection. Most probably, this is because the chosen coefficients are describing locally in time the presence of the evoked response. Nevertheless, we have assumed that physiology has no or small support on the selected coefficients. This is an approximation that may be improved by a more accurate description of the physiology. The continuous wavelet model presented here is not the appropriate framework to model such a component because of the large number of coefficients to be considered. In this case, a discrete wavelet approach is appropriate [30].

 

Fig. 11. Estimation errors for HbO (red) and HbR (blue) responses for random stimuli. Two methods are compared: the General Linear Model (GLM) and the Complex Wavelet technique (CPLX). In each case, we indicate the t-score in parenthesis. For the complex cased, the error mean is closer to zero showing that the bias has been diminished. The t-score values reflects this.

Download Full Size | PDF

5. Hemodynamic evoked optical response: real data

Up until now, we considered real data. In the previous section we argued that even with small signals the phase synchrony may contain information that is of interest to detect the response to paradigm. We propose here to provide two explicit examples of its usage.

5.1. Phase-Locking and HRF estimation: the response in a visual experiment

Seven healthy volunteers were recruited for the study, six males and one female, all right handed (mean age: 28, range: 23–38 years) and had no history of any major psychiatric disorder. This study was approved by the institutional review board of the Centre de Recherche de l’Institut Universitaire de Gériatrie de Montréal (CRIUGM). Written consent was obtained from all subjects prior to the study and subjects received financial compensation after the study. None of the subjects failed to show significant hemodynamic activity and none repeated the DOI data acquisition. All subjects were right-eye dominant. The visual task chosen was identical to that of [16] where negative BOLD signal was observed and localized by a fMRI study. The participants were seated in front of a computer screen (1074×768 resolution) at a distance of 65 cm. Sources and detectors were positioned over the visual cortex (Fig. 12) located with the MRI of the subject and on-line calibrated with a neuronavigation tool. This montage defines N=24 optodes as shown on Fig. 12. The subject had to look at the stimulation on the screen with the dominant eye (the non-dominant eye was blocked by an eye-piece). The stimulus consisted of an up and down moving target changing direction randomly, either in the right side of the screen or in the left side. The subjects were then asked to count the number of direction changes. Since the stimulation is block-designed, we computed the synchronization indices and considered the group average over the set of subjects and strength of the stimuli. Plots of Fig. 13 display the mean synchrony in HbO and HbR measured on all the optode pairs of the montage. We observe that some optode pairs, mostly in the inferior occipital region have higher synchrony, especially in the HbR channel, as expected. The two series of experiments (stimuli in the left side and then in the right side of the computer screen) exhibit the same amount of synchrony. This is due to the connectivity between the ipsi and contralateral visual region [24] that balances the activity between the two hemispheres. Also, the size of the visual stimuli does not reveal any consequence on the synchrony that occurs independently from the intensity of the response. Optode pairs 6 (left hemisphere) and 18 (right hemisphere) localize the brain activity in terms of synchrony. We observe that the HbR channel is more localized than in the HbO channel. The upper occipital synchronies do not show so much contrast between optode pairs although we observe some asymmetry in HbO in favor of the right hemisphere. In order to see if our reconstruction is coherent, we considered a Principal Component Analysis (PCA) of the reconstructed signals and computed the weight of the first component on each optode pair. As a separate validation, we confirmed that the synchrony indices were following a similar trend as the PCA analysis.

 

Fig. 12. Optodes configuration (left) registered over the cortex of the subject (middle) and visual stimulation (right).

Download Full Size | PDF

5.2. Phase-Locking and contrast: a motor task experiment

The second experiment is a finger tapping experiment involving both hands. Here one of the issues found in previous work was the localization of the tapping zone when imaging both sides of the cortex. Typically, finger tapping increases blood flow over the whole cortex and it has been shown before [15] that doing a PCA of the data enables to isolate and better localize the activation. Here we present an approach based on synchrony that will lead to similar results. The protocol consisted of interlaced right-hand, and left-hand finger tapping (5 taps at 2 Hz) presented randomly at intervals between 15 seconds to 20 seconds. We will refer to left and right stimuli as L5 and R5 respectively. The optodes were positioned on both right and left motor cortex using identical geometry on each side. The precise distinction between left and right signals can be found by averaging and doing a PCA but here we propose to use synchrony by using the usual synchrony index:

 

Fig. 13. Block averaged synchronization index per pair number for HbR (left column) and HbO (right column). Results from two experiments are shown in different colors.

Download Full Size | PDF

 

Fig. 14. Optodes configuration for the left side of the head.

Download Full Size | PDF

From the above indices of both side, we can build a contrast function distinguishing right and left activation (independent of the signal strength on each side of the head) as follows:

These indices are shown in Fig. 15. We observe that without any filtering or PCA analysis, the synchrony index provide a clear localization and distinction of activation sites. Finally, synchrony-based reconstruction of the hemodynamic responses on selected optode pairs (the ones that show a contrast between the two conditions) are shown in Fig. 16. We observe consistency between the HbR contrast and the amplitude of the HbR responses that consistently increase in the ipsilateral side of the stimulation.

 

Fig. 15. Top:Map of the synchrony indices for all optode pairs for the HbR channel (Similar map for the HbO channels does not exhibit such a significative contrast between optode pairs). We indicate in bold face the most contrasted optode pairs, they are those located over the motor area of the brain’s subject.

Download Full Size | PDF

5.3. Application to random stimulus

Our simulations showed that it is possible to use complex wavelets with event-related design paradigms. In this section we provide a simple experimental example. A finger tapping experiment was performed with an event related design. Location of the optode pairs over the left motor cortex were similar to the previous setup (see Fig. 14) but optode-pair labels are different, as seen on the top of Fig. 17 defining the experimental configuration. The stimulus is composed with 80 events of 2 sec duration each, randomly generated with a 4 to 20 seconds inter-event delay. The expected evoked response is illustrated on Fig. 17 (a). The scalogram of the HbR channel (optode pair 6) is represented in Fig. 17 (b). This is where the most activity has been detected using the t-score as defined in (20). The strength of the response on each optode pairs and the corresponding t-score are displayed on (c,d) and (e,f) respectively. Activity is located over the motor region with a spatial profile that modulates with the location. Moreover the statistical significance of the recovered values (t-scores) is high at those same locations for both HbR and HbO. We hope to elaborate this approach and evaluate its performance on a more extensive experimental set in forthcoming work.

 

Fig. 16. Synchrony-based reconstructions of the HRF for the HbR channels over the motor area. The solid lines exhibit the dominant response in case of the Left (red) and Right (blue) stimulus.

Download Full Size | PDF

6. Discussion and conclusion

The present study emphasizes the use of the phase of the complex wavelet coefficients of DOI data, in the detection and estimation of some evoked response to stimuli. Working with periodic block paradigms, the phase allows to point out wavelet coefficients that reproduce the paradigm despite the physiological background. In some respect, our approach mainly relies on the phase resetting synchronization with the periodic stimulus. The case of random event-related stimulus is more challenging, but may also take advantage of the phase information of the wavelet representation of the signal. Very preliminary results indicate that this is possible, but an implementation with discrete wavelets will be necessary to make the drift estimation problem tractable.

On their own, the methods developed here do not include physiological modeling, rather they allow better physiology separation. We have shown the ability of these new bases to separate physiology enables better estimations when compared to similar static techniques. We believe two factors are behind those results, the use of a synchrony measure providing adaptive filtering and the potential presence of 1/f noise. The separation may provide new insights in more accurate models of response estimation whereby physiology is used as an input to perform dynamical modeling [7].

Separately, we developed a new instrument to study neuronal activation, synchrony. Here the ability of analytical wavelets to define an instantaneous phase in a concrete manner opens the door for a new measure, the phase-lock of the signal with itself. The advantage of this measure is that it can be used to weight the reconstruction of the hemodynamic response function based on our knowledge of a given protocol. It thus provides a protocol-based, filtering technique.

 

Fig. 17. Event related finger tapping experiment and the optode pairs setup (TOP). (a) The expected response, (b) the scalogram of the HbR channel for pair 6, (c,d) the estimated values for the (top,bottom) optode-pair rows, HbO in red, HbR in blue. (e,f) the correponding t-statistics for the same estimates.

Download Full Size | PDF

Finally, the use of wavelets in general can be of significant benefits when noise has 1/f behavior. It is known that the wavelet basis whitens such correlated noise [26–30] and if present, we could use such basis to perform a general linear model as we discussed in a previous section, to recover event-related activations. In diffuse optical imaging, noise has many sources, some of which are not correlated and the presence of 1/f noise and its influence have not been quantified. Quantifying its presence as well as bias in associated estimators will be the subject of further work.