2.1.

Overview

The image analysis approach was developed to classify MSC phenotype using phase-contrast micrographs of monolayer culture. Figure 1 shows the overall flow of the method; each step is detailed in the subsequent sections. Following preprocessing, the algorithm estimates cell density. Morphological operations and thresholding detect regions of the image that contain cells. Candidate markers are localized within these regions to identify if the cell region is a single cell or cell cluster. A cell cluster is further segmented to identify individual cells. After the segmentation of each cell, the algorithm extracts several human-engineered morphometric and textural features. Cell segmentation and feature extraction algorithms are built using the comprehensive set of reference-standard algorithms provided by the Image Processing Toolbox in MATLAB 9.5 (R2018b). Machine learning classifiers are trained using these features to distinguish between the RS and SR phenotype. Classification models were developed in Jupyter Notebook 6.0.1 using Python 3.5.6 libraries.

Fig. 1

Pipeline of algorithm to classify mesenchymal stem cells in phase-contrast micrographs. Example images illustrate (a) low and (b) moderate density. Cell ROI in image (b) include (c) individual cells and (d) clusters of cells, differentiated by the number of candidate markers (blue) inside RoIs (white).

2.2.

Dataset

Human bone marrow-derived MSCs were seeded at $100\text{cells}/{\mathrm{cm}}^2$ under standard conditions of expansion⁴⁰ and imaged on the second and fourth day after the culture was prepared to capture variation in phenotype as cells proliferate. For the MSCs employed in this study, a density of $\sim 1000\text{cells}/{\mathrm{cm}}^2$ is expected at day 2 and $\sim 6500\mathrm{per}{\mathrm{cm}}^2$ at day 4.¹² A Motic AE31 phase-contrast microscope with a $10\times$ objective and Moticam 1SP 1.0 MP camera was used to acquire culture images. All the images collected for this study have a size of $1280\times 1024\text{pixels}$ and a resolution of $1.56\text{pixels}/\mu \mathrm{m}$. Cell culture and image capture were repeated three times to generate the dataset for training and testing of the algorithm. Using Adobe Photoshop and Microsoft Paint software, cells were manually segmented and labeled as RS or SR phenotype by an individual with more than 15 years of experience in culturing MSCs. Images from two cultures served as ground truth for training, and images from the third culture were used for independent testing as detailed in Table 1. The training dataset of 71 images consisted of 472 cells with 307 cells labeled as RS and 165 cells as SR. The algorithm was validated with 36 phase-contrast micrographs having 186 cells with 121 RS cells and 65 SR cells. Each cell from the segmented ground truth was also characterized as RS or SR by 20 people trained to visually identify MSC phenotype to further analyze the generalizability of the method.

Table 1

Mesenchymal stem cell culture dataset.

2.3.

Image Preprocessing and Area Fraction Estimation

The input RGB phase-contrast micrograph is converted to grayscale ($I_{\text{gray}}$), shown in Fig. 2(a), and preprocessed to reduce the effect of undesired imperfections introduced during imaging. Contrast of $I_{\text{gray}}$ is adjusted to increase the intensity variation between the cells and the substrate, making cell regions more detectable for segmentation. Cell edges are sharpened using unsharp masking and then filtered using anisotropic median-diffusion to remove unwanted artifacts while improving the signal-to-background ratio without distorting edges [Fig. 2(b)].²⁷ The preprocessed image ($I_{\text{preprocessed}}$) is further processed by Sobel filter to identify cell boundaries in the images [Fig. 2(c)]. Sobel operator highlights regions with maximum intensity change, detecting edges above a sensitivity threshold of 1. Once object outlines are obtained, dilation and closing are performed to connect the detected edges. This is followed by flood-fill operation to remove holes from the filled regions.

Fig. 2

Example image of MSCs shows the steps involved in detecting initial cell regions. (a) The input phase-contrast micrograph is converted to grayscale and (b) preprocessed using contrast adjustment, sharpening, and anisotropic filtering. (c) Edges are detected using Sobel operator. (d) Edges are connected and filled using dilation, closing, and flood-fill operation to obtain initial cell regions ($I_{\text{initial}-\text{region}}$).

