Reduced-order modeling of light transport in tissue for real-time monitoring of brain hemodynamics using diffuse optical tomography

Ernesto E. Vidal-Rosas; Stephen Billings; Ying Zheng; John E. Mayhew; David K. Johnston; Aneurin J. Kennerley; Daniel Coca

doi:10.1117/1.JBO.19.2.026008

13 February 2014 Reduced-order modeling of light transport in tissue for real-time monitoring of brain hemodynamics using diffuse optical tomography

Ernesto E. Vidal-Rosas, Stephen Billings, Ying Zheng, John E. Mayhew, David K. Johnston, Aneurin J. Kennerley, Daniel Coca

Author Affiliations +

Journal of Biomedical Optics, Vol. 19, Issue 2, 026008 (February 2014). https://doi.org/10.1117/1.JBO.19.2.026008

Abstract

This paper proposes a new reconstruction method for diffuse optical tomography using reduced-order models of light transport in tissue. The models, which directly map optical tissue parameters to optical flux measurements at the detector locations, are derived based on data generated by numerical simulation of a reference model. The reconstruction algorithm based on the reduced-order models is a few orders of magnitude faster than the one based on a finite element approximation on a fine mesh incorporating 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.¹^,² 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.³

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.¹ 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.⁴^,⁵ 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.⁶ 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⁷ 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.⁸ 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,⁹ which produce sparse stiffness matrices, can reportedly reduce the time required to solve the forward problem ( $\sim 36 times$ faster for a 3-D geometry). To speed up parameter estimation, Eames et al.¹⁰ 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.¹¹^,¹² However, this choice may lead to erroneous results if the models are not coregistered.⁴ Other strategies employed to reduce the complexity of the reconstruction algorithms include the use of nonlinear multigrid optimization¹³^,¹⁴ and adaptive mesh techniques.¹³^,¹⁴ Solving the inverse problem using analytic methods can dramatically reduce the reconstruction time,¹⁵ but the application of these methods is limited because analytic solutions are only available for some cases, such as the slab or the cylinder.¹⁶

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,¹⁷ 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.¹⁸ 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,¹⁹ tissue oxygenation,²⁰ and hemodynamics.²¹^–²³ This experiment focuses on the hemodynamic response under hypercapnia.

The animal used was a female hooded Lister rat weighing $\sim 400 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 g / 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 mg / 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% $O_{2}$ , 80% $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 Hz$ . A more detailed description of the scanning device can be found in Schmitz et al.²⁴ 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 $Ω \subset {IR}^{3}$ with boundary $\partial Ω$ , characterized by absorption and scattering coefficients $μ_{a} (r)$ , $μ_{s} (r)$ , $r \in Ω$ , 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 Ω$ , for each point source $q_{j}$ of the array of sources $q = [q_{1}, \dots, q_{s}]$ distributed on $\partial Ω$ .

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

Eq. (1)

y (j) = P_{j} u (r), j = 1, \dots, s,

from the space of parameter functions

