2.1.

Cell Culture Techniques

NMuMG cells (provided by Dr. Theresa Vincent) were grown in Dulbecco’s modified Eagle’s medium (Sigma-Aldrich, Czech Republic) supplemented with GlutaMAX™ (Life Technologies, Czech Republic), 10% fetal bovine serum (Sigma-Aldrich, Czech Republic), $100\mathrm{U}/\mathrm{ml}$ penicillin, and $0.1\mathrm{mg}/\mathrm{ml}$ streptomycin (Life Technologies, Czech Republic). The cells were harvested by trypsinization and transferred into 10 sterilized observation chambers $\mu$-Slide I (Ibidi GmbH, Germany) at $50{\text{cells}/\mathrm{mm}}^2$ to ensure a low density for segmentation of individual cells. Once the cells were transferred to the observation chambers, they were kept in the incubator and imaged the day after. The 10 chambers were divided into control and $\mathrm{TGF}\beta$-treated ($10\mathrm{ng}/\mathrm{ml}$) with imaging beginning immediately after treatment.

2.2.

CCHM

Quantitative phase imaging of cells was performed by CCHM,^14,24,25 now also available as Q-Phase (TESCAN ORSAY HOLDING, a.s., Brno, Czech Republic). The optical setup of the microscope is based on Mach–Zehnder type interferometer modified for incoherent, off-axis holographic microscopy as shown in Fig. 1(b). The illumination is formed by a low-coherence source (halogen lamp) while the beam is split into two separated optical pathways—reference and object arm. Both arms contain matching condensers, objectives, and tube lenses. In the reference arm, the diffraction grating is placed in order to ensure the achromatic formation of the interference pattern (hologram) in the output plane. The hologram is recorded by the CCD camera and numerically reconstructed using a Q–Phase software (TESCAN ORSAY HOLDING, a. s., Brno, Czech Republic). The numerical reconstruction of the image is based on carrier removal in the Fourier plane.²⁶ The hologram is Fourier transformed using the 2-D fast Fourier transform (FFT) algorithm. The image spectrum in extracted by a windowing operation, whereas the window is centered at the carrier frequency. The frequency origin is translated to the center of the window, and the 2-D inverse FFT is applied to obtain the complex amplitude. Amplitude and phase are derived from the complex amplitude as modulus and argument, respectively. Since the values in the raw phase image are wrapped on the interval ($-\pi ,\pi$), the phase unwrapping algorithm²⁷ is applied. After the reconstruction, the image can still be burdened by the optical aberrations of the imaging system, imperfect adjustment of the microscope, or possibly by surrounding temperature changes. This issue is solved by the subtraction of the compensation surface described in detail in Ref. 28. In this way, a final unwrapped and compensated phase image is obtained. Such reconstructed QPI is proportional to the optical path difference of the two arms according to the following equation:

Fig. 1

Imaging in CCHM. (a) Model of an adhered cell in the observation chamber imaged by CCHM. (b) Optical setup of CCHM: light source (S), relay lens (L), beamsplitters (BS), condensers (C), specimen (SP), reference object (RO), microobjectives (O), tube lenses (TL), diffraction grating (DG), output lenses (OL), output plane (OP), and detector (D).

The use of incoherent illumination enables strong suppression of coherent noise and parasitic interferences. Moreover, the low illumination power of the source ($0.2\mu \mathrm{W}/{\mathrm{cm}}^2$) is not likely to influence the physiological functions of imaged cells, making CCHM very convenient for live cell imaging.

2.3.

Image Acquisition

The NMuMG cells were imaged by CCHM. During the experiment, the samples were illuminated with the halogen lamp through the interference filter ($\lambda =650\mathrm{nm}$, 10 nm FWHM). Microscope objectives (Nikon Plan Fluor $20\times /0.5$) were utilized for the imaging, providing the field of view $140\mu \mathrm{m}$. For the purpose of classification, it was essential to acquire a reasonably large number of cells undergoing EMT. Therefore, six fields of view were imaged with a 5-min interval, each $\sim 1\mathrm{mm}$ apart from each other. Each chamber was imaged for 48 h in the presence or absence of $\mathrm{TGF}\beta$ to obtain the time-lapse QPI for the classification. The media were not changed during the imaging period and conditions within the microscope mimicked that of the cell incubator (temperature 37°C) to ensure the cells were not subjected to stress.

