Detection and Measurement of the Intracellular Calcium Variation in Follicular Cells

This work presents a new method for measuring the variation of intracellular calcium in follicular cells. The proposal consists in two stages: (i) the detection of the cell's nuclei and (ii) the analysis of the fluorescence variations. The first stage is performed via watershed modified transformation, where the process of labeling is controlled. The detection process uses the contours of the cells as descriptors, where they are enhanced with a morphological filter that homogenizes the luminance variation of the image. In the second stage, the fluorescence variations are modeled as an exponential decreasing function, where the fluorescence variations are highly correlated with the changes of intracellular free Ca2+. Additionally, it is introduced a new morphological called medium reconstruction process, which helps to enhance the data for the modeling process. This filter exploits the undermodeling and overmodeling properties of reconstruction operators, such that it preserves the structure of the original signal. Finally, an experimental process shows evidence of the capabilities of the proposal.


Introduction
Calcium (Ca 2+ ) is an ubiquitous intracellular ion signaling responsible for controlling many cellular processes [1,2]. Ca 2+ acts as second messenger triggering pathological events, such as cells injury and death, as well as participating in pathological conditions, such as hypertension, cardiac arrhythmia, hematological problems, muscular diseases, and hormonal disorders, among others [1,3,4]. The role of Ca 2+ in these diseases has just begun to be understood. Some Ca 2+ that induce pathological conditions have been found with the help of substances that interfere with the movement or activation of Ca 2+ [3,5]. Due to the importance of Ca 2+ , there are several numbers of optical and nonoptical techniques, which have been developed for analyzing Ca 2+ either dynamics or concentration. Fluorescent microscopy techniques are frequently used to observe the variation of intercellular Ca 2+ concentration applying chemical fluorescents indicator as markers. Those indicators stimulate the cells causing a fluorescence effect [6]. The fluorescence is detected by microscopy and CCD sensors. To compute the cells with the greatest fluorescence variations, the user must select them manually and analyze all the images in the sequence in order to determine the changes of fluorescence over time. However, it consumes time and human resources being susceptible to error measuring. By such circumstance, the process for segmenting the cells and the study of dynamics of intracellular Ca 2+ represent the main objective in this paper.
An important step in the study of intracellular Ca 2+ consists of the segmentation of each individual cell. The image segmentation results are difficult because the environmental changing conditions are usually uncontrollable. In the literature typically depending on the properties of the cells, the segmentation method is proposed. For instance, some works use neural networks approaches [7], hierarchical threshold [8], or multiscale morphology [9], just to mention a few.
Watershed-plus-marker approach, on the other hand, is the traditional image segmentation method based on morphological mathematical [10,11]. The success of this method depends mainly on the correct detection of the image's markers. The markers can be detected manually or automatically. Automatic approaches help the specialist to save time and resources. However, there are factors that affect the performance of an automatic detection of markers such as noise, cells occlusions, and abrupt changes in the images. Those factors can lead these algorithms to overand undersegmentation, that is, create regions containing partial or multiple cells. In this sense, several approaches have been developed for improving cell segmentation. In this study, we introduce a method to analyze automatically the intracellular calcium variation. This approach consists in two stages: (1) the image enhancement and cell segmentation and (2) the calcium variation modeling. The image enhancement is carried out with a top-hat filter, which homogenizes the luminance conditions. The process of cells segmentation is performed using the marked-controlled watershed transform and filters by reconstruction, which is used to detect markers efficiently, after the homotopy of the gradient image was computed. To measure the calcium intracellular variations, the volume of each marked cell is computed. The procedure consists of estimating the volume using the luminance intensities for each marked cell along the time. After, least squares fitting (LSF) method is applied to create a model of the behavior of the variation of the fluorescence. The behavior model created is considered as an exponential decreasing function. To enhance the development of the model of the behavior of calcium, a novel morphological filter named as medium filter is introduced. This filter smoothes the fluorescence measures, exploiting the under-and-overmodeling of reconstruction operator, preserving the information structure of original signal.
The paper is organized as follows. In the next section, a review of some morphological filters is presented. In Section 3, a method based on marked controlled watershed transform to detect automatically the cells is shown. In Section 4, we show the procedure to estimate the volume of each marked cell starting from the fluorescence intensity along time; after applied least squares fitting and morphological medium filter, the model of fluorescence behavior is created. Finally, results and conclusions can be found in the last section.

