6.1. The problem of saturation

The saturation phenomenon is based upon the fluorescence detection system. The fluorescence detection system utilized in experiments presented in this report is a photomultiplier detector tube (PMT). This detector does not count individual photons; rather it requires a certain minimum number of photons to activate an electrode which will emit a small number of electrons which are subsequently amplified in a stepwise fashion. The readout of the PMT is on a 256 grey scale level. Occasionally the amplified signal can reach saturation if the fluorescence output of the calcium signal being detected is very high. Normally, the settings of the PMT are adjusted so as not to reach saturation, however, detection of the low end of the fluorescence signal is very important.

Notice that the grayscale values of some of the pixels that represent cell 2, shown in , reach a ceiling of 255; see . This is especially noticeable after the cell received the oxytocin stimulus around 1 minute into the experiment. Because the individual pixel values reflect the Ca^{2+} level in the cell, Ca^{2+} summary measures will undoubtedly be affected if the pixels reach the ceiling of 255. Also notice the variability in individual pixel values. The bottom panel of shows the intensity of 20 pixels over time and it is clear that some may reach maximum intensity values that are larger than 255 and some at much lower values. We do not model the behavior of individual pixels in this work but it can certainly be considered in the future.

Two questions are immediate. First, how is the EigenSignal affected when pixel values reach the saturation level of 255? Second, how should one process the Ca

^{2+} signal once pixel saturation has been detected? We implement the algorithm introduced by

Gabriel and Zamir (1979) and let the saturated pixels be missing data to address this issue. To our knowledge, there are no methods in the literature available to deal with the clarification of Ca

^{2+} signal curves and our attempt is the first of its kind.

6.2. The weighted SVD

Although the first EigenSignal is a reasonable measure to use when summarizing the Ca^{2+} signals of the pixel-wise matrices, there are drawbacks if used without adjustment. If there are too many pixels that reach the saturation point, the signal can be under or over-estimated at different time points in the experiment and a misinterpretation of the signal amplitude can occur. It is intuitive to understand how the signal can be under-estimated due to saturation, but over-estimation of the summary signal is certainly an unexpected phenomenon which we will explain. Because pixels reach a saturation, the signal summary is undoubtedly affected by this lack of information on maximal values attained by such pixels. In the course of the time series where many pixels reach saturation, the signal is under-estimated around the peaks, but the effect of this under-estimation results in an over-estimation during a time where no pixels reached saturation. and show this phenomenon. The over-estimation is due to the normalization requirement of the singular vectors and the right skewness of the cross-sectional intensity distribution.

To correct these over- and under-estimation effects, we must remove the effect of the saturated pixels and recalculate the Ca

^{2+} signal without their influence. One approach is to simply remove every row of the pixel-wise data matrix which contains a saturated pixel and recompute the Ca

^{2+} signal using the resulting matrix; however, this could lead to the removal of a significant number of rows from the data matrix. Instead, we propose to implement the weighted SVD, WSVD, using the low rank matrix approximation of

Gabriel and Zamir (1979) where we introduce the use of indicators in the weights, as in

Beckers and Rixen (2003), to treat the saturated pixels as missing data and use a clever choice of weights that allows for accurate recovery of the original signal.

It is important to note that our “missing data” is not really missing, we know that the saturated pixels must at least attain a value of 255. Hence, if we observe values that are below this threshold in our imputation, we would certainly know that we’ve made an error. We implement a check in our algorithm that gives us a flag if a value that is initially saturated falls below its saturation point.

Although an EM type method could certainly be applied in this context, we choose to use the iterative algorithm of

Gabriel and Zamir (1979) because of its speed in convergence and because we do not wish to make distributional assumptions about the data. Now, because the signal variance in our data follows the behavior of the signal itself, we opt to use a variance weight in the imputation scheme. Of course initialization of the algorithm is tricky, in particular, when

*w*_{ij} = 0, but our use of the first EigenPixel and EigenSignal to initialize the algorithm proves to work well; see

Gabriel and Zamir (1979) for more discussion on initialization.

The premise of our approach is that each cell has a “true” Ca^{2+} signal and we are not able to observe that signal because there are only a finite number of pixel values available to capture it. The details of our implementation are provided below.

Let

**u** and

**v** be the first EigenPixel and EigenSignal associated with the second SVD used to extract the putative Ca

^{2+} signal, which includes saturated pixels, so that

**u** and

**v** comprise most of the pixel and signal information of some cell. Continuing with the example from the previous section, the matrix of interest is

. Let the dimensions of the

be

*n* ×

*m*; because most of the variation is explained by the first component in the SVD, the rank one approximation can be obtained by minimizing the error

with respect to

**u** and

**v**. We also wish to weight each term in the double summation so that it removes the influence of saturated pixels and takes into account the appropriate variation. We let the weights be

*w*_{ij} = I

_{ij}/(

*u*_{i}υ

_{j})

^{2}, where I

_{ij} = 0 when

, that is, pixel is saturated, and I

_{ij} = 1, otherwise.

Beckers and Rixen (2003) proposed using an indicator to deal with missing data. We supplement this approach by using (

*u*_{i}υ

_{j})

^{2} to scale the term in the summation of (

6.1) so that the variance is no larger than 1. Now our new minimization problem becomes