All time-lapse images of cells were gathered in the database. The database consisted of six 48-h-long records. Since none of the cells remained in the field of view for the whole imaging due to migration, 150 min (30 time-lapse images with interval 5 min) were determined as an optimal length of the time-lapse record for one cell. Hundred and eighty cells were chosen for the monitoring. Based on their morphology, the cells were labeled by the expert biologist as either epithelial (95 cells) or mesenchymal (85 cells). The cells with uncertain class membership were not considered and were excluded from the database. The two types of classified cell morphologies are shown in Fig. 2.

Fig. 2

Examples of segmented QPIs of (a) epithelial and (b) mesenchymal phenotype gained by CCHM 2 and 30 h after the application of $\mathrm{TGF}\beta$, respectively. During EMT, the cell morphology changed from rounded to elongated, with the cell mass distributing relatively equally over the cell area, while the cell area increased significantly. QPIs are shown in grayscale in units of $pg/\mu {\mathrm{m}}^2$ recalculated from phase (in radians) according to Davies.³¹

2.4.

Classification

Classification, as a category of supervised machine learning, aims to build a model that makes predictions based on self-learning procedure on known labeled data. In the case of cell classification, the algorithm identifies patterns in the input images and trains a model based on labels assigned to the cells by an expert biologist. Such trained model is able to classify cells in new, previously unseen images. The essential requirement for successful classification is a sufficiently large database of labeled cell images, in which the classifier is trained. The overview of the classification process based on time-lapse QPI is shown in Fig. 3 and is described in more detail in the following paragraphs.

Fig. 3

Overview of the proposed feature-based classification process based on time-lapse QPI. First, image preprocessing is carried out. The cells in the image are segmented from the background and identified as ROIs. Cell features are extracted for every ROI. Feature values in several time instants create a time series. Dynamic features are extracted from the time-series while creating the feature vectors representing behavior of cells. The data are split into training and testing set. Both training and testing data are labeled by the expert biologist. After the feature selection, testing data form input for the classifier. The classifier is trained on the training data and prepared to perform the classification on testing data.

2.5.

Image Preprocessing and Feature Extraction

Before the classification itself, the cells were first segmented from the background in the time-lapse images by marker-controlled watershed segmentation approach.³³ The individual cells were tracked using the cell tracking algorithm scripted in MATLAB (MathWorks, Inc.). The algorithm performs cell tracking by linking every segmented cell in the given frame to the nearest cell in the next frame. Only cells remaining in the field of view throughout the specified time were considered for assessment. Further, highly overlapping cells where the segmentation was not clear were not included, nor were cells located on the border of the image. The included cells were identified as separate regions of interest (ROIs), where each ROI was represented by a set of cell features—a procedure referred to as feature extraction. Two types of features were extracted from each cell: morphometric and QPI features. Morphometric features mostly reflect the shape of the cell. These features involved footprint area, perimeter of the footprint area, convex area, perimeter of the convex area, solidity, roundness, indentation, eccentricity, extent, and centroid of the cell. QPI features are extracted from the phase values of the cell in QPI and therefore, contain quantitative information about the dry mass density distribution within the cell. QPI features were composed of the total phase of the cell, average phase, median, variance, standard deviation, skewness, and kurtosis of the phase values. Both types of features are described in more detail in our previous work.¹³ In the next step, all extracted features undergo normalization in order to scale the feature values to a fixed range from 0 to 1.

Each cell in a time instant therefore is represented by a feature vector composed of the cell features captured. Since every cell was recorded in time, each cell feature provides a univariate time series composed of the values of cell features over time. Accordingly, the consideration of all cell features, gives rise to a multivariate time series. There are two possible ways for the representation of time series. In the first, the values of time series itself represent the input for the classification, which will be referred to as a value-based approach. In the second, the feature-based approach, the time series is further represented by the newly defined time-lapse features, which subsequently form the time-lapse feature vector.

In order to explain the formation of the final time-lapse feature vector in the feature-based approach, the brief notation will be introduced. Let $\mathbf{X}=\{X_1,X_2,\dots ,X_Q\}$ represent a collection of $Q$ multivariate time series, where $Q$ is the number of cells in the experiment. Each multivariate time series $X_i$ is formed by $n$ observations ($n$ is the number of time points) and $d$-dimensional variable ($d$ is the number of cell features) as shown in Fig. 4(a). The multivariate time series $X_i$ can be written as

Fig. 4

