## 1.

## Introduction

In recent years there has been an emergence of an assortment of imaging modalities for noninvasively studying the brain. Among these, functional magnetic resonance imaging (fMRI)^{1, 2, 3, 4} and diffuse optical tomography (DOT)^{5, 6, 7, 8} are two techniques that have been developed largely in parallel to study cerebral functional hemodynamic responses. While both of these technologies are being applied successfully to a wide range of similar neuroscience and clinical topics, there are intrinsic limitations to each method, which are imposed by the governing physics of each technology (reviewed in Refs. 9, 10). For example, while fMRI techniques such as blood oxygen level dependent (BOLD) can provide a measurement of blood oxygen saturation changes with fairly high spatial resolution (typically
$2\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}4\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
for functional studies), these signals are physiologically ambiguous, owing to the indirect relationship between changes in the transverse relaxation rate of hydrogen nuclei
$\left(\Delta {R}_{2}^{*}\right)$
and physiological hemodynamic parameters (i.e., deoxyhemoglobin and blood oxygen saturation) (reviewed in Ref. 11). Although such ambiguity does not impede the use of BOLD for mapping the spatial patterns of evoked changes, this does limit the use of BOLD to directly relate physiological parameters between subjects without additional calibration methods. Calibration of the BOLD signal is possible by inducing isometabolic changes in cerebral blood flow using hypercapnia or similar vasoactive agents.
^{12, 13, 14, 15, 16, 17, 18} However, these hypercapnic-calibration methods require the subject to inhale increased levels of carbon dioxide gas for prolonged periods of time (up to several minutes). This procedure is both technically challenging and subject to several possible sources of systematic error^{19} that may render the technique difficult to translate to clinical applications. While the use of hypercapnia-calibrated fMRI techniques to provide quantitative measurements of blood oxygen saturation changes has been important in applying MR techniques to study metabolism, an alternative to hypercapnia calibration is needed to make the estimation of functional
$\mathrm{CMR}{\mathrm{O}}_{2}$
changes more routine. As
$\mathrm{CMR}{\mathrm{O}}_{2}$
is more directly related to neural-metabolic coupling, these measurements could have significant impact in better understanding the connections between neural and hemodynamic function in health and disease (reviewed in Ref. 20).

