A spectrally constrained dual-band normalization technique for protoporphyrin IX quantification in fluorescence-guided surgery.

We report a dual-band normalization technique for in vivo quantification of the metabolic biomarker, protoporphyrin IX (PpIX), during brain tumor resection procedures. The accuracy of the approach was optimized in tissue simulating phantoms with varying absorption and scattering properties, validated with fluorimetric assessments on ex vivo brain tissue, and tested on human data acquired in vivo during fluorescence-guided surgery of brain tumors. The results demonstrate that the dual-band normalization technique allows PpIX concentrations to be accurately quantified by correction with reflectance data recorded and integrated within only two narrow wavelength intervals. The simplicity of the method lends itself to the enticing prospect that the method could be applicable to wide-field applications in quantitative fluorescence imaging and dosimetry in photodynamic therapy.

Multiresolution MR elastography using nonlinear inversion.

PURPOSE: Nonlinear inversion (NLI) in MR elastography requires discretization of the displacement field for a finite element (FE) solution of the "forward problem", and discretization of the unknown mechanical property field for the iterative solution of the "inverse problem". The resolution requirements for these two discretizations are different: the forward problem requires sufficient resolution of the displacement FE mesh to ensure convergence, whereas lowering the mechanical property resolution in the inverse problem stabilizes the mechanical property estimates in the presence of measurement noise. Previous NLI implementations use the same FE mesh to support the displacement and property fields, requiring a trade-off between the competing resolution requirements. METHODS: This work implements and evaluates multiresolution FE meshes for NLI elastography, allowing independent discretizations of the displacements and each mechanical property parameter to be estimated. The displacement resolution can then be selected to ensure mesh convergence, and the resolution of the property meshes can be independently manipulated to control the stability of the inversion. RESULTS: Phantom experiments indicate that eight nodes per wavelength (NPW) are sufficient for accurate mechanical property recovery, whereas mechanical property estimation from 50 Hz in vivo brain data stabilizes once the displacement resolution reaches 1.7 mm (approximately 19 NPW). Viscoelastic mechanical property estimates of in vivo brain tissue show that subsampling the loss modulus while holding the storage modulus resolution constant does not substantially alter the storage modulus images. Controlling the ratio of the number of measurements to unknown mechanical properties by subsampling the mechanical property distributions (relative to the data resolution) improves the repeatability of the property estimates, at a cost of modestly decreased spatial resolution. CONCLUSIONS: Multiresolution NLI elastography provides a more flexible framework for mechanical property estimation compared to previous single mesh implementations.

Whole-brain MR-registered cryo-imaging of a porcine-human glioma model to compare contrast agent biodistributions.

As rapidly accelerating technology, fluorescence guided surgery (FGS) has the potential to place molecular information directly into the surgeon's field of view by imaging administered fluorescent contrast agents in real time, circumnavigating pre-operative MR registration challenges with brain deformation. The most successful implementation of FGS is 5-ALA-PpIX guided glioma resection which has been linked to improved patient outcomes. While FGS may offer direct in-field guidance, fluorescent contrast agent distributions are not as familiar to the surgical community as Gd-MRI uptake, and may provide discordant information from previous Gd-MRI guidance. Thus, a method to assess and validate consistency between fluorescence-labeled tumor regions and Gd-enhanced tumor regions could aid in understanding the correlation between optical agent fluorescence and Gd-enhancement. Herein, we present an approach for comparing whole-brain fluorescence biodistributions with Gd-enhancement patterns on a voxel-by-voxel basis using co-registered fluorescent cryo-volumes and Gd-MRI volumes. In this initial study, a porcine-human glioma xenograft model was administered 5-ALA-PpIX, imaged with MRI, and euthanized 22 hours following 5-ALA administration. Following euthanization, the extracted brain was imaged with the cryo-macrotome system. After image processing steps and non-rigid, point-based registration, the fluorescence cryo-volume and Gd-MRI volume were compared for similarity metrics including: image similarity, tumor shape similarity, and classification similarity. This study serves as a proof-of-principle in validating our screening approach for quantitatively comparing 3D biodistributions between optical agents and Gd-based agents.

Subzone based magnetic resonance elastography using a Rayleigh damped material model.