Overview of the time-lapse feature extraction process. (a) Each cell is represented by the multivariate time series composed of univariate time series (formed by cell feature values obtained within the time period). (b) Time-lapse feature extraction from univariate time series composing a partial time-lapse feature set. Individual segments represent the group of time-lapse features obtained by the extraction technique. The length of the segments indicates the approximate number of extracted time-lapse features for the group. (c) Construction of the final time-lapse feature vector. The final vector representing a single cell is formed by the concatenation of partial time-lapse feature sets belonging to a cell. In addition, the motion and PCA features are added.

We will consider the $j$’th component of the $i$’th time series $X_{ij}=\{X_{ij1},\dots ,X_{ijn}\}$ to be a univariate time series. Therefore, the univariate time series will be composed of the values of one cell feature recorded in time. For each univariate time series $X_{ij}$, a feat time-lapse feature vector $M=(m_1,m_2,\dots ,m_L)$ is formed, where each $m$ is a time-lapse feature extracted from the time series and $L$ is the number of time-lapse features. In this way, each time series $X_{ij}$ is transformed into a partial time-lapse feature vector $M_{ij}$. Each multivariate time series is therefore transformed into $d$ time-lapse feature vectors. The vectors are then concatenated into a final time-lapse feature vector of $d\times L$ dimensions.

The idea of the transformation is extracting information which would otherwise not be obvious as well as reducing the number of features compared to the value-based approach and therefore lowering computational time. The latter is obvious when encountering longer time series. There are several possible methods used for dealing with the feature-based representation of the time series. The employed feature extraction techniques are briefly described in the following paragraph.

Statistical features carry global information about the time series. The following metrics were chosen in order to statistically represent the structure of the time series: mean value, median value, standard deviation, minimum value, maximum value, skewness, and kurtosis. Fourier transform features are formed by the most significant coefficients gained by FFT algorithm. Wavelet transformation features are formed by detail and approximation coefficients computed by discrete wavelet transform algorithm. The trend is represented by the coefficients obtained by the linear least squares fitting of the time series and characterizes a long-term change in the mean value of a cell feature. Approximate entropy quantifies the unpredictability of fluctuations in the time series. The presence of repetitive patterns of fluctuation in a time series renders it more predictable and leads to a relatively small approximate entropy. Symbolic aggregate approximation (SAX) features are gained by SAX method,³⁴ which is composed of two steps: piecewise aggregate approximation (PAA)³⁵ and the conversion of a PAA sequence into a string composed of letters, where the original time series is converted to a symbol string.

All of the above time-lapse features were extracted from each of the univariate time series and created a partial time-lapse feature vector as shown in Fig. 4(b). Subsequently, the partial time-lapse feature vectors obtained from each univariate time series were concatenated into a final time-lapse feature vector, whereas other extracted time-lapse features (principal components analysis and motion features) were added on the tail as shown in Fig. 4(c).

Principal components analysis (PCA) features were gained by applying PCA³⁶ on the whole multivariate time series, while mapping the multivariate data into a lower dimensional space. Motion features are composed of accumulated distance (overall distance travelled by the cell between the initial and the end point during the time interval), Euclidean distance (length of the straight line between the cell’s starting and end point reached during the time of monitoring), velocity (overall distance travelled by the cell over the elapsed time), and directionality (ratio of the Euclidian and accumulated distance). The position of the cell was determined by the cell centroids for all calculated motion features. Both PCA and motion features were added into the final time-lapse feature vector.

In the value-based approach, the extraction of time-lapse features is omitted, since the final time-lapse feature vector is composed of the raw data (values in each time point) contained in the multivariate time series. The final time-lapse feature vector is created by concatenating the univariate time series behind each other.

In both approaches, the final time-lapse feature vector represents a unique behavioral pattern of a cell. Before passing the vectors to the classification algorithms, the time-lapse feature values are scaled to a fixed range from 0 to 1. The example of a set of final time-lapse feature vectors gained by feature-based approach can be seen in Fig. 5, where the first 32 rows represent feature vectors extracted from epithelial cells and the other 35 rows from mesenchymal cells with the columns representing individual time-lapse feature values. The data are further split into training and testing set and are labeled by the expert biologist. Since the final time-lapse feature vectors are of substantial size, the next step is the selection of features with the highest potential to distinguish between the given classes, which would then form input for the machine learning classification algorithms.

Fig. 5

