*a priori*anatomical information acquired by magnetic resonance imaging. We demonstrate the accuracy and speed of the approach using a phantom experiment and through numerical simulation of brain activation in a rat’s head. The applicability of the approach for real-time monitoring of brain hemodynamics is demonstrated through a hypercapnic experiment. We show that our results agree with the expected physiological changes and with results of a similar experimental study. However, by using our approach, a three-dimensional tomographic reconstruction can be performed in ∼

**3 s**per time point instead of the 1 to 2 h it takes when using the conventional finite element modeling approach.

## 1.

## Introduction

Diffuse optical tomography (DOT) is a noninvasive imaging modality that employs near-infrared light to interrogate optical properties of biological tissue.^{1}^{,}^{2} Compared with alternative imaging modalities, such as functional magnetic resonance imaging (fMRI), DOT has several advantages, including simplicity, portability, reduced cost, and faster acquisition speeds. It also has disadvantages, particularly a lower spatial resolution than MRI.^{3}

Reconstructing three-dimensional images of optical properties of tissue using DOT involves solving two distinct problems, the forward and the inverse problem. The first problem is to predict the distribution of light at the detectors using a model of light propagation to tissue. The second problem is to estimate the optical parameter functions describing the optical properties of the tissue, which minimize the difference between the experimental and model-predicted measurements.

Due to the highly scattering nature of tissue at near-infrared wavelengths, DOT requires accurate computational modeling of light propagation to tissue.^{1} For this reason, anatomically realistic three-dimensional (3-D) tissue models, based on *a priori* structural information provided by an alternative imaging modality such as MRI, are increasingly used to maximize the quality of the image reconstruction.^{4}^{,}^{5} On the other hand, as the approximation improves, the complexity of the forward and inverse problems increases, and the fact that the estimation of parameter functions is an infinite-dimensional optimization problem makes DOT a computationally demanding imaging problem.^{6} This approximation-versus-accuracy dichotomy is one of the computational bottlenecks in DOT, which becomes highly relevant for the 3-D case. As a consequence, DOT is not used routinely for bedside monitoring or outside the hospital environment.

Previous efforts to address the computational bottleneck have focused on either reducing the complexity of the forward model while maintaining accuracy (model order reduction) or accelerating the convergence of the numerical reconstruction algorithms. Zhai and Cummer^{7} used model order reduction techniques, common in structural dynamics, to decrease the complexity of the system of equations derived from finite element method (FEM); an adjustable parameter controls the level of the approximation. This solution decreased by an order of magnitude the time that the conventional approach can take. Kolehmainen et al.^{8} employed a Bayesian approximation error model to compensate the loss of accuracy when using a coarse mesh; however, reconstruction times are still far from achieving real-time continuous monitoring.

Alternative approximation schemes for the diffusion equations (DE) using, for example, orthonormal compactly supported wavelets,^{9} which produce sparse stiffness matrices, can reportedly reduce the time required to solve the forward problem ($\sim 36\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{times}$ faster for a 3-D geometry). To speed up parameter estimation, Eames et al.^{10} proposed a sensitivity-based Jacobian reduction technique, which delivers a 14-fold speed increase. Including hard-priors in the reconstruction helps reduce computation time by reducing the dimension of the parameter set to be estimated.^{11}^{,}^{12} However, this choice may lead to erroneous results if the models are not coregistered.^{4} Other strategies employed to reduce the complexity of the reconstruction algorithms include the use of nonlinear multigrid optimization^{13}^{,}^{14} and adaptive mesh techniques.^{13}^{,}^{14} Solving the inverse problem using analytic methods can dramatically reduce the reconstruction time,^{15} but the application of these methods is limited because analytic solutions are only available for some cases, such as the slab or the cylinder.^{16}

This paper proposes an efficient practical alternative to the analytic approach to solving the 3-D optical tomographic reconstruction problem for continuous wave (CW) DOT systems, which takes advantage of *a priori* anatomical information extracted from MRI scans and can be used without restrictions for 3-D domains with complex geometries. The approach overcomes the computational bottlenecks of existing iterative image reconstruction methods, making it possible to use DOT in real time. The proposed method employs forward models of the propagation of light in tissue, which are constructed as parsimonious radial basis function (RBF) maps,^{17} which relate directly the optical parameters of the medium to the outward flux measurements at detector locations. These highly optimized, efficiently computable models dramatically accelerate the image reconstruction algorithms. In a preliminary study, we evaluated an approach based on polynomial models using numerical simulation.^{18} Here, we demonstrate comprehensively the applicability and performance of the superior RBF-based approach using numerical simulation, a phantom study, and by performing 3-D reconstruction of hemodynamics in a rat’s brain during a hypercapnia experiment using an anatomically correct finite element model of the rat’s head, derived from MRI scans. We evaluate both the speed and quality of reconstruction and compare the results with those obtained by a standard FEM approach.

The approach advocated in this paper makes it possible to perform fast high-density tomographic reconstruction with modest computational power, compared to current approaches, but preserving the high accuracy of a fine anatomically correct finite element mesh. The algorithms could be run on a handheld device, thus enabling continuous, real-time monitoring of patients over extended periods at the bedside, inside as well as outside the hospital.

## 2.

## Methods

## 2.1.

### Animal Preparation and Experimental Procedure