PURPOSE: Recently, the attenuating behavior of soft tissue has been addressed in magnetic resonance elastography by the inclusion of a damping mechanism in the methods used to reconstruct the resulting mechanical property image. To date, this mechanism has been based on a viscoelastic model for material behavior. Rayleigh, or proportional, damping provides a more generalized model for elastic energy attenuation that uses two parameters to characterize contributions proportional to elastic and inertial forces. In the case of time-harmonic vibration, these two parameters lead to both the elastic modulus and the density being complex valued (as opposed to the case of pure viscoelasticity, where only the elastic modulus is complex valued). METHODS: This article presents a description of Rayleigh damping in the time-harmonic case, discussing the differences between this model and the viscoelastic damping models. In addition, the results from a subzone based Rayleigh damped elastography study of gelatin and tofu phantoms are discussed, along with preliminary results from in vivo breast data. RESULTS: Both the phantom and the tissue studies presented here indicate a change in the Rayleigh damping structure, described as Rayleigh composition, between different material types, with tofu and healthy tissue showing lower Rayleigh composition values than gelatin or cancerous tissue. CONCLUSIONS: It is possible that Rayleigh damping elastography and the concomitant Rayleigh composition images provide a mechanism for differentiating tissue structure in addition to measuring elastic stiffness and attenuation. Such information could be valuable in the use of Rayleigh damped magnetic resonance elastography as a diagnostic imaging tool.

Monitoring of hemodynamic changes induced in the healthy breast through inspired gas stimuli with MR-guided diffuse optical imaging.

PURPOSE: The modulation of tissue hemodynamics has important clinical value in medicine for both tumor diagnosis and therapy. As an oncological tool, increasing tissue oxygenation via modulation of inspired gas has been proposed as a method to improve cancer therapy and determine radiation sensitivity. As a radiological tool, inducing changes in tissue total hemoglobin may provide a means to detect and characterize malignant tumors by providing information about tissue vascular function. The ability to change and measure tissue hemoglobin and oxygenation concentrations in the healthy breast during administration of three different types of modulated gas stimuli (oxygen/ carbogen, air/carbogen, and air/oxygen) was investigated. METHODS: Subjects breathed combinations of gases which were modulated in time. MR-guided diffuse optical tomography measured total hemoglobin and oxygen saturation in the breast every 30 s during the 16 min breathing stimulus. Metrics of maximum correlation and phase lag were calculated by cross correlating the measured hemodynamics with the stimulus. These results were compared to an air/air control to determine the hemodynamic changes compared to the baseline physiology. RESULTS: This study demonstrated that a gas stimulus consisting of alternating oxygen/carbogen induced the largest and most robust hemodynamic response in healthy breast parenchyma relative to the changes that occurred during the breathing of room air. This stimulus caused increases in total hemoglobin and oxygen saturation during the carbogen phase of gas inhalation, and decreases during the oxygen phase. These findings are consistent with the theory that oxygen acts as a vasoconstrictor, while carbogen acts as a vasodilator. However, difficulties in inducing a consistent change in tissue hemoglobin and oxygenation were observed because of variability in intersubject physiology, especially during the air/oxygen or air/carbogen modulated breathing protocols. CONCLUSIONS: MR-guided diffuse optical imaging is a unique tool that can measure tissue hemodynamics in the breast during modulated breathing. This technique may have utility in determining the therapeutic potential of pretreatment tissue oxygenation or in investigating vascular function. Future gas modulation studies in the breast should use a combination of oxygen and carbogen as the functional stimulus. Additionally, control measures of subject physiology during air breathing are critical for robust measurements.

Considerations for NIR-I and short-wave infrared (SWIR) fluorescence imaging within a clinical operating room.

Short-wave infrared (SWIR/NIR-II) fluorescence imaging has received increased attention for use in fluorescence-guided surgery (FGS) due to the potential for higher resolution imaging of subsurface structures and reduced autofluorescence compared to conventional NIR-I imaging. As with any fluorescence imaging modality introduced in the operating room, an appropriate accounting of contaminating background signal from other light sources in the operating room is an important step. Herein, we report the background signals in the SWIR and NIR-I emitted from commonly-used equipment in the OR, such as ambient and operating lights, LCD screens and surgical guidance systems. These results can guide implementation of protocols to reduce background signal.

Sensitivity study and optimization of a 3D electric impedance tomography prostate probe.

