1.

Introduction

Breast cancer is the most frequent cancer in women in industrialized countries. More than 200,000 new cases were expected in the United States in 2003.¹ In the process of diagnosis in the United States 1.2 million women a year undergo biopsies, with between 75 and 80 of those tests proving benign.² This means that just in the United States about 900,000 biopsies are unnecessary each year. Therefore, the reduction of such biopsies is an urgent need. It would increase patient comfort and decrease costs.

The number of unnecessary biopsies could be reduced by decreasing the number of equivocal findings of x-ray mammography through enhanced sensitivity and specificity. A possible way to achieve this goal is to combine the actual gold standard method in mammography, x-ray mammography, with alternative diagnostic procedures. Optical mammography could be such an adjunct method. It started in the 1920s when Cutler transilluminated a female breast with visible light in a darkened room.³ Today, near-infrared (NIR) light in the wavelength range between 650 and 1200 nm is used, since it penetrates biological tissue several centimeters. Further advantages of NIR optical imaging are the possible differentiation of soft tissues owing to their different absorption or scatter at NIR wavelengths, and its ability to yield functional information. This is related to the specific absorption of NIR light by natural chromophores as a result of changes in blood oxygenation.⁴ Based on these physical properties, NIR devices have been developed and are available as products.^{5 6}

The clinical use of NIR breast scanners, which were developed in the 1990s, has been limited. Specificity was too low,⁷ which means that too many women with benign tissue alterations had been sent for a diagnostic workup. Interest in optical mammography has increased again with the advent of fluorescent contrast agents.^{8 9} Recently tumor-specific contrast agents (so-called smart contrast agents) have been reported.⁸ These agents consist of self-quenched fluorophore dyes (e.g., Cy5.5) mounted on larger molecules (backbones) via specific peptides. These peptides are recognized and cleaved by specific proteases (e.g., those expressed in tumors), thus releasing the fluorophores, which are now unquenched and free to fluoresce in the NIR wavelength range when activated with light at the excitation range of the fluorophore used.

The possible adjunctive use of NIR fluorescence imaging for detecting breast cancer requires the fusion of geometric information about lesions from NIR measurements with morphological images from other modalities. This data fusion is facilitated by using an NIR mammographylike measurement device. It is described in Sec. 2 of this paper. A mathematical method to obtain geometric information from NIR data is the reconstruction of three-dimensional fluorescence distributions from data measured at the body’s surface.^{10 11 12} In the literature, various reconstruction methods have been reported. However, these algorithms, especially iterative ones, are often very time-consuming.

To avoid these time-consuming reconstructions, we investigated the extraction of geometric information about fluorophore-tagged focal breast lesions from NIR surface data by means of a localization algorithm, namely, the recently presented space-space multiple signal classification (MUSIC) algorithm.¹³ MUSIC was originally invented for the analysis of space-time data in radar technology.¹⁴ Meanwhile it has been successfully transferred to analyze space-time data in biomagnetism^{15 16 17} and functional magnetic resonance,¹⁸ and space-frequency data in electrical imaging of the breast.²⁰

The MUSIC method requires a data model for the “interpretation” of the measured data. This model has to describe the geometric and physical properties of signal generators and the volume containing these generators. In our NIR optical application, the signal generators are focal fluorescence spots, i.e., fluorophore-tagged malignant lesions, at different positions below a measurement array at the body region’s surface. They give rise to different signal behavior owing to different distances of exciting laser sources from the fluorescence spots in the medium. This different signal behavior is mathematically required in order that multiple spots can be localized with the MUSIC algorithm. To obtain a patient-independent model geometry, the breast containing the spots is viewed as an unbounded optical medium. The signal generators have been simulated as pointlike spots. In reality, the spots are spatially extended. The geometric shape of spots can in general be incorporated by a multipole expansion of an agglomeration of elemental pointlike spots. Such an incorporation of the shape of signal generators has been shown to be successful in electrical diagnosis of the breast.²⁰ In principle, MUSIC also allows us to calculate the strength of the signal generators in a second step. This is, however, beyond the scope of this paper. The term “space-space” MUSIC is chosen since the data are measured at different surface positions, owing to different positions of the excitation lasers.