Small animals provide excellent models to investigate ischemia,^{19} tissue oxygenation,^{20} and hemodynamics.^{21}22.^{–}^{23} This experiment focuses on the hemodynamic response under hypercapnia.

The animal used was a female hooded Lister rat weighing $\sim 400\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{g}$, kept in a 12 h dark/light cycle environment at a constant temperature of 22°C. Food and water were provided *ad libitum*. The rat was anesthetized with urethane ($1.25\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{g}/\mathrm{kg}$ i.p.), and atropine (0.1 ml) was administered to reduce mucous secretions throughout the course of surgery. Rectal temperature was maintained at 37°C during the surgery and experimental procedures by means of a homoeothermic blanket (Harvard apparatus). The animal was tracheotomized to allow artificial ventilation. Ventilation parameters were adjusted to maintain blood gas measurements within physiological limits. Measurement of mean arterial blood pressure (MABP) was performed after the left and right femoral veins were cannulated; at this stage, phenylephrine (0.13 to $0.26\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mg}/\mathrm{h}$) was infused to keep the blood pressure between physiological limits (MABP: 100 to 110 mmHg).

The skin and underlying muscle were removed from the left side of the rat’s head; the skull was thinned with a drill to achieve better contact with the optode holder, which later was fixed with dental cement. Once the surgery was finished, the animal was placed on a platform in the prone position, fixed via two ear bars and a bite bar.

The optodes were adjusted in a plastic honeycomb structure with 12 holes. The centre-to-centre distance between holes was 4.2 mm.

The rat was ventilated artificially with normal gas mixture (20% ${\mathrm{O}}_{2}$, 80% ${\mathrm{N}}_{2}$). During the hypercapnia experiment, a step increase of 10% carbon dioxide was used, while the oxygen and nitrogen ratio was kept constant. The experiment lasted 300 s: 60 s of baseline followed by an induced hypercapnia period of 120 s and the remaining time the rest period.

## 2.2.

### Instrumentation

The device used to perform the measurements was the dynamic near-infrared optical tomography (DYNOT) apparatus manufactured by NIRx Medical Technologies, (Los Angeles), which operates in continuous mode. The device can acquire measurements in parallel at four wavelengths: 725, 760, 810, and 830 nm at a sampling frequency of $\sim 4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{Hz}$. A more detailed description of the scanning device can be found in Schmitz et al.^{24} The reconstruction algorithms were written in MATLAB® and image reconstruction was performed in MATLAB® R2007a running on a standard PC with a single core 3 GHz Intel Pentium microprocessor and 1 GB RAM.

## 2.3.

### Image Reconstruction Algorithm

## 2.3.1.

#### Forward problem

For a medium $\mathrm{\Omega}\subset {\mathrm{IR}}^{3}$ with boundary $\partial \mathrm{\Omega}$, characterized by absorption and scattering coefficients ${\mu}_{a}(r)$, ${\mu}_{s}(r)$, $r\in \mathrm{\Omega}$, solving the forward problem involves numerical simulation of photon fluence rate measurements ${\{y(j)\}}_{j=1,s}$ from $d$ detectors $y(j)=[{y}_{1}(j),\dots ,{y}_{d}(j)]$ on $\partial \mathrm{\Omega}$, for each point source ${q}_{j}$ of the array of sources $q=[{q}_{1},\dots ,{q}_{s}]$ distributed on $\partial \mathrm{\Omega}$.

The forward problem for the source ${q}_{j}$ is described by the following parameters-to-output mapping:

from the space of parameter functions $u(r)=[{\mu}_{a}(r),{\mu}_{s}^{\prime}(r)]$ to the space of measurements $y(j)$. The forward mapping is obtained by combining a light propagation model (the forward model) with a model of the measurement system.A model of light propagation through tissue that is commonly used in applications involving continuous wave DOT systems is the diffusion approximation of the equation of radiative transfer.^{2}

## (2)

$$-\nabla \xb7D(r)\nabla {\varphi}_{j}(r)+{\mu}_{a}(r){\varphi}_{j}(r)={q}_{j}(r),\phantom{\rule[-0.0ex]{2em}{0.0ex}}r\in \mathrm{\Omega},$$## (3)

$$D(\xi )\frac{\partial {\varphi}_{j}(\xi )}{\partial n}+\frac{1}{2A}{\varphi}_{j}(\xi )=0,\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}\xi \in \partial \mathrm{\Omega},$$^{25}

For any given source ${q}_{j}$, the variable measured by the detector located at ${\xi}_{i}\in \partial \mathrm{\Omega}$ is the outward flux ${\gamma}_{j}({\xi}_{i})$. The corresponding measurement equations are given by

## (4)

$${\gamma}_{j}(\xi )=-D(\xi )\overrightarrow{n}(\xi )\xb7\nabla {\varphi}_{j}(\xi ),\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\xi \in \mathrm{\Omega},\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}{y}_{i}(j)={\gamma}_{j}({\xi}_{i}),\phantom{\rule{0ex}{0ex}}j=1,\dots ,s;\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}i=1,\dots ,d.$$A numerical solution to the DE [Eq. (2)] can be obtained for arbitrary tissue geometries using a finite element.^{25}^{,}^{26} Solving the forward problem is computationally the most expensive step of the reconstruction algorithm for 3-D problems.

## 2.3.2.

#### Generation of the 3-D finite element mesh of the rat’s head