In current clinical practice, the primary diagnostic method for testing for prostate cancer is ultrasound-guided biopsy. In this paper, we consider using a sonolucent array of electrodes, printed on a thin Kapton layer and positioned on the imaging window of a transrectal ultrasound probe, as a method for providing coregistered electrical and ultrasound imaging of the prostate. As the electrical properties of malignant tissues have been shown to differ significantly from benign tissues, the estimation of the electrical properties is expected to be helpful in distinguishing certain beginning pathologies from cancer and in improving the detection rate that current biopsy methods provide. One of the main difficulties in estimating electrical properties of tissues with this electrode configuration is the rapid decay of the sensitivity with distance from the sensing array. In order to partially overcome this difficulty, we propose to use prior information from the ultrasound (US). Specifically we intend to delineate the boundaries of the prostate from the US, to subdivide the organ into a small number of voxels and to estimate the conductivity as constant on each of these subvolumes. We use a 3D forward model based on the finite element method for studying the sensitivity of a simulated segmented prostate for three different electrode array designs. The three designs present different electrode areas and inter-electrode gaps. Larger electrodes are desirable as they present a better contact, but we show that as they result in smaller inter-electrode gaps, shunting currents can be significant and the sensitivity is reduced. Because our clinical measurement system employs a single current source, we consider tetrapolar measurement patterns for evaluating these electrode configurations. Optimal measurement patterns are well defined for adaptive systems, where multiple currents are injected at the same time. For the electrode array designs we consider, which are three dimensional, there are no established systematic methods for forming sets of linearly independent tetrapolar measurement patterns. We develop a novel method for automatically computing a full set of independent tetrapolar measurement patterns that maximizes the sensitivity in a region of interest (ROI). We use these patterns in the forward modeling and sensitivity studies. In addition to the electrode arrays on the probe, we study the use of a further configuration, where a distal electrode is positioned on the exterior of the body and used for current injection.

Hyperspectral data processing improves PpIX contrast during fluorescence guided surgery of human brain tumors.

Fluorescence guided surgery (FGS) using aminolevulinic-acid (ALA) induced protoporphyrin IX (PpIX) provides intraoperative visual contrast between normal and malignant tissue during resection of high grade gliomas. However, maps of the PpIX biodistribution within the surgical field based on either visual perception or the raw fluorescence emissions can be masked by background signals or distorted by variations in tissue optical properties. This study evaluates the impact of algorithmic processing of hyperspectral imaging acquisitions on the sensitivity and contrast of PpIX maps. Measurements in tissue-simulating phantoms showed that (I) spectral fitting enhanced PpIX sensitivity compared with visible or integrated fluorescence, (II) confidence-filtering automatically determined the lower limit of detection based on the strength of the PpIX spectral signature in the collected emission spectrum (0.014-0.041 µg/ml in phantoms), and (III) optical-property corrected PpIX estimates were more highly correlated with independent probe measurements (r = 0.98) than with spectral fitting alone (r = 0.91) or integrated fluorescence (r = 0.82). Application to in vivo case examples from clinical neurosurgeries revealed changes to the localization and contrast of PpIX maps, making concentrations accessible that were not visually apparent. Adoption of these methods has the potential to maintain sensitive and accurate visualization of PpIX contrast over the course of surgery.

Viscoelastic power law parameters of in vivo human brain estimated by MR elastography.

The noninvasive imaging technique of magnetic resonance elastography (MRE) was used to estimate the power law behavior of the viscoelastic properties of the human brain in vivo. The mechanical properties for four volunteers are investigated using shear waves induced over a frequency range of 10-50Hz to produce a displacement field measured by magnetic resonance motion-encoding gradients. The average storage modulus (µR) reconstructed with non-linear inversion (NLI) increased from approximately 0.95 to 2.58kPa over the 10-50Hz span; the average loss modulus (µI) also increased from 0.29 to 1.25kPa over the range. These increases were modeled by independent power law (PL) relations for µR and µI returning whole brain, group mean exponent values of 0.88 and 1.07 respectively. Investigation of these exponents also showed distinctly different behavior in the region of the cerebral falx compared to other brain structures.

Suitability of poroelastic and viscoelastic mechanical models for high and low frequency MR elastography.

