|Home | About | Journals | Submit | Contact Us | Français|
Image-guided near infrared spectroscopy (IG-NIRS) can provide high-resolution vascular, metabolic and molecular characterization of localized tissue volumes in-vivo. The approach for IG-NIRS uses hybrid systems where the spatial anatomical structure of tissue obtained from standard imaging modalities (such as MRI) is combined with tissue information from diffuse optical imaging spectroscopy. There is need to optimize these hybrid systems for large-scale clinical trials anticipated in the near future in order to evaluate the feasibility of this technology across a larger population. However, existing computational methods such as the finite element method mesh arbitrary image volumes, which inhibit automation, especially with large numbers of datasets. Circumventing this issue, a boundary element method (BEM) for IG-NIRS systems in 3–D is presented here using only surface rendering and discretization. The process of surface creation and meshing is faster, more reliable, and is easily generated automatically as compared to full volume meshing. The proposed method has been implemented here for multi-spectral non-invasive characterization of tissue. In phantom experiments, 3–D spectral BEM-based spectroscopy recovered the oxygen dissociation curve with mean error of 6.6% and tracked variation in total hemoglobin linearly.
The field of diffuse optical tomography (DOT) has evolved significantly over the past two decades as a method for probing biological tissue in-vivo to gain a better understanding of vascular and metabolic progressions in disease. DOT has undergone several preliminary clinical trials [1–3] by demonstrating the ability to study physiology of tissue through unique mechanisms not available through existing modalities. These include molecular, structural and vascular characterization of tissue using intrinsic contrast processes. Tumor to background contrasts are typically 100% or more [4,5] and use of extrinsic contrasts has further increased the biological and molecular contrasts available in-vivo [6,7]. A major focus of on-going studies has been to increase the spatial resolution and create complementary imaging systems by combining different modalities, such as MRI-DOT , X Ray-DOT , and Ultrasound-DOT . These hybrid systems have bypassed the resolution limit of DOT (3–5 mm) in favor of the sub-millimeter resolution possible from MRI, X-Ray and Ultrasound by utilizing a-priori structural information to guide the DOT image recovery. This approach merges the strengths of structural information provided by a high resolution modality with the strengths of functional content provided by optical spectroscopy Researchers are now at the threshold of larger-scale multi-institutional clinical trials, especially in radiology (for diagnosis of breast cancer ) and oncology (for neo-adjuvant treatment monitoring ) across patient populations and geographies. As hybrid DOT systems make their way into the clinic through large-scale trials, algorithms behind spectral image recovery must be optimized for imaging and processing large datasets. This is complicated by the challenge of optical imaging in 3D, which is mandatory for clinicians who are accustomed to viewing volumetric images from the existing standard imaging modalities. It is also well-known that imaging in 3–D is more accurate than 2–D for DOT systems . As such, this study focuses on providing a solution to the proposed problem of image-guided spectroscopy of tissue.
To date, most 3D optical imaging studies have been limited to stand-alone DOT systems. 3–D image reconstruction was employed by Enfield et al  in time-resolved optical mammography without priors. In this work, the authors used a hemispherical fluid-filled cup as the breast interface, creating hemispherical meshes for computation on clinical data.
Eppstein et al  showed the feasibility of 3–D optical imaging by using simulated data from 3–D slab-shaped volumes which were approximated by a domain decomposition technique that used subdomains. Srinivasan et al  used a data-subset algorithm to minimize the computational burden by dynamically choosing optimal data-sets which carried information regarding major features in the domain. Their work presented simulations and experiments also in slab shapes. Dehghani et al used an underdetermined version of the Moore-Penrose generalized inverse to compute 3-D DOT estimates for a clinical subject using a conical shaped mesh. All of these studies have been for stand-alone DOT systems and using regular shapes (slab or cylinder) with numerical volume-based methods such as finite element or finite difference.
Hybrid systems combining MRI,×ray or ultrasound with NIR systems provide superior spatial resolution at millimeter scales. Several algorithms have been developed to incorporate the anatomical tissue structure from an alternate modality into the NIR image recovery. These studies have predominantly been in 2–D exploring Bayesian  and regularized methods [18,19] to improve image resolution. In 3-D, Dehghani et al  used exact knowledge of boundaries and assumption of homogeneous piece-wise regions in a cylindrical shape in 3–D and showed nearly 100% recovery in experiments. Further work by Yalavarthy et al incorporated structural information in a generalized framework for optical imaging [21,22]. Azar et al  have demonstrated co-registration of non-concurrent NIR and MRI images and segmentation of MR images to enable priors to be applied in 3-D and have applied this to phantom models and three clinical cases. None of these have obtained the 3–D imaging volume itself from MRI, instead using assumptions about the shape based on a 2–D slice in the plane of the optical detectors to reconstruct. Hence these do not address the issues of volumetric meshing of complex arbitrary shapes that arise from brain and breast MRI segmentation that is a necessary step to compute on MRI-derived volumes.
In our experience, one of the key challenges faced by hybrid systems towards routine clinical use lies is in image-segmentation and volume rendering of complex arbitrary shapes. Volume meshing is not a simple task. It is a complex procedure and can involve several manual interventions on a case-by-case basis, which severely limit the clinical workflow. It is practically impossible, using current resources, to obtain volume breast meshes with interior boundaries preserved, in a regular manner. With a pressing need to optimize IG-NIRS for large-scale clinical trials, our motivation here is to evaluate the potential of an alternate numerical method that can circumvent the issues of volume meshing and enable processing of large data-sets in the clinic.
Our previous work in combining MRI priors into DOT recovery has indicated that the use of hard priors in image reconstruction provides the most accurate results . Hard priors impose the tissue boundaries obtained from MRI into the image recovery process with the assumption of piece-wise constant tissue types and can be described as Image-guided NIR spectroscopy (IG-NIRS) allowing region-based localized and non-invasive characterization of tissue types. This use of hard priors allows us to investigate alternate numerical methods for IG-NIRS that may be more efficient. One such method is the boundary element method (BEM)  which presents an attractive alternate to the existing finite-element method by using only surface description to model light propagation in tissue. It is well known that surfaces can be generated easily and reliably as compared to volume rendering. Most commercial software packages automatically create surfaces with minimal user intervention.
The flowchart in Figure 1 shows an illustration of the basic steps involved in obtaining surface and volume meshes for computation. Both require segmentation of the MR images followed by surface discretization. BEM can be used directly following the surface tessellation whereas FEM requires additional steps to build a volume mesh suitably tagged with boundaries delineating the volume into different tissue types. This has implication in the clinical workflow especially when processing large number of data sets. A drawback in BEM is the use of full matrices as compared to sparse matrices in FEM. However, the faster mesh generation and lower mesh sizes allow for faster computational speed.
BEM was applied in simulations using 3–D MRI segmented head models by Zacharopoulos et al  to reconstruct for shape and optical properties. In the work here, we have assumed the shape of heterogeneities to be known from MRI and reconstructed directly for chromophore concentrations (oxy and deoxy-hemoglobin and water) and scattering parameters. Results from phantom experiments show accurate monitoring of oxygenation as described by the oxygen dissociation curve, and a linear increase was observed with variation in total hemoglobin.
The phantom experiments were carried out in a stand-alone DOT system developed at Dartmouth, described in detail in Mcbride et al . Briefly, the automated NIR tomography device was constructed to perform cross-sectional imaging of the breast. Optical fibers arranged in a circular array (with both radial and vertical degrees of freedom) have the capability to collect data from sets of 16 source and 15 detector locations per source. Intensity modulated light (at 100MHz) at six discrete wavelengths between 660 and 850 nm is detected by photomultiplier tubes. The detected light is electrically mixed with a reference signal to yield a low-frequency (500 Hz) signal from which measurements of amplitude and phase shift were recorded.
The measured data from the tissue/phantom was first calibrated to account for small offsets due to source–detector fiber transmission, alignment characteristics, and errors in discretization or model–data mismatch. This was accomplished using a homogeneous fitting algorithm  which utilized Newton Raphson method for fitting of two parameters: the slope of the phase (log(θ)) with respect to distance from the source location, and the slope of the logarithm of intensity times distance with respect to distance (log(rI)). This fitting procedure used two solutions: (1) analytical solution for infinite medium, followed by (2) numerical solution to the diffusion equation on the homogeneous domain of relevant size. By using the homogeneous fitting algorithm on both reference phantom measurements and unknown tissue measurements, the calibrated measurements and initial estimates for optical properties were obtained. The initial estimates of the optical properties were converted into concentrations of oxyhemoglobin, deoxyhemoglobin and water and scatter using spectral fitting and used as the starting point in the NIRS reconstruction.
Reconstruction in DOT comprises of the forward and the inverse problem. The forward problem is defined as follows: Given the tissue parameter distribution [µ], to find the light fluence distribution [Φ] within the tissue domain Ω. The inverse problem is defined as: Given a distribution of sources [q] and a distribution of measurements [M] giving a measure of [Φ] at the boundary dΩ, finding the tissue parameter distribution [μ] within domain Ω. For the forward problem, we used the frequency domain diffusion equation to model light propagation in tissue. In the diffusion approximation, the assumption is that the diffuse intensity encounters many particles and is scattered almost uniformly in all directions, and therefore its angular distribution is almost uniform isotropic. This equation exists under the assumption that scatter dominates over absorption which is true in the case of several tissue types, including the human breast, in the wavelength region of 650–1350 nm . This differential equation is written as[29,31]:
where Φ(r,ω) is the isotropic fluence at modulation frequency ω and position r, κ (r) is the diffusion coefficient, µa (r) is the absorption coefficient, c is the speed of light in the medium and q0(r,ω) is an isotropic source. The diffusion coefficient can be written as
where µs is the reduced scattering coefficient. Under the assumptions that the internal surfaces of the various tissue types in the domain are known a priori, and the optical properties are piece-wise constant in each tissue type, the equation can be recast into a modified Helmholtz equation given by:
The fundamental Green’s function solution is readily available for this equation and can be used to derive the BEM formulation. The details of this formulation can be found elsewhere . The forward solver was implemented in Matlab using Fortran code for building the internal matrices and parallelized using Matlab Distributive computing toolbox. The computational time for a two region forward problem with a mesh size of 3015 nodes using 4 processors was 51 seconds.
The inverse problem was solved using a gradient-based iterative procedure based on Newton’s method which minimizes the least-squares functional:
where M is the total number of measurements at each wavelength, and are the measured and calculated fluence respectively, at the boundary for each measurement point j. The update in the optical properties µ is given by,
where Φ refers to the change in boundary data, is the Jacobian, the matrix containing the sensitivity of the boundary data to a change in optical property µa and diffusion coefficient κ given by and α is the regularization parameter chosen empirically based on phantom experiments. The Jacobian was calculated using a perturbation approach to approximate the required derivatives by perturbing µa and κ in each region in turn and calculating the resulting change in the boundary measurements.
By coupling the measurements obtained at all the wavelengths together and applying spectral relationships defining absorption characteristics of oxyhemoglobin, deoxyhemoglobin and water as well as scattering relationships defined empirically, into the reconstruction, it is possible to obtain the NIRS parameters directly . This use of spectral constraints is known to improve accuracy of recovered spectroscopic parameters and reduce cross-talk [35–37]. We implemented this into the BEM-based reconstruction so that the update in the chromophore concentrations and scatter could be recovered directly. The structure of the Jacobian is changed due to the addition of spectral features and is detailed in Srinivasan et al . The spectral BEM recovery, implemented for the first time in DOT, was used to obtain IG-NIRS estimates from phantoms and human subjects.
The phantom experiments presented in this paper were carried out using liquid homogeneous phantoms (for obtaining oxygen dissociation curve) and gelatin heterogeneous phantoms (for characterizing changes in total hemoglobin). For the first experiment, a solution with 1% Intralipid and 1% whole blood was placed in a plastic cup of 70 mm diameter. The pO2 was measured independently using a chemical microelectrode, after calibration of the electrode overnight in saline solution. The 1% whole blood was found to have 18µM hemoglobin, and the oxygenation of the solution was reduced by varying the pO2 values from 150mm to 0mm Hg by addition of yeast. Using the spectral BEM toolbox, we reconstructed values of oxygen saturation, total hemoglobin, water concentration and scatter.
For the second experiment, a phantom using gelatin with whole blood added for absorption and Titanium dioxide for scatter, was made. This phantom of diameter 82 mm had a 2.5cm hole drilled 1cm from the boundary which was filled with a saline solution containing blood with 0.75% Intralipid for scattering. The inclusion’s blood concentration was varied systematically from 20 to 44µM and measurements were taken. The dataset was used to reconstruct spectral BEM estimates of the inclusion and background.
A surface mesh was created for the tissue cup phantom comprising of 2595 nodes from 5186 triangular elements using Netgen . The BEM model was used to obtain homogeneous estimates of NIR optical properties using the homogeneous fitting procedure described in Section 2.2 using Newton-Raphson method. This minimization procedure successfully updated the optical properties which were then used with spectral fitting to obtain homogeneous values for total hemoglobin, oxygen saturation, water content and scatter parameters. The results are shown in Figure 2 for the NIRS estimates. The recovered oxygen saturation followed the expected values from literature with a mean error of 6.6%. The total hemoglobin and water stayed constant with the variation in pO2 as expected with some crosstalk between them at low pO2 values (below 8mm Hg). The average for recovered total hemoglobin was 17µM with a standard deviation of 0.9µM compared to expected value of 18µM measured for the whole blood using a clinical co-oximeter. The water recovered average was 98% with standard deviation of 3% comparable to expected 100% for liquid solution. The reduced scattering coefficient at 785nm stayed constant with average of 1.29 mm−1 and standard deviation of 0.03 mm-1.
A surface mesh comprising of two surfaces (outer and inclusion) for the gelatin phantom was created using Netgen with mesh size of 2154 nodes and 4300 triangles. The datasets obtained with varying hemoglobin were calibrated using the homogeneous fitting algorithm and used to reconstruct for the NIRS parameters using the 3–D BEM spectral reconstruction for the two region problem. The recovered values in the inclusion and the background are shown in Figure 3(a). As expected, the total hemoglobin in the inclusion increased linearly with change in concentration while the background stayed fixed. The other NIRS parameters (oxygen saturation, water and scattering) stayed constant with this variation in hemoglobin as expected, as shown in Figure 3(b).
A key advantage to using BEM is the ease in meshing since only surface discretization is necessary. Volume meshing is a demanding task that takes up to 10–100 times the computational effort required to create surface meshes. This has implication in the clinical work-flow especially when processing large number of data sets. As IG-NIRS systems move into clinical trials, optimization for large scale trials become especially important and the use of BEM facilitates this need. The use of spectral constraints ensures accuracy in the recovered chromophores and minimal cross talk between the parameters.
We have validated this approach in phantom experiments by showing linear tracking of changes in total hemoglobin and recovering the hill-curve response of oxygen saturation with changes in partial pressure of oxygen. The next stage of demonstration of this approach will likely involve processing a complete cohort of data to analyze the biological parameter, rather than the computational parameters, as was the focus here.
This work was supported by NIH grants PO1CA80139 and R01EB007966-01A1.