A 3-D finite element mesh that incorporates high-resolution anatomical prior information was derived from pixel images obtained by a 7-Tesla high field animal magnet (Bruker BioSpin, Billerica, Massachusetts). The honeycomb optode holder was filled with clear liquid to generate contrast in the images at the optodes location. A sample of the pixel images is given in Fig. 1(a), where the location of the optodes is clearly visible.

Each image was segmented into skin, skull, muscle, and brain using image processing techniques, and then all slices were stacked together to build the 3-D model of the rat’s head displayed in Fig. 1(b), which later was converted into a tetrahedral mesh consisting of 169,793 elements and 33,944 nodes [Fig. 1(c)]. These steps were accomplished using ScanIP&ScanFE® (Simpleware, LTD, Exeter, UK). The element connectivity matrix for each part and the nodes coordinates were imported into MATLAB®. Each compartment was assigned specific absorption and scattering values: ${\mu}_{a}=0.02\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for skin, ${\mu}_{a}=0.005\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for skull, ${\mu}_{a}=0.015\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for brain, and ${\mu}_{a}=0.022\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for muscle; similarly, ${\mu}_{s}^{\prime}=0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for skin, ${\mu}_{s}^{\prime}=1.63\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for skull, ${\mu}_{s}^{\prime}=1.63\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for brain, and ${\mu}_{s}^{\prime}=1\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ for muscle.^{27}

Spatial information can also be incorporated to guide the inverse solver by constraining the solution to lie within the brain and not in the overlapping tissues (skull, muscle, and scalp). Boas et al.^{5} reconstructed absorption changes only for nodes lying within the brain, and for the remaining nodes, the absorption change was set to zero. Bluestone et al.^{19} constrained the reconstruction algorithm by averaging the gradient at every node within the scalp, muscle, and skull but varying independently the absorption values for the brain. The former procedure was selected given that it can be easily implemented and also produces smaller Jacobian matrices. The region of interest (ROI) defined comprised the area of the brain directly under the location of the optode holder, starting from the cortex up to 7 mm inside the brain. The ROI is displayed in Fig. 1(d) and consists of 7366 tetrahedral elements and 1882 nodes.

## 2.3.3.

#### Reduced-order forward model

The approach proposed to speed up the reconstruction process involves constructing a sparse approximation of the nonlinear mapping [Eq. (1)] using a data set generated by numerical simulation of the forward model given in Eqs. (2) to (4), using a finite element approximation on a fine mesh. The FEM-based forward model is called the complete model.

The reduced-order model (ROM) of Eq. (1) can be expressed in its component form as

## (5)

$${\widehat{y}}_{i,j}(t)={f}_{i,j}[u(t)]+{e}_{i,j}(t),\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}j=1,\dots ,s;\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}i=1,\dots ,{M}_{j};\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}i\ne j,$$To generate the data needed to estimate and validate the reduced forward model, ${N}_{D}={N}_{E}+{N}_{V}$ (typically ${N}_{D}=\phantom{\rule{0ex}{0ex}}1000$, ${N}_{E}={N}_{V}$) uniformly distributed random absorption values are used to compute detector measurements using the FEM approximation of the forward model in Eqs. (2) to (4). The random absorption coefficients were limited to the range $-50$ to 200% the value of the background absorption coefficient. ${N}_{E}$ data points are used to estimate the model, and the remaining ${N}_{V}$ data points are used to validate the model. For each input $\mathbf{u}(t)=[{\mu}_{a}({r}_{1},t),\dots ,{\mu}_{a}({r}_{K},t)]$, the forward model is simulated to generate the corresponding measurements that are used to estimate the model.

The model [Eq. (5)] is assumed to be of the form

## (6)

$${\widehat{y}}_{i,j}(t)=\sum _{l=1}^{{N}_{i,j}}{\theta}_{l}^{i,j}{p}_{l}^{i,j}[\mathbf{u}(t)]+{e}_{i,j},$$## (7)

$${p}_{k}^{i,j}(\mathbf{u})=\phi (\mathbf{u},{\mathbf{c}}_{k})={r}_{k}^{2}\text{\hspace{0.17em}}\mathrm{log}({r}_{k}),$$^{28}provide approximation rates nearly as good as the best-in-class multiquadric methods,

^{29}and have the advantage of not requiring the optimization of an additional shape parameter.

^{30}

Given the input-output data set ${D}_{i,j}^{E}={\{\mathbf{u}(t),{y}_{i,j}(t)\}}_{t=1,{N}_{E}}$, the mapping in Eq. (5) is derived such that it provides a sparse approximation of the underlying nonlinear function using ${N}_{ij}$ terms (${N}_{i,j}\ll {N}_{D}$). Moreover, for every source-detector combination, only the absorption values corresponding to a small number of nodes in the mesh will affect the measurements. Therefore, in practice, the number of inputs in Eq. (5) will be much smaller than the number of nodes, ${n}_{K}^{i,j}\ll K$.

As a result, the overall combination of source-detector models provides a reduced-order representation compared to the full finite element approximation. In simulation studies we have performed, we found that finite-dimensional models derived from data in this manner are as accurate as the standard finite element approximation but requires far less parameters.

## 2.3.4.

#### Model estimation

For every source-detector pair $i,j$, we solve the nonparametric regression problem involving two steps.