The details of our MUSIC algorithm are described in Sec. 3. The results from simulated data that are presented in Sec. 4 are discussed in Sec. 5. The paper ends with conclusions and an outlook on further work in Sec. 6.

2.

Measurement Geometry

The measurement system that should allow us to record spatially distributed surface data that are due to different laser source positions is a mammographylike device. It is assumed to have two measurement planes between which the breast is fixed and compressed. One of the planes contains laser excitations that activate the fluorescence of the fluorophore-marked tumor from different positions. This approach allows us to construct a rather inexpensive device for which no special optical coupling is needed. In addition, the images of the fluorescent spots can be merged directly with the corresponding x-ray mammography images. Figure 1 shows the geometric setup of optical measurement sensors and excitation lasers for a potential optical mammography device. These could also be integrated into the compression plates of an x-ray mammography system.

Figure 1

Mammographylike measurement and excitation geometry.

3.

Space-Space MUSIC

The problem to be solved is the 3-D localization of fluorescence spots in biological tissue. The mathematical structure of the localization problem is the same as that in estimating the direction of arrival of wavefronts impinging on a sensor array in radar measurements,¹⁴ or in estimating the three-dimensional positions of focal brain activities in magnetoencephalographic measurements,^{15 16 17} or estimating the three-dimensional positions of electrically polarized focal lesions in bioimmittance measurements.²⁰ In these applications, either space-time or space-frequency data are recorded with a sensor array at different time instants or at different excitation frequencies, respectively. Data of these types allow derivation of a signal subspace of the respective data space. At each point within the volume of interest, a distance measure between the signal subspace and the application-specific model subspace is calculated. The locations of minimum distances give rise to peaks in a cost function that are identified as positions of the signal generators. It is straightforward to transfer this localization method to NIR data from a sensor array. Instead of time or frequency, the second physical data parameter is assumed to be the position of the laser source that excites the fluorophore-tagged lesions.

3.2.

Localization with MUSIC

Space-space MUSIC analyzes two-dimensional spatial data, in this case photon intensities J_n(r _m), measured at M surface positions r _m (m=1,…,M). The fluorophore is assumed to be excited by external light sources (i.e., laser sources) at a series of N positions, r _n (n=1,…,N), where these sources are turned on one after the other. The M-dimensional data vectors of photon intensities J_n(r _m) thus obtained:

Eq. (5)

$${\mathbf{J}}_{\mathit{n}}={[{\mathit{J}}_{\mathit{n}}({\mathbf{r}}_1),\mathrm{\dots },{\mathit{J}}_{\mathit{n}}({\mathbf{r}}_{\mathit{M}})]}^{\mathit{T}}$$

are considered as column vectors of an M×N intensity data matrix J:

Eq. (6)

$$\mathbf{J}=({\mathbf{J}}_1,\mathrm{...},{\mathbf{J}}_{\mathit{N}}),$$

Note that for the MUSIC localization we do not need to know the positions of the excitations (since they are not incorporated in the model). The index n is simply used to denote the different datasets.

The signal subspace is obtained from the data by performing a singular-value decomposition (SVD) of J. In tensor notation, the SVD is given by

Eq. (7)

$$\mathbf{J}={\displaystyle \sum _{\mathit{k}=1}^{\min (\mathit{M},\mathit{N})}}{\mathit{S}}_{\mathit{k}}{\mathbf{U}}_{\mathit{k}}\otimes {\mathbf{v}}_{\mathit{K}}^{\mathit{T}},$$

where s_k are the singular values, u _k are M-dimensional (M-D) and v _k are N-D orthonormal singular vectors. The vectors u _k depend only on the sensor indices, whereas the vectors v _k depend only on the indices denoting the positions of the excitation light sources. The vectors u _k can be considered as excitation-independent intensity distributions related to the two measurement arrays. Thus they are regarded as basis maps. To visualize these basis maps, the vectors have to be decomposed into an upper and a lower array part, and then reshaped according to the corresponding sensor array. Thus, two 2-D maps are related to a singular vector u _k.

