2.1.

Notation

In the following descriptions, we use the notations superscript $T$ for the transpose operator, ⊗ for the Kronecker tensor product, and $\mathbf{I}$ for the identity matrix. In addition, for modality specific operators, a subscript will be used to reference the modality.

2.2.

General Model Description

The functional contrast underlying both the BOLD and DOT signals derives from similar changes in hemoglobin concentrations and blood oxygen saturation. However, the details of the relationships between the measurements and underlying physiology are based on vastly different biophysical principles in the two modalities. These principles must be properly considered to incorporate both MRI and optical datasets into a single model. In Fig. 1 , we present a schematic outline of our state model, which outlines how the measurement models for both MRI and optical are connected to a common model of physiological contrast. The framework for this model is described by three components, which are described in the following sections: 1. a multidimensional linear model by which a set of spatial-temporal basis functions is used to describe functional and systemic components of oxy- and deoxyhemoglobin concentration changes; 2. a set of measurement model equations describing the measurement biophysics for each modality and connecting the observations and underlying physiology; and 3. the least-squares minimization routine in which the experimental observations are fused to create a joint estimate of the underlying physiology.

Fig. 1

Schematic outline of model. The multimodal fusion model is based on a state-space approach and consists of three components. The set of coefficients to be estimated $\left(\beta \right)$ multiplies a convolution kernel of spatial and temporal basis functions in the multidimensional linear model to model the volumetric changes in oxy- and deoxyhemoglobin due to evoked activation and systemic fluctuations. Linear observation models connect the underlying changes in these hemodynamic variables to expected measurements by both DOT and BOLD technologies. The states are finally estimated by minimizing the error with respect to the experimental data using a Bayesian formulation of the linear inversion operation.

2.3.

Multidimensional Linear Model

To describe the underlying physiological changes within the brain, we assume that changes in oxy- and deoxyhemoglobin can be described as linear combinations of a set of spatial and temporal basis functions designated to capture the functional and physiological hemodynamic fluctuations. In principle, this approach is similar to the general linear model (GLM) which has been previously introduced for fMRI^{36, 37} and optical^{38, 39} analysis, but has been modified here to include both spatial and temporal basis supports. Motivated by the anatomy of the brain and head and physiology of hemodynamic fluctuations, we introduced four distinct pairings of spatial and temporal basis functions in our model. These groups were; 1. the functional brain elements (denoted as subscript “functional” or $F$ ), 2. a global brain physiology basis (subscript “brain global” or $B$ ), 3. the cerebral spinal fluid (CSF) layer surrounding the brain (subscript CSF or $C$ ), and 4. the superficial skin layer (subscript “skin” or $S$ ). These four groups of spatial basis functions are diagrammed in Fig. 2 .

Fig. 2

Basis functions in linear model. A set of spatial and temporal basis functions are used to regularize the model estimates of evoked and systemic hemodynamic changes. To model the different aspects of the physiology, four sets of basis sets are used. In the volume corresponding to either gray or white brain matter, a gamma-variant function (and derivative) was used as a temporal basis (subplot C). These allow the modeling of the differing temporal dynamics of both oxy-and deoxyhemoglobin changes. A spatial basis of overlapping Gaussian spheres is used to reduce the dimensionality of the state space and replace a smoothing operation on the data (subplot B). To model the systemic contributions, the volume corresponding to skin, skull, and cerebral spinal fluid (CSF) layers are grouped into basis functions and given sine and cosine temporal functions to model their dynamics.

The first group of canonic functions representing the individual functional brain elements was derived from a spatial basis set composed of overlapping Gaussian spheres positioned on a hexagonal grid (as shown in Fig. 2). This basis is equivalent to a spatial smoothing kernel and is used in place of a separate smoothing operation applied to the data, as is usually typical of fMRI analysis. This basis also serves to reduce the unknown degrees of freedom of the model and makes the model more computationally tractable. We used a $6\text{-}\mathrm{mm}$ standard deviation Gaussian kernel, which generated between 200 and 250 independent spatial functions (depending on the exact anatomy of the subject’s brain). These basis functions are restricted to those voxels that had been identified as either gray or white cortical matter by a tissue segmentation algorithm using the MRI anatomical images (see Methods in Sec. 3). The Gaussian spatial kernels were truncated at the anatomical boundaries between the brain and the CSF layers to impose a cortical constraint on the reconstruction of functional activation. The basis functions and reconstructions were limited to the contralateral hemisphere (opposite to hand movement) where optical measurements were recorded. The matrices representing each of the four types of spatial basis sets (denoted $\mathbf{G}$ in our model) have matrix dimensions of the number of independent basis functions (e.g., 200 to 250 for the functionally associated regions) by the number of voxels in the volume used for the optical and MRI forward models. When analyzing the experimental data, there are 28,672 columns in each of these matrixes (the fMRI images had $64\times 64$ in-plane resolution and seven axial slices). To model the temporal dynamics of evoked functional changes within each of these Gaussian bases, we used a linear combination of a two-parameter ( $\sigma$ and $\tau$ ) modified gamma function and its derivative (dispersion function), as given by the equations