1. Determine a minimal set of relevant input nodes (the RBF input layer) ${\mathbf{u}}_{i,j}={\{{u}_{l}^{i,j}\}}_{l\in {I}_{i,j}}$, where ${I}_{i,j}$ is a subset of the full set of nodes $\{1,\dots ,K\}$, with dimension ${n}_{K}^{i,j}\ll K$.

2. Given a dictionary of candidate RBF basis functions, determine a minimal subset of regressors ${\{{p}_{k}\}}_{k=1,{N}_{i,j}}$, ${N}_{i,j}\ll {N}_{D}$ and the associated parameter vector ${\{{\theta}_{k}\}}_{k=1,{N}_{i,j}}$, which minimize the prediction error over the data set ${D}_{i,j}$.

Step 1. The input layer of the RBF network was selected using

*a priori*information derived from photon measurement density functions (PMDF),^{31}which characterize the sensitivity of measurements to changes of the optical parameters inside the medium. For each source-detector pair, the PMDF function was computed based on the FEM approximation over a fine mesh,^{32}and only the nodes within a region where values of PMDF are above a given threshold (1 to 5%.) are used to form the input layer of the RBF approximation. For polynomial approximations,^{18}^{,}^{33}it can be shown that all the terms selected in the model by alternative algorithms correspond to nodes within the high-sensitivity area.^{10}The union of nodes selected for all source-detector combinations ${I}^{\mathrm{ROI}}=\phantom{\rule{0ex}{0ex}}\underset{i,j}{\cup}{I}_{i,j}$ defines a region where changes of optical properties can be reconstructed for the given source-detector positions. In most applications, prior information from a secondary imaging modality or human anatomy can be used to specify an ROI and the placement of optodes.^{5}Step 2. Given the selected input vector, a set of candidate RBF functions was constructed by considering each data point as the center of an RBF function, i.e., ${\mathbf{c}}_{k}={\mathbf{u}}_{i,j}(k)$, $k=1,\dots ,\phantom{\rule{0ex}{0ex}}{N}_{E}$. The minimal set of basis functions in Eq. (6) is determined using a greedy orthogonal regression algorithm

where ${\mathrm{err}}_{i}$ is the % contribution of the $i$’th selected term and ${C}_{d}$ is a predefined threshold that can be used to control the trade-off between precision and computational speed. In order to avoid overfitting, a standard cross-validation method was implemented to fine tune the number of terms (basis functions) in the models, based on a separate validation data set.^{34}from an initial set of candidate RBF functions. The algorithm is based on the forward regression principle, which involves adding RBF functions to the model one at the time. At each step, each candidate basis function that is not already in the model is evaluated for inclusion in the model. The candidate regressors are ranked according to their contribution to reducing the variance of the error as measured by the error reduction ratio (ERR) index.^{34}The selection process is stopped when

## 2.3.5.

#### Reduced-order forward model of the rat’s head

For each node within the ROI, 1000 uniform distributed random absorption values within the range ${\mu}_{a}\in [0.0075,0.0045]\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ were generated. The synthetic measurements due to these changes were obtained by solving the DE for each sample.

For each of the 132 source-detector pairs, a reduced-order RBF (thin-plate spline) model was estimated using the first 500 samples; the remaining data were used for validation. The root mean squared error (RMSE), calculated for each source-detector pair $i-j$ as

## 2.3.6.

#### Image reconstruction

The absorption coefficients were obtained using the modified version^{19} of the normalized difference method,^{35} which involves solving the following optimization problem:

## (8)

$${\mathbf{u}}_{\mathrm{ROI}}^{*}=\underset{{\mathbf{u}}_{\text{ROI}}}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}\text{\hspace{0.17em}}\sum _{j=1}^{s}\sum _{j=1}^{{M}_{j}}{\left[\frac{\frac{{y}_{i,j}(t)}{{y}_{i,j}^{0}}{\widehat{y}}_{i,j}(0)-{f}_{i,j}({\mathbf{u}}_{\text{ROI}})}{\frac{{y}_{i,j}(t)}{{y}_{i,j}^{0}}{\widehat{y}}_{i,j}(0)}\right]}^{2},$$^{36}using spatial constraints.

^{5}

## 2.4.

### Calculation of Hemoglobin Concentration

The absorption coefficient of the tissue is wavelength dependent. The DYNOT instrument allows the parallel acquisition of absorption changes at four wavelengths. In this study, the reconstruction of hemodynamic changes was performed using only two near-infrared wavelengths of 760 and 830 nm. At these wavelengths, the absorption coefficient at a given location is assumed to be a linear combination of local oxyhemoglobin and deoxyhemoglobin concentrations ($\mathrm{\Delta}[{\mathrm{HbO}}_{2}]$ and $\mathrm{\Delta}[\mathrm{HbR}]$).^{19}

## (9)

$$\mathrm{\Delta}{\mu}_{a}(\lambda )={\epsilon}_{{\mathrm{HbO}}_{2}}(\lambda )\mathrm{\Delta}[{\mathrm{HbO}}_{2}]+{\epsilon}_{\mathrm{HbR}}(\lambda )\mathrm{\Delta}[\mathrm{HbR}],$$^{5}