Using the same structuring element sizes for morphological operations to process images with a high number of cells and images with a low number of cells leads to poor segmentation performance. Apart from the number of cells and clusters, the variances in their size and shape also become significant as the culture grows over time. Using area fraction (AF) as a deciding factor for cell density addresses this problem as different parameters are used for low and moderate density levels to identify cell regions and the markers inside them. Also, every image is evaluated based on their density estimate rather than through the potentially erroneous assumption that duration in culture is a robust predictor of density. The cell density-based criterion automatically triggers optimizations in the algorithm so as to ensure comparable performance at a range of cell density levels.

AF of binary image ($I_{\text{initial}-\text{region}}$) is used to obtain an estimate of the input micrograph’s cell density. It is calculated by the algorithm as the percentage of white pixels in the image as given as

2.4.

Candidate Cell Region Detection

Cell region detection is conducted as a semantic segmentation⁴¹ step to identify pixels belonging to cells as defined in the truth, and hence, regions detected may contain more than one cell. Figure 3 shows the process of candidate cell region detection. For less dense images, $I_{\text{initial}-\text{region}}$ is used directly to define regions with potential cell objects. Preprocessing is designed to detect objects with high sensitivity [Fig. 3(a)]; therefore, each object is evaluated to remove image artifacts and identify candidate cell regions using thresholding [Fig. 3(b)] and morphological operations [Fig. 3(c)]. For size thresholding, the detected object is removed if the area is less than a threshold value determined by the minimum, maximum, mean, or the standard deviation of the area of all foreground objects in the image. This adaptive approach ensures that the threshold values are not overly biased toward the training set as it takes relative object sizes in each image into account to understand if it is likely to be a cell or not. Similarly, the intensity-based thresholds are calculated using the maximum and minimum intensity values inside that object obtained from the pixel positions in $I_{\text{preprocessed}}$. Objects labeled as cells in the training set contain bright pixels in the cytoplasm and/or dark pixels inside the nucleus after contrast adjustment. Lack of both indicates that the object is not a cell as it has a relatively uniform intensity range similar to substrate. For shape thresholding, circularity and ellipticity features of the object are calculated. From the training data, circularity (mean ± std. dev.) of MSCs and phase imaging artifacts have been found to be $0.43\pm 0.18$ and $0.82\pm 0.05$, respectively, where 1 represents a perfect circle. Ellipticity is measured as

Fig. 3

Illustration of the steps involved in identifying candidate cell regions. (a) $I_{\text{initial}-\text{region}}$ obtained after edge detection and morphological operations are used as the input. (b) Thresholding using size, intensity, and shape criteria removes detected objects that are not cells. (c) Opening and erosion optimize the shape of cell areas. (d) Clearing the image border removes incomplete cells resulting in the final candidate cell regions. (Borders of cell regions are highlighted in blue.)

For moderately dense images, the edge detection step is performed with a reduced sensitivity threshold of 0.5 to cover all cell edges. Dilation and closing are performed with different structuring element sizes to get new potential cell regions. The same thresholding methods as less dense images are carried out on these objects with different threshold values to keep them relevant for images with more objects. Thresholding is also repeated more times compared to low cell density processing as a lower edge detection threshold can cause more false detections. Thresholding is followed by morphological opening, closing, and border clearing to obtain the $I_{\text{cell}-\text{region}}$ for images having moderate cell density.

2.5.

Candidate Marker Detection

In phase-contrast micrographs of cells, intensity is brightest at the cell boundaries where the phase shift is maximum due to the optical path difference (refractive index and thickness) between cells and substrate, and is darker within the cell due to relative uniformity within the cell.⁴² The darkest regions inside cells are taken as candidate cell markers to identify individual cells and to segment cells in clusters since each cell contains one prominent regional minimum. The image $I_{\text{preprocessed}}$ is further processed using Gaussian filtering and median filtering to remove unnecessary noise and false local minima that do not belong to the marker. Then, contrast-limited adaptive histogram equalization improves the contrast of the regional minimum. Finally, morphological reconstruction is performed using the histogram equalized image along with $I_{\text{cell}-\text{region}}$ as the mask to obtain $I_{\text{marker}-\text{processed}}$, shown in Fig. 4(a).