Continuous wave (cw)-based DOT has several complementary features to fMRI methods, including the ability to record a spectroscopic measurement of both oxygenated (oxyhemoglobin,
$\mathrm{Hb}{\mathrm{O}}_{2}$
) and deoxygenated (deoxyhemoglobin, HbR) forms of hemoglobin. In comparison to fMRI, optical methods generally have very high temporal resolutions, with acquisition rates capable of more than several hundred hertz. This resolution is much faster than needed to capture the typical slow evoked responses and fast enough to prevent aliasing of systemic physiological signals, such as cardiac pulsation and other physiology, which can be a major source of noise in fMRI studies due to undersampling.^{21, 22} A drawback of the DOT technology is its lower spatial resolution, which is intrinsically limited by the propagation of photons through highly scattering biological tissue (reviewed in Ref. 23) and by the typically low number of optical measurement pairs recorded. Although DOT has the theoretical potential to provide quantitatively accurate measurements of hemoglobin concentration changes in the brain, in practice this can seldom be achieved because of the partial-volume effects introduced by the low spatial resolution and depth sensitivity of this method. In addition, the tomographic reconstruction of hemoglobin changes from optical measurements is generally an underdetermined and ill-posed inverse problem.^{7} Tomographic images can be improved with a greater number of measurement combinations, including overlapping measurements to provide more uniform sensitivities;^{24, 25} however, regularization schemes must still be used to constrain the image reconstructions of the underlying absorption changes. In general, the accuracy of these reconstructions depends on the method and amount of regularization applied. In recent years, a great deal of attention has been given to this topic (reviewed in Ref. 26); however, more work is still needed. One promising approach—the incorporation of prior knowledge of the spatial location of the hemodynamic change by either anatomical-based^{8, 27, 28, 29} or functionally-based priors^{30}—improves the quantitative ability of DOT by constraining the solutions to the image reconstruction problem, and thus minimizing the errors introduced by partial-volume effects. With respect to optical imaging of the brain, the use of functional MRI data as such a statistical prior for the location of brain activation area has been suggested to improve DOT reconstructions.^{31} While it is believed that the introduction of statistical priors from structural or functional MRI may improve the localization of the optical signal, the implementation of such methods still has several unresolved issues. In particular, regularized reconstructions require a choice for the proper weight of the prior, as recently reviewed in Gibson, Hebden, and Acridge.^{26} In one extreme, the use of a strict (hard) prior (e.g., Ref. 32) will produce images with identical spatial resolution as the original prior (e.g., the functional MR image). However, this constraint assumes that the value, location, and boundaries of the prior have negligible uncertainty. Although the signal quality of fMRI images has greatly improved in recent years due to advancements in pulse-sequence design, RF coil design, and a move to higher magnetic field strengths, background physiology, intertrial variability, and other subject-related factors are still non-negligible sources of error in these measurements and will contribute to uncertainty in a fMRI-based prior. On the other hand, the use of a statistical prior (e.g., Refs. 30, 33), while favorable in respect to the inclusion of the statistical uncertainty about the prior MRI information, requires knowledge of the proper statistical weight for the constraint. The optimal choice of this weighting depends on the relative measurement noise in both the fMRI and optical signals, and requires a proper statistical model of measurement noise. Concurrent multimodal measurements are unique in that physiological noise (for example, intertrial variability of the evoked response) is simultaneously recorded by each modality, while measurement noise is usually independent between instruments. This property of concurrent measurements provides an opportunity to use mutual information within multimodal measurements to help define the optimal statistical weighting of each modality in a joint image reconstruction. This concept of a bottom-up data fusion model has been previously introduced for neural imaging methods such as multimodal electroencephalography (EEG) and magnetoencephalography (MEG),^{34, 35} but has not yet been demonstrated for multimodal hemodynamic measurements or optical methods. In this work, we describe a new analysis method for fusion of simultaneously acquired DOT and BOLD data that provides a joint estimate of the underlying physiological contrast giving rise to the concurrent measurements from both modalities. This approach makes use of the statistical properties of concurrent measurements and the commonality of the underlying physiology and fluctuations giving rise to these measurements. We use a Bayesian framework to jointly estimate brain activation changes from MR and optical using a single image reconstruction step. This approach enables us estimate oxy- and deoxyhemoglobin changes in the brain, with better spatial accuracy than DOT image reconstructions alone through the incorporation of time-varying spatial information from BOLD observations. Because the fMRI information constrains the spatial extent of the reconstruction, this helps to correct partial volume errors associated with optical reconstructions alone. Likewise, the spectroscopic information of the optical data defines the deoxyhemoglobin calibration of the BOLD signal. We find that the resulting fusion images contain quantitative information about micromolar changes in hemoglobin based on the cross-calibration of these two modalities.

We first present numerical simulations to examine the quantitative accuracy of hemoglobin estimates by our data fusion methods. Next, we apply the model to experimental data recorded simultaneously with DOT and BOLD imaging during a $2\text{-}\mathrm{s}$ duration finger-walking task in five subjects.

## 2.

## Theory

## 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.

## 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
.

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

## Eq. 1

$${t}_{2}\left(t\right)=\frac{2}{{e}^{-1}}[\left(\frac{t-\tau}{\sigma}\right)-{\left(\frac{t-\tau}{\sigma}\right)}^{3}]\cdot \mathrm{exp}\left[\frac{-{(t-\tau )}^{2}}{{\sigma}^{2}}\right].$$^{40}or statistical parametric mapping (SPM)