The concentration changes in oxyhemoglobin (${\mathrm{HbO}}_{2}$) and deoxyhemoglobin (HbR) were calculated by solving a linear system of equations, using absorption coefficients estimated from measurements at ${\lambda}_{1}=760$ and ${\lambda}_{2}=830\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, as follows:^{19}

## (10)

$$\mathrm{\Delta}[\mathrm{HbR}]=\frac{{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{1})\mathrm{\Delta}{\mu}_{a}({\lambda}_{2})-{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{2})\mathrm{\Delta}{\mu}_{a}({\lambda}_{1})}{{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{1}){\epsilon}_{\mathrm{HbR}}({\lambda}_{2})-{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{2}){\epsilon}_{\mathrm{HbR}}({\lambda}_{1})}\phantom{\rule{0ex}{0ex}}\mathrm{\Delta}[{\mathrm{HbO}}_{2}]=\frac{{\epsilon}_{\mathrm{HbR}}({\lambda}_{2})\mathrm{\Delta}{\mu}_{a}({\lambda}_{1})-{\epsilon}_{\mathrm{HbR}}({\lambda}_{1})\mathrm{\Delta}{\mu}_{a}({\lambda}_{2})}{{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{1}){\epsilon}_{\mathrm{HbR}}({\lambda}_{2})-{\epsilon}_{{\mathrm{HbO}}_{2}}({\lambda}_{2}){\epsilon}_{\mathrm{HbR}}({\lambda}_{1})}.$$## 3.

## Results

## 3.1.

### Algorithm Evaluation Using a Phantom

## 3.1.1.

#### Experiment setup

To validate the reconstruction approach, a phantom study was carried out. We used a cylindrical container made of antireflective plastic (diameter 80 mm, height 300 mm, wall thickness 0.10 mm) full of skimmed milk (0.1% fat). A perturbation was added [Fig. 2(a)] to the media in the form of another inner plastic cylinder (diameter 20 mm, height 300 mm) full of milk at a different fat concentration (0.81%). The optical properties of the background media^{37} at 810 nm are ${\mu}_{a}=0.0024\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. Similarly, for the perturbation, ${\mu}_{a}=0.024\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and ${\mu}_{s}^{\prime}=0.4\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$.

Previously, milk has been used experimentally as strong scattering media^{38}39.40.^{–}^{41} because its optical absorption and scattering parameters are close to those of most living tissue,^{42} where $0.001<{\mu}_{a}<0.01\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and $0.5<{\mu}_{s}^{\prime}<10\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$.

The experiment consisted of reconstructing the size and location of the perturbation using differential measurements, that is, measurements taken before and after the inner cylinder was introduced into the container. Sixteen colocated sources and detectors were equally distributed around the cylinder at a height of 150 mm [Fig. 2(a)] providing $240\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{source}/\phantom{\rule{0ex}{0ex}}\text{detector}$ channels.

## 3.1.2.

#### Estimation and validation of reduced-order model

The geometry was discretized using 1344 triangular elements and 715 nodes. An ROI was defined to provide *a priori* information to the reconstruction algorithm; this area is indicated with the dotted line in Fig. 2(a). To generate the data needed for model estimation and validation, 1000 uniformly distributed random absorption values for each node were used to compute detector measurements using the FEM approximation of the DE. The data spanned from ${\mu}_{a}=0.0020$ to $0.03\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. The first 500 records were used as estimation data and the remaining records as validation data.

A TPS-RBF was used for model representation. The first step is the determination of the input layer; for this purpose, the Jacobian was calculated for each source-detector pair, and those nodes lying within a 5% threshold of the Jacobian and also located inside the ROI were selected as the input nodes of the network [Fig. 2(b)].

The minimal set of basis functions was selected using the orthogonal regression algorithm described in Sec. 2.3.4. The evolution of the ERR and the RMSE relative to the number of model terms are shown in Figs. 3(a) and 3(b), respectively. To show the performance of the ROM approach in the calculation of outward flux for source-detector pair 3-7, measurement predictions calculated using both the DE and the ROM are displayed in the upper panel of Fig. 4(a).

## 3.1.3.

#### Performance on the reduced-order model in the reconstruction of absorption changes in highly scattering media

Four different TPS-RBF models with ${C}_{d}=0.05$, 0.1, 0.2, and 0.3 were tested. Based on the image correlation coefficient (ICC), no significant variations of image quality were found, and therefore, the model with ${C}_{d}=0.3$ was selected for achieving faster image reconstructions than the other models (less terms are evaluated). Considering that the location of the perturbation for this experiment is known, it is possible to assess the qualitative spatial accuracy of reconstructed images. The ICC is given by^{35}

## (11)

$$\mathrm{ICC}(A,B)=\frac{1}{N-1}\frac{\sum _{i}({x}_{A}^{i}-{\overline{x}}_{A}^{i})({x}_{B}^{i}-{\overline{x}}_{B}^{i})}{\sqrt{\sum _{i}{({x}_{A}^{i}-{\overline{x}}_{A}^{i})}^{2}}\sqrt{\sum _{i}{({x}_{B}^{i}-{\overline{x}}_{B}^{i})}^{2}}},$$Examples of the recovered absorption change obtained using the complete and the $\mathrm{RBF}({C}_{d}=0.3)$ models are shown in Fig. 5; the exact location of the perturbation is indicated with a discontinuous line. For this reconstruction, the complete model required 100 iterations using the conjugate gradient descent (CGD) method and took $\sim 80\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$; the ROM required only six iterations and took $\sim 2.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$.

