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

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.


Introduction
Diffuse optical tomography (DOT) is a noninvasive imaging modality that employs near-infrared light to interrogate optical properties of biological tissue. 1,2Compared 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,5On 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. 6This approximation-versus-accuracy dichotomy is one of the computational bottlenecks in DOT, which becomes highly relevant for the 3-D case.As a consequence, DOT is not used routinely for bedside monitoring or outside the hospital environment.
Previous efforts to address the computational bottleneck have focused on either reducing the complexity of the forward model while maintaining accuracy (model order reduction) or accelerating the convergence of the numerical reconstruction algorithms.Zhai and Cummer 7 used model order reduction techniques, common in structural dynamics, to decrease the complexity of the system of equations derived from finite element method (FEM); an adjustable parameter controls the level of the approximation.This solution decreased by an order of magnitude the time that the conventional approach can take.Kolehmainen et al. 8 employed a Bayesian approximation error model to compensate the loss of accuracy when using a coarse mesh; however, reconstruction times are still far from achieving real-time continuous monitoring.
Alternative approximation schemes for the diffusion equations (DE) using, for example, orthonormal compactly supported wavelets, 9 which produce sparse stiffness matrices, can reportedly reduce the time required to solve the forward problem (∼36 times faster for a 3-D geometry).To speed up parameter estimation, Eames et al. 10 proposed a sensitivitybased 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,12However, this choice may lead to erroneous results if the models are not coregistered. 4ther strategies employed to reduce the complexity of the reconstruction algorithms include the use of nonlinear multigrid optimization 13,14 and adaptive mesh techniques. 13,14Solving 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. 16his 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. 18Here, 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
2][23] This experiment focuses on the hemodynamic response under hypercapnia.
The animal used was a female hooded Lister rat weighing ∼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.

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 ∼4 Hz.A more detailed description of the scanning device can be found in Schmitz et al. 24 The reconstruction algorithms were written in MATLAB® and image reconstruction was performed in MATLAB® R2007a running on a standard PC with a single core 3 GHz Intel Pentium microprocessor and 1 GB RAM.