^{41}software packages). Since the temporal dynamics of the oxy- and deoxyhemoglobin responses are known to differ, we used separate timing parameters ( $\sigma $ and $\tau $ ) for each hemoglobin species and denote these with subscript $\mathrm{Hb}{\mathrm{O}}_{2}$ and HbR, respectively. These timing parameters were empirically estimated by a nonlinear fit to the group average of the region-of-interest average of the DOT time courses as a preprocessing step of this analysis. The empirical values of $\tau $ and $\sigma $ used in the model were ( $0.1\phantom{\rule{0.3em}{0ex}}\mathrm{s}\pm 0.4\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , $6.7\phantom{\rule{0.3em}{0ex}}\mathrm{s}\pm 0.3\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ ) and ( $1.8\phantom{\rule{0.3em}{0ex}}\mathrm{s}\pm 0.4\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ , $6.7\phantom{\rule{0.3em}{0ex}}\mathrm{s}\pm 0.3\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ ) for oxy- and deoxyhemoglobin, respectively (standard errors estimated from the five individual subjects). We note that the inclusion of the second dispersion support in the linear model allows for sufficient flexibility in the temporal shape of the estimated response to model each of the individual subject’s data, as demonstrated in previous MRI studies (e.g., Ref. 42). For example, in this study, we found that a linear combination of these two temporal functions can account for most of the evoked responses in each of five subjects ( ${R}^{2}\mid \mathrm{Hb}{\mathrm{O}}_{2}=0.81\pm 0.08$ and ${R}^{2}\mid \mathrm{HbR}=0.86\pm 0.03$ ; $p<0.001$ for all; for the average of five $\text{subjects}\pm \mathrm{StdErr}$ ).

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]:

## Eq. 2

$${\mathbf{T}}_{\text{Functional},\mathrm{Hb}{\mathrm{O}}_{2}}=\left[\begin{array}{c}{\mathbf{t}}_{1,\mathrm{Hb}{\mathrm{O}}_{2}}\\ {\mathbf{t}}_{2,\mathrm{Hb}{\mathrm{O}}_{2}}\end{array}\right]\otimes \mathbf{U}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{0.3em}{0ex}}{\mathbf{T}}_{\text{Functional},\mathrm{Hb}\mathrm{R}}=\left[\begin{array}{c}{\mathbf{t}}_{1,\mathrm{Hb}\mathrm{R}}\\ {\mathbf{t}}_{2,\mathrm{Hb}\mathrm{R}}\end{array}\right]\otimes \mathbf{U},$$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.^{43} For each of these basis groups, we paired the spatial basis with a series of sine and cosine functions (
$1\u221520\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}1.0\phantom{\rule{0.3em}{0ex}}\mathrm{Hz}$
in
$1\u221520\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:

## Eq. 3

$$\mathbf{G}\mathbf{T}=\left[\begin{array}{cccc}({G}_{F,\mathrm{Hb}{\mathrm{O}}_{2}}\otimes {\mathbf{T}}_{F,\mathrm{Hb}{\mathrm{O}}_{2}})& ({G}_{B,\mathrm{Hb}{\mathrm{O}}_{2}}\otimes {T}_{B,\mathrm{Hb}{\mathrm{O}}_{2}})& ({G}_{C,\mathrm{Hb}{\mathrm{O}}_{2}}\otimes {T}_{C,\mathrm{Hb}{\mathrm{O}}_{2}})& ({G}_{S,\mathrm{Hb}{\mathrm{O}}_{2}}\otimes {T}_{S,\mathrm{Hb}{\mathrm{O}}_{2}})\\ ({G}_{F,\mathrm{Hb}\mathrm{R}}\otimes {\mathbf{T}}_{F,\mathrm{Hb}\mathrm{R}})& ({G}_{B,\mathrm{Hb}\mathrm{R}}\otimes {T}_{B,\mathrm{Hb}\mathrm{R}})& ({G}_{C,\mathrm{Hb}\mathrm{R}}\otimes {T}_{C,\mathrm{Hb}\mathrm{R}})& ({G}_{S,\mathrm{Hb}\mathrm{R}}\otimes {T}_{S,\mathrm{Hb}\mathrm{R}})\end{array}\right],$$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:

## Eq. 4

$$\left[\begin{array}{c}\Delta \mathrm{Hb}{\mathrm{O}}_{2\mid 1}\\ \vdots \\ \Delta \mathrm{Hb}{\mathrm{O}}_{2\mid k}\\ \Delta \mathrm{Hb}{\mathrm{R}}_{1}\\ \vdots \\ \Delta \mathrm{Hb}{\mathrm{R}}_{k}\end{array}\right]=\left[\begin{array}{c}\Delta \mathrm{Hb}{\mathrm{O}}_{2}\\ \Delta \mathrm{Hb}\mathrm{R}\end{array}\right]=\mathbf{G}\mathbf{T}\cdot \beta ,$$## 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.1.

