Research Papers: Imaging

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

[+] Author Affiliations
Ernesto E. Vidal-Rosas

University of Sheffield, Department of Automatic Control and Systems Engineering, Mappin Street, Sheffield, S1 3JD, United Kingdom

Stephen A. Billings

University of Sheffield, Department of Automatic Control and Systems Engineering, Mappin Street, Sheffield, S1 3JD, United Kingdom

Ying Zheng

University of Sheffield, Centre for Signal Processing in Neuroimaging and Systems Neuroscience, Department of Psychology, Western Bank, Sheffield, S10 2TN, United Kingdom

John E. Mayhew, David Johnston

University of Sheffield, Centre for Signal Processing in Neuroimaging and Systems Neuroscience, Department of Psychology, Western Bank, Sheffield, S10 2TN, United Kingdom

Aneurin J. Kennerley

University of Sheffield, Centre for Signal Processing in Neuroimaging and Systems Neuroscience, Department of Psychology, Western Bank, Sheffield, S10 2TN, United Kingdom

Daniel Coca

University of Sheffield, Department of Automatic Control and Systems Engineering, Mappin Street, Sheffield, S1 3JD, United Kingdom

J. Biomed. Opt. 19(2), 026008 (Feb 13, 2014). doi:10.1117/1.JBO.19.2.026008
History: Received June 21, 2013; Revised January 8, 2014; Accepted January 10, 2014
Text Size: A A A

Open Access Open Access

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 3s per time point instead of the 1 to 2 h it takes when using the conventional finite element modeling approach.

Figures in this Article

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 Cummer7 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 (36times 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 optimization13,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.

Animal Preparation and Experimental Procedure

Small animals provide excellent models to investigate ischemia,19 tissue oxygenation,20 and hemodynamics.2123 This experiment focuses on the hemodynamic response under hypercapnia.

The animal used was a female hooded Lister rat weighing 400g, 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.25g/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.26mg/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% O2, 80% N2). 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.

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 4Hz. 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.

Image Reconstruction Algorithm
Forward problem

For a medium ΩIR3 with boundary Ω, characterized by absorption and scattering coefficients μa(r), μs(r), rΩ, solving the forward problem involves numerical simulation of photon fluence rate measurements {y(j)}j=1,s from d detectors y(j)=[y1(j),,yd(j)] on Ω, for each point source qj of the array of sources q=[q1,,qs] distributed on Ω.

The forward problem for the source qj is described by the following parameters-to-output mapping: Display Formula

y(j)=Pju(r),j=1,,s,(1)
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.2Display Formula

·D(r)ϕj(r)+μa(r)ϕj(r)=qj(r),rΩ,(2)
where ϕj(r) is the spatially varying diffuse photon density at position r given the source qj, μ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Ω is usually represented by an isotropic point source qj(r)=δ(rrj), where rj 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 n(ξj). The boundary condition usually employed is of the Robin type. Display Formula
D(ξ)ϕj(ξ)n+12Aϕj(ξ)=0,ξΩ,(3)
where the term A accounts for the refractive index boundary mismatch between the interior and exterior medium.25

For any given source qj, the variable measured by the detector located at ξiΩ is the outward flux γj(ξi). The corresponding measurement equations are given by Display Formula

γj(ξ)=D(ξ)n(ξ)·ϕj(ξ),ξΩ,yi(j)=γj(ξi),j=1,,s;i=1,,d.(4)

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.

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.

Graphic Jump LocationF1 :

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.02mm1 for skin, μa=0.005mm1 for skull, μa=0.015mm1 for brain, and μa=0.022mm1 for muscle; similarly, μs=0.5mm1 for skin, μs=1.63mm1 for skull, μs=1.63mm1 for brain, and μs=1mm1 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.

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 Display Formula