PURPOSE: Descriptions of the structure of brain tissue as a porous cellular matrix support application of a poroelastic (PE) mechanical model which includes both solid and fluid phases. However, the majority of brain magnetic resonance elastography (MRE) studies use a single phase viscoelastic (VE) model to describe brain tissue behavior, in part due to availability of relatively simple direct inversion strategies for mechanical property estimation. A notable exception is low frequency intrinsic actuation MRE, where PE mechanical properties are imaged with a nonlinear inversion algorithm. METHODS: This paper investigates the effect of model choice at each end of the spectrum of in vivo human brain actuation frequencies. Repeat MRE examinations of the brains of healthy volunteers were used to compare image quality and repeatability for each inversion model for both 50 Hz externally produced motion and ≈1 Hz intrinsic motions. Additionally, realistic simulated MRE data were generated with both VE and PE finite element solvers to investigate the effect of inappropriate model choice for ideal VE and PE materials. RESULTS: In vivo, MRE data revealed that VE inversions appear more representative of anatomical structure and quantitatively repeatable for 50 Hz induced motions, whereas PE inversion produces better results at 1 Hz. Reasonable VE approximations of PE materials can be derived by equating the equivalent wave velocities for the two models, provided that the timescale of fluid equilibration is not similar to the period of actuation. An approximation of the equilibration time for human brain reveals that this condition is violated at 1 Hz but not at 50 Hz. Additionally, simulation experiments when using the "wrong" model for the inversion demonstrated reasonable shear modulus reconstructions at 50 Hz, whereas cross-model inversions at 1 Hz were poor quality. Attenuation parameters were sensitive to changes in the forward model at both frequencies, however, no spatial information was recovered because the mechanisms of VE and PE attenuation are different. CONCLUSIONS: VE inversions are simpler with fewer unknown properties and may be sufficient to capture the mechanical behavior of PE brain tissue at higher actuation frequencies. However, accurate modeling of the fluid phase is required to produce useful mechanical property images at the lower frequencies of intrinsic brain motions.

Comparative theoretical performance for two types of regional hyperthermia systems.

Regional hyperthermia systems have drawn attention because of their potential for depositing power noninvasively in deep-seated tumors. Two such systems that have received clinical attention because of their ability to deposit significant amounts of power in tissue are magnetic induction devices and annular phased array applicators. In this paper, theoretical calculations for the specific absorption rate (SAR) and the resulting temperature distributions for these systems are compared. The finite element method is used in the formulation of both the electromagnetic and thermal boundary value problems. Six detailed patient models based on CT-scan data from the pelvic, visceral, and thoracic regions are generated to simulate a variety of tumor locations. In general, the annular phased array deposited more power within the tumor and produced better temperature distributions than the magnetic induction device. However, the ratio of the maximum power absorbed by the tumor to the maximum power absorbed in normal tissue does not appear to be high enough for either device to heat significant portions of perfused tumors to therapeutic temperatures under a wide range of physiological conditions. The results contained herein should aid the physician in comparative treatment planning with existing regional hyperthermia systems.

Optimization of the absorbed power distribution for an annular phased array hyperthermia system.

One of the systems under investigation for producing hyperthermia noninvasively for treating deep-seated tumors is the annular phased array. This device consists of two rings of eight electromagnetic apertures that are placed concentrically about the long axis of the patient and radiate energy toward the center. Previous theoretical and clinical studies have concentrated primarily on systems where the amplitude and phase of the signal applied to each aperture were the same, and these studies have shown that the system is capable of depositing power deep within the patient. Nevertheless, in many situations the system was not capable of producing desirable temperature distributions in the tumor and normal tissue. In this paper we report on a 2-dimensional theoretical investigation where an optimization routine was used to select the amplitude and phases of each of eight apertures. The optimization procedure and resulting calculations were based on CT scans of patients with tumors. The electrical and thermal properties of the different organs and tissues were taken into account. The optimization routine tried to achieve uniform absorbed power in the tumor region with zero absorbed power outside. Using the optimized amplitudes and phases, the SAR (specific absorption rate, W/kg) was calculated for the array. The results show that in general the optimization procedure was successful in that the power deposited within the tumor volume was increased with less power deposited into normal tissue when compared to the equal amplitude and phase case. This SAR data was then used as the input to a program based on the bioheat transfer equation, which calculated the temperature distribution in the patient model for an assumed set of blood perfusion rates. Depending on the location, size of the tumor, and blood perfusion rates, the improvement in the percentage of the tumor brought to therapeutic temperature varied from 0% to as much as 80%.

Theoretical temperature profiles for concentric coil induction heating devices in a two-dimensional, axi-asymmetric, inhomogeneous patient model.