#### Blood oxygen level dependent measurement model

The BOLD signal has a complex origin arising from both intra- and extravascular tissue.^{4, 48, 49} These signals depend on not only the user acquisition parameters, such as echo time, magnetic field strength, and imaging echo type, but on features of the subject anatomy as well, such as vascular architecture and the orientation between blood vessels and the imaging fields.^{50} In this work, we simplify the BOLD measurement model by considering only the extravascular (EV) signal contribution. The EV signal is believed to compose the majority component of the BOLD signal at the 3 tesla (T) field strength.^{49} The EV signal is linear to changes in the concentration of deoxyhemoglobin per volume tissue and is given by the equation^{4, 48, 51}

## Eq. 5

$$\frac{\Delta S\left(t\right)}{{S}_{0}}=\text{BOLD}\left(t\right)\approx -4.3\cdot {\upsilon}_{o}\cdot {V}_{o}\cdot {E}_{o}\cdot TE\cdot \frac{\Delta \left[\mathrm{Hb}\mathrm{R}\left(t\right)\right]}{\left[\mathrm{Hb}{\mathrm{R}}_{0}\right]},$$^{52}). ${E}_{o}$ is the resting oxygen extraction fraction, and ${V}_{o}$ is the baseline blood volume to tissue fraction. $TE$ is the echo time used in the MR pulse sequence ( $30\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$ in this study). $\left[{\mathrm{HbR}}_{0}\right]$ is the baseline deoxyhemoglobin concentration. In general, the coefficients of this equation are not known with the certainty required for quantitative imaging, and moreover, may vary between subjects according to the baseline state or spatial region.

^{53}For this reason, empirically calibrated fMRI techniques have been proposed to remove the influences of these parameters by using ratios of the blood flow and oxygen saturation (i.e., BOLD) signals.

^{12}

To account for these unknown calibration parameters in the BOLD measurement equation, these factors are lumped into a single parameter $(\alpha \equiv -4.3\cdot {\upsilon}_{o}\cdot {V}_{o}\cdot {E}_{o}\cdot TE\u2215\left[{\mathrm{HbR}}_{0}\right])$ . In our method, we use the information in the optical data to help determine the BOLD calibration factor, as we describe in Methods in Sec. 3. In our measurement model, the BOLD signal is a projection of the deoxyhemoglobin component of the model plus an additive noise term $\left({v}_{\text{BOLD}}\right)$ by the equation

## Eq. 6