The overall temporal model of the functional component of the hemodynamic signals can be expressed as the convolution of the experimental stimulus timing $\left(\mathbf{U}\right)$ and the functional impulse response functions [ $t_1$ and $t_2$ in Eq. 1 for either oxy- or deoxyhemoglobin]:

In addition to the first pairing of spatial-temporal basis functions, which is used to model the evoked functional response, the global brain, CSF, and skin groups of temporal basis functions are used as systemic regressors of background physiology. These were included to model nuisance physiological contributions to the BOLD and/or DOT measurements by using larger (“super-voxel”) representations of the brain, CSF, and skin layers, as described in Fig. 2 and similar to the methods previously introduced to model systemic contributions to DOT signals.⁴³ For each of these basis groups, we paired the spatial basis with a series of sine and cosine functions ( $1∕20\text{to}1.0\mathrm{Hz}$ in $1∕20\text{-}\mathrm{Hz}$ steps) to describe the temporally oscillating systemic physiology. The skin-confined basis group models the systemic contributions to predominantly the optical signals and the brain-confined nuisance regressor models physiology common to both modalities.

In total, the set of all spatial-temporal bases is combined into a single convolution operator $\left(\mathbf{G}\mathbf{T}\right)$ , giving a total of eight unique basis groups (four groups for both $\mathrm{Hb}{\mathrm{O}}_2$ and HbR). This matrix is formed by the convolution of the individual spatial $\left(G\right)$ and temporal $\left(T\right)$ basis functions at each time instance for each of the eight groups. The matrix $\mathbf{G}\mathbf{T}$ is given by the equation:

The basis function matrix $\left(\mathbf{G}\mathbf{T}\right)$ describes the spatial-temporal supports for the image reconstruction. Following the standard notation used previously to describe the general linear model in fMRI analysis (i.e., Ref. 44), the spatial and temporal dynamics of the underlying hemoglobin changes are described by a linear sum of the spatial-temporal basis function weighted by a vector of unknown coefficients denoted $\beta$ . The vector containing the modeled hemoglobin changes (both systemic and evoked) at each volume element and time point is thus given by the equation:

2.4.

Observation Models

The second component of our fusion model (as shown in Fig. 1) incorporates the measurement processes that relate the underlying physiological changes [oxy- and deoxyhemoglobin given by Eq. 4] to the observations of each modality. The measurement equations for DOT and fMRI describe the biophysics of each instrument measurement. In this work, we have approximated this biophysics using linear approximations to these measurement equations. The validity of these linear simplifications has been discussed by a number of researchers that have examined the consistency of these DOT and fMRI modalities with these hypothesized underlying biophysical theories by correlating the temporal (reviewed in Ref. 10) and spatial^{45, 46, 47} components of measurements. These works have suggested that linearity between the change in deoxyhemoglobin and the BOLD signal can be justified. Note that although our current model has been formulated using these linear approximations, nonlinear extensions of both the optical and fMRI models are theoretically possible using similar methodologies (although the size of the model and computational memory requirements may currently limit this implementation in practice).

2.4.2.

Diffuse optical tomography measurement model

In the DOT technique, near-infrared light is introduced into the head and propagates through the dense scattering layers of the scalp and skull into the brain. A small fraction of the light introduced eventually exits the head and is recorded a distance away from the source position. Hemodynamic changes in oxy- and deoxyhemoglobin concentrations affect the absorption properties of the brain and thus result in changes in the intensity of light as it migrates along a diffuse trajectory through the head. Using multiple measurements taken by an array of light source and detector positions spaced several centimeters apart (shown in Fig. 3 ), the DOT method can be used to spatially resolve these absorption changes and relate these to changes in oxy- and deoxyhemoglobin. The DOT measurement equation is described by the spatial profile of the light propagation through the head, which determines the sensitivity of these measurements. To derive the distribution of photons in a medium with a complex distribution of absorption and reduced scattering coefficients [ ${\mu }_a\left(r\right)$ and ${\mu }_s^{\prime }\left(r\right)$ , respectively], such as the human head, the photon density must be modeled by empirical means using computerized simulations. In Monte-Carlo-based modeling, the distribution of these photons is statistically modeled based on the probability of an absorption or scattering event at each region of space as described by ${\mu }_a\left(r\right)$ and ${\mu }_s^{\prime }\left(r\right)$ .^{54, 55}

Fig. 3

Diffuse optical tomography. Diffuse optical tomography uses spectroscopic measurements of optical absorption changes to record hemoglobin concentration changes. (a) The optical probe was placed over the subject’s primary motor area. (d) The probe contained four source and eight detector positions spaced $2.9\mathrm{cm}$ apart. (b) and (c) Fiducial markers on the probe were visible in the anatomical MR images. (e) The position of the probe and a segmented layer model for each subject were used to generate the optical sensitivity profiles using Monte Carlo methods.