In this paper we report on theoretical calculations for the temperature distributions produced by an rf magnetic induction device that is placed concentrically about the long axis of the patient. A two-dimensional, axi-asymmetric, inhomogeneous patient model was used in conjunction with a numerical moment method for calculating the electric fields in the tissues of the model and a numerical finite element method for calculating the resulting temperature distributions. The electric fields and the absorbed power per unit volume of tissue were calculated for both a thorax and viscera model, each of which included a tumor volume. The absorbed power values were input into the bioheat transfer equation and the temperature distributions were calculated for a wide range of blood flow rates. Based on the steady-state and transient results, our computer simulations predict poor therapeutic temperature profiles for tumors embedded deeply in the thorax and viscera. This heating technique appears to produce significant therapeutic volumes in superficial tumors located not greater than 7 cm in depth. These theoretical calculations should aid the clinician in the evaluation of induction heating devices for their effectiveness in heating deep-seated and superficial tumors.

In vivo electrical impedance spectroscopic monitoring of the progression of radiation-induced tissue injury.

This study evaluates the potential of electrical impedance spectroscopy (EIS) as a noninvasive technique for tracking the progression of radiation-induced damage in normal muscle tissue. Male Sprague-Dawley rats were irradiated locally to the gastrocnemius and biceps femoris muscle. Single doses were administered using a procedure that spares skin and bone. Complex impedance spectral measurements (taken at 50 frequency points between 1 kHz and 1 MHz) were made at monthly intervals using recessed disk electrodes applied to the skin. A histological scoring scheme was developed for evaluation of injury. A strong dose-dependent progression of injury evident in both spectral measurements and histological scoring has been observed. Latent time also appears to be dependent on dose with changes induced by 70 Gy evident by 2 months, changes induced by 90 Gy observed by 1 month, and dramatic changes found within 3 weeks at 150 Gy. Injury was morphologically comparable to the type of damage that occurs in response to small, fractionated doses, but on a much shorter time scale. Increased spectral shift was a consistent indicator of the extent of tissue injury at the time of measurement. The use of a large single dose resulted in an excellent model in terms of inducing a significant progression in tissue injury over a short post-treatment follow-up period in the muscle mass while also providing a consistent location for in vivo electrical impedance measurements. The results show that EIS can follow radiation-induced tissue change, suggesting that EIS has the potential to monitor the types of injury observed in late radiation damage of muscle tissue noninvasively.

Theoretical temperature distributions produced by an annular phased array-type system in CT-based patient models.

Theoretical calculations for the specific absorption rate (SAR) and the resulting temperature distributions produced by an annular phased array (APA)-type system are made. The finite element numerical method is used in the formulation of both the electromagnetic (EM) and thermal boundary value problems. A number of detailed two-dimensional patient models based on CT-scan data from the pelvic, visceral, and thoracic regions are generated to simulate a variety of tumor locations and surrounding normal tissues. The SAR values from the EM solution are put into the bioheat transfer equation, and steady-state temperature distributions are calculated for a wide range of blood flow rates. Based on our theoretical modeling, the APA shows no preferential heating of superficial over deep-seated tumors. However, in most cases for all three regions of the human trunk only fair thermal profiles (therapeutic area near 60%) are obtained in tumors with little or no blood flow and poor temperature patterns (therapeutic area less than 50%) are found in tumors with moderate to high perfusion rates. These theoretical calculations should aid the clinician in the evaluation of the effectiveness of APA-type devices in heating tumors located in the trunk region.

Estimation of oxygen distribution in RIF-1 tumors by diffusion model-based interpretation of pimonidazole hypoxia and eppendorf measurements.

Numerical simulations of oxygen diffusion from the capillaries in tumor tissue were used to predict the capillary oxygen supply within and near hypoxic regions of the RIF-1 tumor. A finite element method to simulate the oxygen distribution from a histology section is presented, along with a method to iteratively estimate capillary oxygen concentrations. Pathological structural data for these simulations came from sections of the tumor stained with hematoxylin and eosin and were used to define the capillary positions and shapes, while overlapping regions of low oxygen concentration were defined by the hypoxia marker pimonidazole. These simulations were used to calculate spatial maps of the oxygen concentration and were tested for their ability to reproduce Eppendorf pO(2) histograms from the same tumor line. This simulation study predicted that capillary oxygen concentrations ranged from zero to above 20 microM, with a dominant peak in the hypoxic regions showing 78% of capillaries with less than 1 microM oxygen concentration, compared to only 12% in the non-hypoxic regions. The results were not highly sensitive to the metabolic oxygen consumption rate, within the range of 2 to 16 microM/s. This numerical method for oxygen capillary simulation is readily adaptable to histology sections and provides a method to examine the heterogeneity of oxygen within the capillaries and throughout the tumor tissue section being examined.