Basic Notions of Morphological Filtering.
Mathematical morphology is mainly based on the so-called increasing transformations [12][13][14]. A transformation is increasing if for all pairs of functions and , with ≤ ⇒ ( ) ≤ ( ). In other words, increasing transformations preserve the order of the relation. A second property is the idempotence; that is, a transformation is idempotent if and only if ( ( )) = ( ). The basic morphological filters are the morphological opening and the morphological closing with a given structuring element, where represents the elementary structuring element (3 × 3 pixels, e.g.) containing its origin and is an homothetic parameter. Thus, the morphological opening and closing are given, respectively, by where the morphological erosion ( ( )) and dilation ( ( )) are expressed as = ⊖ : → inf ( + ) and ( ) = ⊕ : → sup ( − ). Henceforth, the set will be suppressed rendering, the expressions and are equivalent ( = ). When the parameter is equal to one, all parameters are suppressed ( = ).

Opening (Closing) by
Reconstruction. The notion of reconstruction is a very useful concept provided by MM. Reconstruction transformations are built by means of the geodesic transformations. In a gray-level case, the geodesic dilation 1 ( ) (resp., the geodesic erosion 1 ( )) with ≤ (resp., ≥ ) of size 1 given by 1 ( ) = ∧ ( ) (resp., 1 ( ) = ∨ ( )) is iterated until idempotence. Consider two functions and , with ≥ ( ≤ ). Reconstruction transformations of the marker function in , using geodesic dilations and erosions, expressed by ( , ) and * ( , ), respectively, are defined by When the marker function is equal to the erosion or the dilation of the original function in (2), the opening and the closing by reconstruction are obtained:

Automatic Detection of Cells
The watershed-plus-marker approach transformation is a traditional image segmentation method based in mathematical morphology [10,13]. However, this transformation makes use of an extensive set of morphological filters. This transform is used for segmenting images avoiding the oversegmentation [11]. The oversegmentation criterion consists of setting an upper limit in the number of minima regions detected. This process is performed with the minima impositions over the markers, exploiting homotopy property of the operators.
However, it needs to make some assumptions to use this approach. The most important assumption consists of the fact that that the minima represent the center of the object (in the cell case, the core of it); a second assumption consists of the gradient estimation, such that the fact that it could be performed with morphological operators.

Marker Detection.
Due to the features of the images, the nucleus of each cell is used as regional minimum. As matter of fact, a regional minimum of a gray-scale image is a connected component of pixels with uniform altitude without lower neighbors. Before computing the minima, the nucleus of the cells is mainly dark, surrounded by brighter region composed by the cytoplasm. However, the local luminance conditions of each cell differ singly to each other, affecting the detection of the nucleus. For the homogenization of luminance conditions the top-hat transform is used as a local contrast correction filter.
For the { } be sequence of image. The top-hat transformation is defined as follows: where the dimensions of structuring element are related with the luminance conditions of the scenario; in such a way, the luminance distribution in the image may be approximated with the morphological opening . Whenever the dimension of structuring element becomes proportionally similar to the image dimension, it would represent global luminance sources affectations; on the other hand, small dimensions represent local luminance variations effects. This process is illustrated in Figure 1, where each process step is showed and in certain steps the image has been coded in pseudocolor to point out the effects involved at each step. The direct appliance of the maxima transform detects all the maxima, including noise data, as it is appreciated in Figure 2(b). To avoid the extra minima detection, the image is enhanced with a closing by reconstruction operator. This operator allow grouping by all connected local minima, discarding the majority of noise effects. The closing by reconstruction uses a unitary structuring to approximate the original image. For illustration proposes in Figure 2(c) is appreciated the effect of apply the filter in the minima detection process. Observe that some cell nucleus are well detected, but others are omitted; due to the acquisition process the cytoplasm is not completely closed. This situation can be fixed using subsequent images, where additional minima are correctly detected. Thus, in order to obtain the maximum number of minima, a function that captures the occurrence of the minima in the sequence is constructed. Let { } ∈ and { } ∈ be the images of the sequence and the images containing the detected minima, respectively. ( ) is a binary image such that it takes 1 value if the point belongs to a regional minimum and 0 values otherwise. After using the subsequent frames, summation is built as follows: The surface is drawn in Figure 2 They are discarded by a thresholding criterion. In the case of study, each cell is typically about four pixel of radius, which can be discovered using an opening by reconstruction filter. As complementary, a closing operator with a structure of 3 pixels of dimension is used to connect insolated regions. The process is illustrated in Figure 2. First, a morphological closing of size 3 is applied to fill the small holes; the results can be appreciated in Figure 2(c) (before closing filter) and Figure 2(e) (after closing filter). Next, applying the summation (Figure 2(d)), the minima are found, and finally the regions with big and small areas are discarded, denoting the cells.