For brain activation, the changes in the optical properties are generally small (on the order of a few percent), and thus, a linear approximation is believed to be sufficient for predicting the changes in optical measurements produced by localized changes in the optical properties (reviewed in Ref. 23). The spectroscopic forward model for such optical measurements is of the form^{56, 57, 58}

Eq. 7

$$\left[\begin{array}{c}{\mathbf{Y}}_{\mathrm{DOT}}^{830\mathrm{nm}}\left(t\right)\\ {\mathbf{Y}}_{\mathrm{DOT}}^{690\mathrm{nm}}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{\epsilon }_{\mathrm{Hb}{\mathrm{O}}_2}^{830\mathrm{nm}}\cdot {\mathbf{A}}^{830\mathrm{nm}}& {\epsilon }_{\mathrm{Hb}\mathrm{R}}^{830\mathrm{nm}}\cdot {\mathbf{A}}^{830\mathrm{nm}}\\ {\epsilon }_{\mathrm{Hb}{\mathrm{O}}_2}^{690\mathrm{nm}}\cdot {\mathbf{A}}^{690\mathrm{nm}}& {\epsilon }_{\mathrm{Hb}\mathrm{R}}^{690\mathrm{nm}}\cdot {\mathbf{A}}^{690\mathrm{nm}}\end{array}\right]\cdot \left\{\begin{array}{c}\Delta \left[\mathrm{Hb}{\mathrm{O}}_2\left(t\right)\right]\\ \Delta \left[\mathrm{Hb}\mathrm{R}\left(t\right)\right]\end{array}\right\}+\left[\begin{array}{c}{\upsilon }_{\mathrm{DOT}}^{830\mathrm{nm}}\left(t\right)\\ {\upsilon }_{\mathrm{DOT}}^{690\mathrm{nm}}\left(t\right)\end{array}\right],$$

where $\mathbf{Y}$ is the vector of measured optical signal changes for each source-detector pair and ${\mathbf{A}}^{\lambda }$ is a matrix describing the linear projection from absorption changes from within the volume to optical signal changes between each measurement pair. $\lambda$ indicates wavelength dependence (690- and $830\text{-}\mathrm{nm}$ wavelengths were used in this study). Each row of the matrix $\mathbf{A}$ describes the light propagation between a particular optical source and detector pair (shown in Fig. 3), which describes the spatial sensitivity for that measurement. The ${\mathbf{A}}^{\lambda }$ matrix is a projection operator, which integrates absorption changes over the volume to predict the measurements between sources and detectors. It has been shown that this linear operator is approximated by the adjunct product between the photon density distributions for each given source and each given detector involved in each optical measurement.⁷ In Eq. 6, the projection operators are combined with the spectral extinction coefficients $\left({\epsilon }^{\lambda }\right)$ to give a direct forward operator between the concentration changes in $\mathrm{Hb}{\mathrm{O}}_2$ and HbR (per volume element) and the measured optical density changes at multiple wavelengths. In our model, $v^{\lambda }$ is an additive measurement noise term specific to each measurement channel and wavelength.

2.5.

Data Fusion Model

Rather than treating the DOT and BOLD measurements independently or first computing an MR-based spatial map to be used as a reconstruction prior, our model tries to preserve the mutual information in the multimodal data by simultaneously considering measurements from both instruments. We concatenate the two measurement equations [Eqs. 6, 7] into a single joint-observation operator

2.6.

Bayesian Inversion Routine

To estimate the coefficients for the spatial-temporal basis sets describing changes in $\mathrm{Hb}{\mathrm{O}}_2$ and HbR, the linear model [Eq. 10] is solved using a Bayesian formulation of the linear inverse operator⁵⁹

In Eq. 12, ${\mathbf{R}}^{-1}$ and ${\Gamma }^2{\mathbf{Q}}^{-1}$ represent the inverses of the covariance for the observation noise and state noise models, respectively. We have chosen this formulation over the Tikhonov regularization (Moore-Penrose generalized inverse), which is more commonly used in DOT, because the Bayesian formulation separates the regularization parameters into both an observation noise $\left({\mathbf{R}}^{-1}\right)$ and a state noise $\left({\Gamma }^2{\mathbf{Q}}^{-1}\right)$ preconditioning term. A similar regularized inverse scheme was recently described by Yalavarthy ⁶⁰ This weighted least-squares method provides a better framework to introduce differential observation noise for each modality. From a statistical standpoint, ${\Gamma }^2{\mathbf{Q}}^{-1}$ is a covariance prior on the states $\beta$ . This definition offers some intuitive guidance on how to tune the inversion, since we expect physiological fluctuations to be approximately on the order of micromolar magnitudes. In practice, however, the multidimensional linear model is only an approximation of the actual physical system, and the assumed Gaussian distributed prior on $\beta$ is imperfect. These approximations effectively mean that the Bayesian inverse is a regularization scheme that must be tuned, and we have introduced a scalar parameter $\Gamma$ for that purpose. In addition, we define the variance of the observation (measurement) noise model $\left(\upsilon \right)$ as a diagonal matrix formed from the variance [var( )] of each measurement, i.e.,