Fig. 4

Pipeline for segmentation of markers inside the RoIs (shown in yellow) from $I_{\text{cell}-\text{region}}$. (a) $I_{\text{preprocessed}}$ is processed to get an image ($I_{\text{marker}-\text{processed}}$) suitable for detecting markers. (b) Markers are obtained from $I_{\text{marker}-\text{processed}}$ using H-minima transform with higher threshold A. (c) Markers are obtained from $I_{\text{marker}-\text{processed}}$ using H-minima transform with lower threshold B. (d) Markers from both the thresholds are merged for ROIs satisfying the perimeter threshold criterion and also for ROIs with zero markers from threshold A.

Markers are segmented using two different threshold values (A and B) for H-minima transform.^43–45 Used independently, the higher threshold value of A results in under-detection of markers, and the lower threshold value of B results in false positives (FPs); therefore, the two are combined. First, a very high value is used as threshold A for H-minima transform to localize potential markers inside the ROIs from $I_{\text{cell}-\text{region}}$. Morphological opening and binary area opening are done to remove objects that aren’t the regional minimum [Fig. 4(b)]. Then, minima are obtained with a lower threshold, value B. These minima outputs are dilated and closed to get the candidate markers [Fig. 4(c)]. In the case of zero markers from threshold A in any ROI, markers detected using B are added to that region. Markers from both the thresholds are also merged for potential cluster ROIs using perimeter as a criterion [Fig. 4(d)].

Over-detection error is managed for cell regions with more than one marker using distance thresholding along with area thresholding and morphological operations. The Euclidean distance between the centroid of the markers within a cell region is computed, and if the distance is too small, the regional maximum is over-segmented. Thus, one of the two close markers with the smaller area is removed. Finally, dilation and erosion are performed to get rid of over-segmentations.

2.6.

Cell Segmentation and Validation

ROIs and markers are combined as shown in Fig. 5(a) for the final instance segmentation⁴¹ step to detect and delineate each cell in the image. A region with no marker is not considered a cell, only one marker is labeled as a single cell [Fig. 5(b)], and a region with more than one marker is treated as a cluster of cells [Fig. 5(c)]. The marker count within each cluster indicates the number of cells in the cluster. Each cluster region is segmented into individual cells [Fig. 5(d)] using marker-controlled watershed, which overcomes the limitations of standard watershed technique, such as over and under-segmentation, using markers.^46,47

Fig. 5

Representation of the cell segmentation steps. (a) $I_{\text{cell}-\text{region}}$ is integrated with the results from marker detection (Markers shown in blue). (b) ROIs with single markers are identified as individual cells, and (c) ROIs with more than one marker are identified as clusters of cells. (d) Marker-controlled watershed is carried out to segment individual cells within clusters (watershed ridge lines shown in red). (e) Results from the steps (b) and (d) are combined to get the segmentation output.

The algorithm’s performance for cell detection and segmentation was validated using the expert-defined truth. For cell detection, the true positives (TP) are given by the number of individual cells and cells inside clusters that are correctly detected by the algorithm. The false negatives (FNs) are of two types, namely, missed detections (FN1), where individual cells are not detected, and under-detections (FN2), where the number of detected cells within a cluster is less than the actual number. Similarly, the two types of FPs are false detections (FP1) and over-detections (FP2). FP1 are debris or image artifacts that are erroneously detected as individual cell objects, and FP2 is the number of extra objects detected within a cluster by the algorithm. Sensitivity (S) and precision (P) were computed using