Gradient
Operator. The watershed-plus-marker approach makes use of the gradient operator to impose the markers. In this sense, the morphological gradient can be used as contrast detector. Let ( ) be a function defined in Z 2 and the basic structuring element of 3 × 3 dimensions, centered at point . Then, the transformation is defined as follows for a discrete space: There are other two versions of gradients in mathematical morphology, the internal and the external gradients defined, respectively, as follows: Experimentally it is appreciated that external gradient offers smoother and thick borders (see Figure 3(c)), instead of internal gradient offering defined and clear borders as it is illustrated in Figure 3(c).

Imposed of Minima by Reconstruction.
Once the markers cells signals are detected, these are imposed with their minima on the gradient image. To carry out this task the following procedure is performed. Let and be the set of markers computed as commented above and the gradient image, respectively. After, two news functions are built: the first one consists of a thresholding function ( ), which is defined as ( ) = 255 if ∉ and ( ) = 0 if ∈ , while the second one is built through the gradient image as ( ) = ( ) if ∉ and ( ) = 0 if ∈ . Furthermore, the dual morphological reconstruction of ( ) inside of ( ) is made, denoted by * ( , ). The function * ( , ) only has the minima of , such that the watershed transformation is applied. Figure 3(d) illustrates the results getting after applying watershed segmentation.

Modeling the Intracellular Calcium Dynamic
In this section, we deal with the problem of modeling the intracellular calcium dynamic. The procedure consists in three parts: the estimation of the calcium volume, the fitting of an exponential curve, and the calculus of the error.

Modeling Intracellular Calcium Variations.
As it is appreciated in Figure 4, the dynamics of the calcium stimulus has an exponential behavior. Then, the purpose consists of creating a model of the decreasing behavior stimuli suffered by each cell, where the region of interest is located among the global maxima and the end of the signal. However, due to the noise, it is not possible to detect easily the maximum. To attenuate this inconvenience, an automatic process is performed detecting the maxima for the function { ( )} ∈ . The process consists of a sequential alternating filter in onedimensional scenario. The alternating filter is constituted by a sequence of one closing by reconstruction followed by one opening by reconstructioñ(̃( ))( ) where the size of is varied into the interval [0, ]. The filter undermodels the original signal, smoothing the signal wave and allowing the detection of the global maxima efficiently. Figure 5 illustrates the detection of a representative maximum detected that corresponds to a connected element in one dimension space. The center of the connected element represents the maxima location, such that it is estimated with the mean of the connected elements; that is, modeled as a polynomial decay time-decreasing function as follows: where are the polynomial parameters, where the data used is taken from the maxima to the end of the data. The parameter estimation is performed by least squares the follow expression [

⋅⋅⋅
] and represent the polynomial order for fitting. For illustration purposes in Figure 6, a fitting sample is showed. The exponential help to model and analyze the decrease of the intensity registered in each cell.

