Automatic Quantitative Segmentation of Myotubes Reveals Single-cell Dynamics of S6 Kinase Activation

Abstract

Automatic cell segmentation is a powerful method for quantifying signaling dynamics at single-cell resolution in live cell fluorescence imaging. Segmentation methods for mononuclear and round shape cells have been developed extensively. However, a segmentation method for elongated polynuclear cells, such as differentiated C2C12 myotubes, has yet to be developed. In addition, myotubes are surrounded by undifferentiated reserve cells, making it difficult to identify background regions and subsequent quantification. Here we developed an automatic quantitative segmentation method for myotubes using watershed segmentation of summed binary images and a two-component Gaussian mixture model. We used time-lapse fluorescence images of differentiated C2C12 cells stably expressing Eevee-S6K, a fluorescence resonance energy transfer (FRET) biosensor of S6 kinase (S6K). Summation of binary images enhanced the contrast between myotubes and reserve cells, permitting detection of a myotube and a myotube center. Using a myotube center instead of a nucleus, individual myotubes could be detected automatically by watershed segmentation. In addition, a background correction using the two-component Gaussian mixture model permitted automatic signal intensity quantification in individual myotubes. Thus, we provide an automatic quantitative segmentation method by combining automatic myotube detection and background correction. Furthermore, this method allowed us to quantify S6K activity in individual myotubes, demonstrating that some of the temporal properties of S6K activity such as peak time and half-life of adaptation show different dose-dependent changes of insulin between cell population and individuals.

Key words: time lapse images, cell segmentation, fluorescence resonance energy transfer, C2C12, myotube

Introduction

Live-cell fluorescence imaging provides the advantage of observing a biological system function with spatiotemporal resolution. With recent advances in imaging techniques and fluorescence biosensors, it has become clear that variability in intracellular states induces heterogeneity in the dynamic behavior of individual cells exposed to a uniform stimulus (Hughey et al., 2015; Snijder and Pelkmans, 2011; Tomida et al., 2015; Yao et al., 2016). Therefore, single-cell studies have uncovered the importance of understanding cellular dynamics and heterogeneity.

For single-cell studies, there is considerable need to develop algorithms that enable automatic segmentation and signal intensity quantification in individual cells from a large number of time-lapse images. Automatic segmentation methods have thus far been developed mainly for round shape and mononuclear cells (Bajcsy et al., 2015; Kodiha et al., 2011). For example, the marker based watershed segmentation is widely studied and used for efficient cell segmentation (Chalfoun et al., 2014; Hodneland et al., 2016; Wong et al., 2010). Typical use of the watershed segmentation is to flood a region from the stained nuclei as an initial flooding source until it touches the neighbors or cellular boundaries. However, such a method that requires a nucleus as a marker cannot be used directly for polynuclear cells, such as differentiated C2C12 myotubes.

C2C12 cells, derived from mice myoblasts, have been used as model cell lines for the study of cell differentiation and muscle functions in vitro (Fodor et al., 2017; Rahar et al., 2018; Rapizzi et al., 2008; Soltow et al., 2013). By differentiation induction of C2C12 cells, a large subpopulation of myoblasts forms elongated polynuclear myotubes, while a subpopulation of the cells, called reserve cells, remains undifferentiated (Fig. 1) (Yoshida et al., 1998). When combined with heterogeneity in the cell population, the elongated polynuclear form of the myotubes makes it difficult to perform automatic quantitative segmentation of individual myotubes using conventional methods developed for round shapes and mononuclear cells without a heterogeneous cell population. The elongated polynuclear form property of myotubes leads to over-segmentation by direct use of watershed segmentation. In addition, reserve cells that occupy spaces between the myotubes, make it difficult to automatically identify background regions.

Fig. 1

Differentiated C2C12 myotubes stably expressing Eevee-S6K. Differentiated C2C12 cells stably expressing Eevee-S6K (Yellow). Nuclei were stained with DAPI (Blue). Because myotubes are much thicker than reserved cells, myotubes (red arrows) were much brighter than reserve cells (green arrows). Because of the relatively weaker fluorescence signal in reserve cells, only nuclear signals, but not cell body signals, were visible.

Here we developed an automatic quantitative segmentation method of a myotube using a watershed segmentation of summed binary images and two-component Gaussian mixture model (GMM). We established a C2C12 cell line that stably expressed Eevee-S6K, a fluorescence resonance energy transfer (FRET) biosensor, which monitor S6 kinase (S6K) activity, and acquired time-lapse fluorescence images both of cyan fluorescent protein (CFP) and FRET-induced yellow fluorescent protein (FRET-YFP) in differentiated C2C12 myotubes (Aoki and Matsuda, 2009). We made each image of FRET-YFP binary and summed the images. This procedure resulted in an enhanced contrast between myotubes and reserve cells, thereby permitting selective detection of myotubes. In addition, to identify individual myotubes, we used a myotube center, rather than a nucleus, as a marker of individual myotubes. To identify the myotube center, we applied a distance transform to the binary images of FRET-YFP. A distance transform is an operation that converts a binary image to a distance map image where all pixels have a value corresponding to the Euclidean distance to the nearest boundary pixel (Leymarie and Levine, 1992). We generated binary distance map images and summed the images, permitting enhancement of myotube centers. Using a myotube center instead of a nucleus, individual myotubes were segmented automatically by watershed segmentation.

Removing the noise signal from background regions is an essential prerequisite for obtaining quantitative measurements. However, it is difficult to identify background regions because of the presence of reserve cells that occupy the spaces between the myotubes. We used a two-component GMM that fit a fluorescence intensity histogram and estimated background intensity without manual selection of the background regions. The use of the two-component GMM ensures quantitative and objective measurement of signal intensity.

Thus, in this study, we provide a framework for automatic quantitative segmentation of individual myotubes that combines watershed segmentation using automatic myotube segmentation and background correction. Furthermore, this method allowed us to quantify S6K activity in individual myotubes, demonstrating that some of the temporal properties of S6K activity such as peak time and half-life of adaptation show different dose-dependent changes of insulin between cell population and individuals.

Results

Automatic quantitative segmentation of elongated polynuclear myotubes using time lapse images

We developed an automatic quantitative segmentation method of a myotube from fluorescence time lapse images of differentiated C2C12 cells stably expressing Eevee-S6K (Aoki and Matsuda, 2009). The automatic quantitative segmentation method consists of three steps, Step I, II and III. Step I is pre-processing, where the use of median filter and white top-hat filter reduce variation of fluorescence intensity of the images. Step II is myotube segmentation, where the use of summed binary time-lapse images enhance myotubes and myotube centers, rather than reserve cells. Step III is background correction where we used a two-component GMM that fit a histogram of fluorescence intensity and estimated mean background intensity without manual selection of the background regions. The parameters used in each Steps were summarized in Table I.