## 3.2.

### Simulated 3-D Reconstruction Using a Rat’s Head Model

The proposed approach was validated further through numerical simulation studies using the model of a rat’s head described in Sec. 2.3.2.

The simulation studies involved 12 optodes arranged in a honeycomb pattern, which are located at the top of the head as shown in Fig. 1(b). This results in 132 source-detector combinations. The optode configuration is the same as that used in real experiments.

A perturbation was specified in the form of an inclusion embedded in the brain [Fig. 6(a)] with time-varying optical properties. In the first simulation study, the absorption values corresponding to nodes lying inside the inclusion were obtained by sampling a quasiperiodic function.^{35} The full 3-D FEM model given by Eqs. (2) to (4) was simulated numerically to generate synthetic optode measurements that were used to reconstruct the dynamics of the absorption parameters.

The 3-D reconstructions of the inclusion for the first time point using the FEM and the ROM (Sec. 2.3.5) are displayed in Figs. 6(b) and 6(c), respectively. The isosurfaces correspond to absorption thresholds of 10 and 30% of their maximum estimated variation, $\mathrm{max}{(\mathrm{\Delta}{\mu}_{a})}_{\mathrm{FEM}}=0.0013\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$ and $\mathrm{max}{(\mathrm{\Delta}{\mu}_{a})}_{\mathrm{ROM}}=\phantom{\rule{0ex}{0ex}}0.0025\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{-1}$. To evaluate the quality of the reconstructions, ICC was calculated in each iteration [Fig. 7(a)]. Both algorithms achieved similar quality, but the reduction in time is extremely significant, as speed ratio is approximately three orders of magnitude. Specifically, the ROM approach achieves the maximum ICC value only after four iterations and $\sim 2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, while the FEM-based method achieved similar quality after five iterations, but at a considerably larger amount of time ($\sim 1\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{h}$).

The FEM approach achieves the best ICC score (a 3% improvement over the ROM approach) after 20 iterations and more than 4 h (ROM speed-up factor is $\sim 1200$). In both cases, the image reconstruction times for the FEM approach are prohibitively large, making DOT impractical in cases that require rapid clinical diagnosis. In practice, accuracy can be improved by using more accurate ROMs.

A point at the centre of the inclusion was selected to show the reconstruction of the dynamic changes using both methods. The estimated changes of optical absorption at this point displayed in Fig. 7(b) along with the original perturbation signal show that the ROM approach provides a better estimate of the dynamic changes in the absorption coefficient. To demonstrate the consistency of image reconstruction and speed-up, in Fig. 8 we show the ICC and reconstruction times for both approaches for the entire simulation.

## 3.3.

### Simulated 3-D Reconstruction of Hemodynamic Response to Hypercapnia in the Superior Sagittal Sinus of a Rat’s Head

In this simulation study, a perturbation was specified in the form of an inclusion located in the sagittal sinus of the rat cortex, as shown in Figs. 9(a) and 9(b). The shape of the inclusion resembles a section of a blood vessel. We generated changes of absorption coefficients for mesh nodes located inside the inclusion using profiles of hemodynamic responses to hypercapnia [Fig. 9(c)] in the sagittal sinus recorded in previous experimental studies^{43} using fMRI and optical imaging spectroscopy. The experiments followed the same protocol as detailed in Sec. 2.1.

Absorption changes were computed by substituting the concentration profiles for oxyhemoglobin and deoxyhemoglobin, shown in Fig. 9(c), in Eq. (9), using extinction coefficients at wavelengths of 760 and 830 nm, for each of the two chromophores. At the baseline, we assumed that the total hemoglobin concentration was 100 *μ*M and the hemoglobin oxygen saturation was 50%. The finite element approximation of the diffusion model [Eqs. (2) and (4)] and the measurement equation [Eq. (4)], computed for the 3-D mesh, was simulated to generate synthetic measurements of the outward flux at each optode location for each time point. The synthetic measurements were subsequently used in reconstruction. Absorption changes were reconstructed using the complete and approximated models. The number of iterations was limited to 10. The ROM was selected using ${C}_{d}=1$ as the stopping criteria. The FEM-based approach took 80 min to run 10 iterations for each time point. In contrast, the ROM approach took $\sim 6\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ to compute absorption changes for each time point.

Once the absorption changes were reconstructed at each wavelength, hemoglobin concentrations were derived as before using Eq. (10). Figure 10(c) shows coronal sections of changes in total hemoglobin $\mathrm{\Delta}[\mathrm{HbT}](x,y,z,t)$ for $y=-4$, $-1$, and 2 mm at $t=120\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, reconstructed using the ROM approach.

The estimated changes in oxyhemoglobin (${\mathrm{HbO}}_{2}$), deoxyhemoglobin (HbR), and total hemoglobin (HbT) for the nodes lying within the inclusion were averaged and the values are displayed in Fig. 10(a). Clearly, the ROM approach provides a much better estimate of the magnitude of the hemodynamic changes. To further improve the accuracy, we derived linear calibration functions for the FEM- and ROM-based reconstruction algorithms, which later were used to calibrate the time series of hemodynamic changes in the rat barrel cortex reconstructed using experimental data, following the hypercapnia challenge detailed in Sec. 2.1.

The calibration functions are as follows:

The calibrated changes of ${\mathrm{HbO}}_{2}$, HbR, and HbT reconstructed using the FEM and ROM approaches are shown in Fig. 10(b). The reconstruction errors for each method are shown in Fig. 10(c). The errors can be reduced further by fitting higher-order polynomials, but this may result in overfitting the data at a particular location. We found that the linear calibration functions given in Eqs. (12) and (13) could be used to calibrate hemodynamic changes reconstructed at other locations.

## 3.4.

### Experimental 3-D Reconstruction of Hemodynamic Response to Hypercapnia in the Barrel Cortex of a Rat’s Head

Measurements were obtained using the DYNOT instrument at sampling frequency of $\sim 3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{Hz}$. Each sample for wavelengths at 760 and 830 nm was processed with both the FEM- and the ROM-based image reconstruction algorithms, and a full tomography reconstruction of absorption changes was obtained. Reconstruction of absorption changes was achieved after 5 to 10 iterations using either the ROM or FEM approaches. However, our approach took on average $\sim 3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, while the latter took between 1 and 2 h per time point. Absorption changes were later converted into hemoglobin changes using Eq. (10).

To validate the results of our experimental study, we compared the reconstructed hemodynamic changes during hypercapnia with MRI measurements of changes (average across seven animals) in the cerebral blood volume (CBV-fMRI), obtained by using a paramagnetic contrast agent (MION),^{43} following identical hypercapnic perturbations (onset time of the 10% step increase of ${\mathrm{CO}}_{2}$ was at ${t}_{\text{on}}={t}_{0}+60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, ${t}_{0}=-60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, with the offset at ${t}_{\text{off}}={t}_{\text{on}}+120\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$).

Figure 11(a) shows a representative coronal section CBV-MRI statistical parameter map superimposed on a $256\times 256$ structural scan. A coronal section of the change in total hemoglobin $\mathrm{\Delta}[\mathrm{HbT}](x,y,z,t)$ for $y=2$, obtained using the ROM approach at $t=120\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$, is shown coregistered with the anatomical MRI map of the same in rat [Fig. 11(b)]. This coronal slice is centered to the area directly under the optode holder. Coregistration was easily achieved given that the finite element mesh was derived directly from the MRI data; this minimized the influence of systematic errors due to geometry mismatch and optode positioning, which affect reconstructed images.^{44} To visualize the time evolution of hemodynamic changes during the experiment, we positioned the two ROIs [two rectangles in Fig. 11(b)] in the rat barrel cortex to correspond to a superficial (0 to 1 mm) region and a deep (1 to 2 mm) cortical region, which were used in Ref. 43 to measure changes in blood volume following hypercapnia and whisker stimulation.

The changes in oxyhemoglobin (${\mathrm{HbO}}_{2}$) and deoxyhemoglobin (HbR) for the nodes lying within the two ROIs were averaged, and the final hemodynamic responses [Fig. 11(c)] were inferred using the calibration functions given in Eqs. (12) and (13). The changes in HbT for the deep and superficial cortical layers, reconstructed using the ROM and FEM approaches, are shown in Fig. 11(d) superimposed on the CBV-weighted fMRI reconstructions obtained for the corresponding superficial and deep regions.

We found that the CBV response in the superficial cortical layer ${\mathrm{CBV}}_{\mathrm{sup}}$ is commensurate with the fractional changes in HbT (superficial) reconstructed from DOT measurements using ROM. The time series obtained using the ROM approach show very similar onset gradients and peak values (26 *μ*M versus 24 *μ*M), but ${\mathrm{CBV}}_{\mathrm{sup}}$ starts to increase earlier and returns faster to the baseline. It is very likely that this reflects interanimal variability as the DOT and CBV measurements were not measured concurrently in the same animal. The CBV response in the deep cortical layer ${\mathrm{CBV}}_{\text{deep}}$ peaks at 33.8 *μ*M, while HbT (deep) time series reconstructed from DOT measurements using the ROM approach peak at $\sim 30.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{M}$. This reflects the reduced depth sensitivity of DOT compared with fMRI.

Figure 12 displays the evolution of $\mathrm{\Delta}[{\mathrm{HbO}}_{2}](x,y,z,t)$, $\mathrm{\Delta}[\mathrm{HbR}](x,y,z,t)$, and $\mathrm{\Delta}[\mathrm{HbT}](x,y,z,t)$ for $y=2$ at $t=\phantom{\rule{0ex}{0ex}}\{\mathrm{0,60,120,180}\}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$. The reconstruction at $t=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ shows minimal activity during the baseline period. Reconstructions at $t=60$ and 120 s show an increase in oxyhemoglobin and a decrease in deoxyhemoglobin, which is the typical response to an increase of ${\mathrm{CO}}_{2}$. The reconstructed spatial maps agree with similar DOT and fMRI studies.^{43}^{,}^{44}

## 4.

## Discussion

The results demonstrate that our fast image reconstruction approach for DOT, which uses highly optimized maps between the space of optical parameters and the space of measurements, generates 3-D images in a fraction of the time taken by standard FEM-based algorithms, while preserving image quality.