$${\mathbf{Y}}_{\text{BOLD}}\left(t\right)=[0\phantom{\rule{0.3em}{0ex}}\alpha \cdot \mathbf{I}]\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\}+{\upsilon}_{\text{BOLD}}\left(t\right).$$## 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}

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\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(t\right)\\ {\mathbf{Y}}_{\mathrm{DOT}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(t\right)\\ {\upsilon}_{\mathrm{DOT}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(t\right)\end{array}\right],$$^{7}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

## Eq. 8

$$\mathbf{L}=\left[\begin{array}{cc}{\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ 0& \alpha \cdot \mathbf{I}\end{array}\right].$$## Eq. 9

$$\left[\begin{array}{c}{\mathbf{Y}}_{\mathrm{DOT}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\mathbf{Y}}_{\mathrm{DOT}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\mathbf{Y}}_{\text{BOLD}}\end{array}\right]=(\mathbf{L}\otimes {I}_{n})\cdot \mathbf{G}\mathbf{T}\cdot \beta +\left[\begin{array}{c}{\upsilon}_{\mathrm{DOT}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\upsilon}_{\mathrm{DOT}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\upsilon}_{\text{BOLD}}\end{array}\right],$$## Eq. 11

$$\mathbf{Y}={[{\mathbf{y}}_{1}^{T}\phantom{\rule{0.3em}{0ex}}{\mathbf{y}}_{2}^{T}\phantom{\rule{0.3em}{0ex}}\cdots \phantom{\rule{0.3em}{0ex}}{\mathbf{y}}_{n}^{T}]}^{T}.$$## 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^{59}

## Eq. 12

$$\widehat{\beta}={({\mathbf{H}}^{T}\cdot {\mathbf{R}}^{-1}\cdot \mathbf{H}+{\Gamma}^{2}\cdot {\mathbf{O}}^{-1})}^{-1}\cdot {\mathbf{H}}^{T}\cdot {\mathbf{R}}^{-1}\cdot \mathbf{Y}.$$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
^{60} 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.,

## Eq. 13

$$\mathbf{R}=\left[\begin{array}{ccc}\mathrm{var}\left({\upsilon}_{\mathrm{DOT}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\right)& 0& 0\\ 0& \mathrm{var}\left({\upsilon}_{\mathrm{DOT}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\right)& 0\\ 0& 0& \mathrm{var}\left({\upsilon}_{\text{BOLD}}\right)\end{array}\right].$$## 3.

## Methods

## 3.1.

### Experimental Methods

The experimental protocol used in this study has been previously described.^{61} 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.^{62} The DOT imager has 18 lasers—nine lasers at
$690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
$\left(18\phantom{\rule{0.3em}{0ex}}\mathrm{mW}\right)$
and nine at
$830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
$\left(7\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
(18 and
$7\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\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)=-\mathrm{log}[{I}^{\lambda}\left(t\right)\u2215{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,^{45} 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\phantom{\rule{0.3em}{0ex}}\text{tesla}$
Allegra scanner (Siemens Medical Systems, Erlangen Germany). fMRI data were collected using a gradient echo planar imaging (EPI) sequence
$[\mathrm{TR}\u2215\mathrm{TE}\u2215\alpha =500\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u221530\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u221590\phantom{\rule{0.3em}{0ex}}\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}\u2215\mathrm{TI}\u2215\mathrm{TE}\u2215\alpha =2530\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22151100\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22153.25\phantom{\rule{0.3em}{0ex}}\mathrm{ms}\u22157\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
]. The BOLD time series was preprocessed using a motion-correction algorithm^{64} 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\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ and deep $(\sim 25\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ from the surface of the head (shown in Fig. 4 ). The hemodynamic response was simulated with a maximum peak value of $8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ and $-2.5\phantom{\rule{0.3em}{0ex}}\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.

## 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.

## Eq. 14

$$\left[\begin{array}{c}\Delta O{D}_{\text{BOLD}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(\alpha \right)\\ \Delta O{D}_{\text{BOLD}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\left(\alpha \right)\end{array}\right]=\left[\begin{array}{cc}{\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{830\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\\ {\epsilon}_{\mathrm{Hb}{\mathrm{O}}_{2}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}& {\epsilon}_{\mathrm{Hb}\mathrm{R}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\cdot {\mathbf{A}}^{690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}}\end{array}\right]\cdot \left[\begin{array}{c}0\\ {\widehat{\mathbf{Y}}}_{\text{BOLD}}\u2215\alpha \end{array}\right].$$^{65, 66}and compared to the measured optical data to find the calibration factor $\left(\alpha \right)$ using the equation

## 15.

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
$+\u2215-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.

## Results

## 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\u2215\Gamma =\left(1\phantom{\rule{0.3em}{0ex}}\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,^{61} 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\u2215\Gamma )$
of
$10\phantom{\rule{0.3em}{0ex}}\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).

## 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 .

DOT-only | Data fusion | BOLD calibration (α) | |||
---|---|---|---|---|---|

Subject | $\Delta \mathrm{Hb}{\mathrm{O}}_{2}$ | $\Delta \mathrm{Hb}\mathrm{R}$ | $\Delta \mathrm{Hb}{\mathrm{O}}_{2}$ | $\Delta \mathrm{Hb}\mathrm{R}$ | |

A | $5.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-2.9\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $6.1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.2\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.41\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{M}$ |

B | $0.3\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $3.9\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.52\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{M}$ |

C | $3.1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $4.4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.7\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.24\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{M}$ |

D | $4.9\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-1.0\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $4.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.3\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.19\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{M}$ |

E | $7.3\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-3.0\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $1.6\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.39\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{M}$ |

Group | $4.2\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-2.0\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $4.1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.9\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$ | $-0.55\phantom{\rule{0.3em}{0ex}}\%\text{-}\text{BOLD}\u2215\mu \mathrm{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\u2215\Gamma =10\phantom{\rule{0.3em}{0ex}}\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
^{43} using optical measurements and Glover, Li, and Ress^{67} 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}\u2215\mu \mathrm{M}\pm 0.40\%\text{-}\mathrm{BOLD}\u2215\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\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}100\phantom{\rule{0.3em}{0ex}}\mu \mathrm{M}$
,^{70} 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\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}-0.9\%\text{-}\mathrm{BOLD}\u2215\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
^{12} 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}\u2215\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.

## Discussion

Concurrent multimodal measurements are observations of common underlying changes in underlying functional contrast. Our fusion model combines this mutual information from optical and fMRI modalities into a joint estimate of underlying hemoglobin changes. This approach provides a framework to combine the advantages of the high spatial and temporal resolutions contained between both modalities. In this work, we have developed a model that incorporates the biophysical principles that describe the relationships between the optical and fMRI measurements and underlying cerebral physiology. Using a single image reconstruction step, we can obtain direct estimates of hemoglobin changes that are simultaneously consistent with all sets of observations. This allows us to use the high spatial resolution of fMRI as a spatial prior to improve the optical reconstruction, while at the same time, to use high temporal and spectroscopic information from the DOT as priors on the reconstruction of the BOLD functional changes. The fusion of the higher spatial information from fMRI measurements and the spectroscopic information of the optical technique provided cross-calibrated estimates of hemoglobin changes. The Bayesian framework used in this model allows us to optimally incorporate data from different sensors based on knowledge of the statistical errors in each measurement type. In comparison to methods using fMRI as a statistical prior for the optical reconstruction, our approach incorporates multimodal information based on the statistical properties of concurrent measurements. We believe that this approach provides a framework to more efficiently consider concurrent multimodal datasets and to utilize the mutual information in the signals to improve the accuracy of estimates of the functional response. The advantages of this approach have been discussed for data fusion of magnetoencephalography and electroencephalography data using similar mutual information models (e.g., Refs. 34, 35). We believe that future extensions of our method will offer similar utilities for hemodynamic imaging.

## 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^{73} and blood flow,^{74} 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\phantom{\rule{0.3em}{0ex}}\text{tesla}$
.^{49} 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\phantom{\rule{0.3em}{0ex}}\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^{77} or Windkessel^{82} models will enable the incorporation of more modalities into the model, such as blood flow.^{83} 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
^{79} proposed a similar state-space model for the time-series analysis of simultaneous BOLD and EEG, which could directly extend the work preformed here.

## 6.

## Conclusions

Concurrent multimodal measurements have unique statistical properties, and while an increasing number of publications have implored multimodal acquisition methods, further advancements need to be made to improve specific analysis techniques optimized for this unique form of data. In this work, we describe a fusion model that allows the direct reconstruction of hemoglobin changes from simultaneous optical and fMRI data. This model combines the advantages of both optical and fMRI methods, and allows the estimation of hemoglobin changes with improved temporal, spatial, and quantitative accuracy. Using concurrent optical and fMRI measurements, we estimate the calibration of the 3T BOLD signal to be $-0.55\%\pm 0.40\%$ signal change per micromolar change of deoxyhemoglobin. This is the first demonstration of a multimodal-based calibration of the BOLD signal.

## Acknowledgments

Huppert was funded by the Howard Hughes Medical Institute predoctorial fellowship program. This work was supported by the National Institutes of Health (R01-EB002482, RO1-EB001954, R01-EB006385, T32-CA09502, and P41-RR14075). The authors would like to thank Christiana Andre, Danny Joseph, Div Bolar, and Joe Mandeville, Rick Hoge, and Maria Franceschini for their help in collecting this data. The authors would also like to thank Ville Kolehmainen, Jari Kaipio, and Simon Arridge for useful discussions of this work.