Example of the final time-lapse feature vectors concatenated into matrix. Elements of the matrix contain the (normalized) time-lapse feature values and are visualized using color: from blue (low values) to yellow (high values). First 32 rows represent time-lapse feature vectors extracted from epithelial cells and the other 35 rows from mesenchymal cells.

2.6.

Feature Selection

Since the time-lapse feature vectors are composed of a high number of features, feature selection is performed in order to reduce the dimensionality of the data, which leads to lower computation complexity and makes the training of the classification algorithms less time-consuming. Moreover, in this case, when the number of observations is limited in comparison to the large number of features, the limited observations may lead the learning algorithm to overfit to the noise. Reducing the number of features is therefore, in this case, an essential step before the classification.

We applied the filter approach for the feature selection.³⁷ First, the $t$-test was applied to each feature and the $p$-value for each feature was compared as a measure of the feature’s ability to discriminate between the two classes. To estimate the order of class separation by the features, the empirical cumulative distribution function (CDF) of the $p$-values was computed. There were $\sim 15\%$ of features, which have the $p$-values close to zero and 30% of features having the $p$-values smaller than 0.05. It can be concluded that there are roughly 200 features in the original time-lapse feature set, which have the potential to separate the two cell classes. In the value-based approach, there are $\sim 18\%$ of features, which have the $p$-values close to zero and 30% of features having the $p$-values smaller than 0.05. CDF of the $p$-values showed that there are roughly 150 features from the original time-lapse feature set having rather high discriminative power.

The features were subsequently ordered by their $p$-values. In order to define the appropriate number of features to be selected, the classification error (the number of misclassified observations divided by the number of observations) as a function of the number of features was plotted. To obtain the classification error, several classification algorithms were employed. The results of the classification error in feature-based approach is shown in Fig. 6. The classification error was computed for different numbers of features between 2 and 25. The final number of selected features was determined as the mean value of the results produced by employing different classification algorithms. In the feature-based approach, the filter feature selection method obtains the smallest classification error when 10 features are engaged. Only these 10 features with the highest discriminative power are kept in the reduced time-lapse feature vectors used for the classification. In the value-based approach, 12 features were determined as optimal.

Fig. 6

Classification error as a function of the number of features in feature-based approach.

2.7.

Classification Algorithms

After the features with the highest potential to distinguish between the epithelial and mesenchymal cell classes were selected, they create the input for the classification algorithms. Since the performance of the classification is highly dependent on the selection of the classification algorithm, we employed several supervised machine learning algorithms to correctly compare the performance of the classification based on single-time-point and time-lapse QPI. The following algorithms were tested in several possible variations with differently set parameters: decision trees (complex, medium, and simple tree, with defined maximum number of splits 100, 20, and 4, respectively), discriminant analysis (linear and quadratic discriminant), support vector machines (with linear, quadratic, cubic, and Gaussian kernel), $k$-nearest neighbor classifier (KNN), (fine KNN with $k=1$ and Euclidean distance, medium KNN with $k=10$ and Euclidean distance, cosine KNN with $k=10$ and cosine distance, cubic KNN with $k=10$ and cubic distance, weighted KNN with $k=10$, and weighted by the inverse square of the Euclidean distance) ensemble classifiers (bagged trees, boosted trees, subspace discriminant, and subspace KNN) and artificial neural network (feed-forward backpropagation neural network with one hidden layer containing 10 hidden neurons).

Several performance measures were calculated from the confusion matrix for each classification algorithm: accuracy, precision, recall, and $F$-score. Classification accuracy of a classifier is calculated as the ratio of the sum of the principal diagonal values to the sum of all values in the confusion matrix. It expresses the ratio of correctly classified examples by the classifier. Precision is the ratio of correctly classified positive examples to the total number of positive examples, while recall is the ratio of correctly classified positive examples to the all examples in actual class. $F$-score can be interpreted as a harmonic mean of precision and recall. Fivefold cross validation was used to evaluate the performance of the classification algorithms. The data were partitioned into five randomly chosen subsets of roughly equal size. One subset (testing set) was used for testing of the classifier, which had been trained on the remaining subsets (training set). This process was repeated five times, such that each subset was used for the validation. Since cross validation does not use all the data for training, it is a commonly used method to avoid overfitting. The overall performance of the classification algorithm was determined as the mean of performance measure values reached in the iterations. The whole classification procedure was performed in MATLAB.