The presented algorithm was also assessed by comparing its cell detection and segmentation results with the conventional U-Net,⁴⁹ a deep learning technique for biomedical image segmentation. The U-Net architecture was modeled based on the standard implementation proposed by Ronneberger et al. The conventional U-Net was trained and tested on the same datasets to avoid any bias. Sensitivity, precision, and DICE evaluation metrics were used to compare performance of the algorithm reported here to the standard U-Net.

2.7.

Feature Extraction

FP objects (FP1 and FP2) and all cells belonging to under-detected (FN2) clusters were excluded from the dataset to avoid training or validating classification models using wrong objects. The final dataset used for feature extraction and classification is detailed in Table 2.

Table 2

Mesenchymal stem cell culture dataset used for feature extraction and classification.

Human engineered descriptors of the correctly detected cells were automatically extracted for the task of classifying each cell into RS or SR phenotype. A total of 30 features, consisting of a combination of size, shape, and first-order and second-order statistical texture measures, were computed.^50–52 The names of the features are listed in Table 3. The morphometric features were extracted to distinguish spindle-shaped RS cells from flattened SR cells, whereas texture-based features were extracted to use spatial distribution of intensity for differentiating between SR cells that are flattened and RS cells that have a more prominent phase-contrast halo around their cell body.^16,40 Each first-order feature was calculated for segmented cell regions in $I_{\text{gray}}$, $I_{\text{preprocessed}}$, and $I_{\text{reconstructed}}$ giving three measures for each first-order feature. The second-order features were measured for the gray-level co-occurrence matrix (GLCM) of $I_{\text{gray}}$, $I_{\text{preprocessed}}$, and $I_{\text{reconstructed}}$. Also, each GLCM feature was computed in 24 different orientations yielding 72 measures for each second-order feature. Only one measure out of the three and 72 first and second-order feature measures, respectively, was selected by finding the measure with the highest area under the curve (AUC) of the receiver operating characteristic (ROC) curve^53,54 for distinguishing between RS and SR cells.

Table 3

Human-engineered features extracted for the segmented MSCs.

The 30 features were sorted in descending order of their AUC value before computing the correlation matrix to ensure that features with higher AUCs were retained. Features with correlation $>0.8$ were removed to reduce redundancy and optimize the computation for the classifier.

2.8.

Cell Classification and Validation

The selected features were used to train linear and non-linear classifiers to find the most suitable model for our application. Features were transformed to have zero mean and unit variance before training, and the validation/test data were scaled using the training parameters. As seen in Table 1, the dataset is imbalanced with RS cells being more prevalent than the SR cells. To avoid poor classification of SR cells due to its low prevalence, synthetic minority oversampling technique⁵⁵ was used during training.

Linear kernel support vector machine (LSVM), radial-basis kernel support vector machine (RSVM), linear discriminant analysis (LDA), $K$-nearest neighbor (KNN), and logistic regression (LR) models were trained to classify MSCs as RS or SR phenotype. The models were trained using features from both day 2 and day 4 together, and day 2 and day 4 individually. As the feature correlation and relevance for data from day 2, day 4, and combination of day 2 and 4 would be different, features were selected for each of them individually. It was observed that the classifier could learn better from features of cells whose images were acquired on the same day rather than the combination of features of cells from two different days. AUC values were obtained using ROCKIT software,⁵⁶ and all the models were compared using the average AUC value from five-fold cross-validation.⁵⁷ Further, ensemble classifiers based on soft-voting method were also trained with the top two classifiers for day 2 and day 4 separately.⁵⁸ Finally, the performance of the selected models for the independent test dataset was evaluated using AUC, sensitivity, and specificity. Sensitivity and specificity were determined by selecting a threshold that minimizes $(1-\text{sensitivity}{)}^2+{(1-\text{specificity})}^2$.⁵⁸

3.1.

Cell Detection and Segmentation

The algorithm’s ability to accurately locate MSCs was evaluated using sensitivity and precision. Table 4 shows the breakdown of correctly detected cells and undetected cells for RS and SR cells in the training and independent testing sets. FNs are further defined as missed detections (FN1) or under-detections (FN2) as described in Sec. 2.6. The algorithm correctly detected cells with sensitivity $>0.95$ for both RS and SR cell phenotypes in the training set. A sensitivity over 0.8 in the independent testing set indicates the robustness of the algorithm.