y^i,j(t)=fi,j[u(t)]+ei,j(t),j=1,,s;i=1,,Mj;ij,(5)
where 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 qj. fi,j is a nonlinear function for source-detector pair i-j, u(t)=[u1(t),,uK(t)], uk(t)=μa(rk,t) is the absorption value for the k’th node at a given time t, K is the total number of nodes, and ei,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, ND=NE+NV (typically ND=1000, NE=NV) 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. NE data points are used to estimate the model, and the remaining NV data points are used to validate the model. For each input u(t)=[μa(r1,t),,μa(rK,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 Display Formula

y^i,j(t)=l=1Ni,jθli,jpli,j[u(t)]+ei,j,(6)
where pki,j(x) is an RBF. The RBF chosen in this study is the thin-plate spline. Display Formula
pki,j(u)=φ(u,ck)=rk2log(rk),(7)
where ck denotes the center of the RBF and rk=uck2 is the Euclidian norm. Thin-plate splines are conditionally positive kernels that are invariant under scaling,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 Di,jE={u(t),yi,j(t)}t=1,NE, the mapping in Eq. (5) is derived such that it provides a sparse approximation of the underlying nonlinear function using Nij terms (Ni,jND). 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, nKi,jK.

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.

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) ui,j={uli,j}lIi,j, where Ii,j is a subset of the full set of nodes {1,,K}, with dimension nKi,jK.
  2. Given a dictionary of candidate RBF basis functions, determine a minimal subset of regressors {pk}k=1,Ni,j, Ni,jND and the associated parameter vector {θk}k=1,Ni,j, which minimize the prediction error over the data set Di,j.

  • 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 IROI=i,jIi,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
  • 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., ck=ui,j(k), k=1,,NE. The minimal set of basis functions in Eq. (6) is determined using a greedy orthogonal regression algorithm34 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 Display Formula
    100i=1nerriCd,
    where erri is the % contribution of the i’th selected term and Cd 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.

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[0.0075,0.0045]mm1 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 ij as Display Formula

RMSEij=t=1N[yi,j(t)y^i,j(t)]2N
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 Cd values. The thin-plane-spline (TPS)-RBF models used in reconstruction (estimated for Cd=0.3) had between 8 and 399 inputs and between 11 and 92 RBFs.

Image reconstruction

The absorption coefficients were obtained using the modified version19 of the normalized difference method,35 which involves solving the following optimization problem: Display Formula

uROI*=argminuROIj=1sj=1Mj[yi,j(t)yi,j0y^i,j(0)fi,j(uROI)yi,j(t)yi,j0y^i,j(0)]2,(8)
where yi,j(t) is a measurement taken at time t, yi,j0 is the predefined reference state, y^i,j(0)=fi,j(uROI0) is the predicted measurement corresponding to the reference medium, and uROI0 is the initial vector of the optical parameters specified in Sec. 2.3.2. The optimization problem was solved using an iterative, conjugate-gradient algorithm36 using spatial constraints.5

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 (Δ[HbO2] and Δ[HbR]).19Display Formula

Δμa(λ)=εHbO2(λ)Δ[HbO2]+εHbR(λ)Δ[HbR],(9)
where εHbO2(λ) and εHbR(λ) are the extinction coefficients for oxyhemoglobin and deoxyhemoglobin, respectively, and λ is the photon wavelength.5

The concentration changes in oxyhemoglobin (HbO2) and deoxyhemoglobin (HbR) were calculated by solving a linear system of equations, using absorption coefficients estimated from measurements at λ1=760 and λ2=830nm, as follows:19Display Formula

Δ[HbR]=εHbO2(λ1)Δμa(λ2)εHbO2(λ2)Δμa(λ1)εHbO2(λ1)εHbR(λ2)εHbO2(λ2)εHbR(λ1)Δ[HbO2]=εHbR(λ2)Δμa(λ1)εHbR(λ1)Δμa(λ2)εHbO2(λ1)εHbR(λ2)εHbO2(λ2)εHbR(λ1).(10)

Algorithm Evaluation Using a Phantom
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 media37 at 810 nm are μa=0.0024mm1 and μs=0.5mm1. Similarly, for the perturbation, μa=0.024mm1 and μs=0.4mm1.

Graphic Jump LocationF2 :

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 media3841 because its optical absorption and scattering parameters are close to those of most living tissue,42 where 0.001<μa<0.01mm1 and 0.5<μs<10mm1.

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 240source/detector channels.

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.03mm1. 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).

