Diffuse Optical Tomography Using Bayesian Filtering in the Human Brain

: The present work describes noninvasive diffuse optical tomography (DOT), a technology for measuring hemodynamic changes in the brain. These changes provide relevant information that helps us to understand the basis of neurophysiology in the human brain. Advantages, such as portability, direct measurements of hemoglobin state, temporal resolution, and the lack of need to restrict movements, as is necessary in magnetic resonance imaging (MRI) devices, means that DOT technology can be used both in research and clinically. Here, we describe the use of Bayesian methods to filter raw DOT data as an alternative to the linear filters widely used in signal processing. Common problems, such as filter selection or a false interpretation of the results, which is sometimes caused by the interference of background physiological noise with neural activity, can be avoided with this new method.


Introduction
Functional brain imaging has provided substantial information regarding how dynamic neural processes are distributed in space and time. Some imaging modalities to study brain function use functional magnetic resonance imaging (fMRI) and, although these measurements can cover the entire brain, they require costly infrastructure. Technical limitations in fMRI devices, such as a fixed scanner, contraindication with metal implants, scanner noise, and stress associated with fear, are avoided or reduced using optical imaging devices. In addition, optical imaging techniques, such as functional near infrared spectroscopy (fNIRS), can measure changes in oxyhemoglobin (HbO) and deoxyhemoglobin (HbR) at a much higher sampling rate than fMRI devices. Therefore, contrary to positron emission tomography (PET) [1] or X-ray computed tomography [2], optical imaging does not require a contrast agent, which, hence, avoids issues, such as limited doses for infants or anaphylactic reactions. Moreover, optical imaging techniques are low-cost when compared to other neuroimaging techniques, such as fMRI or magnetoencephalography (MEG). These circumstances have potentiated the optical imaging techniques development in recent years in research, diagnosis, and prognostic studies.
The first biomedical application of NIR light in the human tissue was published in 1977 by Jöbsis [3], showing blood oxygenation measurements while using NIR light of approx. 700-1300 nm wavelength. This is known as the optical window, because NIR light at this wavelength can pass through the tissue, to a depth of around 3-4 cm [4]. Because of the presence of the natural chromophores in the tissue (HbO & HbR) that absorb NIR light, both are used as cerebral activation markers [5].
NIR light is applied on the subject head through LEDs or optical fibers, which are combined as a source-detector pair, recording the changes in the optical density generated by changes in the optical properties in each tissue layer of the human head. Diffuse optical tomography (DOT) devices measure hemodynamic changes that are correlated with neural activity, caused by dynamics in blood volume, blood flow and blood oxygenation. Certain physiological signals, such as heart rate or ventilation rate, are involved in systemic blood oxygen and cerebral hemodynamics and, hence, influence the scalp layer [6]. During functional DOT neuroimaging experiments, these physiological signals create background noise. This noise is stronger than those in anatomical regions and mixed with absorbed signals from the cortex, thereby generating short-term variability that involves spatial and temporal changes throughout the brain and scalp [7]. Some researchers have devised methods to filter out the background physiological noise from the signals generated by neural activation. The most common methods are linear filtering, such as the application of band pass filter [8], low pass filter [9], high pass filter [10], adaptive filter [11], or principal component analysis [12], depending on the task or study area. A common problem during the processing of DOT data, which generally occurs during signal filtering, is that the judgment of the researcher might lead to an error in the cutoff frequencies selected. Because of this, a more robust and independent procedure is necessary. Using the Bayesian algorithm application as a filter method on the raw DOT data will allow for an independent procedure, without the influence of the researcher, paradigm, or cerebral area for study, which is more accurate than the linear methods used so far. Here, we describe a Bayesian method to model physiological data and separate it from the neural response, which avoids the common filtering controversies during the DOT data processing.

