In-vivo molecular imaging is a rapidly growing research area that provides a translational bridge from in-vitro research to clinical applications, and serves as an important tool for drug discovery and development.1, 2, 3, 4 Despite the limited penetration depth in tissue, optical methods are a powerful technique for labeling and tracking cells and/or disease in tissues, and offer advantages in throughput, ease of use, temporal resolution, and cost effectiveness over other imaging modalities like magnetic resonance imaging, thus finding increased use on small animals. Bioluminescent and fluorescent proteins are widely used for tracking cell growth and monitory gene expression in small animals. A broad range of targeted and activatable fluorescent probes are being developed for use in multiple disease areas, showing promise for moving into clinical settings. Fluorescent proteins, dyes, and nanoparticles are all being engineered toward emission in the red to near-infrared wavelength range where tissue absorption is lower and endogenous tissue autofluorescence is reduced. Distinct spectral responses and a large number of wavelength choices allow multiplexing of several probes.
These recent developments are driving imaging technology toward multispectral fluorescence detection capability, along with a concurrent evolution of software tools for separating, or unmixing, multiple probes based on spectral signature. Both bioluminescent and fluorescent reporters can be unmixed based on spectral features.5, 6 Fluorescence imaging has unique challenges due to the presence of tissue autofluorescence generated by excitation light. Spectral unmixing can improve contrast and sensitivity by separating a fluorescent probe of interest from background autofluorescence.7, 8 By scanning the excitation or emission filters, images collected at different wavelengths can be processed to separate multiple fluorophores, including the autofluorescence component, based on their distinct spectral signatures.
Spectral unmixing is a well established technique that is utilized in many fields such as remote sensing,9 chemometrics,10, 11 microarray technology,12 etc. It has also been applied to fluorescence microscopy, which uses acquisition systems similar to in-vivo fluorescence imaging but works on a microscopic scale.
In general, techniques available for fluorescence microscopy can be divided into two categories:13 1. classification-based algorithms, including supervised classification analysis, unsupervised clustering analysis, and principal component analysis;14 and 2. linear unmixing.5, 15, 16, 17, 18, 19 Classification-based methods assume each pixel only contains one label and thus are not suitable for the quantitative analysis of pixels with possibly colocalized labels. Linear unmixing is the most commonly used method to analyze the image with mixed contributions from colocalized labels. A fundamental assumption for linear unmixing is that the spectroscopic measurement in a pixel is the linear combination of individual fluorophores’ pure spectra scaled by the concentration in that pixel plus some residue , i.e., .
Consider a multispectral imaging problem that can be discretized to a matrix formis the spectral data matrix with rows corresponding to pixels in the image and columns corresponding to different wavelengths (i.e., bandpass filters); is the by concentration profile where each row represents the composition of the mixed fluorophores at the corresponding pixel; is the by spectra matrix where each row is the pure spectrum of the corresponding fluorophore measured at these wavelengths; and is the by matrix of residual not explained by the model, and ideally should be close to the experimental error. Equation 1 is often recognized as a bilinear model, where is a linear function of both and . To accurately solve for the concentration profile, several criteria have to be met. First, fluorophores considered in the model have to account for all the significant signals present in the scene. Secondly, the number of wavelengths should at least equal the number of fluorophores required to resolve, i.e., . Lastly, the pure spectra of those fluorophores should be a-priori determined, either measured individually or interactively picked from a pure pixel with signal from only one fluorophore. Finally, can be determined by solving a least-square problem that minimizes:
While this approach works well for microscopy, is not ideal for in-vivo fluorescence imaging due to the inherent differences in the two imaging techniques. In fluorescence microscopy, samples are diluted in a thin transparent solution or collected directly from a very thin focus plane; in-vivo imaging measures the diffuse radiance from the surface of scattering animal tissue, and thus a fluorophore might make a significant contribution to a cloud of surrounding pixels instead of only the pixel geometrically associated with it. Thus, the quantity in Eq. 1 has to be considered as the signal distribution of fluorophores instead of the concentrations of fluorophores in the pixel. Furthermore, photons from a fluorophore located inside the animal may be heavily attenuated before reaching the surface due to light scattering and absorption. Such attenuation varies with depth, tissue type, and wavelength, so the emitted fluorescence spectrum can be considerably shifted relative to the in-vitro spectrum. Exact modeling of the spatial distribution of fluorophores and the spectral response of light in tissue requires a rigorous 3-D diffusion model, but our purpose in this work is to develop a significantly improved spectral unmixing relative to using the measured in-vitro spectra based on a first-order 2-D approximation to the complex 3-D problem. Such an algorithm should not rely on a-priori knowledge of in-vitro spectra, but ideally would automatically solve for the in-vivo spectral components as well as the signal distribution for each fluorophore.
Spectral unmixing methods have been successfully used for in-vivo fluorescence imaging6 by Mansfield using liquid crystal tunable filters, but details of the unmixing algorithm are unpublished to the best of our knowledge. Recent progress in chemometrics suggests that a method called the multivariate curve resolution alternating least-square method (MCR-ALS) would be suitable for our purpose. It was developed to extract simultaneously pure concentration and spectral profiles from spectroscopic data gathered during chemical process monitoring,10, 11 and has been successfully applied to other applications that are based on the same bilinear model described in Eq. 1, such as microarray technology12 and imaging spectroscopy.20
We present our spectral unmixing algorithm for in-vivo imaging based on MCR-ALS, including the general MCR-ALS framework and the necessary adaptation to fit the in-vivo imaging application. The image acquisition system used for the study is described, which operates in both reflectance and transillumination geometries. Experimental results are presented to examine the performance of the algorithm. Data were collected on a mouse phantom and live mouse models with varying depth and concentration of fluorescent labels such as dyes and nanoparticles.
Methods and Materials
Multispectral Imaging System
The images in this study were acquired using an IVIS® Spectrum Imaging System (Caliper Life Sciences/Xenogen, Alameda, California) as shown in Fig. 1 . It supports in-vivo spectral imaging for both fluorescent and bioluminescent probes. This instrument uses a , -wide, high-sensitivity back-thinned charge coupled device (CCD) camera thermoelectrically cooled to to minimize the dark current and the associated thermal noise. The CCD has high quantum efficiency in the range of and thus covers the visible (VIS) and near-infrared (NIR) range of interest of optical imaging. Light collection is accomplished with an f/1 lens with zoom capability. The minimum detectable radiance for the system is specified at less than , enabling detection of weak signals from deep tissue.
For fluorescence imaging, the instrument uses a quartz tungsten halogen lamp with broadband emission from . The lamp output is delivered to a 12-position excitation filter wheel using a fiber optic bundle. The excitation wheel includes ten bandpass filters spaced every from . The emission light is focused on the CCD through the lens system with an integrated 24-position emission filter wheel. The emission wheel is installed with 18 bandpass filters spaced every from . The spacing and bandwidth of excitation and emission filters are chosen at a medium level to achieve a reasonable spectral resolvability without sacrificing signal per channel and spectral range. More than one hundred excitation and emission filter combinations are available in the VIS and NIR range for spectral unmixing.
For fluorescence imaging, the instrument can switch between reflectance and transillumination modes by rotating a mirror, as illustrated in Fig. 1. Reflectance mode refers to the conventional modality (epifluorescence mode) that the excitation light is projected to the top surface of the object, where the emitted fluorescence light is collected. Reflectance mode favors the reporters near the surface, since the relative strong autofluorescence background generated from the top surface can overwhelm the fluorescent signals of the deeper reporters. In transillumination mode, the excitation light is delivered to an translation assembly under the stage and focused to the bottom side of the object. An anodized aluminum plate with a grid of holes is used to hold the object and guide the light through. The autofluorescence contribution is often reduced, since most autofluorescence is generated from the bottom surface and is absorbed as it travels to the top surface. When imaging a fluorophore located at the center of the subject, illuminating from the bottom would create the same amount of fluorophore signal to the detector as illuminating from the top, generating much less autofluorescence contribution. Therefore, the transillumination mode offers improved access to the deep fluorophores during in-vivo imaging.
To provide a consistent measurement against different imaging conditions so that signals are traceable over time and on different instruments, the imaging system is absolutely calibrated for spectral radiance traceable to the National Institute of Standards and Technology (NIST).7, 21 This calibration provides the conversion of CCD electron counts to isotropic radiance on the object surface by taking into account losses through the optics and aperture (f-stop), field of view (FOV), and camera settings. For epifluorescence imaging particularly, the power and distribution of the excitation light are recorded so that measurements can be normalized to the incident excitation intensity profile to report how efficiently the target yields the fluorescence. Instrument calibration is critical to the quantification of images with different imaging conditions, as is typically encountered during longitudinal studies.
Multivariate Curve Resolution Algorithm
General Multivariate Curve Resolution Alternating Least-Square Scheme
The basic framework of MCR-ALS is shown in Fig. 2 . Initially the number of spectral components (number of fluorophores) is specified to determine the size of the bilinear problem, and an initial estimate for either or in Eq. 1 is made. We then employ an iterative method for solving two alternating least-square problems, i.e., minimization of Eq. 2 over for given , and minimization of Eq. 2 over for given . To help converge to a unique and realistic solution, additional constraints can be added during iterations. These constraints come from knowledge about the spectra such as nonnegativity or low and high pass filters to force the spectra to zero at certain wavelengths. Constraints reduce the solution space by enforcing some boundaries and help solve the rotation and scaling ambiguity inherent in the bilinear model. Here, the rotation and scaling ambiguity refers to the fact that any arbitrary by orthogonal matrix ( , where is identity matrix) can result in another suitable solution and for the bilinear model. A special case of such an orthogonal matrix is a diagonal matrix often known as the intensity ambiguity. Usually this cannot be solved by setting the constraints, and thus a normalization procedure should be enforced on either or throughout the iteration. The convergence is achieved when the absolute change of the residual norm is below a threshold (typically 0.1%). After convergence, both and are determined.
In some cases, the number of spectral components is not known. To determine the number of spectral components, one approach is based on the principal component analysis (PCA), i.e., examining the variance of data explained by a selected number of principal components.11, 22 Since principal components are sorted in terms of the variance they explain, when an additional principal component only affects a small margin of the total explained variance, the selected principal components have accounted for the majority signals and the rest mainly contribute to the random noise. Principal components are orthogonal to each other, and though they have no biological meanings, they imply the independent components present in the data.
The selection of initial estimate may influence the final solution, mainly because there are multiple local minimums in the solution space for problems with many degrees of freedom.12 There are several approaches to choosing the initial estimate. Using uniformly distributed positive random numbers or random samples of the raw data to build is the simple and less biased method, but the resulting estimate is either too far away from the desired solution or numerically unstable, and therefore provides inconsistent solutions.12 Some systematic methods such as evolving factor analysis (EFA)11 and simple-to-use self-modeling mixture analysis(SIMPLISMA)20 have been developed to extract the pure component spectra, but they are generally more suitable for data containing the pure state of each component than the highly mixed in-vivo fluorescence data.
To make the best initial estimate, our approach is based on statistical data analysis using the -means clustering method. The idea is to partition the dataset into groups on some measure of distance between samples, and then use the centroids of these groups as the initial estimate. An important step is to select a distance measure, which will determine how the similarity of two samples is calculated. The correlation distance measure is chosen because it emphasizes the shape of the spectrum and is not affected by the scaling due to the concentration. Alternatively, the initial estimate can be loaded a priori from a spectral library or selected spectra from a region of interest (ROI) on the image. This option gives more control of the unmixing algorithm, especially for in-vitro imaging where pure components are relatively easy to determine.
Constraints during the ALS are critical for MCR-ALS to generate a meaningful solution. In chemometrics, several types of constraints have been introduced,11, 23 and they are non-negativity (both concentrations and spectra cannot have negative values), unimodality (concentration of temporal profiles or pure spectra only have one peak/maxima), and equality (the possibility to fix known values in the concentration profiles or in the spectra during optimization). These constraints need to be adapted and translated to use in in-vivo fluorescence imaging.
Non-negativity constraints are very straightforward and are applied after each calculation for or using least-squares method or . Instead of setting negative elements to zero, a more rigorous method of non-negative least squares (NNLS)24 is used to recalculate the solution at the columns of or the rows of that comprise negative elements rather than the whole matrix. Unimodality constraints are suitable for fluorescence spectra since fluorophores usually only have single excitation or emission peaks but are not applicable to the signal distribution matrix (concentration profiles in chemometrics). To ensure only one peak in the spectrum during the ALS optimization, the secondary peak is either cut vertically or horizontally,25 or corrected by averaging similarly to unimodal least-squares algorithms.26
Equality constraints refer to the linear constraints that can be added to the least-square problems based on the information about the spectra or concentration profiles. The full matrices discussed earlier in Eq. 2 are actually performed one wavelength at a time or one pixel at a time. Incorporating the equality constraints, the ALS problem for spectral unmixing can be rewritten as two constrained least-squares problems:represents the ’th wavelength, represents the ’th pixel, and are the ’th column and the ’th row of the spectral data matrix , is the ’th column of the spectral matrix , is the ’th row of the concentration matrix , is the constraint matrix that selects the constrained components in , is the vector that contains the fixed spectra for the selected components at the ’th wavelength, is the constraint matrix that selects the constrained components in , and is the vector that contains the fixed concentration for the selected components at the ’th pixel. Equality constraints can be used to improve the wellness of least-squares problem and help reduce the cross talk between components. In chemometrics, equality constraints are used to incorporate information such as closure (the sum of the concentration is constant), zero concentration (set the concentration to zero for species known to be absent in a mixture), and selectivity (specify a portion of the data that only contains contribution from a single component).27 For in-vivo fluorescence imaging, similar information can be incorporated as equality constraints as well. In some situations when the probe being used is known, a portion of its spectrum may be predefined. For example, the emission spectrum of most fluorophores has a sharp increase from the shorter wavelength range to its emission peak, and then decreases relatively slowly toward longer wavelengths. Such characteristics are observed for both in-vitro and in-vivo imaging. If the rising edge of one fluorophore is known, it is relatively safe to set its spectra to zero at wavelengths lower than this edge. We call this equality constraint a high pass filter because it is similar to adding a high pass filter to filter out the unwanted spectrum. Similarly, the low pass filter can be used to define the falling edge to eliminate the unwanted spectrum at the longer wavelengths. High pass and low pass constraints are effective to separate colocalized fluorophores that have different emission peaks. Similar to the concepts of zero concentration and selectivity, pure mask and ROI mask can be used to define equality constraints for the distribution maps. Pure mask assigned to a component determines a region in the image that only contains this component, in other words, the contribution of other components in this mask is zero. ROI mask associated with a component determines a region in the image where this component might exist, in other words, outside of this mask the contribution of this component is zero.
An underlying question rarely addressed is how to apply the variety of constraints to the intended component. For example, assuming there are two spectral components, one is the tissue autofluorescence with a broad spectrum and the other is a red fluorescent dye with a sharp spectrum, if we want to apply the unimodality constraint or high pass constraint to the fluorophore, how do we ensure the constraint is applied to the dye not the autofluorescence? This reflects the component-wise rotation ambiguity and the mismatch between the initial estimate and its constraints. To disentangle such ambiguity, spectra are sorted in ascending order of their center wavelengths (intensity weighed average), i.e., by the color of the fluorophores. Therefore in this case, the red fluorescent dye is tagged next to the yellowish autofluorescence, and the constraints can be applied accordingly.
In addition to the sorting, the topology of the excitation and emission wavelengths is recorded to handle scans of the excitation filters, or emission filters, or both sets. It can also unmix multiple bioluminescent signals, where there is no excitation light. The wavelength topology can guide the algorithm to parse each spectrum to multiple sections so that the spectral constraints can still be applied.
Experiments and Results
For the in-vivo experiments discussed in this section, athymic nude (nu/nu) mice were anesthetized using 1.5 to 2% isofluorane through a nose cone system during image acquisition. The animals were fed with a low-chlorophyll food to reduce autofluorescence in the intestinal region. Images shown in this section are cropped from full size ones to save space.
The ability to separate multiple fluorophores, especially when they have significant spectral overlap, was tested in an in-vitro experiment. Two fluorescent dyes, Alexa Fluor®(AF) 680 and 750 and two quantum dots (QD) 700 and 800 (Invitrogen, Carlsban, California), were diluted and placed in a 96-well plate. Each fluorophore occupied one well. Raw images were taken using ten emission filters from , spaced every , and excited in reflectance mode at . The image overlays where fluorescent radiance images are superimposed on background grayscale photographs using rainbow pseudocolor are shown in Fig. 3a . From the top row to the bottom row, the fluorophores are AF680, QD700, QD800, and AF750, respectively. Four components were chosen to unmix, since the autofluorescence background was at a negligible level and signals were mostly contributed by four fluorophores. Spectral unmixing was performed with default settings that include non-negativity constraints.
Unmixed distribution maps are shown in Fig. 3b, where each column is displayed with its own pseudocolor scale; the unmixed images can be quantified to obtain the signal distribution of the fluorophore in the raw image. A composite image is also created to better visualize unmixed fluorophores in a single RGB color image. In this case, the distribution maps of AF680, QD700, QD800, and AF750 are combined to create a false color image shown as the last image in Fig. 3b by putting them into green, blue, magenta, and red channels, respectively. Figure 3b indicates that four fluorophores are successfully separated from the raw images where signals show a significant overlap. Improved contrast and clarity are obtained in the individual distribution maps. The unmixed spectra are shown in Fig. 3c and are similar to the published spectra for these fluorophores. Quantum dots and Alexa Fluor® dyes can be differentiated even though they have similar peak emission wavelengths.
Multiple Reporters in Vivo
This experiment gives an example of multispectral imaging of a mouse containing subcutaneous injections of of two fluorescence probes, AF680 and AF750, using IVIS Spectrum. AF680 was injected near the neck and AF750 was injected at the lower body. Figure 4a shows the raw spectral images in units of efficiency acquired using ten emission filters from and excited in the reflectance mode at . Clearly tissue autofluorescence contributed substantial signals in addition to the dyes. We assume the tissue autofluorescence is a single spectral component, and therefore selected three components to unmix. The unimodality constraint was used on the AF750 to guide the unmixing. Figure 4b shows the unmixed distribution maps for three components, i.e., autofluorescence, AF680, and AF750. Figure 4c shows the unmixed in-vivo spectra for these components. The composite image is shown as the last image in Fig. 4b. Compared to the in-vitro spectra of AF680 (cyan line) and AF750 (magenta line), the in-vivo AF680 spectrum (blue line) shows more significant red shift than that of AF750, which can be explained by the steeper absorption gradient with respect to wavelength in the shorter wavelength range. To demonstrate the benefit of our approach over the conventional method using the in-vitro spectral library, we used the in-vitro spectra of AF680 and AF750 together with the autofluorescence spectrum to perform the linear decomposition. The results are shown in Fig. 4d. It is clear that directly using the in-vitro spectra leads to severe cross talk among components in this experiment.
In certain animal experiments, reporters can be physically colocalized, so it is important to examine the efficacy of the unmixing algorithm in such a situation. In this experiment, two polyethylene tubes containing Alexa Fluor® dyes (left side: a mixture of AF680 and AF750; right side: AF680) were implanted next to the kidneys in a mouse. Figure 5a shows the raw spectral images in units of efficiency acquired using eight emission filters from using two excitation filters. The first four images were excited at , but the last four images were excited at to improve the AF750 signals. Since tissue autofluorescence was not significant, only two components were selected to unmix. A high pass filter at was used for AF750 to help the unmixing. Figure 5b shows the unmixed distribution maps for AF680 and AF750. The result matched the expectation, in that AF680 was resolved on both sides of the mouse but AF750 only appeared on right side. The yellowish mixture in the composite image corresponds to the mixing of green and red channels representing AF680 and AF750, respectively. Unmixed spectra in Fig. 5c shows well separated two spectral profiles.
This experiment gives an example of the multispectral imaging of a mouse implanted with a polyethylene tube containing quantum dots (QD) 805 in the transillumination mode. The tube was placed near the left kidney and the transillumination source was positioned approximately below the tube. Figure 6a shows the raw spectral images in units of radiance acquired using six emission filters from and excited at . In the first few images when QD800 has low emission, tissue autofluorescence is strong. Figure 6b shows the unmixed distribution maps for two components, i.e., autofluorescence and QD800. The unmixed autofluorescence pattern is similar to the light transmission. After unmixing, QD800 is successfully isolated. Figure 6c shows the unmixed in-vivo spectra for autofluorescence and QD800.
In previous examples, fluorophores of interest are relatively strong compared to tissue autofluorescence. Figure 7 gives an example to test the situation when tissue autofluorescence dominates. AF680 were sealed in the tip of a -diam rod and the rod was placed about deep underneath the low back of a phantom mouse (Caliper Life Sciences, Alameda, California). The phantom mouse was imaged using six emission filters from and excited in the reflectance mode at . The dye barely can be distinguished from the raw images shown in Fig. 7a, but is clearly separated in the unmixed maps shown in Fig. 7b. Here a high pass filter at for AF680 was applied to guide the unmixing algorithm.
Spectral unmixing outputs the distribution maps and the component spectra. To use the unmixing results to quantify the fluorophores, one should first understand the meaning of those outputs. Unmixing mathematically is a decomposition process asdenotes the measured image at the ’th wavelength and represents the decomposed image containing the signal only contributed by the ’th spectral component. With respect to Eq. 2, , where is the distribution map of the ’th component and is the spectral intensity of the ’th component at the ’th wavelength. In our algorithm, since each spectrum is normalized by its peak value, becomes unity at the peak wavelength, and thus indeed represents the measured signal of the component at its peak wavelength, and the spectrum implies the proportional distribution at other wavelengths. As long as the peak wavelength is included in multispectral analysis, the distribution map provides a consistent measure of the fluorophore signal, regardless of how many filters are used to unmix. Alternately, it has been suggested to use the integrated signal over all the wavelengths to quantify13 so that higher readings are achieved, but it makes it difficult to compare results between analyses done with different filter setups.
It also should be recognized that the signal distribution map bears a complex relationship to the fluorophore concentration map, not just because the fluorescence yield varies with excitation and emission wavelength, but also because the light transport is determined by the surrounding tissue. It is possible to correct the quantum yield with calibration, but accounting for the light transport is much more challenging and usually requires solving a complex 3-D model for photon diffusion. Therefore, the distribution map should not be regarded as a measure of absolute fluorophore concentration, but as a tool to compare fluorescence sources at similar imaging conditions and comparable depth.
Scanning Excitation Filters
Most typically, emission filters are used to scan through the emission peaks of multiple fluorophores when performing spectral unmixing. All the previous examples shown in this work utilize emission filter scans primarily. However, the choice of scanning emission filters is relatively arbitrary. There may be situations where scanning excitation filters, or scanning both excitation and emission filters, is desirable. One example would be separating quantum dots from an organic dye with a similar emission spectrum. Since the quantum dots and organic dyes have very different excitation spectra, it would be preferable to scan excitation filters over emission filters. An example of unmixing using excitation filters is shown in Fig. 8, where we unmix the same target used in Fig. 7 and achieve similar results.
In-vivo molecular imaging of bioluminescent or fluorescent reporters and probes has become a standard tool in academic laboratories and drug discovery companies, based on low cost, ease of use, and high throughput. However, due to the scattering and absorption of photons in tissue, and the presence of tissue autofluorescence, interpretation of images can be challenging.
Spectral unmixing is a powerful tool for reducing the effects of background autofluorescence and separating multiple probes or reporters. By separating the different spectral components and visualizing their distribution, users obtain improved clarity in interpreting their experimental results. We introduced a powerful and flexible spectral unmixing algorithm based on the multivariate curve resolution method that supports both epi- and transillumination fluorescence configurations, allows both excitation and emission filter scans, and also enables unmixing of bioluminescent sources. By optionally adding spectral constraints such as a unimodality or low/high-pass filters, the unmixing problem is further constrained, leading to improved robustness. The algorithm also supports the unmixing of unknown spectra, which is encountered frequently in animal models because tissue absorption can shift the spectrum by an unknown amount. Region of interest (ROI) tools are available to quantify unmixed probe distributions, and pseudocolor renderings are available for visualization of distribution maps. The algorithm is successfully demonstrated on numerous subjects including in-vitro well plates, phantom devices, and live animal models.
The authors would like to acknowledge Tamara Troy and Ning Zhang for their generous help in obtaining animal experiment data.