Table I Parameters used for pre-processing

Parameter	Value (pixel)	Description
median filtera	5	side length of filter window
white top hatb	35	side length of filter window
white top hatc	30	side length of filter window
threshold sized	10000	threshold area of myotube

(a) is the common value in detection of myotubes and detection of myotube centers. (b) is used for detection of myotubes. (c) is used for detection of myotube centers. (d) is used after watershed segmentation.

Step I: Pre-processing of time lapse images

Differentiated C2C12 cells are polynuclear and elongated in shape. Therefore, a myotube region in an obtained image shows a non-uniform distribution of fluorescence. For proper segmentation of myotubes, image smoothing of obtained images is essential. We first applied a median filter to remove outliers (Huang et al., 1979) and then applied a white top-hat filter to reduce variation of fluorescence intensity (Fig. 2).

Fig. 2

Pre-processing of fluorescence images. In Step I-i, we applied median filter with 5×5 square window to the raw FRET-YEP time-lapse images. In Step I-ii, we applied white top-hat to images processed by the median filter. Filter window of the white top-hat is 35×35 and 30×30 in detection of myotubes and detection of myotube centers, respectively.

Step II: Segmentation method of myotubes

In Step II, we detected myotubes and myotube centers using thresholding methods including the triangle method and Otsu’s method (Otsu, 1979; Zack et al., 1977) (Fig. 3, Fig. 4), and the individual myotubes were segmented by watershed segmentation using the detected myotubes and myotube centers (Bajcsy et al., 2015) (Fig. 5). Step II consisted of three groups that contained 11 substeps in total; Step II-i to II-iii, detection of myotubes; Step II-iv to II-x, detection of myotube centers; Step II-xi, watershed segmentation.

Fig. 3

Detection of myotubes. (A) Detection of myotubes consisted of three substeps; Step II-i, first binarization; Step II-ii, summation; Step II-iii, second binarization. In step II-ii, the intensity increases from blue to red. (B) FRET-YFP intensity histogram at frame 0. Red line indicates the threshold value, which was determined by the triangle method. (C) Intensity histogram of the summed binary image. The red line indicates the threshold value, as determined by Otsu’s method. (D) Labeled images of detected myotubes. One color corresponds to one continuous myotube region. Black denotes a region of either reserve cells or background.

Fig. 4

Detection of myotube centers. (A) Detection of myotube centers consists of seven substeps; Step II-iv, first binarization of fluorescence time-lapse images; Step II-v, distance transform of the first binarized images; Step II-vi, second binarization of the distance map images; Step II-vii, Summation of the second binary images; Step II-viii, third binarization of the summed images; Step II-ix, labeling of the third binary image; Step II-x, denoising of the labeled image. Note that in Step II-v and II-vii, colors denote intensity in a pixel, whereas, in Step II-ix and II-x, one color corresponds to one continuous region. Black indicates either reserve cells or background regions. (B) Skeletonized image of a labeled image. The red arrow head indicates short length debris. (C) Length histogram of a skeletonaized image. We assumed that length is a number of pixels in a continuous region of a skeletonized image. The red line indicates threshold length determined by the triangle method.

Fig. 5

Watershed segmentation. (A) Myotubes and myotube centers were detected from time-lapse FRET-YFP images. A detected myotube was used as a boundary for watershed segmentation. A detected myotube center was used as a marker for watershed segmentation. Watershed segmentation was performed using detected myotubes and myotube centers, areas that were less than 10000 were removed. (B) Overlays of raw FRET-YFP time-lapse image at frame 0 and the result of watershed segmentation. Yellow lines indicate contours of segmented myotube regions.

Step II-i to II-iii: Detection of myotubes. Detection of myotubes consists of three substeps, Step II-i to II-iii (Fig. 3A); Step II-i, first binarization of FRET-YFP time-lapse images; Step II-ii, summation of the first binarized images; Step II-iii, second binarization of the summed image.

In Step II-i, we detected candidate regions of myotubes in each image by making them binary (Step II-i first binarization in Fig. 3A). Since the fluorescence intensity histogram of the FRET-YFP was skewed to the left, the triangle method was used to determine the threshold for the first binarization (Fig. 3B) (Zack et al., 1977). In addition, in order to prevent fluctuation of the threshold value due to the number of pixels, we standardized the total area of the intensity histogram to be one prior to using the triangle method. However, the binary images included myotubes and fragmented regions of reserve cells. To restrict this to only myotubes, enhanced contrast between myotubes and reserve cells will be needed.

In Step II-ii, we summed the binary images (Step II-ii Summation in Fig. 3A). Because myotubes were thicker and brighter than reserve cells, the summation of binary images allowed the contrast between myotubes and reserve cells to be enhanced.

In Step II-iii, we made the summed image binary to detect myotubes (Step II-iii Second binarization in Fig. 3A). Since the intensity histogram of the summed image was bimodal, Otsu’s method was used to determine the threshold for the second binarization (Fig. 3C). Otsu’s method is a statistical thresholding method to find a threshold that minimizes the within-class variance (Otsu, 1979).

Thus, we detected myotubes by successive thresholding methods, including the triangle method and Otsu’s method, using the summed binary image of FRET-YFP time-lapse images (Fig 3D). However, we still detected some of the myotubes as one continuous region rather than as individual myotubes. Further segmentation will be needed to detect individual myotubes.

Step II-iv to II-x: Detection of myotube centers. For segmentation of individual cells, watershed segmentation has been used most frequently in fluorescence imaging (Bajcsy et al., 2015). In watershed segmentation, a stained cell nucleus often is used as a marker of a cell because the majority of the cell type is mononuclear and round. However, conventional watershed segmentation methods using a nucleus as a marker of an individual cell are not designed to identify polynuclear and elongated shape cells such as myotubes. Therefore, we used a myotube center, instead of a nucleus, as a marker of a myotube for watershed segmentation.

Detection of myotube centers consists of seven substeps, Step II-iv to II-x (Fig. 4A); Step II-iv, first binarization of FRET-YFP time-lapse images; Step II-v, distance transform of the first binary images; Step II-vi, second binarization of the distance map images; Step II-vii, summation of the second binary images; Step II-viii, thrird binarization of the summed images; Step II-ix, labeling of the third binary image; Step II-x, denoising of the labeled image.