Optical Data Acquisition
A DYNOT 232 instrument (NIRx Medizintechnik GmbH, Berlin, Germany) was used to acquire the DOT data. The system performs continuous-wave measurements while using two frequency encoded laser sources at 760 nm and 830 nm with a sampling rate of 1.81 Hz in a time multiplexed scanning fashion. The equipment provides a high dynamic measurement range needed for diffuse tomography multi-distance measurements. NIR light travelled to and from the DOT device through optic fibers (optodes). Figure 1 shows an array of 64 fiber optic probes that were separated by 1 cm, which were used to measure hemodynamic changes in the frontal cortex. The fiber optic probes acted as detectors and 32 of them acted as source (colocalization), thereby providing 2048 optical channels. The optical properties for each tissue layer in the human head were measured according to the distance between the source-detector pairs that were placed on the head surface. At a distance of around 3-4 cm between a source-detector pair, the photons contain information regarding the optical properties from both intracerebral and extracerebral layers, whereas, at a distance of 1-2 cm between the source-detector pair, the photons contain mostly information from extracerebral layers [13]. The photon trajectory follows a banana path (arrows); therefore, detections at a 10-20 mm distance from the source (S) provide information from extracerebral layers, while detections 30-40 mm from the source also provide information from intracerebral layers. Source-detector pairs are separated by 1 cm in frontal and sagittal views. A total of 64 optical fiber forms the matrix, providing 2048 optical channels thanks to colocalization. Note the rigid structure holding the tip of the optical fibers in place in order to ensure the optical fiber-scalp coupling (b).

Physiological Signal Modeling
The Kalman filter [14] can be used to extract the Bayesian distribution of the state in the Gaussian state space model, and it is defined as: where ( ⃗ ) ∈ ℝ is the state at time , ( =0,1,2…), ⃗( ) ∈ ℝ is the measurement at time , is the state transition matrix, is the observation model, ⃗( )~ (0, ) is the Gaussian process noise, ⃗( )~ (0, ) is the gaussian measurement noise, is the covariance of process noise, and is the covariance of measurement noise. However, Kalman filters are based on linear dynamical systems discretized in the time domain. Therefore, at each discrete time increment, a linear operator is applied to the state to generate the new state, including some noise. Physiological signals, such as the respiration rate, contain periodicity and are time-varying; therefore, the Equation (1) must be expressed in a continuous dynamic model following the Wiener process [15]: where models the dynamics of the state and the state-dependent noise; ⃗ ( ) is a white noise process with a given spectral density matrix . If Δ = − , the weak solution for continuoustime stochastic differential equation [16] can be expressed as: Here, if we define = ( Δ ), we arrive at an equivalent model to the one expressed in Equation (1).
When assuming that there is a reference sensor to measure the respiratory rate ( ), the signal can be modeled as a zero-mean periodic signal with period frequency , which can be approximated to expanded Fourier series: As our interest is periodic signal modeling ( ), where the frequency is a function of time, then the Equation (5) can be written as: where refers to the physiological signals in terms of the frequency in time . Neither coefficients a nor b are considered time-varying, but simply amplitude constants.
Equation (5) can be described as a differential equation, where the constant frequency is replaced by a time-varying one in order to avoid the problem of both a and b being constant coefficients, and to maintain the continuity property in the frequency domain. Moreover, physiological signals have changes in the phases and amplitudes of the harmonics, which can be modeled by adding white noise ( ), giving us: This Equation (6) can then be written as a stochastic differential equation [16]: Because we are modeling the frequency trajectory as a constant signal with a noise term, it is not necessary to add a frequency-dependent term to explain an exact state space. Hence, because the modeled physiological signal in the stochastic space model is suitable for Kalman filters when the frequency is known in a time-varying linear space, it can be written as (equivalent to Equation (1)): where , = 1 for = 1, … , and other all values are zeros; = (1 0 1 0 . .1 0).