Forward problem
For a medium Ω ⊂ IR 3 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 fyðjÞg j¼1;s from d detectors yðjÞ ¼ ½y 1 ðjÞ; : : : ; y d ðjÞ on ∂Ω, for each point source q j of the array of sources q ¼ ½q 1 ; : : : ; q s distributed on ∂Ω.
The forward problem for the source q j is described by the following parameters-to-output mapping: yðjÞ ¼ P j uðrÞ; j¼ 1; : : : ; s; from the space of parameter functions uðrÞ ¼ ½μ a ðrÞ; μ 0 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. 2 • DðrÞ∇ϕ j ðrÞ þ μ a ðrÞϕ j ðrÞ ¼ q j ðrÞ; r∈ Ω; where ϕ j ðrÞ is the spatially varying diffuse photon density at position r given the source q j , μ a is the absorption coefficient, is the diffusion coefficient, and μ 0 s is the reduced scattering coefficient.The collimated source incident at ξ j ∈ ∂Ω 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 ñðξ j Þ.The boundary condition usually employed is of the Robin type.
where the term A accounts for the refractive index boundary mismatch between the interior and exterior medium. 25or any given source q j , the variable measured by the detector located at ξ i ∈ ∂Ω is the outward flux γ j ðξ i Þ.The corresponding measurement equations are given by γ j ðξÞ ¼ −DðξÞñðξÞ • ∇ϕ j ðξÞ; ξ ∈ Ω; y i ðjÞ ¼ γ j ðξ i Þ; j ¼ 1; : : : ; s; i ¼ 1; : : : ; d: A numerical solution to the DE [Eq.( 2)] can be obtained for arbitrary tissue geometries using a finite element. 25,26Solving 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.
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, μ 0 s ¼ 0.5 mm −1 for skin, μ 0 s ¼ 1.63 mm −1 for skull, μ 0 s ¼ 1.63 mm −1 for brain, and μ 0 s ¼ 1 mm −1 for muscle. 27patial 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 ŷi;j ðtÞ ¼ f i;j ½uðtÞ þ e i;j ðtÞ; j ¼ 1; : : : ; s; i ¼ 1; : : : ; M j ; i ≠ j; (5 where ŷ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Þ; : : : ; 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, 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Þ; : : : ; μ 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 ŷi;j ðtÞ ¼ where p i;j k ðxÞ is an RBF.The RBF chosen in this study is the thin-plate spline.
where c k denotes the center of the RBF and r k ¼ ku − c k k 2 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. 30iven the input-output data set D E i;j ¼ fuðtÞ; y i;j ðtÞg t¼1;N E , the mapping in Eq. ( 5) is derived such that it provides a sparse approximation of the underlying nonlinear function using N ij terms (N i;j ≪ 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 i;j K ≪ 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.

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 ¼ fu i;j l g l∈I i;j , where I i;j is a subset of the full set of nodes f1; : : : ; Kg, with dimension n i;j K ≪ K. 2. Given a dictionary of candidate RBF basis functions, determine a minimal subset of regressors fp k g k¼1;N i;j , N i;j ≪ N D and the associated parameter vector fθ k g k¼1;N i;j , which minimize the prediction error over the data set D i;j .
Step 1.The input layer of the RBF network was selected using a priori information derived from photon measurement density functions (PMDF), 31 which characterize the sensitivity of measurements to changes of the optical parameters inside the medium.For each source-detector pair, the PMDF function was computed based on the FEM approximation over a fine mesh, 32 and only the nodes within a region where values of PMDF are above a given threshold (1 to 5%.) are used to form the input layer of the RBF approximation.For polynomial approximations, 18,33 it can be shown that all the terms selected in the model by alternative algorithms correspond to nodes within the high-sensitivity area. 10The union of nodes selected for all source-detector combinations I ROI ¼ ∪ i;j 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. 5tep 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; : : : ; N E .The minimal set of basis functions in Eq. ( 6) is determined using a greedy orthogonal regression algorithm 34 from an initial set of candidate RBF functions.The algorithm is based on the forward regression principle, which involves adding RBF functions to the model one at the time.At each step, each candidate basis function that is not already in the model is evaluated for inclusion in the model.The candidate regressors are ranked according to their contribution to reducing the variance of the error as measured by the error reduction ratio (ERR) index. 34The selection process is stopped when 100 − 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.

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 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 3) had between 8 and 399 inputs and between 11 and 92 RBFs.

Image reconstruction
The absorption coefficients were obtained using the modified version 19 of the normalized difference method, 35 which involves solving the following optimization problem: y i;j ðtÞ y 0 i;j ŷi;j ð0Þ − f i;j ðu ROI Þ y i;j ðtÞ y 0 i;j ŷi;j ð0Þ 3 7 5

2
; where y i;j ðtÞ is a measurement taken at time t, y 0 i;j is the predefined reference state, ŷi;j ð0Þ ¼ f i;j ðu 0 ROI Þ is the predicted measurement corresponding to the reference medium, and u 0 ROI 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 36 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 (Δ½HbO 2 and Δ½HbR). 19 a ðλÞ ¼ ε HbO 2 ðλÞΔ½HbO 2 þ ε HbR ðλÞΔ½HbR; (9 where ε HbO 2 ðλÞ and ε HbR ðλÞ are the extinction coefficients for oxyhemoglobin and deoxyhemoglobin, respectively, and λ is the photon wavelength. 5he 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: 19 3 Results

Experiment setup
To validate the reconstruction approach, a phantom study was carried out.We used a cylindrical container made of antireflective plastic (diameter 80 mm, height 300 mm, wall thickness 0.10 mm) full of skimmed milk (0.1% fat).A perturbation was added [Fig.2(a)] to the media in the form of another inner plastic cylinder (diameter 20 mm, height 300 mm) full of milk at a different fat concentration (0.81%).The optical properties of the background media 37 at 810 nm are μ a ¼ 0.0024 mm −1 and μ 0 s ¼ 0.5 mm −1 .Similarly, for the perturbation, μ a ¼ 0.024 mm −1 and μ 0 s ¼ 0.4 mm −1 .Previously, milk has been used experimentally as strong scattering media [38][39][40][41] because its optical absorption and scattering parameters are close to those of most living tissue, 42 where 0.001 < μ a < 0.01 mm −1 and 0.5 < μ 0 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.

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.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).Four different TPS-RBF models with C d ¼ 0.05, 0.1, 0.2, and 0.3 were tested.Based on the image correlation coefficient (ICC), no significant variations of image quality were found, and therefore, the model with C d ¼ 0.3 was selected for achieving faster image reconstructions than the other models (less terms are evaluated).Considering that the location of the perturbation for this experiment is known, it is possible to assess the qualitative spatial accuracy of reconstructed images.The ICC is given by 35 ICCðA; where x i A and x i B denote the intensities of the i'th pixel in images A and B, respectively, and x i A and x i B 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 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 ∼80 s; the ROM required only six iterations and took ∼2.5 s.

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. 35The full 3-D FEM model given by Eqs. ( 2) to ( 4) was simulated numerically to generate synthetic optode measurements that were used to reconstruct the dynamics of the absorption parameters.
The 3-D reconstructions of the inclusion for the first time point using the FEM and the ROM (Sec.2.3.5) are displayed in Figs.6(b) and 6(c), respectively.The isosurfaces correspond to absorption thresholds of 10 and 30% of their maximum estimated variation, 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 ∼2 s, while the FEMbased method achieved similar quality after five iterations, but at a considerably larger amount of time (∼1 h).
The FEM approach achieves the best ICC score (a 3% improvement over the ROM approach) after 20 iterations and more than 4 h (ROM speed-up factor is ∼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.

Simulated 3-D Reconstruction of Hemodynamic
Response to Hypercapnia in the Superior Sagittal Sinus of a Rat's Head In this simulation study, a perturbation was specified in the form of an inclusion located in the sagittal sinus of the rat cortex, as shown in Figs.9(a) and 9(b).The shape of the inclusion resembles a section of a blood vessel.We generated changes of absorption coefficients for mesh nodes located inside the inclusion using profiles of hemodynamic responses to hypercapnia [Fig.9(c)] in the sagittal sinus recorded in previous experimental studies 43 using fMRI and optical imaging spectroscopy.
The experiments followed the same protocol as detailed in Sec.2.1.Absorption changes were computed by substituting the concentration profiles for oxyhemoglobin and deoxyhemoglobin, shown in Fig. 9(c), in Eq. ( 9), using extinction coefficients at wavelengths of 760 and 830 nm, for each of the two chromophores.At the baseline, we assumed that the total hemoglobin concentration was 100 μM and the hemoglobin oxygen saturation was 50%.The finite element approximation of the diffusion model [Eqs.( 2) and ( 4)] and the measurement equation [Eq.( 4)], computed for the 3-D mesh, was simulated to generate synthetic measurements of the outward flux at each optode location for each time point.The synthetic measurements were subsequently used in reconstruction.Absorption changes were reconstructed using the complete and approximated models.The number of iterations was limited to 10.The ROM was selected using C d ¼ 1 as the stopping criteria.The FEMbased approach took 80 min to run 10 iterations for each time point.In contrast, the ROM approach took ∼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.
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: Δ½HbO 2 FEM;cal ¼ 5:26Δ½HbO 2 FEM − 0.94 • ROM approach: Δ½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.

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 ∼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 ∼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), 43 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 × 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. 44To visualize the time evolution of hemodynamic changes during the experiment, we positioned the two ROIs [two rectangles in Fig. 11(b)] in the rat barrel cortex to correspond to a superficial (0 to 1 mm) region and a deep (1 to 2 mm) cortical region, which were used in Ref. 43 to measure changes in blood volume following hypercapnia and whisker stimulation.
The changes in oxyhemoglobin (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 ∼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 ¼ f0;60;120;180g 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. 43,44

Discussion
The results demonstrate that our fast image reconstruction approach for DOT, which uses highly optimized maps between the space of optical parameters and the space of measurements, generates 3-D images in a fraction of the time taken by standard FEM-based algorithms, while preserving image quality.
In this work, we used RBFs 17 to construct our reduced-order forward light models.We found this modeling approach to be more powerful and better suited for this modeling task compared to the polynomial approximation technique. 18The 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 ¼ X q i¼1 n i ; where n i ¼ 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, 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 HbO 2 , HbR, and HbT using the Beer-Lambert law and only performed full tomographic reconstruction at specific time points of interest.In contrast, in our study, we performed 3-D reconstructions at each time point and have shown that following the hypercapnic challenge, the reconstructed changes in total hemoglobin, which are believed to be approximately the CBV changes in the barrel cortex, are in good agreement with blood volume changes measured in a previous study 43 using the CBV-MRI technique. 47,48ven on a modestly powered desktop PC and using code written in MATLAB®, which is significantly slower than compiled code, it took ∼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 49 as well as adults [50][51][52] in real time and at much higher resolution 53 than what is achievable using the standard algorithms.In a clinical context, in the absence of subject-specific anatomical images, atlas-based head models registered on the subject's head, such as those proposed by Ferradal, 54 could be used to obtain ROMs that are suitable for real-time processing.

Fig. 1
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.
− ŷi;j ðtÞ 2 N s 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. 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.

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

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

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

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

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

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

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

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

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