3.1.

Experimental Methods

The experimental protocol used in this study has been previously described.⁶¹ The data used in this analysis have been previously reported in that paper for comparison of the temporal characteristics of optical methods and fMRI. Five healthy, right-handed subjects (4 male, 1 female) were imaged in this experiment. The task consisted of a two-second duration finger-walking on the right (dominant) hand. The Institutional Review Board at Massachusetts General Hospital approved these procedures.

3.2.

Diffuse Optical Tomography Acquisition and Preprocessing

We used a multichannel continuous-wave optical imager (CW4, TechEn Incorporated, Milford, Massachusetts) to obtain the measurements as previously described.⁶² The DOT imager has 18 lasers—nine lasers at $690\mathrm{nm}$ $\left(18\mathrm{mW}\right)$ and nine at $830\mathrm{nm}$ $\left(7\mathrm{mW}\right)$ —and 16 detectors of which only four source positions and eight detectors were used here. The laser wavelengths were 690 and $830\mathrm{nm}$ (18 and $7\mathrm{mW}$ , respectively).

The DOT probe was made from flexible plastic strips with plastic caps inserted in it to hold the ends of the $10\text{-}\mathrm{m}$ -source/detector fiber optic bundles. The probe consisted of two rows of four detector fibers, and one row of four source fibers arranged in a rectangular grid pattern and spaced $2.9\mathrm{cm}$ between nearest neighbor source-detector pairs (shown in Fig. 3). This plastic probe was then secured to the subject’s head centered over the contralateral primary motor cortex (M1) via Velcro straps and foam padding. The $10\text{-}\mathrm{m}$ fibers were run through the magnet bore to the back of the scanner and through a port into the control room, where they were connected to the DOT instrument.

Following collection and separation of source-detector pairs, the timing of the DOT data was synchronized to the MRI images (as described in Ref. 61). The data were down-sampled using a Nyquist filter to the same $2\text{-}\mathrm{Hz}$ sample frequency as the fMRI. The raw light fluence measurements were converted to changes in optical density by the negative log of the normalized incident light intensity $\{\Delta O{D}^{\lambda }\left(t\right)=-\log [I^{\lambda }\left(t\right)∕I^{\lambda }\left(0\right)]\}$ .