In Step II-iv, we made binary images using the triangle method for detection of candidate regions of myotubes.

In Step II-v, we transformed the binary images to distance map images that emphasized myotube centers (Step II-v Distance transform in Fig. 4A). However, some regions of the myotube centers were failed to be emphasized due to the presence of nuclei, and therefore were over-segmented.

In Step II-vi, the distance map images were made binary using the triangle method (Step II-vi Second binarization in Fig. 4A), and then summed to emphasize myotube centers (Step II-vii Summation in Fig. 4A). Because the nuclear positions in a myotube were changed frequently, the summed image attenuated the influence of nuclear existence.

In Step II‑viii, the summed image was made binary using Otsu’s method for detection of myotube centers (Step II-viii Third binarization in Fig. 4A).

In Step II-ix, we labeled the binary image with one continuous region as one myotube center (Step II-ix Labeling in Fig. 4A). Although the labeled image seemed to show a successful detection of myotube centers, it contained lots of small debris ranging in length from one to dozens of pixels (Fig. 4B and C). Therefore, we removed debris by the short length, where a length was defined as a number of pixels of a skeletonized region. Skeletonization is a thinning algorithm that removes outer pixels of a region to find its medial axis. For accurate myotube segmentation, these debris need to be properly removed.

Since the length histogram of the myotube centers was skewed to the left, Step II-x used the triangle method to determine the threshold (Fig. 4C). We detected individual myotube centers by removing regions smaller than the threshold (Step II-x Denoising in Fig. 4A, Table I).

Throughout the first ten substeps (Step II-i to II-x), we obtained an individual myotube region as a boundary for watershed segmentation, and an individual myotube center as a marker for watershed segmentation. Thus, the method using a myotube center instead of a nucleus enabled us to segment individual myotubes by the watershed segmentation.

Step II-xi: Watershed segmentation. In Step II-xi, we applied watershed segmentation to the myotubes detected in Step II-iii using the detected myotube center in Step II-x as a marker. Individual myotubes were segmented by removing regions smaller than a threshold (Fig. 5).

Evaluation of segmentation

In this study, implementation of the triangle method and Otsu’s method enabled an automatic segmentation of myotubes. However, because both the triangle method and Otsu’s method are based on the intensity histogram of an image, the segmentation results could be affected by excitation intensity and number of time-lapse images. Therefore, we compared the segmentation performance of our method to other traditional segmentation methods (Otsu’s method and triangle method) with various transmittance of excitation light and the various number of time-lapse images (Fig. 6). Ground truth of myotubes were generated from the MIP image by manual segmentation (Fig. S1), and Jaccard index (Levandowsky and Winter, 1971) was calculated between ground truth and each method. The Jaccard index is an indicator of the similarity between two regions (Bajcsy et al., 2015). The higher value of the index, the more similarity. Note that, like the ground truth, the Otsu’s method and the triangle method performed segmentation using MIP images. It was an advantageous performance comparison to Otsu method and triangle method. Also note that it is difficult to prepare perfect ground truth for cells with ambiguous boundaries between cells.

Fig. 6

Accuracy of segmentation. (A) Left panel: Representative segmentation results with the indicated transmittance of excitation light (25%, 12%, 6% and 3%). Yellow lines are contours of individual myotubes. Right panel: Representative segmentation results when changing the number of time-lapse images by resampling every two to ten images. Full shows segmentation results using the full set of the images. Every 2, Every 5 and Every 10 indicate segmentation results obtained by resampling the images at every two, five and ten time points, respectively. Transmittance was set to 25%. (B) Jaccard indices in each condition in (A) (median±iqr, n=8 stage positions). Segmentation result of ground truth was used as the reference to calculate the Jaccard index. N. S. (Not significant), p>0.05 (Welch’s t-test). p values were corrected by Bonferroni correction.

Our method showed significantly higher Jaccard index than the traditional segmentation methods with 12 % or less transmittance of excitation light (Fig. 6A left panel, 6B left panel). This result indicates that our method is more robust to change in excitation intensity than the traditional segmentation methods. For reduction of number of images, there is no significant difference under every conditions. (Fig. 6A right panel, 6B right panel). These results suggest that the segmentation performance of our method is comparable to the traditional segmentation method against the number of images used. In addition, our method has an advantage that adjacent myotubes can be separated by implementing detection of myotube center as a marker for watershed segmentation.

Step III: Background correction method

Due to various factors, including the setup of the optical system, properties of the detector, and the fluorescent probe, a fluorescence image often contains background intensity variation in the same field. In addition, live-cell fluorescence imaging often is performed with weak excitation light to prevent photo-bleaching and photo-toxicity, which generates an image with a low signal to noise ratio. In ratiometric data, such as the FRET ratio, failure to identify of background regions leads to significant artifacts in signal intensity quantification. Therefore, proper background correction is needed for quantitative signal acquisition (Zimmermann et al., 2003).

Manual selection of background regions near the regions of interest (ROIs) has been used widely for background correction (Ceelen et al., 2007; Horie et al., 2015; Rahar et al., 2018). However, in the case of differentiated C2C12 cells, it is difficult to manually identify background regions, because reserve cells occupy spaces between myotubes (Fig. 1). Here, we used a two-component Gaussian mixture model (GMM) to estimate background regions.

Background correction consisted of five substeps (Fig. 7); Step III-i, maximum intensity projection (MIP) of FRET-YFP time-lapse images, Step III-ii, binarization of the MIP image; Step III-iii, NOT AND (NAND) operation on the raw time lapse fluorescence images and binary MIP image; Step III-iv, estimation of background intensity using two-component Gaussian mixture model (GMM). Step III-v, signal intensity quantification of individual myotubes.

Fig. 7

Background correction. Background correction consists of five substeps; Step III-i, maximum intensity projection (MIP); Step III-iii, NOT AND (NAND) operation; Step III-iv, two-component Gaussian mixture model (GMM). Step III-v, signal intensity quantification. In Step III-iii, the intensity increases from blue to red. In Step III-iv, a red line indicates estimated background intensity distribution and a red dashed line indicates an average of the distribution as the estimated background intensity. A yellow line indicates the estimated intensity distribution of a region that included reserve cells and dead cells. In Step III-v, each blue line corresponds to the time series of each individual myotube. A red line indicates the mean time series of the FRET ratio of myotubes in each background correction.

MIP in Step III-i is a processing technique that keeps only the pixels of maximum intensity along the z-axis of stack images.

In Step III-ii, binarization of MIP images derived from FRET-YFP time-lapse images enabled reliable separation of myotubes from other regions composed of reserve cells and background regions (Step III-ii in Fig 7).

