2.2. Unmixing

Assuming that several fluorescent species are present in a sample, they can all potentially contribute to the signal collected in a given pixel. In principle, determining the relative abundances of these fluorophore species via bleaching can be achieved by fitting the amplitudes of a multi-exponential decay at each pixel. This basic problem occurs in many arenas from magnetic resonance imaging [22] to fluorescence lifetime imaging [23,24] and nuclear physics [25]. There are a wide range of computational approaches to this challenge, such as maximum likelihood estimation [25], the method of least-squares [26], method of moments [27] and the Gardner Transform [28]. However, our problem is more general as it can include both spectral and bleaching information. Accordingly, we extract the spectral and photobleaching characteristics from the data set itself. This self-calibrated approach avoids using physical models of photobleaching that may not be consistent with real world samples. We will refer to the spectral-bleaching characteristic of a fluorophore as its spectral-bleaching fingerprint. An analogous quantity has recently been used to unmix fluorescent probes based on their fluorescence lifetime and emission spectra [24]. Though the combination of spectral and photostability information is the most general form of our technique, we note that we can also use photostability information by itself, as will be seen in the Results section.

To separate pixel-by-pixel contributions of each fluorescent label, we use the MATLAB non-negative least-squares (NNLS) function lsqnonneg to solve the relevant linear unmixing problem [19]:

where I_k(x, y) is the (measured) intensity at pixel (x, y) for frame index k (spectrally concatenated if spectral information is included), a_i(x, y) ≥ 0 is the relative scalar abundance of fluorophore species i at pixel (x, y) and v_ik is the k^th entry of the T-element (spectral-) bleaching fingerprint for fluorophore type i. Note that the system is overdetermined when T > N (N is the number of fluorescent labels in the sample). The extra information in this overdetermined system is crucial for noise suppression. The non-negativity prior for a_i prevents the unmixed abundances from reaching nonphysical negative values and improves unmixing fidelity. The i^th image - the abundance map of fluorophore type i - is given by a_i(x, y) and ideally contains only signal from the i^th fluorescent species.

If there are regions of the sample where each fluorescent probe exists in isolation, then the (spectral-) bleaching fingerprints v_ik can be manually identified. However, for many samples, this is not the case.

For such situations (typically cellular samples), the bleaching fingerprint of each fluorophore along with their pixel-by-pixel abundances, are simultaneously estimated by non-negative matrix factorization (NMF) [20]. This algorithm attempts to find non-negative bleaching traces and fluorophore abundances that solve the mixing problem of Eq. (1) in the least squared sense, using the alternating least squares (ALS) procedure [20]. The non-negativity constraint restricts possible solutions for the bleaching traces and the abundance maps to those with only positive values. As with NNLS above, this guarantees that the result is consistent with the fact that the intensity is a non-negative quantity. We perform NMF by using the built-in MATLAB non-negative matrix factorization function nnmf. We supply the nnmf function with the principal components of photobleaching curves as the initial estimates of the bleaching characteristics. For added robustness, we take the optimal solution out of three replicates – the first seeded with principal components as initial solutions and the following two with random initial guesses. Typical results required no more than 25 iterations.

For highly multiplexed samples, the NMF approach above may converge to unphysical solutions where bleaching traces increase in intensity over time, yielding incorrect unmixing results. This can be mitigated by imposing additional restrictions to the solution space. To this end, we implement a modification to the standard ALS NMF procedure, which we call non-increasing NMF (NI-NMF). At each iteration, if the value of a bleaching trace estimate B(t) at time t = t_n+1 exceeds the value of the bleaching trace estimate at time t_n, then we set B(t_n+1) = B(t_n). This is modification is performed for each time point in succession so that the value of the bleaching trace is monotonically decreasing. We run NI-NMF for 25 replicates, with the first one being seeded with the principal components, and select the result with the lowest mean-squared error. For all NI-NMF and NMF procedures it was found to be beneficial to exclude dim pixels to reduce the influence of noise and decrease computation time. We typically exclude all pixels with that are dimmer than 1-10% of the brightest pixel in the data set. Once the estimate for the bleaching curves is found, the least squares solution for the abundances at every pixel in the image is found by using the timelapse data and the estimated bleaching curves together with the MATLAB backslash operator (a QR solver), and then setting negative abundances to zero. This is equivalent to including all pixels on the final iteration of the ALS (NI-) NMF algorithm.