Application of linear circuit models to impedance spectra in irradiated muscle.

We have applied a number of modeling schemes to previously reported in vivo electrical impedance measurements on irradiated and normal muscle in the hind legs of rats. Specifically, seven-parameter parallel pathways and embedded membrane circuit models have been fit to group averages of impedance spectra measured at different doses and time points. Correlations between histologically scored tissue sections and model parameters have also been determined. The results show that both models produce good fits to the experimental observations, especially in the case of the irradiated tissues. The correlations between histology scores and circuit parameters were, however, higher with the embedded model. Trends in the spectra and the model parameters were found to agree with the expected changes in tissue pathophysiology associated with the progression of tissue injury from radiation exposure. Quantitative correlations with specific histological criteria were less conclusive, suggesting that more information may be needed to refine the model architecture if model parameters are to be explicitly related to the types and extent of tissue damage induced by radiation treatments.

Calibration of near-infrared frequency-domain tissue spectroscopy for absolute absorption coefficient quantitation in neonatal head-simulating phantoms.

Frequency-domain tissue spectroscopy is a method to measure the absolute absorption coefficient of bulk tissues, assuming that a representative model can be found to recover the optical properties from measurements. While reliable methods exist to calculate absorption coefficients from source-detector measurements less than a few centimeters apart along a flat tissue volume, it is less obvious what methods can be used for transmittance through the larger tissue volumes typically associated with neonatal cerebral monitoring. In this study we compare the use of multiple distance frequency-domain measurements processed with (i) a modified Beer-Lambert law method, (ii) an analytic infinite-medium diffusion theory expression, and (iii) a numerical finite element solution of the diffusion equation, with the goal of recovering the absolute absorption coefficient of the medium. Based upon our observations, the modified Beer-Lambert method provides accurate absolute changes in the absorption coefficient, while analytic infinite-medium diffusion theory solutions or finite element-based numerical solutions can be used to calculate the absolute absorption coefficient, assuming that the data can be measured at multiple source-detector distances. We recommend that the infinite-medium multi-distance method or the finite element method be used across large tissue regions for calculation of the absolute absorption coefficient using frequency-domain near-infrared measurements at multiple positions along the head.

Spatially varying optical property reconstruction using a finite element diffusion equation approximation.

A finite element reconstruction algorithm for optical data based on a diffusion equation approximation is presented. A frequency domain approach is adopted and a unified formulation for three combinations of boundary observables and conditions is described. A multidetector, multisource measurement and excitation strategy is simulated, which includes a distributed model of the light source that illustrates the flexibility of the methodology to modeling adaptations. Simultaneous reconstruction of both absorption and scattering coefficients for a tissue-like medium is achieved for all three boundary data types. The algorithm is found to be computationally practical, and can be implemented without major difficulties in a workstation computing environment. Results using simulated data suggest that qualitative images can be produced that readily highlight the location of absorption and scattering heterogeneities within a circular background region of close to 4 cm in diameter over a range of contrast levels. Absorption images appear to more closely identify the true size of the heterogeneity; however, both the absorption and scattering reconstructions have difficulty with sharp transitions at increasing depth. Quantitatively, the reconstructions are not accurate, suggesting that absolute optical imaging involving simultaneous recovery of both absorption and scattering profiles in multicentimeter tissues geometries may prove to be extremely difficult.

Frequency-domain near-infrared photo diffusion imaging: initial evaluation in multitarget tissuelike phantoms.

In this paper, an initial evaluation of our finite element based frequency-domain image reconstruction algorithm is performed for experiments where multiple millimeter-sized heterogeneities are embedded within a tissue-equivalent (optically) background medium having multicentimeter dimensions. The cases considered consist of several interesting geometry and optical property contrast combinations including (i) two different-sized targets with the same contrast at three different separation distances; (ii) two different-sized targets with different contrasts at two different separation distances; and (iii) three targets with the same and different sizes and contrasts, respectively. The reconstruction algorithm that has been used is an enhanced version of our originally developed regularized least squares approach that now includes total variation minimization, dual meshing, and spatial low-pass filtering. Quantitative measures of image quality including the size, location, and shape of the embedded heterogeneities along with errors in their recovered optical property values are presented. The results show that multiple targets can be clearly detected for all combinations of locations, sizes, and contrast levels considered, but the quantitative nature of this detection is influenced by these parameters.