Graphic Jump LocationF3 :

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.

Graphic Jump LocationF4 :

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.

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

Four different TPS-RBF models with Cd=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 Cd=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 by35Display Formula

ICC(A,B)=1N1i(xAix¯Ai)(xBix¯Bi)i(xAix¯Ai)2i(xBix¯Bi)2,(11)
where xAi and xBi denote the intensities of the i’th pixel in images A and B, respectively, and xAi and xBi 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(Cd=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 80s; the ROM required only six iterations and took 2.5s.

Graphic Jump LocationF5 :

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

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.

Graphic Jump LocationF6 :

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.0013mm1 and max(Δμa)ROM=0.0025mm1. 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 2s, while the FEM-based method achieved similar quality after five iterations, but at a considerably larger amount of time (1h).

Graphic Jump LocationF7 :

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

Graphic Jump LocationF8 :

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.

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

Graphic Jump LocationF9 :

Simulated hypercapnia using the rat head mesh. (a) and (b) Location and size of the simulated absorber. (c) Oxyhemoglobin (HbO2) 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 Cd=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 6s 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=120s, reconstructed using the ROM approach.

Graphic Jump LocationF10 :

Reconstructed changes in HbO2 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=120s. (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 (HbO2), 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: Display Formula
    Δ[HbO2]FEM,cal=5.26Δ[HbO2]FEM0.94Δ[HbR]FEM,cal=5.26Δ[HbR]FEM+0.39.(12)
  • ROM approach: Display Formula
    Δ[HbO2]ROM,cal=0.93Δ[HbO2]ROM2.81Δ[HbR]ROM,cal=0.91Δ[HbR]ROM0.90.(13)

The calibrated changes of HbO2, 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.

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 3Hz. 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 3s, 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 CO2 was at ton=t0+60s, t0=60s, with the offset at toff=ton+120s).

Figure 11(a) shows a representative coronal section CBV-MRI statistical parameter map superimposed on a 256×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=120s, 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 43 to measure changes in blood volume following hypercapnia and whisker stimulation.

Graphic Jump LocationF11 :

Reconstructed changes in HbO2, 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=120s. The two rectangles mark the superficial and deep cortical areas used to compute blood volumes. (c) Reconstructed changes in HbO2 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 (HbO2) 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 CBVsup 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 CBVsup 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 CBVdeep peaks at 33.8 μM, while HbT (deep) time series reconstructed from DOT measurements using the ROM approach peak at 30.3μM. This reflects the reduced depth sensitivity of DOT compared with fMRI.

Figure 12 displays the evolution of Δ[HbO2](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=0s 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 CO2. The reconstructed spatial maps agree with similar DOT and fMRI studies.43,44

Graphic Jump LocationF12 :

Spatial maps of reconstructed changes in HbO2 and HbR in the rat brain based on experimental DOT measurements at different time points during hypercapnia.

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 RBFs17 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 Display Formula

M=i=1qni,whereni=ni1(n+i1)i,n0=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,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 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 HbO2, 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 study43 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 3s 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 newborns49 as well as adults5052 in real time and at much higher resolution53 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.

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

Arridge  S. R., Schotland  J. C., “Optical tomography: forward and inverse problems,” Inv. Prob.. 25, (12 ), 123010  (2009). 0266-5611 CrossRef
Arridge  S. R., “Optical tomography in medical imaging,” Inv. Prob.. 15, (2 ), R41 –R93 (1999). 0266-5611 CrossRef
Gibson  A. P., Hebden  J. C., Arridge  S. R., “Recent advances in diffuse optical imaging,” Phys. Med. Biol.. 50, (4 ), R1 –R43 (2005). 0031-9155 CrossRef
Yalavarthy  P. K. et al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express. 15, (13 ), 8043 –8058 (2007). 1094-4087 CrossRef
Boas  D. A., Dale  A. M., Franceschini  M. A., “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage. 23, , S275 –S288 (2004). 1053-8119 CrossRef
Yalavarthy  P. K. et 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). 0094-2405 CrossRef
Zhai  Y., Cummer  S. A., “Fast tomographic reconstruction strategy for diffuse optical tomography,” Opt. Express. 17, (7 ), 5285 –5297 (2009). 1094-4087 CrossRef
Kolehmainen  V. et al., “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A. 26, (10 ), 2257 –2268 (2009). 0740-3232 CrossRef
Landragin-Frassati  A. et al., “Application of a wavelet-Galerkin method to the forward problem resolution in fluorescence diffuse optical tomography,” Opt. Express. 17, (21 ), 18433 –18448 (2009). 1094-4087 CrossRef
Eames  M. E. et al., “An efficient Jacobian reduction method for diffuse optical image reconstruction,” Opt. Express. 15, (24 ), 15908 –15919 (2007). 1094-4087 CrossRef
Yalavarthy  P. K. et al., “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys.. 34, (6 ), 2085 –2098 (2007). 0094-2405 CrossRef
Kanmani  B., Vasu  R. M., “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). 0031-9155 CrossRef
Guven  M. et al., “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: I,” Inv. Prob.. 23, (3 ), 1115 –1133 (2007). 0266-5611 CrossRef
Ye  J. C. et al., “Nonlinear multigrid algorithms for Bayesian optical diffusion tomography,” IEEE Trans. Image Process.. 10, (6 ), 909 –922 (2001). 1057-7149 CrossRef
Wang  Z. M. et 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). 0146-9592 CrossRef
Markel  V. A., Schotland  J. C., “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E. 70, (5 ), 056616  (2004). 1063-651X CrossRef
Powell  M. J. D., “Radial basis functions for multivariable interpolation: a review,” in Algorithms for Approximation. , pp. 143 –167,  Clarendon Press ,  New York  (1987).
Vidal-Rosas  E. E. et al., “Real-time diffuse optical tomography using reduced-order light propagation models,” in  Proc. of the 8th IASTED Int. Conf. on Biomedical Engineering , pp. 378 –385,  ACTA Press  (2011).
Bluestone  A. Y. et al., “Three-dimensional optical tomographic brain imaging in small animals, Part 1: Hypercapnia,” J. Biomed. Opt.. 9, (5 ), 1046 –1062 (2004). 1083-3668 CrossRef
Dunn  A. K. et al., “Spatial extent of oxygen metabolism and hemodynamic changes during functional activation of the rat somatosensory cortex,” Neuroimage. 27, (2 ), 279 –290 (2005). 1053-8119 CrossRef
Culver  J. P. et al., “Evidence that cerebral blood volume can provide brain activation maps with better spatial resolution than deoxygenated hemoglobin,” Neuroimage. 27, (4 ), 947 –959 (2005). 1053-8119 CrossRef
Culver  J. P. et al., “Volumetric diffuse optical tomography of brain activity,” Opt. Lett.. 28, (21 ), 2061 –2063 (2003). 0146-9592 CrossRef
Siegel  A. M. et 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). 0031-9155 CrossRef
Schmitz  C. H. et al., “Dynamic studies of small animals with a four-color diffuse optical tomography imager,” Rev. Sci. Instrum.. 76, (9 ) (2005). 0034-6748 CrossRef
Schweiger  M. et al., “The finite-element method for the propagation of light in scattering media—boundary and source conditions,” Med. Phys.. 22, (11 ), 1779 –1792 (1995). 0094-2405 CrossRef
Paulsen  K. D., Jiang  H. B., “Spatially varying optical property reconstruction using a finite-element diffusion equation approximation,” Med. Phys.. 22, (6 ), 691 –701 (1995). 0094-2405 CrossRef
Hielscher  A. H., Alcouffe  R. E., Barbour  R. L., “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol.. 43, (5 ), 1285 –1302 (1998). 0031-9155 CrossRef
Wendland  H., “Computational aspects of radial basis function approximation,” Studies Comput. Math.. 12, , 231 –256 (2006).CrossRef
Franke  R., “Smooth interpolation of scattered data by local thin plate splines,” Comput. Math. Appl.. 8, (4 ), 273 –281 (1982). 0898-1221 CrossRef
Buhmann  M. D., Radial Basis Functions: Theory and Implementations. , Cambridge Monographs on Applied and Computational Mathematics ,  Cambridge University Press ,  Cambridge  (2003).
Arridge  S. R., “Photon-measurement density-functions. 1. Analytical forms,” Appl. Opt.. 34, (31 ), 7395 –7409 (1995). 0003-6935 CrossRef
Arridge  S. R., Schweiger  M., “Photon-measurement density-functions. 2. Finite-element-method calculations,” Appl. Opt.. 34, (34 ), 8026 –8037 (1995). 0003-6935 CrossRef
Vidal-Rosas  E. E., “Advanced tomographic reconstruction methods for diffuse optical imaging of cerebral heamodynamic response,” PhD Thesis, University of Sheffield (2011).
Billings  S. A., Chen  S., Korenberg  M. J., “Identification of mimo non-linear systems using a forward-regression orthogonal estimator,” Int. J. Ctrl.. 49, (6 ), 2157 –2189 (1989). 0020-7179 
Graber  H. L., Pei  Y. L., Barbour  R. L., “Imaging of spatiotemporal coincident states by DC optical tomography,” IEEE Trans. Med. Imaging. 21, (8 ), 852 –866 (2002). 0278-0062 CrossRef
Nocedal  J., Wright  S., Numerical Optimization. , Springer Series in Operations Research , 2nd ed.,  Springer ,  New York  (2006).
Qin  J., Lu  R., “Hyperspectral diffuse reflectance for determination of the optical properties of milk and fruit and vegetable juice,” Proc. SPIE. 5996, , 232 –241 (2005).
Gratton  E. et al., “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. B: Biol. Sci.. 352, (1354 ), 727 –735 (1997). 0962-8436 CrossRef
Fishkin  J. B. et 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). 1063-651X CrossRef
Gratton  G. et al., “Feasibility of intracranial near-infrared optical scanning,” Psychophysiology. 31, (2 ), 211 –215 (1994). 0048-5772 CrossRef
Fishkin  J. B., Gratton  E., “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). 0740-3232 CrossRef
Cheong  W. F., Prahl  S. A., Welch  A. J., “A review of the optical-properties of biological tissues,” IEEE J. Quantum Electron.. 26, (12 ), 2166 –2185 (1990). 0018-9197 CrossRef
Kennerley  A. J. et al., “Concurrent fMRI and optical measures for the investigation of the hemodynamic response function,” Magn. Reson. Med.. 54, (2 ), 354 –365 (2005). 0740-3194 CrossRef
Bluestone  A. Y. et al., “Three-dimensional optical tomographic brain imaging in small animals, Part 2: Unilateral carotid occlusion,” J. Biomed. Opt.. 9, (5 ), 1063 –1073 (2004). 1083-3668 CrossRef
Eggebrecht  A. T. et al., “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” Neuroimage. 61, (4 ), 1120 –1128 (2012). 1053-8119 CrossRef
Hielscher  A. H., Klose  A. D., Hanson  K. M., “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging. 18, (3 ), 262 –271 (1999). 0278-0062 CrossRef
Mandeville  J. B. et al., “Dynamic functional imaging of relative cerebral blood volume during rat forepaw stimulation,” Magn. Reson. Med.. 39, (4 ), 615 –624 (1998). 0740-3194 CrossRef
Troprès  I. et al., “Vessel size imaging,” Magn. Reson. Med.. 45, (3 ), 397 –408 (2001). 0740-3194 CrossRef
Liao  S. M. et al., “High-density diffuse optical tomography of term infant visual cortex in the nursery,” J. Biomed. Opt.. 17, (8 ), 081414  (2012). 1083-3668 CrossRef
Joseph  D. K. et 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). 0003-6935 CrossRef
Koch  S. P. et al., “High-resolution optical functional mapping of the human somatosensory cortex,” Front. Neuroenerget.. 2, (5 ), 1 –8 (2010). 1662-6427 CrossRef
Zeff  B. W. et 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). 0027-8424 CrossRef
White  B. R., Culver  J. P., “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt.. 15, (2 ), 026006  (2010). 1083-3668 CrossRef
Ferradal  S. L. et 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). 1053-8119 CrossRef

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.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.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