Forward Model and Image Reconstruction
In functional neuroimaging, the goal is to model the changes in the light attenuation measured from the boundary volume ∂ , generated by changes in the optical properties μ within a volume Ω: where depicts measurements onto the boundary volume, based on optical properties (absorption and scattering coefficients) in a position within . J is the Jacobian operator, which relates changes in the measured light intensity on the boundary volume ∂Ω with changes in internal optical properties , and it can be written as = ( )/( µ). This matrix is constructed from a model that is known as the forward model, which is derived from the Radiative Transport Equation (RTE), the mathematic model that is used in DOT and describes light propagation. From simplifying the RTE, we arrive at the diffusion approximation (DA): where ( , ) is the photon density in time and position , ( , ) is an isotropic source, is the absorption coefficient, c is the speed of the light in the medium, and D is the diffusion coefficient.
The DA equation was described in a finite matrix that can be solved using the finite element model (FEM) [17]. FEM discretizes the continuous domain into a finite number of elements based on basis functions. An important aspect during the light propagation modeling is the anatomy model, which adds more accuracy to the spatial distribution of the optical properties [18]. Therefore, the light propagation modeling in the anatomy model impacts the reconstructed image quality. Anatomy images from MRI were used for light propagation modeling at the layer level. These anatomy images had been segmented into five layers aforehand: scalp, skull, cerebrospinal fluid, and gray and white matter. A mesh was then created from each layer, improving the anatomical details to a submillimeter scale.
During image reconstruction from a turbid medium, the photon density (or light intensity) Φ and optical property µ relationship is nonlinear. In DOT measurements, changes in the light intensity ∆Φ that are generated by changes in the inner optical property ∆ can be assumed as relatively differential changes; therefore, the problem can be considered to be linear [19]. A perturbation method [20] allows for us to construct the changes in the optical properties according to the measured light intensities on the surface. During the image reconstruction, the J matrix is directly inverted, as it is ill-conditioned and ill-posed. To invert , the singular value decomposition (SVD) algorithm was used. The diagonal matrix from SVD decomposition contains singular values providing information on the propagation errors from ∆Φ to ∆µ.
was used as regularization parameter in order to reduce the inverse matrix's dimensionality while using a minimum description length (MDL) index [21].

DOT and fMRI Statistics
A new approach for treating DOT volumes as if they were fMRI volumes was used, using the canonical SPM 8 software to obtain cerebral activation maps, as is performed with fMRI images [21]. fMRI and DOT volumes were both filtered to improve the signal to noise ratio, while using a high pass filter based on a discrete cosine transform, with a 64s cut-off period. The regressors were convolved with the canonical HRF functions, according to the design matrix. After estimation, cerebral activation maps were generated by applying a fixed effect model analysis. The cerebral activation maps that were provided by DOT data processing were compared to the cerebral activation maps generated by the fMRI processing method based on the standard General Linear Model (GLM) [22].

DOT Data Filtering
The optical signals propagating through the brain contain several spontaneous fluctuations originating from cardiac pulsation, respiration, and change of blood pressure, which contaminate the signals measured by DOT and induce spatial and temporal changes that may lead to false interpretation of brain activations. The Bayesian method explained in Section 2.3 was applied to raw DOT data to filter out the physiological noise from the hemodynamic response. Figure 2 shows the time series for DOT data both without filtering and after filtering with the Bayesian algorithm at a wavelength of 760 nm. Moreover, it shows the power spectrum density (PSD) for the respiratory fluctuations extracted from the Bayesian method versus the PSD for unfiltered data. Physiological noise and its harmonics are both mixed with the neural signal and they cannot be distinguished from the respiratory component in the pre-filtered data. However, the signal extracted by the Bayesian filter shows the respiratory component around 0.15 Hz. Since cardiac fluctuations vary in a range of 1-2 Hz and the sampling rate used for this example was 1.81 Hz, these recordings are not represented here. The raw DOT data were filtered while using a band pass filter with a cutoff range between 0.01~0.3 Hz, the most commonly used frequencies in literature [23], in order to test the reliability of Bayesian-filtering in the optical channels. Figure 3 shows both wavelengths (760 nm and 830 nm) unfiltered, Bayesian filtered, and band pass filtered. Reconstructed DOT volumes from linearly filtered data do not show statistical significance during the GLM processing (p < 0.05), therefore only the differences at the wavelength level are shown.