Error Model Fitting.
The correct construction of model over data is defined by introducing two measures of error: BIAS error and RMSE error. The first one is a measure of error modeling. The second measure is a precision error modeling criterion. The bias error provides information about how the model fits real data. Negative bias error means that model is undermodeling the data; that is, the model is a function under real data. Consequently, positive bias error represents overmodeling. Values near to zero mean that the model catches the dynamic of real data. Formally, bias error is defined as Bias( , * ) = ∑ =0 − * , where represent real data and * estimated data. Note that when BIAS is equal to zero it does not mean that the model is correct. It means that the same proportions of measures are below and under for real data. Then to quantify the precision error the RMSE is used. This error is the average of absolute differences among real and modeling data. RMSE is defined as follows: RMSE( , * ) = (1/ ) ∑ =1 ( * − ) 2 where * represent modeling function and real data.

Enhancement of Data.
Although the least squared method offers the optimal model, it depends on the data measurement having normal distribution. Then, by the nature of the model, it results hard in verifying that these measures have a normal distribution. As a consequence it is necessary to enhance the data in order to facilitate the convergence of the approach. For simplicity, it is assumed that any signal ( ), resulted from the calculation of the volume of an cell for an time, is affected by additive noise with zero mean as follows: where * is the signal free noise and is the additive noise with mean zero. In fact has zero mean; the true signal data * is located into min{dom( )} and max{dom( )} values. However, given is a random variable, it locally should not present a zero mean, making it difficult to estimate the * value. To figure it out, it is needed to analyze locally the information, inferring the trend, and make an estimation of the expected value. The proposal consists of exploiting certain properties of operators taken from morphology operators. The reconstruction operators are useful because they approximate a surface by iterating successively a marker, getting the other surface that has similar topological properties. The approximation does not keep the original level of detail of its shape, such that it depends on the form and properties of structural element used. It should be considered inconvenient, but, in practical terms, it is its major advantage in sense and it represents the main trend of the original data, eliminating variations less than the structural element (high frequencies) of the original signal and resulting in a new signal that under-or overmodels the original data.
Considering the basic operators by reconstruction (opening and closing), the property of extensive or antiextensive, respectively, cause the fact that the application of each one over a signal results iñ( ) or̃( ) signals such as under-or overmodeling the original. Both of them remain as the global trend of the topological information of . Consequently, the residual presents important topological information. However, the distribution of the data changes slightly: the shape of the derivative of the original and the approximated signal is different, changing the statistical properties of its PDFs. Figure 7 presents the probability density  function (PDF) approximated via its histogram after application over a signal . The histogram of opening operator presents a negative deviation, which means that the surface approximated is undermodeled. On the other hand, when we apply a closing operator, it overmodels the original signal and its histogram is deviated to the positive side of the range. The proposal consists of mixing both filters, preserving the statistical information of the original signal. Noise effects are represented by the high frequencies. These frequencies must be discarded preserving the global trend of original signal . The discarded frequencies are directly related to the size of the structural element and the sampling process; that is, given a structural element of size , it represents a temporality of , where is the mean frequency of acquisition of . Then the process of filtering ( ) is statistical consistent if and only if * , minus preserve the following equality: This is, the density distribution function of the difference between filtered data and original data is a normal distribution centered at the origin. The development of correct statistical filter must satisfy (9), where it is appreciated that opening and closing reconstruction operators provide negative and positive bias information of the approximated surface. The original signal is enveloped by the opening reconstruction and closing reconstruction, respectively, such that̃( ) ≤ ≤̃( ). Consequently, for estimating , using̃( ) and̃( ) and considering that [{ 1... }] = 0, and approximation to is where 1 and 2 are values between [0, 1] and its sum is the unit. In case that̃( ) and̃( ) use the same structural element, 1 = 2 = 0.5. In other cases these values would vary depending on effects of the geometry in the reconstruction process. The filter described above is denoted as a medium reconstruction filter. An extension of this filter implies a sequential form, where the properties of the structural element used in reconstruction stage should be varied as follows: let ( , ) be a function that returns a structural element with particular properties for instant; sequential version of medium reconstruction filter is defined as Note that function ( , ) would vary the size and the topology of the structural element. The topology and size will affect the model that fits the data. The effect of applying the medium filter by reconstruction is illustrated in Figure 8, where in Figure 8(a) are presented the original data (blue color) and the filtered data (red color). As is appreciated, filtered signal follows the main trend of the original signal, discarding the high frequencies, and always statistical properties are kept as it is appreciated in Figure 8 data fit the exponential decreasing function. For a detailed analysis the bias and RMSE error are showed for the case of filter and nonfilter signal. Bias error behaves to close in both scenarios; but the RMSE is deeply reduced, which means that the fitting process results are better, after filter data (Table 1).