Ernesto E. Vidal-Rosas ; Stephen A. Billings ; Ying Zheng ; John E. Mayhew ; David Johnston, et al.
"Reduced-order modeling of light transport in tissue for real-time monitoring of brain hemodynamics using diffuse optical tomography", J. Biomed. Opt. 19(2), 026008 (Feb 13, 2014). ; http://dx.doi.org/10.1117/1.JBO.19.2.026008


Figures

Graphic Jump LocationF1 :

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.

Graphic Jump LocationF2 :

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.

Graphic Jump LocationF6 :

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.

Graphic Jump LocationF3 :

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.

Graphic Jump LocationF4 :

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.

Graphic Jump LocationF11 :

Reconstructed changes in HbO2, 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=120s. The two rectangles mark the superficial and deep cortical areas used to compute blood volumes. (c) Reconstructed changes in HbO2 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.

Graphic Jump LocationF5 :

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

Graphic Jump LocationF7 :

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.

Graphic Jump LocationF8 :

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.

Graphic Jump LocationF9 :

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

Graphic Jump LocationF10 :

Reconstructed changes in HbO2 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=120s. (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.

Graphic Jump LocationF12 :

Spatial maps of reconstructed changes in HbO2 and HbR in the rat brain based on experimental DOT measurements at different time points during hypercapnia.

Tables

References

Arridge  S. R., Schotland  J. C., “Optical tomography: forward and inverse problems,” Inv. Prob.. 25, (12 ), 123010  (2009). 0266-5611 CrossRef
Arridge  S. R., “Optical tomography in medical imaging,” Inv. Prob.. 15, (2 ), R41 –R93 (1999). 0266-5611 CrossRef
Gibson  A. P., Hebden  J. C., Arridge  S. R., “Recent advances in diffuse optical imaging,” Phys. Med. Biol.. 50, (4 ), R1 –R43 (2005). 0031-9155 CrossRef
Yalavarthy  P. K. et al., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express. 15, (13 ), 8043 –8058 (2007). 1094-4087 CrossRef
Boas  D. A., Dale  A. M., Franceschini  M. A., “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” Neuroimage. 23, , S275 –S288 (2004). 1053-8119 CrossRef
Yalavarthy  P. K. et 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). 0094-2405 CrossRef
Zhai  Y., Cummer  S. A., “Fast tomographic reconstruction strategy for diffuse optical tomography,” Opt. Express. 17, (7 ), 5285 –5297 (2009). 1094-4087 CrossRef
Kolehmainen  V. et al., “Approximation errors and model reduction in three-dimensional diffuse optical tomography,” J. Opt. Soc. Am. A. 26, (10 ), 2257 –2268 (2009). 0740-3232 CrossRef
Landragin-Frassati  A. et al., “Application of a wavelet-Galerkin method to the forward problem resolution in fluorescence diffuse optical tomography,” Opt. Express. 17, (21 ), 18433 –18448 (2009). 1094-4087 CrossRef
Eames  M. E. et al., “An efficient Jacobian reduction method for diffuse optical image reconstruction,” Opt. Express. 15, (24 ), 15908 –15919 (2007). 1094-4087 CrossRef
Yalavarthy  P. K. et al., “Weight-matrix structured regularization provides optimal generalized least-squares estimate in diffuse optical tomography,” Med. Phys.. 34, (6 ), 2085 –2098 (2007). 0094-2405 CrossRef
Kanmani  B., Vasu  R. M., “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). 0031-9155 CrossRef
Guven  M. et al., “Effect of discretization error and adaptive mesh generation in diffuse optical absorption imaging: I,” Inv. Prob.. 23, (3 ), 1115 –1133 (2007). 0266-5611 CrossRef
Ye  J. C. et al., “Nonlinear multigrid algorithms for Bayesian optical diffusion tomography,” IEEE Trans. Image Process.. 10, (6 ), 909 –922 (2001). 1057-7149 CrossRef
Wang  Z. M. et 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). 0146-9592 CrossRef
Markel  V. A., Schotland  J. C., “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E. 70, (5 ), 056616  (2004). 1063-651X CrossRef
Powell  M. J. D., “Radial basis functions for multivariable interpolation: a review,” in Algorithms for Approximation. , pp. 143 –167,  Clarendon Press ,  New York  (1987).
Vidal-Rosas  E. E. et al., “Real-time diffuse optical tomography using reduced-order light propagation models,” in  Proc. of the 8th IASTED Int. Conf. on Biomedical Engineering , pp. 378 –385,  ACTA Press  (2011).
Bluestone  A. Y. et al., “Three-dimensional optical tomographic brain imaging in small animals, Part 1: Hypercapnia,” J. Biomed. Opt.. 9, (5 ), 1046 –1062 (2004). 1083-3668 CrossRef
Dunn  A. K. et al., “Spatial extent of oxygen metabolism and hemodynamic changes during functional activation of the rat somatosensory cortex,” Neuroimage. 27, (2 ), 279 –290 (2005). 1053-8119 CrossRef
Culver  J. P. et al., “Evidence that cerebral blood volume can provide brain activation maps with better spatial resolution than deoxygenated hemoglobin,” Neuroimage. 27, (4 ), 947 –959 (2005). 1053-8119 CrossRef
Culver  J. P. et al., “Volumetric diffuse optical tomography of brain activity,” Opt. Lett.. 28, (21 ), 2061 –2063 (2003). 0146-9592 CrossRef
Siegel  A. M. et 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). 0031-9155 CrossRef
Schmitz  C. H. et al., “Dynamic studies of small animals with a four-color diffuse optical tomography imager,” Rev. Sci. Instrum.. 76, (9 ) (2005). 0034-6748 CrossRef
Schweiger  M. et al., “The finite-element method for the propagation of light in scattering media—boundary and source conditions,” Med. Phys.. 22, (11 ), 1779 –1792 (1995). 0094-2405 CrossRef
Paulsen  K. D., Jiang  H. B., “Spatially varying optical property reconstruction using a finite-element diffusion equation approximation,” Med. Phys.. 22, (6 ), 691 –701 (1995). 0094-2405 CrossRef
Hielscher  A. H., Alcouffe  R. E., Barbour  R. L., “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol.. 43, (5 ), 1285 –1302 (1998). 0031-9155 CrossRef
Wendland  H., “Computational aspects of radial basis function approximation,” Studies Comput. Math.. 12, , 231 –256 (2006).CrossRef
Franke  R., “Smooth interpolation of scattered data by local thin plate splines,” Comput. Math. Appl.. 8, (4 ), 273 –281 (1982). 0898-1221 CrossRef
Buhmann  M. D., Radial Basis Functions: Theory and Implementations. , Cambridge Monographs on Applied and Computational Mathematics ,  Cambridge University Press ,  Cambridge  (2003).
Arridge  S. R., “Photon-measurement density-functions. 1. Analytical forms,” Appl. Opt.. 34, (31 ), 7395 –7409 (1995). 0003-6935 CrossRef
Arridge  S. R., Schweiger  M., “Photon-measurement density-functions. 2. Finite-element-method calculations,” Appl. Opt.. 34, (34 ), 8026 –8037 (1995). 0003-6935 CrossRef
Vidal-Rosas  E. E., “Advanced tomographic reconstruction methods for diffuse optical imaging of cerebral heamodynamic response,” PhD Thesis, University of Sheffield (2011).
Billings  S. A., Chen  S., Korenberg  M. J., “Identification of mimo non-linear systems using a forward-regression orthogonal estimator,” Int. J. Ctrl.. 49, (6 ), 2157 –2189 (1989). 0020-7179 
Graber  H. L., Pei  Y. L., Barbour  R. L., “Imaging of spatiotemporal coincident states by DC optical tomography,” IEEE Trans. Med. Imaging. 21, (8 ), 852 –866 (2002). 0278-0062 CrossRef
Nocedal  J., Wright  S., Numerical Optimization. , Springer Series in Operations Research , 2nd ed.,  Springer ,  New York  (2006).
Qin  J., Lu  R., “Hyperspectral diffuse reflectance for determination of the optical properties of milk and fruit and vegetable juice,” Proc. SPIE. 5996, , 232 –241 (2005).
Gratton  E. et al., “Measurements of scattering and absorption changes in muscle and brain,” Philos. Trans. R. Soc. B: Biol. Sci.. 352, (1354 ), 727 –735 (1997). 0962-8436 CrossRef
Fishkin  J. B. et 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). 1063-651X CrossRef
Gratton  G. et al., “Feasibility of intracranial near-infrared optical scanning,” Psychophysiology. 31, (2 ), 211 –215 (1994). 0048-5772 CrossRef
Fishkin  J. B., Gratton  E., “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). 0740-3232 CrossRef
Cheong  W. F., Prahl  S. A., Welch  A. J., “A review of the optical-properties of biological tissues,” IEEE J. Quantum Electron.. 26, (12 ), 2166 –2185 (1990). 0018-9197 CrossRef
Kennerley  A. J. et al., “Concurrent fMRI and optical measures for the investigation of the hemodynamic response function,” Magn. Reson. Med.. 54, (2 ), 354 –365 (2005). 0740-3194 CrossRef
Bluestone  A. Y. et al., “Three-dimensional optical tomographic brain imaging in small animals, Part 2: Unilateral carotid occlusion,” J. Biomed. Opt.. 9, (5 ), 1063 –1073 (2004). 1083-3668 CrossRef
Eggebrecht  A. T. et al., “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” Neuroimage. 61, (4 ), 1120 –1128 (2012). 1053-8119 CrossRef
Hielscher  A. H., Klose  A. D., Hanson  K. M., “Gradient-based iterative image reconstruction scheme for time-resolved optical tomography,” IEEE Trans. Med. Imaging. 18, (3 ), 262 –271 (1999). 0278-0062 CrossRef
Mandeville  J. B. et al., “Dynamic functional imaging of relative cerebral blood volume during rat forepaw stimulation,” Magn. Reson. Med.. 39, (4 ), 615 –624 (1998). 0740-3194 CrossRef
Troprès  I. et al., “Vessel size imaging,” Magn. Reson. Med.. 45, (3 ), 397 –408 (2001). 0740-3194 CrossRef
Liao  S. M. et al., “High-density diffuse optical tomography of term infant visual cortex in the nursery,” J. Biomed. Opt.. 17, (8 ), 081414  (2012). 1083-3668 CrossRef
Joseph  D. K. et 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). 0003-6935 CrossRef
Koch  S. P. et al., “High-resolution optical functional mapping of the human somatosensory cortex,” Front. Neuroenerget.. 2, (5 ), 1 –8 (2010). 1662-6427 CrossRef
Zeff  B. W. et 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). 0027-8424 CrossRef
White  B. R., Culver  J. P., “Quantitative evaluation of high-density diffuse optical tomography: in vivo resolution and mapping performance,” J. Biomed. Opt.. 15, (2 ), 026006  (2010). 1083-3668 CrossRef
Ferradal  S. L. et 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). 1053-8119 CrossRef

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging & repositioning the boxes below.

Related Book Chapters

Topic Collections

PubMed Articles
Advertisement
  • Don't have an account?
  • Subscribe to the SPIE Digital Library
  • Create a FREE account to sign up for Digital Library content alerts and gain access to institutional subscriptions remotely.
Access This Article
Sign in or Create a personal account to Buy this article ($20 for members, $25 for non-members).
Access This Proceeding
Sign in or Create a personal account to Buy this article ($15 for members, $18 for non-members).
Access This Chapter

Access to SPIE eBooks is limited to subscribing institutions and is not available as part of a personal subscription. Print or electronic versions of individual SPIE books may be purchased via SPIE.org.