In Step III-iii, by performing the NAND operation on the raw time-lapse fluorescence images and the binary MIP image, we extracted a region composed of reserve cells and background regions from each time lapse image (Step III-iii in Fig. 7).

In Step III-iv, we estimated background intensity both of CFP and FRET-YFP time-lapse images using two-component GMM (Step III-iv in Fig. 7). In the two-component GMM, the histogram was divided into two components of Gaussian mixture distributions which corresponded to signal distribution of a region that included reserve cells and background regions, and background intensity was estimated as an average of the lower component.

In Step III‑v, we subtracted the estimated background intensity from each time point of the CFP and FRET-YFP time-lapse images, quantified signal intensities of CFP and FRET-YFP in individual myotubes (Fig. S2, right panel), and acquired a time series of the FRET ratio (FRET-YFP/CFP) (Step III-v in Fig. 7).

Evaluation of background correction

For evaluating accuracy of background correction using two-component GMM, we compared the accuracy of background correction using two-component GMM with other automatic background correction methods using raw histogram (RAW) and kernel density estimation (KDE) (Fig. 8). In a background correction using RAW, background intensity was estimated as the mode intensity of the histogram. In background correction using KDE, background intensity was estimated as the mode intensity of the histogram, approximated by kernel density estimation.

Fig. 8

Comparison of accuracy of background corrections using RAW, KDE and two-component GMM. (A) Upper panels: time series of the FRET-ratio (FRET-YFP/CFP) of individual myotubes in RAW, KDE, and GMM, respectively (n=96, eight stage positions). Here, we used the same time lapse fluorescence images for all background corrections using RAW, KDE, and two-component GMM. One blue line corresponds to the time series of one myotube. The red line indicates the mean time series of FRET ratio of myotubes in each background correction. Lower panels: absolute first order differential time-series of the FRET ratio of individual myotubes in each background correction. The red line indicates the mean time series of the absolute first order differential time series. (B) AUC distributions of the absolute first order differential time-series of myotubes in each background correction. *, p<0.05 (Steel-Dwass test). In each violin plot, box plots are shown in the inset, and a white dot denotes the median.

Steps III‑i to III‑iii are common steps for all background corrections using RAW, KDE and two-component GMM. In Step III-iv, background intensity was estimated using RAW, KDE and two-component GMM, respectively. In Step III-v, time series of the FRET ratio (FRET-YFP/CFP) was acquired using each background correction method (Fig. 8A).

For evaluating accuracy of background corrections using RAW, KDE and two-component GMM, we calculated a quantification noise of the time series of FRET ratio in a myotube at a frame index i as d_i and a total quantification noise of a myotube as AUC, described by

  

$$d_i=\left|y_i-y_{i-1}\right|,$$(1)

  

$$AUC=\sum _{i=1}^{N-1}\frac{\text{Δ}t}{2}\left(d_i+d_{i+1}\right),$$(2)

where ${\left\{y_i\right\}}_{i=1}^N$ is the time series of the FRET ratio at a frame index i, N is the total number of frames, Δt is the time interval between frames, $\{{d_i\}}_{i=1}^{N-1}$ is the absolute first order differential time series between y_i and y_i–1 (Fig. 8A, lower panel), and AUC is the area under the curve of d_i in Eq (1).

The median of the AUC in background correction using two-component GMM was significantly smaller than those using RAW or KDE (Fig. 8B), indicating that the two-component GMM gives the lowest quantification noise among the three methods. In addition, GMM has advantages, in that the estimated distribution was continuous and the parameters including the mean and the variance in the Gaussian component were determined automatically.

We compared the accuracy of background corrections for manual estimation and two-component GMM (Fig. S1 and Fig. 9). In manual estimation, we generated the MIP image from FRET-YFP time-lapse images and selected the regions of myotubes and background visually (Fig. S1). For each selected region, we quantified signal intensities of CFP and FRET-YFP, and calculated a time series of the FRET ratio (Manual in Fig. 9A). In two-component GMM, the FRET ratio was calculated from intensities of CFP and FRET-YFP in the manually selected regions (Manual+two-component GMM in Fig. 9A). There was no significant difference in distributions of the AUCs among the three background correction methods, indicating that our proposed method using two-component GMM enabled provided objective and reasonable quantification compared with manual estimation (Fig 9B).

Fig. 9

Comparison of accuracy of background correction using Manual and Manual+two-component GMM. (A) Upper panels: time series of the FRET ratio (FRET-YFP/CFP) of individual myotubes in Manual and Manual+two-component GMM, respectively (n=80, eight stage positions). The result of background correction using automatic myotube segmentation and two-component GMM in Fig. 6 was shown (two-component GMM). One blue line corresponds to the time series of one myotube. The red line indicates the mean time series of the FRET ratio of myotubes in each background correction. Lower panels: absolute first order differential time series of the FRET ratio of individual myotubes in each background correction. The red line indicates the mean time series of the absolute 1st order differential time series. (B) AUC distributions of the absolute first order differential time-series of myotubes in each background correction. N. S. (Not significant), p>0.05 (Steel-Dwass test). In each violin plot, box plots are shown in the inset, and a white dot denotes the median.

In an intramolecular FRET biosensor such as Eevee-S6K, CFP and YFP are linked by a linker domain, and the total noises of the time series of CFP and FRET-YFP should show a high correlation in a steady state. Therefore, we examined the coefficient of determination of the AUCs between CFP and FRET-YFP (Fig. 10). The coefficients of determination of the AUCs in background correction using RAW, KDE and two-component GMM were 0.822, 0.854, and 0.924, respectively (Fig. 10, upper panels). Similarly, the coefficients of determination of AUCs in Manual and Manual+two-component GMM were 0.899 and 0.892, respectively (Fig. 10, lower panels). Thus, background correction using two-component GMM shows the highest coefficient of determination of the AUCs of the absolute first order differential time series between CFP and FRET-YFP, indicating that two-component GMM is the most reasonable method for background correction.

Fig. 10

Coefficients of determination between time-lapse images of CFP and FRET-YFP in each quantification method. One dot corresponds one myotube. The black line is the regression of AUCs of CFP and FRET-YFP.

