Abstract. Analysis of biomolecules in cells essentially relies on fluorescence microscopy. In combination with fully automatic image analysis it allows for insights into biological processes on the sub-cellular level and thus provides valuable information for systems biology studies. In this paper we present two new techniques for automatic segmentation of cell areas and included sub-cellular particles. A new cascaded and intensity-adaptive segmentation scheme based on coupled active contours is used to segment cell areas. Structures on the sub-cellular level, i.e. stress granules and processing bodies, are detected applying a scale-adaptive wavelet-based detection technique. Combining these results yields fully automated analyses of biological processes, and allows for new insights into interactions between different cellular structures and their distributions among different cells. We present an experimental evaluation based on ground-truth data that confirms the high-quality of our segmentation results regarding these aims and opens perspectives towards deeper insights into biological systems for other problems from systems biology.
Keywords: Segmentation, Wavelets, A contrario, Active contours, Fluorescent microscope imaging, Cells, Granules.
Introduction
The analysis of biomolecules in cells essentially relies on fluorescence microscopy imaging and provides valuable information in biomedical as well as systems biology studies. However, image analysis is frequently biased by subjective evaluation or limited in terms of quantitative assessment. One main aspect in this respect is the quantitative analysis of distinct sub-cellular particles in eucaryotic cells. Two classes of such sub-cellular particles are processing bodies (PBs) and stress granules (SGs) (e.g. [2,5,1]), sometimes also subsumed under the global name granules. It is presumed that these play an essential role in RNA metabolism. Processing bodies are characterized by an accumulation of various proteins involved in mRNA degradation [2]. However, if mRNAs are degraded in PBs is still controversially discussed as reviewed for instance in [4]. During cellular stress another kind of RNA-protein granules appears in tissue cultured cells, the so called stress granules (SGs).
Fig. 1. Fluorescently labeled nuclei and sub-cellular particles: (a) two cell nuclei, (b) a cell with labeled PBs, and (c) one with labeled SGs
In contrast to processing bodies, poly-adenylated mRNAs were shown to accumulate in SGs together with various RNA binding proteins whereas factors involved in mRNA degradation are barely observed in these granules [24]. Intriguingly, transient association of SGs and PBs was observed [18]. Together this led to the view that while mRNAs are likely to be degraded within PBs or by factors associated with these, RNAs are rather stabilized upon their association with SGs [16,17]. Due to the transient association of both granule classes it was suggested that mRNAs could thus be directed for degradation upon their transfer from stress granules to processing bodies. To investigate the role of SGs and PBs in controlling mRNA fate it is essential to monitor the subcellular distribution, number, size and contact events of both granule classes in a quantitative manner. Therefore it is instrumental to develop analysis tools which allow the automated and quantitative assessment of these granules in segmented cell areas on the basis of fluorescence microscopy.
In cell biological analyses sub-cellular structures are usually labeled on the single molecule level by direct labeling, indirect immune-staining or by genetic approaches. The latter allows to detect biomolecules directly by the fusion of fluorescent proteins as well as indirectly by an association of the fluorescent fusion proteins with other biomolecules, for instance aptamer-containing exogenous RNAs (see for instance [22]). All these techniques lead to the labeling of biomolecules by fluorochromes with distinguishable spectral properties and thus allow multi-channel imaging. Unfortunately explicit labeling of the complete cell area is typically not possible and enforces to extract it from one of the available channels originally intended for detection of other particles. In the here described experimental set-up, three distinguishable channels were used for the labeling of biological samples. Nuclei were labeled by DAPI staining and were characterized as more or less homogeneously textured bright blobs. Processing bodies (PBs) as well as stress granules (SGs) were labeled by indirect immune-staining using granule specific primary antibodies and subsequently secondary antibodies linked to distinct fluorochromes. This labeling strategy identifies PBs as small bright spots with a small size variance while most SGs appear significantly larger than PBs and show a pronounced variance in size (cf. Fig. 1).
Due to the significant variation in appearance of the different sub-cellular particles there is no integrated segmentation approach that allows to detect all kinds of particles and the cells themselves. Accordingly, we apply two different techniques for detecting PBs and SGs on the one hand, and the cell areas on the other. PBs and SGs detection relies on a scale-adaptive wavelet-based detection approach able to cope with the variance in size of these particles [14]. Cell area segmentation is conducted in the fluorescence channel for PBs in each image. To solve the task we adopt coupled active contour models.
One main contribution of this paper is the scale-adaptive extension of a wavelet-based particle detection algorithm. As the second main contribution we present an extension of active contour models towards a new cascaded segmentation scheme that yields larger flexibility in segmenting objects with specific intensity distributions. In detail, since the target area of a cell does not show significant discontinuities along its border and also cannot be modeled by one homogeneous intensity distribution, customary active contour energy models are not appropriate to solve the given task. Our new cascaded approach overcomes these problems and allows to cope with objects of non-homogeneous intensity distributions and weak discontinuities along their borders by incremental adaptation to the objects’ intensity characteristics.
The remainder of this paper is organized as follows. After reviewing related work in the next section, we introduce our scale-adaptive wavelet-based approach to particle detection in Section 3. In Section 4 the new cascaded technique for intensity-adaptive cell segmentation is presented. Subsequently experimental evaluations based on ground-truth data are discussed (Section 5), and we conclude with some final remarks in Section 6.
Related Work
Detection of cells and of granules in fluorescent microscope images is an instance of the general segmentation problem in image analysis. As these images show special characteristics, adaptations are required and have been proposed.
A level-set based approach for segmentation and tracking of cells is proposed in [10] and [11]. The first frame of a sequence is used for initial segmentation, where the fitting term of the classical Chan-Vese model [6] is replaced with a Gaussian likelihood for the intensity values with unknown variance. Lumped cells are separated using the watershed transform and subsequent region merging. A multi phase level-set technique is employed for tracking, where a coupling term of multiple level-set functions as proposed in [9] is used. Nearby cells are separated via the Radon-Transform in addition to the coupling term.
For the detection of spot-like particles several approaches exist, which often rely on global and local thresholding techniques like Otsu’s global method or the local Niblack operator [23,3]. Further techniques include sampling from an image intensity density estimated via h-dome transform and subsequent clustering of samples [21]. The method in [19,13] is based on wavelet decomposition, but best-suited to detect particles with limited variation in size.
An integrated algorithm for cell area segmentation and sub-cellular particle detection was proposed by [8]. They combine wavelet-based particle detection adopting the approach of [19] with active contour-based cell area segmentation. Unfortunately their snake energy is not transferable to our domain as the cells in our scenario show quite inhomogeneous and grainy intensity distributions which do not allow for easy modeling within a snake energy functional.
Detection of Particles
To detect particles in microscopy images we recently proposed [14] to extend the method of [19,13], because more flexibility with regard to the size of particles to be detected is required by the nature of stress granules.
Particles of a certain size are best represented in a particular interval of scales, therefore one interval is appropriate if all particles in the image share similar size characteristics.
However, if a single scale interval does not properly match the characteristics of the particles to be detected, irrelevant scales are included or important ones are excluded, causing incorrect sizes or shapes in detection, or leading to complete loss of particles. The use of multiple correlation images of – usually overlapping – scale intervals [an, bn] allows the detection of a particle in the scale interval which best suits its characteristics.
However, in many cases the particle is also found in adjacent intervals with incorrect size or shape. This causes spatially overlapping and thus competing particle hypotheses from different scale intervals, which induces the need to select the correct one. To tackle this problem, particle hypotheses which claim the same image regions are organized in a tree with leaves corresponding to the finest scale interval where hypotheses are available.
The amplitude-scale-invariant Bayes estimator [12] is applied to yield estimates of denoised wavelet coefficients. In addition, negative coefficients are set to zero as in our case bright particles on dark background are to be detected. Wavelet coefficients obtained after these processing steps are denoted by WWs (x, y).
The nature of the wavelet transform applied implies high redundancy of coefficients, which is exploited to identify particles using information of adjacent scales. The high correlation across scales motivates the combination of coefficients of adjacent scales s e
[a, b] to wavelet correlation images Coarser scale hypotheses appear closer to the root, because they are assumed to have larger support than the hypotheses of finer scales. This circumstance is illustrated in Fig. 2.
Fig. 2. Hypotheses trees for a clip of two correlation images
The concept of meaningful events [7] is employed to delete nodes of inferior hypotheses from the tree. This meaningfulness can be understood as the significance of an event produced by a background process or under a null hypothesis respectively, and has strong relation to statistical hypothesis testing. The background process H0 employed in our model assumes that no particle exists at the present location and correlation image values are caused by noise. Correlation values are supposed to be pairwise independent, which allows us to represent the probability P (Fi | H0) that particle Fi is produced by the background process H0 as
where
are random variables modeling the correlation image value for interval
observed at pixel position (x, y). We denote with p (Fi) the p-value of Fi, which is the probability to observe correlation values at least as extreme as the values of Fi produced by H0:
The tree of competing hypotheses is reduced comparing p-values. Hypotheses with smaller p-values are kept as they are assumed to be more unlikely caused by noise. To account for the difference in size of the support of hypotheses, the p-value of a node is normalized by the size of its support. Starting from the leaves, a parent’s normalized p-value is compared to the product of the normalized p-values of its children. The children of a hypothesis are removed from the tree, if the following equation is satisfied:
where IFiI denotes the support of hypothesis Fi. Otherwise the parent hypothesis is deleted and its children are connected to the hypothesis in the tree one level above their former parent.
Fig. 3. Examples of fluorescently labeled particles. (a) SG channel. (b) Detected SGs (black) using the scale-adaptive wavelet algorithm. (c) PB channel. (d) Detected PBs (black) using the same algorithm. Image intensity values in (b) and (d) are scaled for better visualization of detected particles.
The described method is able to select detection candidates from the correct scale interval. Examplary detection results for SGs and PBs are shown in Fig. 3 where SGs are hypothesized in intervals [2, 3] and [3,4] and selected from the appropriate interval according to their size.
Cell Segmentation
Active contour models for object segmentation have gained large interest in the computer vision community during the last years as they provide an integrated formal framework for combining image data with a priori knowledge. In practice there are two common ways to model active contours, either implicitly adopting level-set approaches or explicitly in terms of parametric contour models, i.e. snakes. Level-set techniques appear often superior compared to snake-based approaches with regard to the numerical tractability, however, preserving topology is usually easier to guarantee with snakes as topological stability is inherent in the underlying model theory. In our scenario of cell segmentation the number of objects is known in advance due to initial segmentation of nuclei (see Subsection 4.4). Accordingly the topology of the desired segmentation result is fixed and we prefer snakes for this application.
Introduction to Snakes
Explicit contour models rely on the definition of an object contour as a function c : R ^ R2 which maps a parameter value s e [0,1] defined along the contour c to points x in 2D space. In addition, an energy functional E(c(s)) is defined over the contour function c(s) which assesses the contour’s shape and positioning in the image with regard to the given object segmentation task. During the segmentation process the contour is evolved towards a local minimum of the energy, which gives the final segmentation result.
In general the energy functional of a snake model consists of at least two types of energy terms [15]. On the one hand we have an internal energy term Eint(c(s)) solely depending on the contour itself, e.g. on its geometric properties like length and curvature,
where a and β are weighting constants. On the other hand we have an external energy term which is usually derived from the input image to be segmented,
During segmentation the snake is iteratively moved to minimize the energy functional resulting from the sum of both terms:
While in the continous formulation above the energy terms are calculated by integration along the contour, in practice the snake is discritized to a list of points and line segments. Hence, the integrals in equations (6) and (7) transform to sums over all contour points.
The optimization of the snake’s shape and location given the energy functional is accomplished by applying combined implicit and explicit gradient descent techniques. To this end a time t is introduced to model the contour evolution process,
where γ denotes the step-size for the gradient descent step in each iteration. The derivation of the energy functional with regard to the contour function can be calculated applying the corresponding Euler-Lagrange equations for the snake energy functional [15].
Region Homogeneity Critérium
While the internal energy term presented in equation (6) may be seen as the quasistandard for internal snake energy formulation, the adaptation of the model to specific applications is usually accomplished by choosing appropriate external energy terms depending on the characteristics of the target objects. Often gradient-based energies are chosen, however, as the cells in our scenario are not characterized by significant gray-scale discontinuities along their boundaries such energies do not provide sufficient information to properly solve the segmentation task in our case. This is a significant difference from other related papers on cell segmentation for fluorescent microscope images as our cells show specific region intensity profiles (cf. Fig. 4 (a)).
Accordingly, the external energy used in our application is given by a term proposed in [6] which considers image intensity information of the interior region Rin(c(s)) enclosed by the snake and the exterior region Rout (c(s)) lying outside of the snake contour.
Fig. 4. Example image clip showing (a) dark nuclei surrounded by cell area and fluorescently labeled PBs, (b) the mask corresponding to the automatic segmentation result, and (c) the manually labeled ground-truth mask for comparison
In particular, for the both regions Gaussian-like intensity distributions are assumed, i.e. both regions are modeled by constant average intensities denoted μin and μο-ut, and deviations from these intensities are penalized as follows:
where I(x) is the pixel intensity at position x, and Xin and Xout are weighting constants encoding the expected variance within the regions.
For calculating the derivative of this energy term with regard to the contour, equation (10) first needs to be transformed applying Green’s theorem [20, pp. 20] to convert the region integrals into integrals along the contour. Subsequently Euler-Lagrange can again be applied which yields the following result:
where nc is the local normal vector of the contour pointing towards the exterior region. During optimization the term from equation (11) is evaluated for each 2-D point of the discretized snake contour. It pushes the local snake point either towards the interior or exterior region depending on which average intensity μη or μοut is closer to the local intensity value given the expected region variances. Accordingly, the whole snake is subsequently guided towards a location where its interior pixel value distribution is well described by μη, and the remaining exterior pixel value distribution by μοut. Note that μη and μοut are continously updated in each iteration to always provide an optimal approximation of the region intensity distributions with regard to the current segmentation.
Snake Coupling
The typical use case of snake-based segmentation assumes a single snake to be optimized within a given image at a single point in time. In the cell images of our application, however, there are multiple cells to be segmented simultaneously, i.e. information about the contour of one cell directly effects the contour of neighboring cells. In particular, as in reality cells do not overlap the segmentation should also avoid this.
Zimmer et al. proposed in [25] an approach for parallel optimization of several snakes that allows to penalize snake overlap by introducing a coupling energy term Eo between N individual snakes,
where ρ is again a weighting factor that allows to tune the influence of the overlap penalty on the snake energy as a whole. The corresponding derivative with regard to a single snake ci(s) is given by
with nc being again the local normal vector of the contour. tj basically counts the number of regions overlapping with the interior of snake ci at pixel position x.
We adopt the approach for penalizing overlap and include it into our snake energy functional which finally subsumes internal snake energies (Eq. 6), region homogeneity criteria (Eq. 10) for each individual snake and the image background outside of all snake contours, and the overlap penalty from equation (12):
Computing the derivatives of this global energy term with respect to each individual snake results in a set of N evolution equations. These are individually optimized applying gradient descent techniques as given in equation (9).