u (r) = [μ_{a} (r), μ_{s}^{'} (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.²

Eq. (2)

- \nabla \cdot D (r) \nabla ϕ_{j} (r) + μ_{a} (r) ϕ_{j} (r) = q_{j} (r), r \in Ω,

where

ϕ_{j} (r)

is the spatially varying diffuse photon density at position

r

given the source

q_{j}

,

μ_{a}

is the absorption coefficient,

D = {[3 (μ_{a} + μ_{s}^{'})]}^{- 1}

is the diffusion coefficient, and

μ_{s}^{'}

is the reduced scattering coefficient. The collimated source incident at

ξ_{j} \in \partial Ω

is usually represented by an isotropic point source

q_{j} (r) = δ (r - r_{j})

, where

r_{j}

is located at a depth of one scattering length inside the medium, along the direction of the normal vector to the surface at the source location

\vec{n} (ξ_{j})

. The boundary condition usually employed is of the Robin type.

Eq. (3)

D (ξ) \frac{\partial ϕ_{j} (ξ)}{\partial n} + \frac{1}{2 A} ϕ_{j} (ξ) = 0, ξ \in \partial Ω,

where the term

A

accounts for the refractive index boundary mismatch between the interior and exterior medium.²⁵

For any given source $q_{j}$ , the variable measured by the detector located at $ξ_{i} \in \partial Ω$ is the outward flux $γ_{j} (ξ_{i})$ . The corresponding measurement equations are given by

Eq. (4)

γ_{j} (ξ) = - D (ξ) \vec{n} (ξ) \cdot \nabla ϕ_{j} (ξ), ξ \in Ω, y_{i} (j) = γ_{j} (ξ_{i}), j = 1, \dots, s; i = 1, \dots, d .

A numerical solution to the DE [Eq. (2)] can be obtained for arbitrary tissue geometries using a finite element.²⁵^,²⁶ 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.

Fig. 1

Main steps for generating the three-dimensional (3-D) finite element mesh from magnetic resonance imaging (MRI) slices. (a) MRI slice obtained with a 7T MRI scanner. (b) 3-D geometric model of the rat’s head and optodes. (c) Finite element mesh used solved the forward problem. (d) Skull, brain, and region of interest (ROI) used to constrain the inverse solver.

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: $μ_{a} = 0.02 {mm}^{- 1}$ for skin, $μ_{a} = 0.005 {mm}^{- 1}$ for skull, $μ_{a} = 0.015 {mm}^{- 1}$ for brain, and $μ_{a} = 0.022 {mm}^{- 1}$ for muscle; similarly, $μ_{s}^{'} = 0.5 {mm}^{- 1}$ for skin, $μ_{s}^{'} = 1.63 {mm}^{- 1}$ for skull, $μ_{s}^{'} = 1.63 {mm}^{- 1}$ for brain, and $μ_{s}^{'} = 1 {mm}^{- 1}$ for muscle.²⁷

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.⁵ 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.¹⁹ 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

Eq. (5)

{\hat{y}}_{i, j} (t) = f_{i, j} [u (t)] + e_{i, j} (t), j = 1, \dots, s; i = 1, \dots, M_{j}; i \neq j,

where

{\hat{y}}_{i, j} (t)

is the predicted steady-state measurement at the

i

’th optode location at a given time

t

computed using the forward model [Eqs. (2) to (4)] given the source

q_{j}

.

f_{i, j}

is a nonlinear function for source-detector pair

i - j

,

u (t) = [u_{1} (t), \dots, u_{K} (t)]

,

u_{k} (t) = μ_{a} (r_{k}, t)

is the absorption value for the

k

’th node at a given time

t

,

K

is the total number of nodes, and

e_{i, j}

is the approximation error. As the focus is on the estimation of the absorption coefficient, the model [Eq. (5)] assumes that the reduced scattering coefficient of the medium is spatially invariant.

To generate the data needed to estimate and validate the reduced forward model, $N_{D} = N_{E} + N_{V}$ (typically $N_{D} = 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 $u (t) = [μ_{a} (r_{1}, t), \dots, μ_{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

Eq. (6)

{\hat{y}}_{i, j} (t) = \sum_{l = 1}^{N_{i, j}} θ_{l}^{i, j} p_{l}^{i, j} [u (t)] + e_{i, j},

where

p_{k}^{i, j} (x)

is an RBF. The RBF chosen in this study is the thin-plate spline.

Eq. (7)

p_{k}^{i, j} (u) = φ (u, c_{k}) = r_{k}^{2} \log (r_{k}),

where

c_{k}

denotes the center of the RBF and

r_{k} = {‖ u - c_{k} ‖}_{2}

is the Euclidian norm. Thin-plate splines are conditionally positive kernels that are invariant under scaling,²⁸ provide approximation rates nearly as good as the best-in-class multiquadric methods,²⁹ and have the advantage of not requiring the optimization of an additional shape parameter.³⁰

Given the input-output data set $D_{i, j}^{E} = {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_{i j}$ terms ( $N_{i, j} ≪ 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} ≪ 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) $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} ≪ 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} ≪ N_{D}$ and the associated parameter vector ${θ_{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),³¹ 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,³² 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,¹⁸^,³³ it can be shown that all the terms selected in the model by alternative algorithms correspond to nodes within the high-sensitivity area.¹⁰ The union of nodes selected for all source-detector combinations $I^{ROI} = \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.⁵
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., $c_{k} = u_{i, j} (k)$ , $k = 1, \dots, N_{E}$ . The minimal set of basis functions in Eq. (6) is determined using a greedy orthogonal regression algorithm³⁴ 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.³⁴ The selection process is stopped when
$100 - \sum_{i = 1}^{n} {err}_{i} \leq C_{d},$
where ${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.

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 $μ_{a} \in [0.0075, 0.0045] {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

{RMSE}_{i j} = \sqrt{\frac{\sum_{t = 1}^{N} {[y_{i, j} (t) - {\hat{y}}_{i, j} (t)]}^{2}}{N}}

was used to evaluate the performance of each model. In order to investigate the trade-off between accuracy and speed on the calculation of outward flux in comparison to the measurements computed using the DE, four sets of models were estimated for different

C_{d}

values. The thin-plane-spline (TPS)-RBF models used in reconstruction (estimated for

C_{d} = 0.3

) had between 8 and 399 inputs and between 11 and 92 RBFs.

2.3.6.

Image reconstruction

The absorption coefficients were obtained using the modified version¹⁹ of the normalized difference method,³⁵ which involves solving the following optimization problem:

Eq. (8)

u_{ROI}^{*} = \underset{u_{ROI}}{\arg \min} \sum_{j = 1}^{s} \sum_{j = 1}^{M_{j}} {[\frac{\frac{y_{i, j} (t)}{y_{i, j}^{0}} {\hat{y}}_{i, j} (0) - f_{i, j} (u_{ROI})}{\frac{y_{i, j} (t)}{y_{i, j}^{0}} {\hat{y}}_{i, j} (0)}]}^{2},

where

y_{i, j} (t)

is a measurement taken at time

t

,

y_{i, j}^{0}

is the predefined reference state,

{\hat{y}}_{i, j} (0) = f_{i, j} (u_{ROI}^{0})

is the predicted measurement corresponding to the reference medium, and

u_{ROI}^{0}

is the initial vector of the optical parameters specified in Sec. 2.3.2. The optimization problem was solved using an iterative, conjugate-gradient algorithm³⁶ using spatial constraints.⁵

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 ( $Δ [{HbO}_{2}]$ and $Δ [HbR]$ ).¹⁹

Eq. (9)

Δ μ_{a} (λ) = ε_{{HbO}_{2}} (λ) Δ [{HbO}_{2}] + ε_{HbR} (λ) Δ [HbR],

where

ε_{{HbO}_{2}} (λ)

and

ε_{HbR} (λ)

are the extinction coefficients for oxyhemoglobin and deoxyhemoglobin, respectively, and

λ

is the photon wavelength.⁵

The concentration changes in oxyhemoglobin ( ${HbO}_{2}$ ) and deoxyhemoglobin (HbR) were calculated by solving a linear system of equations, using absorption coefficients estimated from measurements at $λ_{1} = 760$ and $λ_{2} = 830 nm$ , as follows:¹⁹

Eq. (10)

Δ [HbR] = \frac{ε_{{HbO}_{2}} (λ_{1}) Δ μ_{a} (λ_{2}) - ε_{{HbO}_{2}} (λ_{2}) Δ μ_{a} (λ_{1})}{ε_{{HbO}_{2}} (λ_{1}) ε_{HbR} (λ_{2}) - ε_{{HbO}_{2}} (λ_{2}) ε_{HbR} (λ_{1})} Δ [{HbO}_{2}] = \frac{ε_{HbR} (λ_{2}) Δ μ_{a} (λ_{1}) - ε_{HbR} (λ_{1}) Δ μ_{a} (λ_{2})}{ε_{{HbO}_{2}} (λ_{1}) ε_{HbR} (λ_{2}) - ε_{{HbO}_{2}} (λ_{2}) ε_{HbR} (λ_{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³⁷ at 810 nm are $μ_{a} = 0.0024 {mm}^{- 1}$ and $μ_{s}^{'} = 0.5 {mm}^{- 1}$ . Similarly, for the perturbation, $μ_{a} = 0.024 {mm}^{- 1}$ and $μ_{s}^{'} = 0.4 {mm}^{- 1}$ .

Fig. 2

Phantom study specifications. (a) Phantom geometry, dimensions, and optode locations. The dotted line indicates the ROI used. (b) Input nodes for the thin-plane-spline radial basis function (TPS-RBF) model corresponding to the source-detector pair 3-7. The selected inputs represent all the nodes inside the ROI for which the Jacobian exceeds the 5% boundary.

Previously, milk has been used experimentally as strong scattering media³⁸^–⁴¹ because its optical absorption and scattering parameters are close to those of most living tissue,⁴² where $0.001 < μ_{a} < 0.01 {mm}^{- 1}$ and $0.5 < μ_{s}^{'} < 10 {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 source / 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 $μ_{a} = 0.0020$ to $0.03 {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).

Fig. 3

Model selection for the phantom study. (a) Unexplained variance as a function of the number of terms of the TPS-RBF model for source-detector pair 3-7. (b) Root mean squared error, as a function of the number of terms, evaluated for the estimation and validation data.

Fig. 4

Comparison of the finite element method (FEM) and reduced-order model (ROM) approaches used in the phantom study. (a) Comparison of the outward flux for source-detector pair 3-7 calculated using both the diffusion equation and the ROM approach. (b) Error (%) between both models. (c) Evolution of image correlation coefficient (ICC) for 100 iterations using the FEM-based (red) and ROM (blue) approaches in the reconstruction algorithm.

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³⁵

Eq. (11)

ICC (A, B) = \frac{1}{N - 1} \frac{\sum_{i} (x_{A}^{i} - {\bar{x}}_{A}^{i}) (x_{B}^{i} - {\bar{x}}_{B}^{i})}{\sqrt{\sum_{i} {(x_{A}^{i} - {\bar{x}}_{A}^{i})}^{2}} \sqrt{\sum_{i} {(x_{B}^{i} - {\bar{x}}_{B}^{i})}^{2}}},

where

x_{A}^{i}

and

x_{B}^{i}

denote the intensities of the

i

’th pixel in images

A

and

B

, respectively, and

x_{A}^{i}

and

x_{B}^{i}

denote the mean intensities of the images. Image

A

corresponds to the original image, and image

B

corresponds to the reconstructed image using the complete or the reduced-order forward models. A perfect match between the original and reconstructed image is achieved when

ICC = 1

. The evolution of the ICC with respect to the number of iterations [Fig. 4(b)] shows that the reconstruction based on the ROM converges much faster.

Examples of the recovered absorption change obtained using the complete and the $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 s$ ; the ROM required only six iterations and took $\sim 2.5 s$ .

Fig. 5

Image reconstruction of the perturbation shown in Fig. 2(a) using (a) FEM-based and (b) ROM approaches.

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.³⁵ 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.

Fig. 6

Simulated 3-D reconstruction using the rat head mesh for a single time point. (a) Location of the simulated perturbation. Reconstructed perturbation using (b) FEM-based approach and (c) ROM approach. The iso-surfaces correspond to the 10% and 30% of the maximum recovered absorption value. The overlapping skin, skull, and muscle tissues are not displayed.

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, $\max {(Δ μ_{a})}_{FEM} = 0.0013 {mm}^{- 1}$ and $\max {(Δ μ_{a})}_{ROM} = 0.0025 {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 s$ , while the FEM-based method achieved similar quality after five iterations, but at a considerably larger amount of time ( $\sim 1 h$ ).

Fig. 7

Comparison of the FEM- and ROM-based reconstruction approaches in the simulated 3-D image reconstruction for a single time point. (a) Evolution of the ICC for the first time sample. (b) Reconstructed time series using both the FEM- and ROM-based methods.

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.

Fig. 8

Comparison of the FEM- and ROM-based reconstruction approaches in the simulated 3-D image reconstruction for all the time points. (a) Evolution of the ICC. The quality of the image decreases for small absorption changes ( $< 3 %$ ). (b) Reconstruction times for the FEM- and ROM-based approaches.

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⁴³ using fMRI and optical imaging spectroscopy. The experiments followed the same protocol as detailed in Sec. 2.1.

Fig. 9

Simulated hypercapnia using the rat head mesh. (a) and (b) Location and size of the simulated absorber. (c) Oxyhemoglobin ( ${HbO}_{2}$ ) and deoxyhemoglobin (HbR) profiles used to generate synthetic optical flux measurements.

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 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 $Δ [HbT] (x, y, z, t)$ for $y = - 4$ , $- 1$ , and 2 mm at $t = 120 s$ , reconstructed using the ROM approach.

Fig. 10

Reconstructed changes in ${HbO}_{2}$ and HbR in the rat brain during hypercapnia based on simulated measurements. (a) Coronal sections of the change in total hemoglobin $Δ [HbT] (x, y, z, t)$ for $y = - 4$ , $- 1$ , 2 mm at $t = 120 s$ . (b) Uncalibrated hemodynamic responses reconstructed using the FEM and ROM approaches superimposed on the original simulated changes. (c) Calibrated hemodynamic responses reconstructed using the FEM and ROM approaches. (d) Reconstruction errors corresponding to calibrated responses.

The estimated changes in oxyhemoglobin ( ${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:

• FEM approach:
Eq. (12)
$Δ {[{HbO}_{2}]}_{FEM, cal} = 5.26 Δ {[{HbO}_{2}]}_{FEM} - 0.94 Δ {[HbR]}_{FEM, cal} = 5.26 Δ {[HbR]}_{FEM} + 0.39 .$
• ROM approach:
Eq. (13)
$Δ {[{HbO}_{2}]}_{ROM, cal} = 0.93 Δ {[{HbO}_{2}]}_{ROM} - 2.81 Δ {[HbR]}_{ROM, cal} = 0.91 Δ {[HbR]}_{ROM} - 0.90 .$

The calibrated changes of ${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 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 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),⁴³ following identical hypercapnic perturbations (onset time of the 10% step increase of ${CO}_{2}$ was at $t_{on} = t_{0} + 60 s$ , $t_{0} = - 60 s$ , with the offset at $t_{off} = t_{on} + 120 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 $Δ [HbT] (x, y, z, t)$ for $y = 2$ , obtained using the ROM approach at $t = 120 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.⁴⁴ 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.

Fig. 11

Reconstructed changes in ${HbO}_{2}$ , HbR, and HbT in the rat brain during hypercapnia based on experimental diffuse optical tomography (DOT) measurements. (a) Spatial map of blood volume changes measured using CBV-MRI. The two rectangles mark the superficial and deep cortical areas used to compute blood volumes. (b) HbT in a coregistered image of MRI and DOT for $t = 120 s$ . The two rectangles mark the superficial and deep cortical areas used to compute blood volumes. (c) Reconstructed changes in ${HbO}_{2}$ and HbR inside the regions defined in (b), using FEM and ROM approaches. (c) Changes in HbT inside the regions defined in (b), reconstructed using FEM and ROM approaches, superimposed on blood volume changes (CBV) inside the regions defined in (a), reconstructed using CBV-MRI.

The changes in oxyhemoglobin ( ${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 ${CBV}_{\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 ${CBV}_{\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 ${CBV}_{deep}$ peaks at 33.8 μM, while HbT (deep) time series reconstructed from DOT measurements using the ROM approach peak at $\sim 30.3 μ M$ . This reflects the reduced depth sensitivity of DOT compared with fMRI.

Figure 12 displays the evolution of $Δ [{HbO}_{2}] (x, y, z, t)$ , $Δ [HbR] (x, y, z, t)$ , and $Δ [HbT] (x, y, z, t)$ for $y = 2$ at $t = {0,60,120,180} s$ . The reconstruction at $t = 0 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 ${CO}_{2}$ . The reconstructed spatial maps agree with similar DOT and fMRI studies.⁴³^,⁴⁴

Fig. 12

Spatial maps of reconstructed changes in ${HbO}_{2}$ and HbR in the rat brain based on experimental DOT measurements at different time points during hypercapnia.

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¹⁷ 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.¹⁸ 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

M = \sum_{i = 1}^{q} n_{i}, where n_{i} = \frac{n_{i - 1} (n + i - 1)}{i}, n_{0} = 1

increases rapidly as the model order and the number of inputs increase. For RBF approximations, the dimension of the candidate set depends mainly on the length of the estimation data, which is especially advantageous when modeling light transport in the 3-D case, where there are hundreds or even thousands of input nodes.

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,⁴⁵ 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.¹⁹ using the model-based iterative image reconstruction,⁴⁶ 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 ${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⁴³ using the CBV-MRI technique.⁴⁷^,⁴⁸

Even on a modestly powered desktop PC and using code written in MATLAB®, which is significantly slower than compiled code, it took $\sim 3 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⁴⁹ as well as adults⁵⁰^–⁵² in real time and at much higher resolution⁵³ 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,⁵⁴ 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

1.

S. R. ArridgeJ. C. Schotland, “Optical tomography: forward and inverse problems,” Inv. Prob., 25 (12), 123010 (2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar

2.

S. R. Arridge, “Optical tomography in medical imaging,” Inv. Prob., 15 (2), R41 –R93 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

3.

A. P. GibsonJ. C. HebdenS. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol., 50 (4), R1 –R43 (2005). http://dx.doi.org/10.1088/0031-9155/50/4/R01 PHMBA7 0031-9155 Google Scholar

4.

P. K. Yalavarthyet al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express, 15 (13), 8043 –8058 (2007). http://dx.doi.org/10.1364/OE.15.008043 OPEXFF 1094-4087 Google Scholar

5.

D. A. BoasA. M. DaleM. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage, 23 S275 –S288 (2004). http://dx.doi.org/10.1016/j.neuroimage.2004.07.011 NEIMEF 1053-8119 Google Scholar

6.

P. K. Yalavarthyet al., “Implementation of a computationally efficient least-squares algorithm for highly under-determined three-dimensional diffuse optical tomography problems,” Med. Phys., 35 (5), 1682 –1697 (2008). http://dx.doi.org/10.1118/1.2889778 MPHYA6 0094-2405 Google Scholar

7.

Y. ZhaiS. A. Cummer, “Fast tomographic reconstruction strategy for diffuse optical tomography,” Opt. Express, 17 (7), 5285 –5297 (2009). http://dx.doi.org/10.1364/OE.17.005285 OPEXFF 1094-4087 Google Scholar

8.

V. Kolehmainenet al., “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A, 26 (10), 2257 –2268 (2009). http://dx.doi.org/10.1364/JOSAA.26.002257 JOAOD6 0740-3232 Google Scholar

9.

A. Landragin-Frassatiet al., “Application of a wavelet-Galerkin method to the forward problem resolution in fluorescence diffuse optical tomography,” Opt. Express, 17 (21), 18433 –18448 (2009). http://dx.doi.org/10.1364/OE.17.018433 OPEXFF 1094-4087 Google Scholar

10.

M. E. Eameset al., “An efficient Jacobian reduction method for diffuse optical image reconstruction,” Opt. Express, 15 (24), 15908 –15919 (2007). http://dx.doi.org/10.1364/OE.15.015908 OPEXFF 1094-4087 Google Scholar

11.

P. K. Yalavarthyet al., “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys., 34 (6), 2085 –2098 (2007). http://dx.doi.org/10.1118/1.2733803 MPHYA6 0094-2405 Google Scholar

12.

B. KanmaniR. M. Vasu, “Diffuse optical tomography using intensity measurements and the a priori acquired regions of interest: theory and simulations,” Phys. Med. Biol., 50 (2), 247 –264 (2005). http://dx.doi.org/10.1088/0031-9155/50/2/005 PHMBA7 0031-9155 Google Scholar

13.

M. Guvenet al., “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: I,” Inv. Prob., 23 (3), 1115 –1133 (2007). http://dx.doi.org/10.1088/0266-5611/23/3/017 INPEEY 0266-5611 Google Scholar

14.

J. C. Yeet al., “Nonlinear multigrid algorithms for Bayesian optical diffusion tomography,” IEEE Trans. Image Process., 10 (6), 909 –922 (2001). http://dx.doi.org/10.1109/83.923278 IIPRE4 1057-7149 Google Scholar

15.

Z. M. Wanget al., “Experimental demonstration of an analytic method for image reconstruction in optical diffusion tomography with large data sets,” Opt. Lett., 30 (24), 3338 –3340 (2005). http://dx.doi.org/10.1364/OL.30.003338 OPLEDP 0146-9592 Google Scholar

16.

V. A. MarkelJ. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E, 70 (5), 056616 (2004). http://dx.doi.org/10.1103/PhysRevE.70.056616 PLEEE8 1063-651X Google Scholar

17.

M. J. D. Powell, “Radial basis functions for multivariable interpolation: a review,” Algorithms for Approximation, 143 –167 Clarendon Press, New York (1987). Google Scholar

18.

E. E. Vidal-Rosaset al., “Real-time diffuse optical tomography using reduced-order light propagation models,” in Proc. of the 8th IASTED Int. Conf. on Biomedical Engineering, 378 –385 (2011). Google Scholar

19.

A. Y. Bluestoneet al., “Three-dimensional optical tomographic brain imaging in small animals, Part 1: Hypercapnia,” J. Biomed. Opt., 9 (5), 1046 –1062 (2004). http://dx.doi.org/10.1117/1.1784471 JBOPFO 1083-3668 Google Scholar

20.

A. K. Dunnet al., “Spatial extent of oxygen metabolism and hemodynamic changes during functional activation of the rat somatosensory cortex,” Neuroimage, 27 (2), 279 –290 (2005). http://dx.doi.org/10.1016/j.neuroimage.2005.04.024 NEIMEF 1053-8119 Google Scholar

21.

J. P. Culveret al., “Evidence that cerebral blood volume can provide brain activation maps with better spatial resolution than deoxygenated hemoglobin,” Neuroimage, 27 (4), 947 –959 (2005). http://dx.doi.org/10.1016/j.neuroimage.2005.05.052 NEIMEF 1053-8119 Google Scholar

22.

J. P. Culveret al., “Volumetric diffuse optical tomography of brain activity,” Opt. Lett., 28 (21), 2061 –2063 (2003). http://dx.doi.org/10.1364/OL.28.002061 OPLEDP 0146-9592 Google Scholar

23.

A. M. Siegelet al., “Temporal comparison of functional brain imaging with diffuse optical tomography and fMRI during rat forepaw stimulation,” Phys. Med. Biol., 48 (10), 1391 –1403 (2003). http://dx.doi.org/10.1088/0031-9155/48/10/311 PHMBA7 0031-9155 Google Scholar

24.

C. H. Schmitzet al., “Dynamic studies of small animals with a four-color diffuse optical tomography imager,” Rev. Sci. Instrum., 76 (9), (2005). http://dx.doi.org/10.1063/1.2038467 RSINAK 0034-6748 Google Scholar

25.

M. Schweigeret al., “The finite-element method for the propagation of light in scattering media—boundary and source conditions,” Med. Phys., 22 (11), 1779 –1792 (1995). http://dx.doi.org/10.1118/1.597634 MPHYA6 0094-2405 Google Scholar

26.

K. D. PaulsenH. B. Jiang, “Spatially varying optical property reconstruction using a finite-element diffusion equation approximation,” Med. Phys., 22 (6), 691 –701 (1995). http://dx.doi.org/10.1118/1.597488 MPHYA6 0094-2405 Google Scholar

27.

A. H. HielscherR. E. AlcouffeR. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol., 43 (5), 1285 –1302 (1998). http://dx.doi.org/10.1088/0031-9155/43/5/017 PHMBA7 0031-9155 Google Scholar

28.

H. Wendland, “Computational aspects of radial basis function approximation,” Studies Comput. Math., 12 231 –256 (2006). http://dx.doi.org/10.1016/S1570-579X(06)80010-8 Google Scholar

29.

R. Franke, “Smooth interpolation of scattered data by local thin plate splines,” Comput. Math. Appl., 8 (4), 273 –281 (1982). http://dx.doi.org/10.1016/0898-1221(82)90009-8 CMAPDK 0898-1221 Google Scholar

30.

M. D. Buhmann, Radial Basis Functions: Theory and Implementations, Cambridge University Press, Cambridge (2003). Google Scholar

31.

S. R. Arridge, “Photon-measurement density-functions. 1. Analytical forms,” Appl. Opt., 34 (31), 7395 –7409 (1995). http://dx.doi.org/10.1364/AO.34.007395 APOPAI 0003-6935 Google Scholar

32.

S. R. ArridgeM. Schweiger, “Photon-measurement density-functions. 2. Finite-element-method calculations,” Appl. Opt., 34 (34), 8026 –8037 (1995). http://dx.doi.org/10.1364/AO.34.008026 APOPAI 0003-6935 Google Scholar

33.

E. E. Vidal-Rosas, “Advanced tomographic reconstruction methods for diffuse optical imaging of cerebral heamodynamic response,” University of Sheffield, (2011). Google Scholar

34.

S. A. BillingsS. ChenM. J. Korenberg, “Identification of mimo non-linear systems using a forward-regression orthogonal estimator,” Int. J. Ctrl., 49 (6), 2157 –2189 (1989). IJCOAZ 0020-7179 Google Scholar

35.

H. L. GraberY. L. PeiR. L. Barbour, “Imaging of spatiotemporal coincident states by DC optical tomography,” IEEE Trans. Med. Imaging, 21 (8), 852 –866 (2002). http://dx.doi.org/10.1109/TMI.2002.801154 ITMID4 0278-0062 Google Scholar

36.

J. NocedalS. Wright, Numerical Optimization, 2nd ed.Springer, New York (2006). Google Scholar

37.

J. QinR. Lu, “Hyperspectral diffuse reflectance for determination of the optical properties of milk and fruit and vegetable juice,” Proc. SPIE, 5996 232 –241 (2005). Google Scholar

38.

E. Grattonet al., “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. B: Biol. Sci., 352 (1354), 727 –735 (1997). http://dx.doi.org/10.1098/rstb.1997.0055 PTRBAE 0962-8436 Google Scholar

39.

J. B. Fishkinet al., “Gigahertz photon density waves in a turbid medium: theory and experiments,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics, 53 (3), 2307 –2319 (1996). http://dx.doi.org/10.1103/PhysRevE.53.2307 PLEEE8 1063-651X Google Scholar

40.

G. Grattonet al., “Feasibility of intracranial near-infrared optical scanning,” Psychophysiology, 31 (2), 211 –215 (1994). http://dx.doi.org/10.1111/psyp.1994.31.issue-2 PSPHAF 0048-5772 Google Scholar

41.

J. B. FishkinE. Gratton, “Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge,” J. Opt. Soc. Am. A, 10 (1), 127 –140 (1993). http://dx.doi.org/10.1364/JOSAA.10.000127 JOAOD6 0740-3232 Google Scholar

42.

W. F. CheongS. A. PrahlA. J. Welch, “A review of the optical-properties of biological tissues,” IEEE J. Quantum Electron., 26 (12), 2166 –2185 (1990). http://dx.doi.org/10.1109/3.64354 IEJQA7 0018-9197 Google Scholar

43.

A. J. Kennerleyet al., “Concurrent fMRI and optical measures for the investigation of the hemodynamic response function,” Magn. Reson. Med., 54 (2), 354 –365 (2005). http://dx.doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

44.

A. Y. Bluestoneet al., “Three-dimensional optical tomographic brain imaging in small animals, Part 2: Unilateral carotid occlusion,” J. Biomed. Opt., 9 (5), 1063 –1073 (2004). http://dx.doi.org/10.1117/1.1784472 JBOPFO 1083-3668 Google Scholar

45.

A. T. Eggebrechtet al., “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” Neuroimage, 61 (4), 1120 –1128 (2012). http://dx.doi.org/10.1016/j.neuroimage.2012.01.124 NEIMEF 1053-8119 Google Scholar

46.

A. H. HielscherA. D. KloseK. M. Hanson, “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging, 18 (3), 262 –271 (1999). http://dx.doi.org/10.1109/42.764902 ITMID4 0278-0062 Google Scholar

47.

J. B. Mandevilleet al., “Dynamic functional imaging of relative cerebral blood volume during rat forepaw stimulation,” Magn. Reson. Med., 39 (4), 615 –624 (1998). http://dx.doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

48.

I. Troprèset al., “Vessel size imaging,” Magn. Reson. Med., 45 (3), 397 –408 (2001). http://dx.doi.org/10.1002/(ISSN)1522-2594 MRMEEN 0740-3194 Google Scholar

49.

S. M. Liaoet al., “High-density diffuse optical tomography of term infant visual cortex in the nursery,” J. Biomed. Opt., 17 (8), 081414 (2012). http://dx.doi.org/10.1117/1.JBO.17.8.081414 JBOPFO 1083-3668 Google Scholar

50.

D. K. Josephet al., “Diffuse optical tomography system to image brain activation with improved spatial resolution and validation with functional magnetic resonance imaging,” Appl. Opt., 45 (31), 8142 –8151 (2006). http://dx.doi.org/10.1364/AO.45.008142 APOPAI 0003-6935 Google Scholar

51.

S. P. Kochet al., “High-resolution optical functional mapping of the human somatosensory cortex,” Front. Neuroenerget., 2 (5), 1 –8 (2010). http://dx.doi.org/10.3389/fnene.2010.00012 FNREJG 1662-6427 Google Scholar

52.

B. W. Zeffet al., “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Natl. Acad. Sci. U. S. A., 104 (29), 12169 –12174 (2007). http://dx.doi.org/10.1073/pnas.0611266104 PNASA6 0027-8424 Google Scholar

53.

B. R. WhiteJ. P. Culver, “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt., 15 (2), 026006 (2010). http://dx.doi.org/10.1117/1.3368999 JBOPFO 1083-3668 Google Scholar

54.

S. L. Ferradalet al., “Atlas-based head modeling and spatial normalization for high-density diffuse optical tomography: in vivo validation against fMRI,” Neuroimage, 85 (Part 1), 117 –126 (2013). http://dx.doi.org/10.1016/j.neuroimage.2013.03.069 NEIMEF 1053-8119 Google Scholar

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.

Biographies for other authors are not available.

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Ernesto E. Vidal-Rosas, Stephen Billings, Ying Zheng, John E. Mayhew, David K. Johnston, Aneurin J. Kennerley, and Daniel Coca "Reduced-order modeling of light transport in tissue for real-time monitoring of brain hemodynamics using diffuse optical tomography," Journal of Biomedical Optics 19(2), 026008 (13 February 2014). https://doi.org/10.1117/1.JBO.19.2.026008

Published: 13 February 2014

Access the abstract

JOURNAL ARTICLE
15 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 3 scholarly publications.

Explore citations on Lens.org

KEYWORDS

3D modeling

Tissue optics

Data modeling

Absorption

Finite element methods

Reconstruction algorithms

Hemodynamics

1.

Introduction

2.

Methods

2.1.

Animal Preparation and Experimental Procedure

2.2.

Instrumentation

2.3.

Image Reconstruction Algorithm

2.3.1.

Forward problem

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

2.3.2.

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

Fig. 1

2.3.3.

Reduced-order forward model

Eq. (5)

Eq. (6)

Eq. (7)

2.3.4.

Model estimation

2.3.5.

Reduced-order forward model of the rat’s head

2.3.6.

Image reconstruction

Eq. (8)

2.4.

Calculation of Hemoglobin Concentration

Eq. (9)

Eq. (10)

3.

Results

3.1.

Algorithm Evaluation Using a Phantom

3.1.1.

Experiment setup

Fig. 2

3.1.2.

Estimation and validation of reduced-order model

Fig. 3

Fig. 4

3.1.3.

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

Eq. (11)

Fig. 5

3.2.

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

Fig. 6

Fig. 7

Fig. 8

3.3.

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

Fig. 9

Fig. 10

Eq. (12)

Eq. (13)

3.4.

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

Fig. 11

Fig. 12

4.

Discussion

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years