We next examined whether our proposed method using two-component GMM can be used for myotubes under other conditions, such as myotubes stimulated with insulin or myotubes expressing an ATP probe stimulated with an electrical pulse stimulation (EPS) (Fig. 11). First, we used insulin-stimulated differentiated C2C12 cells stably expressed Eevee-S6K (Fig. 11A, upper panels). Individual myotubes were identified and quantified by manual estimation or by our proposed method using GMM (Fig. 11B, upper panels). We calculated the correlation coefficient between the time series of the FRET ratio quantified by manual estimation and two-component GMM in each region where the Jaccard index was larger than 0.5 (Fig. 11C, upper panel). Most of the correlation coefficients were larger than 0.98, indicating that our proposed method using two-component GMM is comparable to the manual estimation. In addition, we used EPS-stimulated differentiated C2C12 cells stably expressed mitAT1.03 which is FRET biosensor for monitoring ATP concentration in a mitochondrion (Imamura et al., 2009) (Fig. 11A, lower panel), and individual myotubes were quantified (Fig. 11B, lower panels). Similar to Eevee-S6K, all of the correlation coefficients were larger than 0.98 (Fig. 11C, lower panel), indicating that our proposed method using two-component GMM is comparable to manual estimation.

Fig. 11

Quantification of signal activity of S6K and ATP concentration by applying our proposed method. (A) Time lapse of the FRET ratio image. Differentiated C2C12 cells stably expressing Eevee-S6K were stimulated with 100 nM of insulin at 50 min. Differentiated C2C12 cells stably expressing mitAT1.03 were stimulated with EPS (10 ms with 50 V, 1Hz interval) at 10 min and continued for 15 min. (B) A time series of the FRET-ratio (FRET-YFP/CFP) in response to insulin (n=90, eight stage positions) and EPS (n=10, one stage position) quantified by two-component GMM and quantified by Manual (n=80, eight stage positions, insulin; and n=10, one stage position, EPS). One blue line corresponds one myotube. The red dashed line indicates time points of insulin stimulation. The red filled area indicates a period of EPS. (C) Histograms of Pearson’s correlation coefficients of the time series of the FRET ratio between two-component GMM and Manual in response to insulin (Upper panels) and EPS (Lower panels).

The differences between population and individual myotubes

Our proposed framework aims to quantify time series of fluorescence time lapse images in living individual myotube objectively. Using the framework, we tried to quantify insulin-dependent S6K activation and looked for different properties between cell population and individual myotube.

Myotubes were stimulated with various doses of insulin (0 nM to 100 nM) and time series of S6K activity in individual myotubes were quantified from the obtained fluorescence time lapse images (Fig. 12A). We defined properties of the time series such as Peak, Peak time, AUC, etc., from the time series (Fig. 12A lower right panel, B). For cell population, Peak, AUC, Half-life of adaptation and Intensity at half-life of adaptation increased in a dose-dependent manner of insulin. Variance in Peak time decreased, and Adaptation precision slightly decreased in a dose-dependent manner of insulin.

Fig. 12

Insulin-dependent S6K activation in single myotube and cell population. (A) Time series of the FRET-ratio (FRET-YFP/CFP) in response to various doses of insulin in individual single myotubes. From Control to 100 nM, number of myotubes were n=98, 102, 100, 97, 86, 103 and 98, respectively. One blue line corresponds to the time series of one myotube. The red line indicates the mean series of each time series. Lower right panel is definition of the properties of time series. (B) The cell population response of time series in (A). In each violin plot, box plots are shown in the inset, and a white dot denotes the median. The red line indicates the median series of each properties.

In order to investigate the difference in properties between the cell population and the individuals, we performed correlation analysis with bootstrap subset (Bootstrap) as the cell population and all data (All) as the individuals (Fig. 13). The correlation of Peak time with all other properties were significantly different between the cell population and the individuals. Also, the correlation of Half-life of adaptation with all other properties were significantly different between the cell population and the individuals. These results suggest that Peak time and Half-life of adaptation show different responses between the cell population and individuals. The reason why the Peak time was significantly different in the combinations with all properties may be because the variance of Peak time decreased in a dose-dependent manner of insulin (Fig. 12B). Consistent with this, the correlations between Peak time with all other properties were decreased in a dose-dependent manner of insulin (Fig. S3).

Fig. 13

Correlation between the properties in individual myotubes and population. (A) Spearman’s rank correlation coefficients between the properties. All (n=684) and Bootstrap (n=700) indicate the correlations in all individuals and bootstrap subsets of population of myotubes, respectively. The bootstrap subsets in each dose data were generated by iterating 100 times to randomly sample 10 points and calculate the median. Cyan rectangles; the correlations higher in the cell population than in the individuals, magenta rectangles; the correlation higher in the individuals than in the cell population, yellow rectangles; the reversed correlation between the cell population and individuals. (B) The difference in correlation coefficients between All and Bootstrap. N. S. (Not significant), *p<10^–3, **p<10^–5, ***p<10^–7. z score is a test static of the difference between two correlation coefficients (see Materials and methods).

We further analyzed Half-life of adaptation that shows significant differences from all other properties except Peak time. Differences in properties between the cell population and individuals can be classified into three groups that show higher correlation in the cell population than the individuals, higher correlation in the individuals than the cell population, and reversed relationship between the cell population and individuals (Fig. 13). In the combination of Half-life of adaptation, Peak and AUC showed stronger correlation with Half-life of adaptation in the cell population than in the individuals (Fig. 13A, B). Consistent with this, the distributions of Peak and AUC were more separated in the cell population than the individuals (Fig. S4A). Adaptation precision showed higher correlation with Half-life of adaptation in the individuals than in the cell population (Fig. 13A, B). In addition, the correlation of Half-life of adaptation with Adaptation precision in each dose decreased in a dose-dependent manner of insulin (Fig. S3 bottom center). This suggests that Half-life of adaptation and Adaptation precision in individuals differently responded in a dose-dependent manner of insulin. The correlation between intensity at Half-life of adaptation with Half-life of adaptation in the cell population was positive, whereas that in the individuals was negative (Fig. 13A, B). The distribution of Half -life of adaptation and intensity at Half -life of adaptation for each dose varied in the negative correlation direction in individuals, but not in the cell population (Fig. S4C). Furthermore, the correlations of Half-life of adaptation with intensity at Half-life of adaptation in each dose of insulin were negative in the individuals (Fig. S4 bottom right). These results indicate that individuals possess hidden properties that can not be seen in cell population. Thus, some of the properties of insulin dose-dependent S6K activation differ between the cell population and individuals.

Discussion

Advantages of automatic myotube segmentation