Results and Discuss
The proposal described above is tested under an experimental method that consists of analyzing a sequence of images that contains cells, which are exited applying Flour 4, in order to measure the effect over Ca 2+ belonging to each cell and characterizing its behavior. The process is illustrated in Figure 9. The process diagram sums up the sequence of processing steps done over the sequence. The sequence of images was acquired from biological researchers of the Institute of Neurobiology, Campus UNAM-UAQ. The sequence was obtained from cells of Xenopus laevis frog. The calcium is measured indirectly with its excitation via Fluo-4 (by Molecular Probes). The optical material consists of a microscopy of fluorescence, setting up in an Olympus camera sensor IX71 at 485 to 520 wavelength sensitivity nm (excitationemission, respectably); finally the images were acquired with fast acquisition camera (Evolution QEi Media Cybernetics), at 30 frames per second (Fps) with a resolution of 320 × 240 pixels. Finally, for testing purposes, 1,000 images have been selected, which represents temporary a sequence of 33second length. The cell detection is a tough task because there are many factors which inside directly in the analysis process, as the nonhomogeneity luminance conditions of the images and the conditions that present the cells of interest. After the image acquisition, the following task consists of finding out and segmenting each cell. This process is performed via watershed approach. However, starting from the first frame acquired does not warrant correct cell detection. To make more robust the cell process detection, for each image, the cells are detected, as described in the third th section. Once  the cells are detected, the neighborhood around the cell is considered to analyze the calcium concentration. The calcium concentration is performed by the measure of the luminance of each cell. The relation between the luminance intensity of the cell is highly correlated with the calcium concentration; that is, cells with high luminance have major calcium concentration. Next, the creation of the model of behavior results in a difficult task, because the behavior observed is not linear with the use of autoregression methods being inappropriate [15]. Then, discarding the times where the cell started to become excited, the dynamic of decreasing is modeled with an exponential function via least square method. The selected data includes the maxima location to the vanishing exited behavior. As recreated better the model, before applying least squared, the median filter by reconstruction is applied, which improve the accuracy of the modeling. Finally the results are showed in Figure 10. Observe, in Figure 10(a), that the cells are detected and the dynamic modeled as exponential superposed over measured data is showed (Figure 10(b)). The use of filter dismisses the high frequencies smoothing the behavior of luminance variations. The dismissing of high frequencies adds extra accuracy warranting that the exponential fitting has more significance, although the data are affected by noise effects. Finally, the way of how the cells were segmented represents a framework to analyze the intracellular calcium, which segment automatically the set of cells. This process is convenient in the sense that many of microscopic dynamics could be analyzed efficiently providing better information to the biologists.

Conclusions
In this paper, an automatic method for the study of intracellular calcium based on a marked controlled watershed transform for segmenting stage is presented. A new filter based on reconstruction operators is introduced. Then, having a high precision of cell segmenting and efficient ways to discard the noise measurement result the base for an automatic frameworks analysis as the experimentation shows. Finally, the reconstruction operators applied over one dimension data results usefully in the development of filters that help to create models of the dynamic of the calcium.