Spatio-temporal analysis of the cerebral spontaneous oscillation

Cerebral vasomotion was studied on the human brain in vivo by use of multi-optode frequency domain near-infrared spectroscopy (NIRS). Vasomotion is a spontaneous oscillation with a frequency of 0.1 Hz in the arterial flow. We investigated (1) the fluctuations of cerebral hemodynamics on the dynamical characteristics of cerebral vasomotion and (2) the dynamical coupling between vasomotion in the skin and brain. We found that (1) vasomotion is temporal coherent at least for about 3 min; (2) vasomotion observed from NIRS is low-dimensional chaotic with its fractal dimension of about 4.5; (3) vasomotion is spatially coherent with coherence length of about 1 - 2 cm but cerebral vasomotion is dynamically independent from vasomotion in the skin.


INTRODUCTION
The determination of the optical properties of the brain is important in many fields of medicine and neuroscience, both for monitoring of physiologically important processes and localizing any brain activities. In particular, the measurement of hemoglobin concentrations in the brain can provide important information about the cerebral oxygenation or related to neurovascular coupling. In near-infrared region, the light is highly scattering and low absorbing and the main chromospheres in biological tissue are oxy-and deoxyhemoglobin molecules. The spectral studies of the human brain using the near-infrared spectroscopy (NIRS) have been developed recently [1][2][3] and the baseline values ofthe oxy-and deoxyhemoglobin concentration of the adult human brain have been successfully measured using multi-optode frequency-domain NIRS separated from the scalp noninvasively [4].
The frequency-domain NIRS monitors the dynamical changes in the hemoglobin parameters continuously. Generally, the measured signals are fluctuating a lot with multiple sources and their spectral analysis show rich structure in the frequency domain. The frequency components in the spectral analysis are variable depending on the individual, but mainly we observe primary frequencies around 1, 0.3, 0.1, and 0.04 Hz corresponding to pulse, breathing, vasomotion, and the blood pressure wave, respectively. The rhythms, except 0.1Hz, are systemically driven or controlled by the central nervous system [5,6], whereas the origin of vasomotion is still unknown but it exists in the entire body.
Vasomotion is a rhythmic alternation of vasoconstriction and vasodilation of vessels, which is observed in the diameter measurement of the large and small arteries, and arterioles. It has been showed that 0.1Hz oscillation in NIRS is due to vasomotion by many researchers [7][8][9].
In the NIRS signals, the amplitude ofvasomotion is about 1-10 fold amplitude ofthe pulse and the frequency ofit is variable but narrowly banded within the frequency band of 0. 1±0.04Hz. In the diameter measurements, the patterns of vasomotion are variable depending on the types of arteries and measurement point [10] and the dynamics of vasomotion is determined to be highly nonlinear [11]. The spatial averaging of the NIRS signal is much larger than the spatial unit of the artery network, and therefore it is more reasonable to assume that the slow oscillation is the superposition of numerous vasomotion waves rather than an individual vasomotion at the level of single arteries. Therefore, the 0.1Hz oscillation observed in the NIRS signals is a collective vasomotion wave.
Vasomotion in the brain is particularly interesting with its dynamic coupling with brain activation. Recently, Cimponerin and Kaplan [12 observed from the optical imaging of the cat visual cortex that vasomotion is stimulusdependent. According to their observation, the spatially asynchronous vasomotion became spatially coherent at the start of stimulus as well as phase-locked to the stimulus. This is direct evidence that vasomotion is sensitive to neuronal activation and the phase-locking phenomena suggest that the relation between the neurovascular coupling and vasomotion is nonlinear. In addition, Diehl et al. [13] suggest the phase shift of cerebral vasomotion to the blood pressure wave as a measure of cerebral flow autoregulation, Biswal et al. [14,15] and Lowe et al. [14,15] even try to find a criterion for a functional connectivity map based on the correlation of spontaneous vasomotion. Most of all, the coupling behavior of cerebral vasomotion to the neuronal activation arises as an interesting subject with the increased interest in the area of the functional brain mapping. Usually, the cerebral vasomotion and neurovascular coupling are observed to interact with each other in a complicated manner so that the cerebral vasomotion is introduced as an intrinsic artifact in analyzing signals from neurons and it is challenging to filter.
In this manuscript, we focus on defining the dynamic characteristics of vasomotion in the cerebral hemodynamics assessed by frequency domain NIRS.

Multi-optode frequency domain near-infrared spectroscopy
The propagation of light in tissue is well described by the diffusion equation for semiinfinite [16] or two-layer media where ln[p2M(p)] and ô,,q$ are the slopes ofmodulation and phase shift with respect to p. Using Eqs. (3) and (4), Pa and Ps can be determined directly from the values of 0and M measured at p in the frequency domain. One benefit of these expressions is their independence from the instrumental parameters such as the laser light power, detector sensitivity, and other instrumental factors. However, since the above equations were obtained under the assumption that the amplitude and initial phase are the same for each light source, an instrumental calibration must be performed to correct for any inequalities in the amplitude and phase by calibration with a tissue phantom with known optical properties.
According to Beer-Lambert law, the absorption coefficient can be written as a linear sum of the extinction coefficients multiplied by the concentration of absorber and in the near infrared region, oxyhemoglobin (c0y) and deoxyhemoglobin (cdeoyx) are main absorber. Therefore the absorption coefficient is represented by 2 + CoxyCdexy (5) where g1 is the extinction coefficient of the ith chromophor. Since the p, and are dependent on wavelength of light, we can solve eq. (5) by use oftwo wavelengths of light.

Experimental Methods
The instrumental configuration used in determining the hemoglobin parameters is schematically shown in Figure 1 .We used a multi-channel frequency domain near infrared spectrometer (Imagent, 155, Champaign, IL). In brief, the instrument employs laser diodes (LD) modulated at the frequency, 110MHz and photomultiplier tubes (PMT) whose gain is modulated at a close frequency of 1 10.005MHz to heterodyne the high frequency down to the frequency of 5 kHz. In this study, 32 LDs (16 at 758nm and 16 at 830nm) and our PMTs were used. The LDs were multiplexed so that two LDs with the same wavelength and at the same location were on at a time. The light from LDs was guided by the optical fibers with a diameter of 400 im to the tissue surface and the photons re-emitted from the tissue were collected simultaneously by the fiber bundles with a diameter of 5.6 mm placed several centimeters apart from the source fibers. The collected light was carried to the PMTs and then the signals from the PMTs were digitally processed to yield I, M, and.
We investigated the time series of oxyhemoglobin concentration (c01) as the vasomotion parameter of interest. Data were obtained from an adult forehead with multidistance frequency domain NIRS probe (source-detector distance 1-8cm) at data acquisition frequency of 2.5Hz. Figure 1 Block diagram of frequency-domain near-infrared spectroscopy (Oxy-Imager ISS, Champaign, IL USA) and geometric arrangement of the long-ranged multi-distance probe containing the source and detector fibers on the forehead.

Data analysis 2.3.1. Pattern recognition
A repeatability of a similar pattern in a time series can be characterized by pattern recognition. We use the least squares method to fit a part of interest to the entire time series. To illustrate the linear squares fitting process, suppose one has a test data set {yj} with a size of N that can be fitted to a reference data set {x1} with the same size. We maymodel the reference data by a first degree polynomial such as y = px+q, whereas p and q are fitting coefficients.
We solve the equation for the unknown coefficients p and q with a least squares fitting process by minimizing the normalized summed squares of the residuals j(p,q) ,i.e.

human head
Proc. of SPIE Vol. 5330 31 (yi -(px +q))2 z2(p,q)= '= 2 ' (6) Na where N is the length of the data and c is the sample variance of the data yi. The minimized summed squared residual x2(p,q) is solved by d2 (p, q)/ dp 0 and d2 (p, q)/ dq = 0 . Using this method, one can recognize a particular pattern in a time series by scanning the whole data with the pattern. If the data is the same as the pattern, ,will be zero. On the contrary, if one fails to interpolate the data with the pattern, the minimized j will be a constant function with the mean value of data, therefore, j will be zero by definition.

Cross-correlation study
One effective way to describe the joint properties of two time series is the cross correlation function. It describes the general dependence of one signal to the other. For two time sequences, x(t) at time tand y(t) at time t+'r for an observation period T, a cross correlation function is defined by CCF(r) tim ---Jx(t)y(t + r)dt. (7) T-o 2T T The cross correlation function CCF(-r) indicates the existence of correlation between x(t) and y(t) for a specific time displacement, -r. Therefore, the maximum of CCF(-r) will yield the time delay ofy(t) with respect to x(t).
A cross correlation function within the same signal is called the autocorrelation function. The autocorrelation function for random data provides a general temporal structure of a time series by evaluating the dependence of the value of the time series at one time on the value at another time. Intuitively, it is a measure of how similar a series is with a shifted version of itself. This autocorrelation function is an even function by definition and has a global maximum at the origin, i.e. r = 0. The major application ofan autocorrelation function is to monitor the influence ofa value at any time to values at a future time. For example, a noisy signal without any deterministic nature has a quickly decaying ACF(t) as r increases. Whereas, if the signal is deterministic, ACF('r) will exist over all r and have local peaks at every period of oscillation modes. Hence, ACF(t) is very useful for detecting a characteristic frequency buried in an extraneous noise. If CCF(t) oftwo independent time series is similar to ACF('r) ofone, the signals are coherent to each other.

Fractal dimension
The fractal dimension is a magnification factor of self-similar objects in a fractal geometry that is defined by D = log(number of self-similar unit)/log(magnification factor). Basically, it is a measure of how complicated a self-similar pattern is. This concept may also be applied for determining the dynamic characteristic of a time series. In case of a time series, Df can be the degree of freedom of the dynamical system that the series belongs to and it is a measure of strange attractors ifthe system is deterministic chaos [18]. Ifthe time series is completely random, Djwill be infinite whereas if the time series is a deterministic signals, Djwill be a finite number.
A popular method to calculate a fractal dimension in time series analysis is the correlation dimension D: where C(r) is the number of points of the set located within a smaller distance of the values than a given distance r (correlation integral). The correlation integral C(r) is defined as C(r) lim H(r --x ), (9) i,j=1 i j where H is the Heaviside step function (H(s) = 1 for s>O and H(s) = 0 for s<O). For example, a distance r can be arbitrarily chosen and C(r) can be evaluated by counting the number of points of the neighborhood in the experimental data that falls into the distance r. And the slope in a plot oflog C(r) versus log (r) as r goes to zero will beequal to D.
However, in practice, the experimental data are not continuous but discrete and D will yield different values for different time scale. Therefore, a reliable way to evaluate D is to calculate D for different time scales and find a value that remains constant regardless ofthe change ofthe time scale. The time scale can be easily modified by multiplying an integer m in the time index x1 -x,,,, and usually, m is called an embedding dimension. Then the correlation dimension for this embedding dimension m, C(r,m) is calculated by C(r, m) tim --H(r -Xmj -Xmj ). (10) Ifthe time series is noisy, C(r, m) will increase with increasing m without ever reaching a plateau. But if the time series is deterministic, C(r, m) will reach a plateau for a time scale larger than a certain value of m because of its characteristic of self-similarity. And the plateau value of C(r,m) is the fractal dimension of the time series.

Temporal structure
We first estimated the temporal coherence based on the pattern recognition method. The temporal coherence measures the correlation of two distant events in a time series. We used the nonfiltered data after detrending with a 3rd order polynomial function (Figure 2(a)). Figure 2(b) shows the result from the pattern recognition method. For the pattern recognition analysis, we randomly picked a time window with a size of lOsec and scanned the entire time series and evaluated the summed square of the residuals X2(t). The local minima indicate the appearance of a similar pattern. A detailed analysis of X2(t) shows that local minima tend to be separated by a period of 15 sec and that this behavior We next use the autocorrelation analysis to determine the long range dependence ofthe vasomotion wave. Figure 3 (a) shows ACF of the same time series in Fig. 2(a). Figure 3 (b) shows the same function in a smaller time window. The modulations of ACF with an approximate period of 1 Osec and other oscillations at longer periods are also observed. The main feature of Fig. 3 is that the ACF of the oxyhemoglobin signal decreases to zero slowly in an oscillatory behavior conserving the long range dependence lasting for a significant amount oftime (> 3mm).

Fractal dimension
The fractal dimension of the vasomotion was calculated using the Grassberger-Procaccia method [1 8]. We first chose a time delay of 2.5sec that is approximately ¼ ofthe vasomotion period. The correlation integral C(c,m) was obtained by Eq. (10) by calculating the probability that the pair of neighboring points has a distance smaller than e for an embedding dimension m. The correlation dimension was obtained from the slope by a linear regression in the linear region with distinct slopes of log C( m). The correlation dimension D was plotted with respect to the embedding dimension m in Fig. 4. In order to test the dependency ofthe fractal dimension on the time delay r, we plotted D for different values ofT in Fig. 4. Each plot of D versus m reaches a plateau and its value is similar. The value at the plateau is taken as the estimate ofthe fractal dimension Djand a fractal dimension ofabout 4.5 is suggested from this plots. We analyzed a set of data spatially distributed in the vertical section into the brain to extract information on the dynamical correlation between vasomotion at the skin and the brain surface. The data was obtained from the forehead of an adult with the multidistance probe ranging from 1-8cm of its source-detector distance. Visualizing the correlation between vasomotion is very difficult, when we do not filter the data, because the fluctuation amplitudes of the other physiological signal such as pulse and blood pressure wave are comparable to the vasomotion. Therefore, we filtered the vasomotion with 5th order of Butterworth bandpass filter with passband from 0.06Hz to 0.20Hz and stopband from 0.02Hz to 0.28Hz. Since the phase difference between two parameters is one of our interests, we used zero-phase filtering by processing the input data in both the forward and reverse direction [191. ). : Figure 5 The region of sensitivity for long-ranged multidistance NIRS signal. The horizontal axis is the surface of the forehead and the vercal axis is the depth from the tissue surface. The NIRS signal at each position is obtained from the slopes of . From the given probe, the approximate region of maximal signal sensitivity is shown in Fig. 5 and the penetration I depth z is determined by z = r"2 Pa/s ) " 2 [201. We denote the filtered oxyhemoglobin contents at the th position as c,,(t). n1 indicates the closest point to the skin and as n increases, the measurement point goes deeper. In order to investigate the spatial correlation between vasomotion signals at different depths with the aim of describing the dynamical feature of the propagation of vasomotion waves, we use a cross-correlation technique. The cross-correlation Xn,m between two simultaneously registered time series of c,,(t) and Cm(t) 15 determined by Eq. (7) For all the possible pairs, we draw the cross-correlation map of a test signal c,,(t) with the reference signal .By simple visual inspection of the cross-correlation plots, we joined all the regions that have similar cross-correlation to the autocorrelation at a given location, which indicate that the patterns of slow oscillations are similar at that location. We obtained the map shown in figure 6. A clear pattern of local spatial coherence appears localized into topologically distinct and non-intersecting  x(cm) Figure 6 Topological map of the regions that have similar cross correlation to the autocorrelation at a given location.

DISCUSSION AND CONCLUSION
Both the pattern recognition and autocorrelation analysis show that similar patterns are repeating with the frequency of vasomotion and the correlation dissipates slowly with time. This analysis clearly shows that the long range dependence of the vasomotion is not limited to a few periods of vasomotion but exists for dozens of periods. In other words, the vasomotion wave is a temporally coherent wave even though there is no explicit driving force ofthe fluctuation.
The finite fractal dimension suggests that the dynamics ofvasomotion belong to a low-dimensional dynamical system. If the vasomotion is stochastic, D will increase continuously without reaching a plateau with m. It is interesting that the fractal dimension we obtained is similar to those obtained from diameter measurements. Rosen et al. [21] calculated the fractal dimension ofthe diameter and red blood cell velocity in a single vessel. Fractal dimensions of3-4 were obtained for diameter oscillations and that of 6-7 for velocity oscillations in the same vessel. Slaaf et al. [22] also obtained the fractal dimension of about 4 in the diameter measurement. The signal from NIRS is a volumetric parameter measuring the amount of fresh blood in the arteriole network in a tissue. Therefore, it is reasonable that our result is similar to those for the diameter measurement rather than the flow measurement. However, the vasomotion results from NIRS are statistically averaged over a large scale of arteriole networks, merging the individual vasomotion at a number of arterioles. If any contribution from an arbitrary dynamical factor is added or removed during the signal integration, the g4 . v...