Automatic segmentation of cells in fluorescence images has been developed for mononuclear and round shape cells (Bajcsy et al., 2015; Kodiha et al., 2011). However, an automatic segmentation method for elongated polynuclear cells, such as differentiated C2C12 myotubes has yet to be developed. Indeed, intracellular signal activity and reactive oxygen species in myotubes have been measured in vitro only by manual ROI selection of myotubes (Ceelen et al., 2007; Horie et al., 2015; Rahar et al., 2018). Since our method does not make assumptions on the number of nuclei and cell shape, our method can be used for any type of polynuclear and deformed cells. However, since our method does not implement cell tracking, it is difficult to apply it to moving cells. Using our method for images every fixed period and implementing cell tracking, there is a possibility that our method can be applied to moving cells.

Some recent studies using deep learning showed improvements of accuracy and throughput of segmentation of cells (Van Valen et al., 2016). Deep learning is a powerful tool for segmentation of cells with high accuracy and high throughput, and may be used effectively for elongated polynulcear cells, such as C2C12 myotubes. Despite versatility of deep learning for any type of cells, deep learning requires an enormous training dataset and computational resources in general. Especially, in cells with ambiguous cell boundary such as myotube, deep learning may learn bias of manual segmentation. By contrast, our proposed method does not require any training and few computational resources.

The calculation time of our method was less than 15 min even when processing 131 images including from segmentation to quantification, using AMD Ryzen 7 1800X (Table II). When only the segmentation was performed, our method took about 30 sec for one dataset (131 frames). This is fast enough compared with manual segmentation which took about 15–20 min using Roi tool in Fiji. In background correction, two-component GMM took about 3 sec per image. This is practical time, because our method was fully automated after setting the segmentation parameters.

Table II Computation time

Step	Mean±S.D. (sec)	Max (sec)
Detection of myotubes (Step II-i to II-iii)a	6.13±0.52	6.82
Detection of myotube centers (Step II-iv to II-x)b	24.97±3.72	32.33
Watershed (Step II-xi)c	—	—
Background correction (Step III)d	2.94±0.84	4.73

Because (a), (b) and (c) are batch processes, the computation times were measured for each dataset (1344×1024 pixels, 131 frames) with eight stage positions (n=8). (c) was too fast to measure the computation time. Because of (d) is a sequential process for each image, the computation time was measured in each image (131 frames, eight stage positions, n=1048).

Heterogeneity of cell response

We quantified S6K activity in individual myotubes and demonstrated that some of the properties of S6K activity in a dose-dependent manner of insulin differed between cell population and individuals. This raises the possibility that some of the properties may be underestimated or incorrectly estimated by the cell population analysis. Moreover, single-cell analysis will address whether variation of cell population is derived from intra-cellular variation or inter-cellular variation.

The activity of other signaling molecules is also known to differ between individual cells. For example, in a mitogen-activated protein kinase (MAPK), Asynchronous oscillation of p38 MAPK response to external stimulation was revealed by single-cell analysis (Tomida et al., 2015). In the case of calcium response to ATP stimulation, the variance of calcium response was revealed that caused by internal state of the cells (Yao et al., 2016).

Taken together, single-cell analysis is important to reveal hidden properties of cell system that can not be observed in cell population. Signal transduction in skeletal muscle has thus far been studied by conventional biochemical and molecular biological methods such as western blot, that reflect the activity at cell population level (Areta et al., 2014; Baar and Esser, 1999; Dalle Pezze et al., 2016; Deldicque et al., 2008; West et al., 2016). In contrast, our proposed method can automatically detect and quantify signal activity of single myotube, and will open the door to single-cell analysis in signal transduction of skeletal muscle.

Materials and Methods

Construction of C2C12 cell line stably expressing fluorescent biosensors

C2C12 cells (kindly provided by Takeaki Ozawa, University of Tokyo, Tokyo, Japan) were cultured at 37°C under 5% CO₂ in Dulbecco’s modified Eagle’s medium (DMEM), (High Glucose) with L-Glutamine and Phenol Red (Wako Pure Chemical Industries Limited, Osaka, Japan) supplemented with 10% fetal bovine serum (Nichirei Bioscience Incorporated, Japan). C2C12 cells stably expressing FRET biosensors, Eevee-S6K (plasmid was kindly provided by Kazuhiro Aoki, National Institute for Basic Biology, Aichi, Japan) (Aoki and Matsuda, 2009; Komatsu et al., 2011) and mitAT1.03 (plasmid was kindly provided by Hiromi Imamura, Kyoto University, Kyoto, Japan) (Imamura et al., 2009), were constructed using the PiggyBac Transposase System (System Biosciences, U.S.A.), respectively. Three hundred μL of Opti-MEM (Life technologies, U.S.A.), 4 μL of Lipofectamin 2000 (Invitrogen, U.S.A.), 1.0 μg of PiggyBac transposon vector clone (kindly provided by Kazuhiro Aoki, National Institute for Basic Biology, Aichi, Japan) and 0.2 μg of PiggyBac transposase expression vector (PB210PA-1, Funakoshi, Japan) were mixed and let stand for 5 min. Thereafter, 70% confluent C2C12 cells were plated on a 35 mm dish and transfected with the mixture and incubated for 6 h. For selection of infected cells, the cells were cultured with DMEM (High glucose) with L-Glutamine and Phenol Red containing 10% fetal bovine serum and 20 μg/mL of Blasticidin S Hydrochloride (Wako, Japan). Selected cells were seeded on a Cell Culture Dish 430167 (Corning Incorporated, U.S.A.) and cultured until forming the colonies. The colonies were picked and seeded on a Cell Culture Dish 430167. After seeding and proliferation, the cells were stored at a concentration of 1.0×10⁵ cells/mL with Bambanker (NIPPON Genetics, Japan). Eevee-S6K is localized in the cytosol and mitAT1.03 is localized in the mitochondria.

Differentiation induction of C2C12 cells

For differentiation induction, C2C12 cells were plated on 35 mm/Collagen Coated Dish Collagen type I (IWAKI, Japan) at a concentration of 1.0×10⁵ cells/dish and cultured for two days under conditions described above until confluent, and then confluent cells were cultured in DMEM with L-Glutamine and Phenol Red supplemented with 2% horse serum (Nichirei Bioscience Incorporated) for seven days (Odelberg et al., 2000).

Fluorescence time lapse imaging

C2C12 myotubes were starved in 2 mL of Medium199, Hanks’ Balanced Salts (Life technologies, U.S.A.) overnight and then mineral oil (Sigma Aldrich) was stratified to prevent vaporization of the medium prior to the fluorescence time-lapse imaging. Fluorescence time-lapse imaging was performed with an inverted fluorescence microscope, IX 83 (Olympus, Tokyo, Japan) equipped with a UPLSAPO10X2 objective lens (Olympus), a ORCA-R2 C10600-10B CCD camera (Hamamatsu Photonics), a U-HGLGPS mercury lamp (Molecular Devices, Sunnyvale, CA), XF3075 and XF3079 emission filter for CFP and YFP (Omega Optical), respectively, and an MD-XY30100T-META automatically programmable XY stage (Molecular Devices).