We solve the minimization by alternating between

*u*_{i} and υ

_{j} . Fixing

*j*, we can expand the expression in (

6.2), let

and we get that

solves that portion of the minimization. Similarly, if we fix

*i*,

, where

. The new proposed EigenPixel and EigenSignal vectors are

**u**^{new} =

**u**′/‖

**u**′‖ and

**v**^{new} = v′/‖

**v**′‖ respectively. This gives us a recurrence relation that we can use to obtain a clearer version of the EigenPixel and EigenSignal, where the EigenSignal will represent the clarified Ca

^{2+} signal of interest.

Beckers and Rixen (2003) offer a similar recurrence as a way of imputing missing values in oceanographic data. We change the number of relevant components in the SVD and add a weight that includes a rescaling factor 1/(

*u*_{i}υ

_{j})

^{2}. We include this variance rescaling factor because the variance of the signal and the signal are synchronized and we want to account for that effect. The pseudo code used to program this is shown below:

- Let
**u**^{0} and **v**^{0} be the initial EigenPixel and EigenSignal vectors obtained by taking the SVD of the pixel-wise matrix that comprises all the pixel and signal information about the cell of interest, including saturated values. - The first proposed EigenPixel and EigenSignal are
**u**^{1} = **u**′/‖**u**′‖ and **v**^{1} = **v**′/‖**v**′‖ respectively, where . - The (
*k* + 1)st proposed EigenPixel and EigenSignal are **u**^{k+1} = **u**′/‖**u**′‖ and **v**^{k+1} = **v**′/‖**v**′‖ respectively, where . - Iterate until convergence.

The missing values are imputed by the corresponding *u*_{i}υ_{j} after the convergence of the algorithm. We check to make sure that any imputed value for initially saturated pixels do not fall below its saturated value. In our application of the algorithm to the real data, all imputed values passed this test. Very rare numbers of pixels experience total saturation, that is, I_{ij} = 0 for all *j* = 1,…, 512 given some fixed *i*. This is particularly true if we only consider a subinterval of the 85 minute run. Since it is not possible to impute values for such pixels, they are dropped from our analysis to avoid creating bias.

Consider

, the

*n* ×

*m* pixel-wise matrix of

*n* pixels and

*m* time points. For a fixed pixel

*i*,

has a total of

*m* potential saturated values. If I

_{ij} is the indicator described above and if we let

be the number of nonsaturated values for the

*i*th pixel across time, then we made it a rule to remove the

*i*th pixel if

*p*_{i}/m < 1/8. This means that we remove any pixel row of the matrix

*X* if more than 7/8 of it is saturated.

6.3. Application of the WSVD to simulated data

To evaluate the accuracy of our method, we applied it to a simulated data set where we used sine curves to emulate the behavior of typical cell data as shown in . Our simulated data matrix, and first EigenSignal and EigenPixel are shown in . The data shown in represent the true signal we are trying to recover. Real data, however, have noise and also possess saturated pixels that dampen the signal. To duplicate this behavior, we threshold the data matrix so that everything larger than 0.50 is replaced by 0.50; this mimics a saturation at pixel locations that have values larger than 0.50. To add noise, we introduce realizations from a Gaussian distribution with mean 0 and variance proportional (*u*_{i}υ_{j})^{2}. We introduced this variance into the simulated data because it is consistent with the type of variance observed in the real data and we wanted to emulate that behavior. shows the original and saturated first EigenPixel and EigenSignal curves. The first EigenSignal and EigenPixel vectors from the saturated data are both dampened and exaggerated in different regions.

After applying the weighted SVD to the saturated data matrix, we see that upon convergence of the algorithm the resulting first EigenPixel and EigenSignal both come very close to the original curves. The relative error at every pixel and time point are shown on the right panel of . When we compare the ratio of the error sum of the saturated over the WSVD Eigen vectors, we see an 11 fold difference in the EigenPixel and a 8 fold difference in the EigenSignal. To see the effect of our weight on the results, we removed the weights and repeated the analysis. Comparing the ratio of the error sum of the saturated over the WSVD Eigen vectors with no weights, only a 1 fold difference in the EigenPixel and a 1 fold difference in the EigenSignal are observed.

6.4. Application of the WSVD to actual data

For illustration we apply the weighted SVD to the pixel-wise matrices of cell 2 from the treated group of low hormone level and cell 11 from the control group of high hormone level. shows the first EigenPixel and EigenSignal of both cell 2 and cell 11. We see that in the peak region of both, between 0–4 minutes into the experiment, there is a large difference in the EigenSignal vectors, especially for cell 11. This is not surprising since most of the saturation occurs in the peak region of the experiment, hence, pixel imputation will mostly affect this region. To further explore this phenomenon, we apply the WSVD only in the peak region (0–4 minutes) and results are also shown in . We see that the saturated pixels were dampening the expression in the peak region. This is a key finding since it is believed that Ca^{2+} expression in this peak region could be used to characterize cells studied.

We have shown that the weighted SVD can be used to clarify the Ca^{2+} signal in the cells presented. This is an important step when harnessing Ca^{2+} expression from these cells, especially because Ca^{2+} expression is dampened drastically if we do not take into account the saturation effect.