The number of numerically significant singular values is related to the spots behaving linearly independent with respect to the excitation positions. This behavior is determined by the different distances between the excitation locations and the spots. The basis maps associated with the numerically significant singular values define the signal subspace. These maps show regular structures that are due to the spots in the medium. The residual basis vectors are determined by noise; they define the so-called orthogonal signal subspace. Figure 2 illustrates this for the simulated measurement geometry explained in Sec. 4.1 and two fluorescent tumors with diameters of 6 mm and both 25 mm deep.

Figure 2

Measurement data for a (a) specific excitation and (b) a singular-value spectrum for two tumors and the mammography measurement geometry. (c) and (d) The two basis maps of the signal subspace. (e) Map of the first basis vector of the orthogonal signal subspace. For each map, the contribution of the upper measurement plane (left) and the lower one (right) are shown separately.

The localization of fluorescent spots, which are modeled as pointlike sources, is done by comparing the model subspace defined by Eq. (4) with the signal subspace for different positions of a model spot. The steps of the search procedure are

1. (1) discretization of the search volume,

2. (2) evaluation of a localization function at each grid point, and finally

3. (3) determination of the positions of best fit between the model subspace and the signal subspace from the localization function.

The positions of best fit are interpreted as the centers of gravity of fluorescence spots. Since fluorescence spots give rise to peaks in the intensity distribution measured on the sensor arrays, the search can be restricted to line searches below the peaks of the upper plate, and above the peaks of the lower plate into the respective depth directions, which accelerates the search procedure.

At a search point r the localization function F is chosen as

Eq. (8)

$$F(r)={\left|\left(1-{\displaystyle \sum }_{s=1}^S{u}_s\otimes u_S^T\right)L_{\text{norm}}(r)\right|}^2=|P{L}_{\text{norm}}(r){|}^2,$$

where P is the projector onto the orthogonal signal subspace (where S denotes the number of numerically significant singular values), and

Eq. (9)

$${\mathbf{L}}_{\text{norm}}=\frac{\mathbf{L}}{\left|\mathbf{L}\right|}$$

denotes the normalized basis vector of the 1-D model subspace. The normalization of the model subspace basis vectors compensates for its numerical decrease with increasing depth in order to ensure equal weighting of signal generators from all depths considered.

In the presence of spots, and under the assumption that data errors are tolerable in the model, the localization function exhibits local minima. This indicates the maximum or minimum distance between the model and the orthogonal signal subspace, or the signal subspace, respectively. These local minima are associated with the centers of gravity of the fluorescence spots.

4.

Data

Computer simulations have been carried out to generate data originating from fluorophores in a bounded and optically inhomogeneous medium. This section describes details of the data simulation.

4.1.

Simulation Model

The medium containing fluorescent tagged lesions is assumed to be within a cube of 74×74×56 mm³. Photon diffusion is calculated by numerically solving the stationary fluorescent diffusion equation (see Refs. 10 and 11) using the finite-difference (FD) method. We chose a regular grid of 2 mm in each dimension [giving 37×37×28 cubic (2 mm³) voxels] using a central difference scheme. Robin (or flux) boundary conditions²¹ at the tissue–air (or tissue–device) interface have been implemented.

The optical properties of the medium have been set to μ_a=0.005/mm for the mean absorption coefficient and μ_s ^′=1.1/mm for the mean reduced scattering coefficient. This choice approximates reasonably well the optical properties of female breast tissue.²² The FD simulations allow us to simulate different levels of optical background inhomogeneity. For this purpose, we let the optical properties oscillate sinusoidally around the mean values. The amplitudes of these oscillations are 50 of the mean values. Figure 3 shows the distribution of the optical properties. For our simulations, we assumed a highly specific contrast agent as described in Ref. 8. Thus, we did not introduce background fluorescence.

Figure 3

(a) Inhomogeneous absorption and (b) scattering coefficient distributions incorporated in the simulations.

To simulate a breast compressed between two measurement planes, a planar excitation and measurement array was placed centrally on top of a cube and another detection grid was placed at the bottom. These measurement arrays were assumed to be square grids of 15×15 detectors with 4-mm spacing. The excitation array consists of 5×5 light sources with spacings of 12 mm. The areas of the measurement planes were about 6×6 cm² each. Figure 4 shows the simulated optical medium and the measurement systems, where the two blobs between the measurement planes indicate two fluorescent spots.