All the images of myotubes were 1344×1024 pixels and 0.645 μm/pixel. Time-lapse images of myotubes expressing Eevee-S6K were acquired for 650 min every 5 min (total 131 frames) with eight stage positions. Time-lapse images of myotubes expressing mitAT1.03 were acquired for 45 min every 1 min (total 46 frames) with one stage position. After background subtraction, FRET-YFP/CFP ratio images were created with Meta-Morph software (Universal Imaging, West Chester, PA).

Image analysis

Our proposed method was developed using Python and Python modules, including Mahotas (Coelho, 2013), Scikit-image (van der Walt et al., 2014) and Scikit-learn (Pedregosa et al., 2011) for automatic quantification. Fiji (Schindelin et al., 2012) was used for manual quantification of fluorescence intensity of FRET-YFP and CFP. The FRET-YFP and CFP intensities were averaged over the whole myotube area.

Significance of the difference between two correlation coefficients

Fisher’s z-transformation of correlation coefficients r is given by

  

$$z_i=\frac{1}{2}\ln \frac{1+r_i}{1-r_i},i\in \left\{\mathrm{1,2}\right\},$$(3)

where i is an index of correlation coefficient. A z score of the difference between two correlation coefficients, described by

  

$$z=\frac{z_1-z_2}{\sqrt{\frac{1}{n_1-3}+\frac{1}{n_2-3}}},$$(4)

where n₁ and n₂ are numbers of data corresponding to r₁ and r₂. As the z score follows normal distribution, significance test was performed with significance level α=10^–3, 10^–5, 10^–7.

Availability of data and materials