Table 4

Cell detection sensitivity of the algorithm for training and independent testing.

Precision is the fraction of objects detected by the algorithm that were identified as cells in the ground truth labeling of the dataset. Table 5 shows the number of incorrectly identified objects that are either false detections (FP1) or over-detections (FP2). The algorithm’s precision for cultures used in training was above 0.95. A precision $>0.85$ in the independent testing set demonstrates the generalizability of the algorithm to detect cells with a low FP rate.

Table 5

Cell detection precision of the algorithm for training and independent testing.

The DICE metric was used to evaluate the overlap between the algorithm and the ground truth segmentation of MSCs. DICE scores (mean ± std. dev.) for the training and testing dataset are summarized in Table 6 for all cells outlined in the truth and only for the cells that were correctly detected. The algorithm segmented over 85% and 80% of the manually defined cell areas for training and testing, respectively, regardless of cell phenotype.

Table 6

Cell segmentation DICE score of the algorithm for training and independent testing.

The performance of the algorithm for cell detection and segmentation were analyzed per each cell in Tables 3–5. Sensitivity, precision, and DICE metrics of the algorithm were further examined per image from all cultures and days. The mean ± std. dev of these metrics for the training and testing images are given in Table 7, and it was confirmed that the algorithm could detect and segment cells consistently for each image.

Table 7

Cell detection and segmentation performance of the algorithm per image from training and testing datasets. Mean and standard deviations are calculated for all cells over all images.

Figure 6 shows the ground truth along with algorithm’s segmentation results for two images of both low and moderate density from the test dataset. The difference in performance of the algorithm to localize cells in low and moderately dense images was also reviewed as it was trained with different parameters for these two levels of cell densities. Welch’s $t$-test at 95% confidence level failed to show statistical difference in the values of sensitivity ($p=0.375$), precision ($p=0.191$), and DICE ($p=0.289$) for low and moderately dense cell images.

Fig. 6

Sample images of MSCs from the independent test dataset comparing the ground truth cell outlines (top row) with the segmentation results of the algorithm (bottom row). Columns (a) and (b) are examples of low-density images whose AF values were estimated as 0.04 and 0.08, by the algorithm, respectively. Similarly, examples of moderate density images with AF values of 0.15 and 0.21 are shown in column (c) and (d), respectively. The ground truth images also illustrate that the truncated cells on image borders are not considered for analysis.

The cell detection and segmentation results of the algorithm are compared against the conventional U-Net in Table 8 for both training and testing data. The presented algorithm outperforms the conventional method in all of the evaluation metrics for this task. The low precision of the U-Net was determined to be due to detection of several image artifacts as cells. The algorithm was able to overcome such FP detections through various morphological thresholding steps. Moreover, although the U-Net model used weighted maps to learn the small separation borders, it could not separate cells in clusters as effectively as the marker-controlled watershed technique used in the presented pipeline. Additional work would be needed to modify and optimize U-Net’s architecture for it to be able to localize MSCs more accurately. Thus, it was statically demonstrated using paired two-sample $t$-tests that the algorithm could detect and segment cells significantly better than the state-of-the-art method.

Table 8

Comparison of algorithm’s cell detection and segmentation performance with conventional U-Net.

3.2.

Cell Phenotype Classification