Monte Carlo methods were used to generate the optical sensitivity profiles describing the DOT measurement equation [Eq. 6].^{45, 55, 63} Anatomical MP-RAGE images acquired during the session were segmented into a five-layered head model for each of the subjects,⁴⁵ which was used for these Monte Carlo simulations and also to define the spatial basis functions used in the linear model. From the vitamin E fiducial markers used to mark the optical probe, the locations of the DOT optodes were located. The photon migration paths were sampled from the simulated trajectories of ${10}^8$ photons at both 690- and $830\text{-}\mathrm{nm}$ wavelengths. These simulations were used to calculate the linear optical forward matrices $\left({\mathbf{A}}^{\lambda }\right)$ and substituted into Eq. 7 (see Fig. 3). To create the spectroscopic forward operator, the hemoglobin extinction coefficients were used from the Oregon Medical Laser Center (http://omlc.ogi.edu/).

3.3.

Functional Magnetic Resonance Imaging Acquisition and Preprocessing

BOLD-fMRI measurements were performed using a $3\text{tesla}$ Allegra scanner (Siemens Medical Systems, Erlangen Germany). fMRI data were collected using a gradient echo planar imaging (EPI) sequence $[\mathrm{TR}∕\mathrm{TE}∕\alpha =500\mathrm{ms}∕30\mathrm{ms}∕90\mathrm{deg}]$ with seven $5\text{-}\mathrm{mm}$ oblique orientation slices ( $1\text{-}\mathrm{mm}$ spacing) and $3.75\text{-}\mathrm{mm}$ in-plane spatial resolution. Structural scans were performed using a T1-weighted MPRAGE sequence ( $1\times 1\times 1.33\text{-}\mathrm{mm}$ resolution, $\mathrm{TR}∕\mathrm{TI}∕\mathrm{TE}∕\alpha =2530\mathrm{ms}∕1100\mathrm{ms}∕3.25\mathrm{ms}∕7\mathrm{deg}$ ]. The BOLD time series was preprocessed using a motion-correction algorithm⁶⁴ and mean normalized before being used in the model.

3.4.

Numerical Simulations

We used simulation studies to test the quantitative accuracy of this method. This enables us to examine the ability to quantify hemoglobin changes, since no empirical “gold-standard” method exists that can provide quantitative measurements with which to validate our experimental results. Simulated inclusions were placed shallow $(\sim 18\mathrm{mm})$ and deep $(\sim 25\mathrm{mm})$ from the surface of the head (shown in Fig. 4 ). The hemodynamic response was simulated with a maximum peak value of $8\mu \mathrm{M}$ and $-2.5\mu \mathrm{M}$ for oxy- and deoxyhemoglobin concentration changes, respectively. Varied levels of measurement noise were added to the simulated measurements to create varied levels of contrast-to-noise ratios.

Fig. 4

Characterization of model quantitative accuracy. The quantitative accuracy of the model was examined with observations simulated from a shallow ( $18\mathrm{mm}$ , toprow) and deep ( $25\mathrm{mm}$ , bottomrow) inclusion. The simulated inclusions are shown contoured (black) in each image (axial projections). The boundaries of the cortex are outlined in white. Measurement noise was added to achieve a 5:1, 50:1, and 500:1 contrast-to-noise ratio (CNR). Hemoglobin changes were reconstructed using the fusion model and DOT data alone for comparison. Reconstructions were performed at various regularization amounts to examine the dependence of the recovered magnitude of the response on the regularization parameter. The images are shown for the reconstruction at $1∕\Gamma =1\mu \mathrm{M}$ (indicated by the vertical red line) for the $\mathrm{CNR}=50:1$ case.

3.5.

Optical Calibration of the Blood Oxygen Level Dependent Signal

The BOLD measurement model depends on an unknown calibration factor $\left(\alpha \right)$ , which depicts the several baseline and structural unknown factors within Eq. 5. This calibration factor is determined empirically by comparing the spatial and temporal profiles of the DOT and BOLD data as described in Ref. 45. Using the linear forward operator describing the DOT measurement model [Eq. 7], we project the model estimate of the BOLD signal from its natural volume space representation into the optical measurement (i.e., source-detector) space. By substituting Eq. 6 into Eq. 7, the expected optical signals due to a change in deoxyhemoglobin can be modeled from the model estimated BOLD signal, as discussed in Ref. 45.

In the modified Beer-Lambert equation [top in Eq. 15], $L$ is the linear distance between the position of a source and detector, and $DPF$ is the differential path-length factor, which is a dimensionless coefficient defined as the effective path length through the head divided by the source-detector separation.^{63, 65, 66} This is calculated from the Monte Carlo simulations. Since this projection uses the optical forward model, it accounts for the DOT partial volume errors in the calibration of the BOLD signal. A single calibration factor $\left(\alpha \right)$ was used, which is consistent with findings of the spatial correlation of response amplitudes between the two modalities.^{45, 46} To estimate this parameter in the model, we employ an iteration routine between state and parameter estimation. Alpha $\left(\alpha \right)$ is first estimated via Eq. 15 using the raw optical and MRI data, and then the states $\left(\beta \right)$ are estimated using Eq. 12. The estimated optical and BOLD signals are recovered from the forward model using the estimated value of beta, and the alpha is re-estimated from the modeled BOLD and optical signals using the bottom of Eq. 15. This beta and alpha estimation process was repeated for 20 iterations. Alpha was found to converge below a $+∕-10\%$ variation after a few iterations (approximately five iterations).

3.6.

Estimation of Noise Covariance

To estimate the covariance of the measurement noise $\left(\mathbf{R}\right)$ , we calculated the variance of the residual for each measurement by iteratively solving the state estimate and estimating $\mathbf{R}$ to be the variance of the residual of the measurements. The initial seed of $\mathbf{R}$ was calculated from linear regression of the data with the temporal basis.

An identity matrix was used for the state covariance matrix $\mathbf{Q}$ . The regularization tuning parameter $\left(\Gamma \right)$ was adjusted in the analysis, as noted in Sec. 4. A single regularization parameter was used to scale the variance of both the oxy- and deoxyhemoglobin states.

4.1.

Simulation Studies

To investigate the improvements in the reconstruction accuracy, forward simulations of synthetic data and inverse reconstructions were first performed as described in Sec. 3. The simulated observation vectors were reconstructed using the DOT only and the multimodal fusion (DOT and BOLD) data in the reconstructions. Reconstructed axial slices are shown in Fig. 4 for the shallow and deep simulated inclusions. Using only the DOT data, the image reconstruction of both oxy- and deoxyhemoglobin depend strongly on the regularization parameter used $\left(\Gamma \right)$ . In Fig. 4, reconstructions are represented at regularization values of $1∕\Gamma =\left(1\mu \mathrm{M}\right)$ . To examine the dependence of the reconstruction accuracy on the regularization parameter, we looked at the reconstructed response amplitude and goodness-of-fit of the model for varied levels of regularization. This was repeated for both the shallow and deep inclusions and at several levels of simulated instrument noise (5:1, 50:1, or 500:1 contrast-to-noise ratio). The results, shown in Fig. 4, demonstrate the robustness of the fusion model to the choice of this regularization for deoxyhemoglobin reconstructions. Since the estimate of deoxyhemoglobin change is overdetermined in the presence of both optical and BOLD measurements, accurate estimates of the response are recovered at minimal regularization. When over-regularization is applied, the response is underestimated, as expected. In contrast, reconstructions using only the DOT data show high dependence on the regularization amount. This result is consistent with prior expectations of the behavior of the DOT inverse problem with our measurement geometry. The fusion reconstruction of oxyhemoglobin changes is also dependent on the regularization applied. This is due to the fact that the reconstruction of oxyhemoglobin changes is only indirectly informed by the fMRI data through the off-diagonal terms in the spectroscopic optical forward model [Eq. 7]. Since only the DOT measurements contain direct information about oxyhemoglobin, this component of the inverse problem is still ill posed and requires significant regularization. However, the oxyhemoglobin reconstructions are still modestly improved by the fusion reconstruction over the DOT-only reconstruction.

4.2.

Experimental Studies

Images of evoked hemoglobin changes were reconstructed from the BOLD alone, DOT alone, and multimodal (fusion) experimental datasets for all five subjects, as shown in Figs. 5 and 6 . In Fig. 5, we show an in-depth look at the results for a single subject (subject A), and show the reconstructions of deoxyhemoglobin (upper rows) and oxyhemoglobin (bottom rows). The right column of images in Fig. 5 shows the surface rendering of the activation regions with the central sulcus marked for clarification. In Fig. 6, we show the reconstructed results for the data-fusion method for subjects B through E and comparison to the BOLD model. In agreement with previous analysis of this data,⁶¹ we found localized functional activation in the motor-cortex (M1; precentral gyrus) and primary sensory (S1; post-central gyrus) regions for each of the individual five subjects for the fMRI and fusion reconstructions. The peak amplitudes of the estimated evoked oxy- and deoxyhemoglobin responses are given in Table 1 . In general, the multimodal fusion and BOLD images were similar, although a few notable differences were observed. These differences are most likely the result of differences in the regularization in the two models. A regularization amount $(1∕\Gamma )$ of $10\mu \mathrm{M}$ was used in the DOT and fusion reconstructions, and 10% in the BOLD images, which may account for these differences. In addition, the fusion images are more sensitive than single modality methods to both temporal and spatial registration errors between the two modalities, which can include intrinsic differences due to differential vascular sensitivity. This can introduce disinformation between modalities, which will result in the loss of confidence of an activation event and lower functional effects statistics (e.g., more likely to reject a null hypothesis).

Fig. 5

Direct hemoglobin reconstructions from empirical data. This figure shows (b) the model reconstructions using the DOT alone, BOLD alone, and (a) fusion datasets for deoxy- and oxyhemoglobin changes. Data are represented for subject A. All responses have been normalized for comparison. Reconstructions were performed at a regularization level of $1∕\Gamma =10\mu \mathrm{M}$ for the DOT and fusion models and $1∕\Gamma =10\%$ for the BOLD model. Functional changes are shown masked below the half-max response amplitude. In (a) the approximate location of the optical probe is shown in green. The right-most images show the maximum intensity projections of the activation patterns onto the brain pial surface [surfaces generated using Freesurfer (The Massachusetts General Hospital; see http//surfer.nmr.mgh.harvard.edu)] and displayed using NeuroLens [Université de Montréal; (see http://www.neurolens.org)]. The central sulcus is outlined (green line) in each surface image.

Fig. 6

Comparison of fusion and MRI image reconstructions. In this figure, we show the BOLD (a) and fusion reconstructions (b) for subjects B through E (left to right). Subject A was presented in Fig. 5. Regularization was applied as described for Fig. 5. All images are shown normalized to the maximum response and masked below the half-max amplitude. Axial and coronal projections are shown. The outline of the boundary between the brain and superficial layers is shown in green. The blue rectangle indicates the location of the imaging volume from the fMRI.

Table 1

Comparison of reconstructed hemoglobin amplitudes. The maximum (minimum) response amplitude was calculated for a region of interest in each of the DOT alone and fusion reconstructions for the five subjects. Regions of interest were manually defined in the motor cortex area. The BOLD calibration factor (α) was calculated as described in the text. The values shown are for the reconstructions at a regularization 1∕Γ=10μM .

In the model reconstructions that only used the DOT data, we found that the amplitudes of the estimated hemoglobin changes were dependent on the regularization applied, which was consistent with the simulation results and prior expectations on the behavior of the ill-posed inversion of the optical forward model and minimum norm estimator. The reconstructions shown in Figs. 5 and 6 and values given in Table 1 were obtained at a regularization value of $1∕\Gamma =10\mu \mathrm{M}$ , which provided the best reconstruction from the optical data alone when qualitatively compared to the fusion or BOLD reconstructions. We have chosen to present images using this regularization point to give the best possible representation of optical-only reconstructions for comparison to the data fusion method. With this choice of regularization, the optical-only results were qualitatively close to the fMRI images (as demonstrated in Fig. 5 for subject A), although distinct biases in the spatial locations were notable. Even with the cortical constraint of our model, we noted that our optical-only reconstructions tended to be biased toward the locations of optodes, which, in most cases, displaced the location of the DOT reconstructions relative to the BOLD alone or fusion derived estimates. In contrast to the deoxyhemoglobin reconstructions, the improvements of the fusion model to estimate of oxyhemoglobin changes over the DOT alone model were less dramatic. However, we found that the fusion estimates of oxyhemoglobin were less susceptible to biases toward the optode positions, although they were still superficially biased toward the head’s surface. The optical probe used in this study was based on a simple nearest-neighbor geometry. Recent work has shown that optical-only reconstructions could be further improved using tomographic probe designs, including the use of overlapping measurements (e.g., Refs. 24, 25). These methods would be expected to further improve the accuracy of the reconstructions performed here for both optical-only and fusion methods.

In addition to the functionally evoked responses, our model also includes regressors for the background physiological oscillations modeled as global effects using large spatial basis supports. The inclusion of these regressors accounted for an average of an additional $8.5\%\pm 4.4\%$ of the total model variance based on partial $R^2$ analysis of the model (average of five $\text{subjects}\pm \mathrm{StdErr}$ ; range 3.3%[subject B] to 26.3%[subject E]). As recently demonstrate by Diamond ⁴³ using optical measurements and Glover, Li, and Ress⁶⁷ using fMRI measurements, the use of direct measurements of systemic physiology by pulse-oximtery or noninvasive blood pressure methods as additional model regressors may provide further reduction of these cerebral physiological signals and should be examined in future work.

4.3.

Optical Calibration of Blood Oxygenation Level Dependent

In Table 1, we show the recovered values of the optically calibrated BOLD scaling factor $\left(\alpha \right)$ for each of the five subjects (mean: $-0.55\%\text{-}\mathrm{BOLD}∕\mu \mathrm{M}\pm 0.40\%\text{-}\mathrm{BOLD}∕\mu \mathrm{M}$ ). The high variance in the group estimate may be the result of differences in the baseline volume or oxygen extraction levels between subjects (particularly subject C). In addition, the values of the BOLD model parameters given in Eq. 5 are based on additional assumptions about the size and orientation of blood vessels (i.e., Ref. 50) and may also contribute to the differences observed between the five subjects. We can compare our measured values of alpha with theoretical estimates derived from Eq. 4 and data from previous literature. Assuming a baseline blood volume fraction of 3 to 5%,^{68, 69} a total hemoglobin concentration of $60\text{to}100\mu \mathrm{M}$ ,⁷⁰ and an oxygen extraction fraction of 30 to 40%,^{69, 71} the expected value for this calibration factor can be estimated to be between $-0.3\text{to}-0.9\%\text{-}\mathrm{BOLD}∕\mu \mathrm{M}\text{-}$ -HbR for the MRI acquisition parameters used in this study. This calibration factor has also been determined experimentally by the hypercapnia method of BOLD calibration described by Davis ¹² The hypercapnia calibration method determines the maximal BOLD change possible if all deoxyhemoglobin were displaced from the region (reviewed in Ref. 19). Thus, according to our approximation of a primarily extravascular water contribution to the BOLD signal, as given by Eq. 5, this hypercapnia calibration factor (usually denoted $M$ ) can be related to our alpha parameter by multiplying by the baseline deoxyhemoglobin levels (i.e., $M\approx \alpha \times \left[{\mathrm{HbR}}_0\right]$ ). The empirical value of $M$ has been reported between 7 to 25% as tabulated in Ref. 19, which corresponds to a value of alpha between $-0.4$ and $-1.4\%\text{-}\mathrm{BOLD}∕\mu \mathrm{M}$ using the ranges of baseline total hemoglobin and oxygen extraction cited before. Indeed, both the empirical and theoretical estimates of this calibration factor are consistent with our values calculated for the five subjects, as shown in Table 1. In future work, it will be necessary to validate the optical-calibration approach against other calibration methods like the hypercapnic method.

5.1.

Comparison of Reconstruction Methods

In the comparisons between the optical only, BOLD only, and fusion reconstructions of the simulated data, we found that fusion methods produced the most accurate estimates of hemoglobin changes both in terms of spatial localization and quantitative accuracy. In particular, we found that for the fusion reconstructions of deoxyhemoglobin, the amplitude of these changes was relatively independent of the magnitude of the regularization applied and were fairly robust to errors in underestimation of the regularization parameter $\left(\Gamma \right)$ , as demonstrated in Fig. 4. In the fusion model, enough information is available to make the deoxyhemoglobin reconstruction problem overdetermined. Even at a minimal regularization level, the quantitatively accurate magnitudes of the deoxyhemoglobin responses were still recovered for both the deep and shallow deoxyhemoglobin inclusions. As expected, an over-regularization of the linear model resulted in underestimation of the hemodynamic response in all models. In contrast to deoxyhemoglobin, which is directly informed by the BOLD model, we found that reconstructions of oxyhemoglobin changes in the fusion were still dependent on the regularization. This is because oxyhemoglobin is only indirectly informed by the BOLD signal and still represents an underdetermined problem.

In comparison, we found that our DOT alone model was generally more superficially biased than the fusion model at recovering the spatial profile of both targets, and that the amplitude was very sensitive to the amount of regularization applied. This bias will result in an underestimation of the response amplitudes, since the optical sensitivity profiles fall exponentially with increasing depth. This finding is consistent with prior expectations concerning the accuracy of DOT reconstructions. Although the spatial basis functions used in this work ensured a cortical constraint to functional activations, the reconstructed amplitudes were underestimated by up to several orders of magnitude, in particular for the deep simulated inclusion, and the magnitude of the DOT estimated hemoglobin changes were highly dependent on the regularization parameters used.

The DOT alone and fusion reconstructions of the empirical data were consistent with the findings from the simulations. The spatial locations of the fusion reconstructions were more consistent with the profiles from the BOLD alone. Although the magnitudes of the changes were comparable in our final images, the DOT reconstructions were highly dependent on the regularization applied, which was chosen here to provide the best possible DOT reconstructions based on the BOLD result. Thus, this may be misleading to the quality of the reconstructions that can be obtained routinely by DOT alone.

5.2.

Optically Calibrated Blood Oxygen Level Dependent Signals

In addition to the Bayesian fusion model that we have introduced in this work, we have also presented a method for using the spectroscopic information of the optical data to provide insight into the BOLD signal, allowing for optically calibrated BOLD imaging. In this work, we have assumed that a single calibration factor can be applied to calibrate the BOLD signal. In both the simulation and experimental results, we found that the estimate of this factor converged after only a few iterations of the model. In the simulation results, we found that we could recover this calibration factor. In the experimental results, we found approximate agreement between our estimates of this factor and the expectation from theoretical work approximated at 3T and the empirical hypercapnia-based BOLD calibration.

5.3.

Model Limitations and Future Extensions

While our data fusion method greatly improves the reconstruction and quantitative accuracy of deoxyhemoglobin changes, the improvements to oxyhemoglobin quantification were more modest. In principle, the state covariance matrix $\left(\mathbf{Q}\right)$ used in the Bayesian pseudoinverse model [Eq. 12] can be used to introduce a statistical prior between deoxy- and oxyhemoglobin maps with the inclusion of off-diagonal elements connecting the two matrix quadrants for oxy- and deoxyhemoglobin. However, in the work described here, we have purposefully not included these off-diagonal terms based on recent work that has suggested that underlying vascular structures may displace the locations of the two chromophores due to differential arterial versus venous weightings.^{45, 46, 47, 72} In support of this, several recent fMRI-based studies have also shown spatial displacements between the BOLD signal and MR measures of blood volume⁷³ and blood flow,⁷⁴ which corroborate a spatial displacement of arterial and venous-weighted measurements. In future studies, we suggest that the incorporation of an MR measure of blood volume changes^{73, 75, 76} into this data fusion model may be more appropriate to constrain total-hemoglobin changes and consequently the quantification of oxyhemoglobin.

An assumption made in this model was that we only examined the extravascular contribution to the BOLD signal. This assumption ignores the contribution from the water molecules within the blood vessels that interact directly with the deoxyhemoglobin heme group.^{1, 4, 48, 77} This assumption is supported by recent work by Lu and van Zijl determined that the extravascular signal composes approximately 67% of the intrinsic relaxation $\left(\Delta R_2^{*}\right)$ at $3\text{tesla}$ .⁴⁹ Further evidence of the dominance of the extravascular signal has also been demonstrated in the empirical comparison of optical and BOLD signals, which have shown a strong temporal correlation between the BOLD signal and the optical measure of deoxyhemoglobin (reviewed by Ref. 10). Thus, the source of BOLD functional contrast at $3\text{tesla}$ is expected to be predominantly extravascular, which we believe justifies the assumption in this model in this current work. However, at lower field strengths, alternative MR acquisition protocols, or for more quantitatively accurate results, the intravascular component may become necessary. The inclusion of the intervascular component of the fMRI signal creates a nonlinearity in the BOLD measurement equation, which additionally depends on changes in the venous blood volume that determine oxygen saturation changes (i.e., Ref. 48). The future extension of this model to incorporate this term can be achieved by the replacement of the linear (extravascular BOLD) measurement model in Eq. 5 with a more detailed state-space model of the vascular physiology (i.e., Refs. 78, 79, 80, 81). In future work, vascular models such as the Balloon⁷⁷ or Windkessel⁸² models will enable the incorporation of more modalities into the model, such as blood flow.⁸³ This nonlinear extension of this bottom-up model could allow the direct estimation of cerebral metabolism changes and provide a means to fuse multimodal data from a broader scope of methods. Recent work by Riera ⁷⁹ proposed a similar state-space model for the time-series analysis of simultaneous BOLD and EEG, which could directly extend the work preformed here.