Imaging Results
DOT systems are characterized by high temporal resolution (~500 ms) as compared to fMRI systems (~2000 ms), meaning that slow physiological oscillations are rapidly sampled, causing nearby time points in the DOT signal to be highly correlated. Serial correlations violate common statistical assumptions of the independence of repeated measurements over time and contribute to inaccurate estimations of type I errors, making a statistical analysis that is based on a GLM impossible. Applying the Bayesian filtering on raw DOT data reduces serial correlations between consecutive points over time, hence providing data independence, which will allow for the use of a GLM that is commonly used in the neuroimaging field for mapping cerebral activations. Figure 4 shows T-maps of brain activation for both the blood-oxygen level-dependent (BOLD) signal and the HbO signal (p < 0.001, corrected false discovery rate, FDR). The figure depicts an increased statistical significance in Bayesian filtered DOT images. Brain activation maps generated from unfiltered or band pass filtered (0.01~0.3 Hz) wavelengths have not been represented in this work, because they did not show significant changes (T: 1.282, no significant voxels, p ≤ 0.05, uncorrected; and, T: 1.645, no significant voxels, p ≤ 0.05, uncorrected, respectively). Figure 3. Representation of the average of the optical channels (2048) over time. The blue plot shows raw DOT data, the red plot shows DOT data filtered with a band pass filter, and the black plot shows DOT data after applying the Bayesian filter algorithm for 760 nm (a) and 830 nm (b). Abscissas axes correspond to normalized arbitrary units and the ordinate axis represents the experimental time in frames.
A conjunction analysis was performed for the example shown in Figure 4. The spatial conjunction that is shown in Table 1 depicts common anatomical regions between the hemoglobin molecule HbO and BOLD signal measured by DOT and fMRI equipment, respectively. As a result, 296 voxels with a size of 2 × 2 × 2 are activated by the BOLD signal, where 174 voxels are common with the HbO signal (~59%). . T-maps of brain activation measured by DOT and functional magnetic resonance imaging (fMRI) devices in the same subject. All results were mapped onto individual anatomy. Threshold p < 0.001 corrected false discovery rate (FDR), at the voxel level for BOLD signal (a) and HbO signal (b) in a healthy subject. Note that cerebral activation maps provided by DOT data filtered with a Bayesian method follow a similar spatial distribution to maps provided by fMRI, using the general Linear Model (GLM) as a statistical method for both neuroimaging modalities.
As fMRI has been the most commonly used technique in neuroimaging studies, we assume that the information that is provided by fMRI is reliable. Both the literature [24][25][26] as well as T-maps of fMRI help to corroborate the implication of frontal areas for the cognitive task, allowing for the comparison of T-maps from both fMRI and DOT modalities.
The results show that DOT detects the same areas as fMRI at a functional level. These results can be used to corroborate the reliability of Bayesian filtering applied on the frontal cortex through a cognitive task crossing barrier, such as the presence of frontal sinus where scalp-brain distance vary across the subject [13]. Additionally, cognitive signals generate more subtle relative activation changes, unlike motor or visual tasks. Even so, we still obtain increased statistical significance using cognitive changes in a strict sense, as has been shown in prior works [21]. Our results show that band pass filtering alone does not provide statistical significance in the T-maps; therefore, the comparison with fMRI cannot do it. Previous studies have used this linear filtering to detect brain activations in the frontal lobe [8,27,28]. However, they use a simpler topography approach with fewer sourcedetector pairs (~16) and without complex mathematic models to reconstruct three-dimensional (3D) images, in contrast to DOT. This is a clear example of how data filtering affects image reconstruction.
The high fidelity of Bayesian filtering on raw DOT data has been tested using motor tasks to generate cerebral activity in temporal lobes [29], even during transcranial magnetic stimulation [13]. Some studies have used the Bayesian method in different steps during the DOT reconstruction, such as in the inverse problem [30], in order to improve the depth accuracy for the forward model [31], or to remove the scalp artefact [32], but most of them have tested the method in phantoms or simulation models. Table 1. Representation of number of voxels common by anatomical area between HbO-BOLD (T: 2.33, p < 0.05) Here, we present the new approach on the raw DOT data recorded in the human brain. Other authors have modeled and compared with the BOLD signal as well [33]; however, they have applied the method in other steps during the DOT image reconstruction and have used DOT equipment with different characteristics than in the present work. The physiological changes (cardiac and respiratory) must be recorded during the experimental period in order to apply the Bayesian filter, because they are used as regressors and removed from the raw data. Future studies could record other variables such as the motion artifact time variable, which could also be used as a regressor to improve the quality of the results. Finally, we are inclined to think that the results would be similar even if the physiological modeling approach used phases instead of frequencies.

Conclusions
The problem of filter selection during DOT data processing is solved here with the use of Bayesian filtering, thereby standardizing the data filtering as an independent procedure without the influence of the user, paradigm or cerebral area under study. Hence, false interpretation of the results, which sometimes occurs due to the interference of background physiological noise with neural activity, is avoided. Here we describe the Bayesian modeling using a cognitive example, which shows a 59% overlap between both image modalities, considering that each one has a different sampling rate and measures different molecules present in the blood.