The datasets of time-lapse images, segmentation results and segmentation programs are available from our database (http://kurodalab.org/info/Inoue).

Acknowledgments.

We thank Kazuhiro Aoki for kindly providing FRET biosensor, Eevee-S6K.

We thank Shinsuke Uda for helpful discussions and technical advice. We thank our laboratory members Masashi Fuji and Atsushi Hatano for their critical reading of this manuscript, helpful discussions, and technical assistance with the experiments.

This work was supported by the Creation of Fundamental Technologies for Understanding and Control of Biosystem Dynamics, CREST, of the Japan Science and Technology Agency (JST)(#JPMJCR12W3) and by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number (17H06300, 17H06299, 18H03979). K. Kunida receives funding from a Grant-in-Aid for Young Scientists (B) (#16K19028).

References

•  Aoki,  K. and  Matsuda,  M. 2009. Visualization of small GTPase activity with fluorescence resonance energy transfer-based biosensors. Nat. Protoc., 4: 1623–1631.
•  Areta,  J.L.,  Hawley,  J.A.,  Ye,  J.-M.,  Chan,  M.H.S., and  Coffey,  V.G. 2014. Increasing leucine concentration stimulates mechanistic target of rapamycin signaling and cell growth in C2C12 skeletal muscle cells. Nutr. Res., 34: 1000–1007.
•  Baar,  K. and  Esser,  K. 1999. Phosphorylation of p70(S6k) correlates with increased skeletal muscle mass following resistance exercise. Am. J. Physiol., 276: C120–7.
•  Bajcsy,  P.,  Cardone,  A.,  Chalfoun,  J.,  Halter,  M.,  Juba,  D.,  Kociolek,  M.,  Majurski,  M.,  Peskin,  A.,  Simon,  C.,  Simon,  M.,  Vandecreme,  A., and  Brady,  M. 2015. Survey statistics of automated segmentations applied to optical imaging of mammalian cells. BMC Bioinformatics, 16: 330.
•  Ceelen,  K.K.,  Oomens,  C.W.J.,  Stekelenburg,  A.,  Bader,  D.L., and  Baaijens,  F.P.T. 2007. Changes in Intracellular Calcium during Compression of C2C12 Myotubes. Exp. Mech., 49: 25–33.
•  Chalfoun,  J.,  Majurski,  M.,  Dima,  A.,  Stuelten,  C.,  Peskin,  A., and  Brady,  M. 2014. FogBank: a single cell segmentation across multiple cell lines and image modalities. BMC Bioinformatics, 15: 6598.
•  Coelho,  L. 2013. Mahotas: Open source software for scriptable computer vision. J. Open Res. Softw., 1: e3.
•  Dalle Pezze,  P.,  Ruf,  S.,  Sonntag,  A.G.,  Langelaar-Makkinje,  M.,  Hall,  P.,  Heberle,  A.M.,  Razquin Navas,  P.,  van Eunen,  K.,  Tölle,  R.C.,  Schwarz,  J.J.,  Wiese,  H.,  Warscheid,  B.,  Deitersen,  J.,  Stork,  B.,  Fäßler,  E.,  Schäuble,  S.,  Hahn,  U.,  Horvatovich,  P.,  Shanley,  D.P., et al. 2016. A systems study reveals concurrent activation of AMPK and mTOR by amino acids. Nat. Commun., 7: 13254.
•  Deldicque,  L.,  Sanchez Canedo,  C.,  Horman,  S.,  De Potter,  I.,  Bertrand,  L.,  Hue,  L., and  Francaux,  M. 2008. Antagonistic effects of leucine and glutamine on the mTOR pathway in myogenic C2C12 cells. Amino Acids, 35: 147–155.
•  Fodor,  J.,  Gomba-Tóth,  A.,  Oláh,  T.,  Almássy,  J.,  Zádor,  E., and  Csernoch,  L. 2017. Follistatin treatment suppresses SERCA1b levels independently of other players of calcium homeostasis in C2C12 myotubes. J. Muscle Res. Cell Motil., 38: 215–229.
•  Hodneland,  E.,  Tai,  X.C., and  Kalisch,  H. 2016. PDE Based Algorithms for Smooth Watersheds. IEEE Trans. Med. Imaging, 35: 957–966.
•  Horie,  M.,  Warabi,  E.,  Komine,  S.,  Oh,  S., and  Shoda,  J. 2015. Cytoprotective role of Nrf2 in electrical pulse stimulated C2C12 myotube. PLoS One, 10: 1–13.
•  Huang,  T.,  Yang,  G., and  Tang,  G. 1979. A fast two-dimensional median filtering algorithm. IEEE Trans. Acoust., 27: 13–18.
•  Hughey,  J.J.,  Gutschow,  M.V.,  Bajar,  B.T., and  Covert,  M.W. 2015. Single-cell variation leads to population invariance in NF-B signaling dynamics. Mol. Biol. Cell, 26: 583–590.
•  Imamura,  H.,  Nhat,  K.P.,  Togawa,  H.,  Saito,  K.,  Iino,  R.,  Kato-Yamada,  Y.,  Nagai,  T., and  Noji,  H. 2009. Visualization of ATP levels inside single living cells with fluorescence resonance energy transfer-based genetically encoded indicators. Proc. Natl. Acad. Sci. USA, 106: 15651–15656.
•  Kodiha,  M.,  Bański,  P., and  Stochaj,  U. 2011. Computer-based fluorescence quantification: a novel approach to study nucleolar biology. BMC Cell Biol., 12: 25.
•  Komatsu,  N.,  Aoki,  K.,  Yamada,  M.,  Yukinaga,  H.,  Fujita,  Y.,  Kamioka,  Y., and  Matsuda,  M. 2011. Development of an optimized backbone of FRET biosensors for kinases and GTPases. Mol. Biol. Cell, 22: 4647–4656.
•  Levandowsky,  M. and  Winter,  D. 1971. Distance between Sets. Nature, 234: 34–35.
•  Leymarie,  F. and  Levine,  M.D. 1992. Fast raster scan distance propagation on the discrete rectangular lattice. CVGIP Image Underst., 55: 84–94.
•  Odelberg,  S.J.,  Kollhoff,  A., and  Keating,  M.T. 2000. Dedifferentation of mammalian myotubes induced by msx1. Cell, 103: 1099–1109.
•  Otsu,  N. 1979. A Threshold Selection Method from Gray-Level Histograms. IEEE Trans. Syst. Man. Cybern., 9: 62–66.
•  Pedregosa,  F.,  Varoquaux,  G.,  Gramfort,  A.,  Michel,  V.,  Thirion,  B.,  Grisel,  O.,  Blondel,  M.,  Müller,  A.,  Nothman,  J.,  Louppe,  G.,  Prettenhofer,  P.,  Weiss,  R.,  Dubourg,  V.,  Vanderplas,  J.,  Passos,  A.,  Cournapeau,  D.,  Brucher,  M.,  Perrot,  M., and  Duchesnay,  É. 2011. Scikit-learn: Machine Learning in Python. J. Mach. Learn. Res., 12: 2825–2830.
•  Rahar,  B.,  Chawla,  S.,  Pandey,  S.,  Bhatt,  A.N., and  Saxena,  S. 2018. Sphingosine-1-phosphate pretreatment amends hypoxia-induced metabolic dysfunction and impairment of myogenic potential in differentiating C2C12 myoblasts by stimulating viability, calcium homeostasis and energy generation. J. Physiol. Sci., 68: 137–151.
•  Rapizzi,  E.,  Donati,  C.,  Cencetti,  F.,  Nincheri,  P., and  Bruni,  P. 2008. Sphingosine 1-phosphate differentially regulates proliferation of C2C12 reserve cells and myoblasts. Mol. Cell. Biochem., 314: 193–199.
•  Schindelin,  J.,  Arganda-Carreras,  I.,  Frise,  E.,  Kaynig,  V.,  Longair,  M.,  Pietzsch,  T.,  Preibisch,  S.,  Rueden,  C.,  Saalfeld,  S.,  Schmid,  B.,  Tinevez,  J.-Y.,  White,  D.J.,  Hartenstein,  V.,  Eliceiri,  K.,  Tomancak,  P., and  Cardona,  A. 2012. Fiji: an open-source platform for biological-image analysis. Nat. Methods, 9: 676–682.
•  Snijder,  B. and  Pelkmans,  L. 2011. Origins of regulated cell-to-cell variability. Nat. Rev. Mol. Cell Biol., 12: 119–125.
•  Soltow,  Q.A.,  Zeanah,  E.H.,  Lira,  V.A., and  Criswell,  D.S. 2013. Cessation of cyclic stretch induces atrophy of C2C12 myotubes. Biochem. Biophys. Res. Commun., 434: 316–321.
•  Tomida,  T.,  Takekawa,  M., and  Saito,  H. 2015. Oscillation of p38 activity controls efficient pro-inflammatory gene expression. Nat. Commun., 6: 8350.
•  van der Walt,  S.,  Schönberger,  J.L.,  Nunez-Iglesias,  J.,  Boulogne,  F.,  Warner,  J.D.,  Yager,  N.,  Gouillart,  E., and  Yu,  T. 2014. scikit-image: image processing in Python. PeerJ, 2: e453.
•  Van Valen,  D.A.,  Kudo,  T.,  Lane,  K.M.,  Macklin,  D.N.,  Quach,  N.T.,  DeFelice,  M.M.,  Maayan,  I.,  Tanouchi,  Y.,  Ashley,  E.A., and  Covert,  M.W. 2016. Deep Learning Automates the Quantitative Analysis of Individual Cells in Live-Cell Imaging Experiments. PLOS Comput. Biol., 12: e1005177.
•  West,  D.W.D.,  Baehr,  L.M.,  Marcotte,  G.R.,  Chason,  C.M.,  Tolento,  L., Gomes, A.V,  Bodine,  S.C., and  Baar,  K. 2016. Acute resistance exercise activates rapamycin-sensitive and -insensitive mechanisms that control translational activity and capacity in skeletal muscle. J. Physiol., 594: 453–468.
•  Wong,  L.C.,  Lu,  B.,  Tan,  K.W., and  Fivaz,  M. 2010. Fully-automated image processing software to analyze calcium traces in populations of single cells. Cell Calcium, 48: 270–274.
•  Yao,  J.,  Pilko,  A., and  Wollman,  R. 2016. Distinct cellular states determine calcium signaling response. Mol. Syst. Biol., 12: 894.
•  Yoshida,  N.,  Yoshida,  S.,  Koishi,  K.,  Masuda,  K., and  Nabeshima,  Y. 1998. Cell heterogeneity upon myogenic differentiation: down-regulation of MyoD and Myf-5 generates “reserve cells”. J. Cell Sci., 111(Pt 6): 769–779.
•  Zack,  G.W.,  Rogers,  W.E., and  Latt,  S.A. 1977. Automatic measurement of sister chromatid exchange frequency. J. Histochem. Cytochem., 25: 741–753.
•  Zimmermann,  T.,  Rietdorf,  J., and  Pepperkok,  R. 2003. Spectral imaging and its applications in live cell microscopy. FEBS Lett., 546: 87–92.