In this work, we used RBFs^{17} to construct our reduced-order forward light models. We found this modeling approach to be more powerful and better suited for this modeling task compared to the polynomial approximation technique.^{18} The models are developed based on numerical simulation data generated by an FEM approximation of the diffusion equation over a fine mesh. A greedy algorithm is used to select a minimal number of RBF basis functions from a library of candidate regressors. An important advantage of RBF over polynomial approximations is that it produces a much smaller set of candidate model terms that have to be evaluated.

The number of candidate terms for a polynomial model of order $q$ with $n$ inputs, given by

The RBF models can be constructed to approximate with desired accuracy the predictions of the full forward model, which was used to generate training data. In practice, however, reducing the error below a certain threshold rapidly increases model complexity, which in turn increase reconstruction times, while the quality of the reconstructed image does not improve significantly. For a given geometry and ROI, simulation studies can be carried out to determine the ROM, which offers the optimal trade-off between reconstruction accuracy and speed. Even when more complex forward RBF light models are used, however, the speed-up is still considerable.

If a high-fidelity subject-specific head model is available,^{45} our approach can be used to derive much simpler individualized forward light models covering the relevant field of view, which can be used to perform fast reconstructions.

This study demonstrated the applicability of the proposed DOT reconstruction approach for real-time monitoring of brain hemodynamics through a hypercapnia experiment. The results agree with the expected physiological changes and with the results of a similar experimental study carried out by Bluestone et al.^{19} using the model-based iterative image reconstruction,^{46} which reported typical reconstruction times of $\sim 2$ to 4 h per time point. To circumvent the computational bottleneck, the authors of the previous study used measurements corresponding to a single source-detector pair to calculate the time course of ${\mathrm{HbO}}_{2}$, HbR, and HbT using the Beer-Lambert law and only performed full tomographic reconstruction at specific time points of interest. In contrast, in our study, we performed 3-D reconstructions at each time point and have shown that following the hypercapnic challenge, the reconstructed changes in total hemoglobin, which are believed to be approximately the CBV changes in the barrel cortex, are in good agreement with blood volume changes measured in a previous study^{43} using the CBV-MRI technique.^{47}^{,}^{48}

Even on a modestly powered desktop PC and using code written in MATLAB®, which is significantly slower than compiled code, it took $\sim 3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{s}$ to compute a 3-D reconstruction at each time point. Using compiled C/C++ code, the time to compute a reconstruction could be reduced below the acquisition time of the instrument. The low computational cost of the method means that DOT reconstruction algorithms can be run on a mobile device. Given the lower cost and increased portability of DOT instrumentation, this opens up the possibility to use DOT routinely for continuous monitoring brain function of newborns^{49} as well as adults^{50}51.^{–}^{52} in real time and at much higher resolution^{53} than what is achievable using the standard algorithms. In a clinical context, in the absence of subject-specific anatomical images, atlas-based head models registered on the subject’s head, such as those proposed by Ferradal,^{54} could be used to obtain ROMs that are suitable for real-time processing.

## Acknowledgments

The authors gratefully acknowledge that this work was supported in part by EPSRC grants EP/H00453X/1 (S.A.B. and D.C.) and GR/T01006/01 (D.C.), by BBSRC grant BB/K010123/1 (Y.Z., S.A.B. and D.C.) and by an ERC Advanced Investigator Grant (S.A.B.). E.E.V.R. gratefully acknowledges the support from a grant of the Mexican National Research Council for Science and Technology (CONACYT).

## References

## Biography

**Ernesto E. Vidal-Rosas** is a research associate at the University of Sheffield. He received his BEng and MSc degrees in electronic engineering from the Orizaba Institute of Technology and the National Centre for Research and Technology Development Cuernavaca in 2003 and 2006, respectively, and his PhD degree from the University of Sheffield in 2011. His main research interests include the development and application of novel methods for image reconstruction and analysis for diffuse optical tomography.

**Stephen A. Billings** is a professor in the Department of Automatic Control and Systems Engineering at the University of Sheffield and leads the Signal Processing and Complex Systems research group. His research interests include system identification and information processing for nonlinear systems, narmax methods, model validation, prediction, spectral analysis, adaptive systems, nonlinear systems analysis, wavelets, machine vision, cellular automata, spatiotemporal systems, EEG, fMRI and optical imagery of the brain, synthetic biology, crystal growth, and related fields.

**Ying Zheng** is a professor at the University of Reading. She has been working in the interface of systems engineering and neuroscience for over 10 years. A graduate of the University of Sheffield, she joined the University of Reading in August 2013. In close collaboration with researchers at Sheffield, she has developed mathematical models of neural and hemodynamic signals and their coupling to elucidate on the functions of the brain in normal and diseased states.

**Aneurin J. Kennerley** is MR physicist at the University of Sheffield. He received his MPhys degree in experimental physics from the University of Newcastle upon Tyne in 2002 and his PhD degree in neuroscience from the University of Sheffield in 2006. His research concerns the biophysical modeling of the BOLD fMRI signal using concurrent fMRI and optical imaging spectroscopy. His research won the AstraZeneca prize (2010) and Peter Mansfield prize (2013) for innovative *in vivo* MRI research.

**Daniel Coca** is a professor at the University of Sheffield. He received his MEng degree in electrical engineering from Transilvania University of Brasov in 1993 and his PhD degree in automatic control and systems engineering from the University of Sheffield in 1997. His research interests include modeling, analysis, and control of complex systems, inverse problems, diffuse optical tomography, and biological data analysis using reconfigurable computers.