Figure 4

Excitation configuration on top of the simulated cube and detectors on the top and bottom.

We further assumed that our excitations are limited to approximately 300 mJ/cm² by the laser safety laws²³ for human applications. Given an input fiber area of 3 mm² (i.e., 2-mm diameter) we assumed a light power of 40 mW and an illumination time in continuous-wave mode of 0.25 s/fiber. With 25 fibers turned on sequentially, this gave a total illumination time of about 6 s. For real experiments, this time might increase slightly as a result of switching times and other possible delays.

5.

Results

In this section we present localization results for the data simulated in Sec. 4. The data model underlying the localization algorithm was described in Sec. 3.1. Again, we would like to stress its simplicity, i.e., its unboundedness and its optical homogeneity via a uniform diffusion parameter in Eq. (3). These properties are especially useful in clinical applications. The mean values of the inhomogeneous background optical properties were assumed to be known by some estimation procedure. Also, localizations where the optical background parameters have been under- or overestimated (90 or 110 of the actual parameters, respectively) have been performed.

Figure 5 summarizes the localization results where the frequency of the normalized localization error ε

Figure 5

Summary of the localization results from 96 configurations: frequency of localization events versus the true depth of the tumor and versus the normalized localization error (see text).

For single tumors, the localization only fails in a single case. A 37-mm deep, 2-mm small tumor was only missed by 2 mm in the case of underestimating the mean optical background properties by 10. In the residual cases, the tumors were localized perfectly well. The normalized localization errors were mostly zero (ε=0).

The localizations for the configurations with two tumors only fail if their distances in depth differ too much (e.g., 13 versus 37 mm). In these cases, the tumor that is farther from the plane containing the excitation sources is mislocalized, whereas the tumor closer to the excitation plane is localized correctly. In all other cases (regardless of whether the optical properties were under- or over- or correctly estimated), we got |ε|⩽1, which means that the localization was correct.

6.

Conclusions

The space-space MUSIC method described and tested for a single measurement plane¹³ was transferred and adapted to the analysis of multisensor optical fluorescence data from a mammographylike measurement geometry, i.e., a system with two measurement planes on top and bottom of the medium under investigation. It has been shown that this application of the MUSIC algorithm has allowed successful localization of simulated fluorescence spots.

The effectiveness and the robustness of the algorithm have been demonstrated with the successful localization of spots from simulated noisy data within a bounded and optically inhomogeneous medium. The localization results are, as expected, better than those obtained from data acquired with a “hand-held” single measurement probe.¹³

The approach presented here for the detection and analysis of fluorescence data has several attractive features. The data model is extremely simple: boundary effects are neglected. The localization results were barely influenced by assumed wrong optical background parameters of the infinite medium. In short, discrepancies between the data model incorporated into the MUSIC algorithm and the data simulation model do not affect the quality of the localization results. Therefore, these results indicate that the proposed space-space MUSIC algorithm is patient independent, which would be extremely useful for clinical applications. Of course, this has to be confirmed in further studies.

Another attractive feature of the MUSIC algorithm is that it yields results in real time, i.e., within subseconds. This performance is due to the simplicity of the underlying data model and the evaluation of the cost function along a search path, which has fewer grid points than the entire discretized volume under study.

The space-space MUSIC algorithm as presented can still be extended, and the results are expected to become even better. In view of clinical data, a possible improvement would come from the inclusion of the shape and orientation of a spot through higher-order multipole sources.²⁰

The localization of signal generators is only the first of two steps of a MUSIC algorithm. In a second step, the localized positions of the signal generators (i.e., of the spots) can be used to determine their signal strengths (i.e., the amount of fluorescing molecules) by solving a set of linear equations relating these signal strengths at the found positions and the light intensities measured. Another field of activity is an improved estimate of the optical background parameters. Such an estimate could be obtained by analyzing reflectance data while recording fluorescence data.

In summary, the localization results presented in this paper indicate that the proposed space-space MUSIC algorithm seems to be a promising method for analyzing NIR fluorescence data.