As mentioned previously, only correctly detected individual cells and cells in clusters with correct cell count were used for training and testing the machine learning models for the task of classifying each cell phenotype as RS or SR. LSVM, RSVM, LDA, KNN, and LR were trained and compared using the AUC metric from five-fold cross-validation. These five classifiers trained using object features from both day 2 and day 4 cultures were validated by testing their performance in cross-validation for a combination of cells from “day 2 + day 4,” and day 2 and day 4 independently. It was observed that day 2 + day 4 models did not perform as well for day 2 as they did for day 4 (Table 9). This may be due to models’ bias toward a higher prevalence of day 4 MSC features. Training the classifiers with features from day 2 and day 4 individually improved their performance for both days (Table 9). The outputs from the top two classifiers (RSVM and LR for day 2; LDA and KNN for day 4) based on cross-validation AUC were combined for ensemble classification. Figure 7 shows the classifier agreement plot between the models used for ensemble classification. The $x$ axis in Fig. 7(b) is not continuous due to the discrete probability distribution of the KNN classifier. The disagreement between the top 2 classifiers for each day may be because they learn differently from the same features. “RSVM + LR” and “LDA + KNN” models being a combination of linear and non-linear approaches are able to harness the potential of both the models to make more accurate predictions. This likely explains why fusion classifiers perform better than the individual classifiers during five-fold cross-validation, and hence, they are selected for the image analysis pipeline.

Table 9

Five-fold cross-validation to compare performance of models for distinguishing between RS and SR cells. AUC values are in bold for Day 2 and Day 4 classifiers that perform best during cross-validation.

Fig. 7

Diagonal classifier agreement plots between (a) logistic regression classifier ($x$ axis) and radial kernel support vector machine ($y$ axis) for day 2 culture and (b) KNN ($x$ axis) and linear discriminant analysis ($y$ axis) for day 4 culture. Each point represents a cell whose phenotype is predicted as RS (green) or SR (red) by each classifier.

The selected ensemble classifiers for day 2 and day 4 are further evaluated using the independent test dataset. The algorithm was able to correctly identify cell phenotypes with an AUC of 0.82 for day 2 and 0.79 for day 4. The classification models for both days have been statistically proven to perform better than random chance using the 95% confidence interval for AUCs given in Table 10. These inferences were not corrected for multiple comparisons as only a single statistical test was performed for each day. The fitted ROC curves obtained for the day 2 and day 4 ensemble classifiers are shown in Fig. 8. The algorithm achieves a sensitivity and specificity of over 0.75 for both days, further demonstrating its effectiveness in classifying MSCs based on their phenotype.

Table 10

Performance of the algorithm for classifying MSCs as RS or SR for independent testing.

Fig. 8

Fitted binormal ROC curves illustrating the performances of day 2 and day 4 ensemble classifiers for the task of distinguishing between RS and SR cells during independent testing. The dashed blue line represents the RSVM + LR day 2 classifier and the solid orange line represents the LDA + KNN day 4 classifier.

3.3.

Generalizability of Phenotype Classification

The cells in the training and test dataset were also labeled by a group of 20 individuals with varied level of expertise in working with MSCs. The group comprised of nine people with one to three years of experience, six people with four to nine years of experience, and five people with ten or more years of experience. This analysis was performed to study the generalizability of the developed algorithm as phenotype assessment of MSCs is highly subjective. Each cell was classified as either RS (0) or SR (1) by all 20 individuals, and its average phenotype score (APS) was computed. Based on the range of average scores, the cells were split into categories such as highly RS, moderately RS, uncertain, moderately SR, and highly SR (Table 11). The confidence of the group in identifying a cell as either RS or SR was represented by these five categories; where the highly RS/SR categories meant highest confidence in labeling that cell and uncertain category meant that the phenotype of that cell was almost indeterminate.

Table 11

Classification sensitivity of the algorithm for the five cell phenotype categories.

The output of the ensemble classifiers for day 2 and day 4 from training and testing set was studied to understand which category of the cells were most wrongly predicted by the algorithm. It was noticed that the algorithm had the highest classification sensitivity ($>0.9$) for highly RS and highly SR categories and the least sensitivity (as low as 0.4) for the uncertain category (Table 11). This enabled us to understand that while our image analysis pipeline has the potential to predict the cell phenotype with a sensitivity very close to human interpretation, it did exhibit a similar weakness to trained humans in that it had lower sensitivity in categorizing marginal morphological characteristics. If these limitations are addressed in future iterations of the algorithm, this approach has the potential to out-perform human observers with decades of experience.