Data-driven Maturity Level Evaluation for Cardiomyocytes derived from Human Pluripotent Stem Cells (Invited Paper)

1. Introduction

According to the world health organization report, cardiovascular diseases (CVDs) are the leading cause of death worldwide [1]. An estimated 17.9 million people died from CVDs in 2016, representing 31% of all global deaths [2]. In addition, statistical data have shown that one person dies every 36 seconds in the United States from cardiovascular diseases, and heart disease costs the united states about $219 billion each year from 2014 to 2015 [3]. Therefore, it is urgent to find a definitive treatment for heart disease.

The human pluripotent stem cell-derived cardiomyocytes (hPSC-CMs) represent a potential unlimited cell supply for the cell-based cardiac regeneration therapy, drug toxicity screening [4], and cardiovascular disease modeling [5]. However, current hPSC-CMs in vitro resemble a neonatal (immature) stage in cardiac markers, action potential, and morphology, and lag far behind the standard of mature CMs. To improve hPSC-CM maturity and function, researchers have made major efforts to promote the maturation of hPSC-CMs, including mechanical stress stimulation, electrical stimulation, biochemical cues, 3D cardiac tissue remodeling, modification of the substrate stiffness, and combinatorial co-culture to promote paracrine effects [6,7]. Recently, researchers have demonstrated that a focus on mitochondrial proteomics could help develop metabolically matured hPSC-CMs [8].

Since the maturation process of hPSC-CMs is influenced by many factors, developing efficient maturation strategies for hPSC-CMs requires a deep understanding of how each factor or combined factors direct the maturation process. However, the maturity level of hPSC-CMs is currently evaluated by various biological methods spanning different aspects, including cell morphology and structure, electrophysiology, calcium handling, and gene expression [9,10]. The lack of standardized metrics and protocols to quantify the overall maturity level of the cell in vitro impedes an understanding of the cell maturation process and the development of efficient cell maturation strategies.

To address the above challenges, we propose to evaluate the hPSC-CM maturity level with a data-driven approach. A data-driven model is a technology based on the analysis of the data about a specific system that can find the relationship between variables in the system without explicit knowledge of the physical behavior of the system [11,12,13]. With 48804 gene probe data at multiple in vitro culture time points [14], the data-driven model is a promising method for the quantitative evaluation of the hPSC-CM maturity level.

Based on the gene profile data, which span different stages of cardiac development, we have verified that hPSC-CMs can continually mature toward more adult-like hPSC-CMs up to 120 days. Then, by adopting in vitro culture duration as the proxy to estimate the maturity level of hPSC-CMs, these gene probes were ranked based on their importance to predict the culture duration. Five different predictive methods have been proposed and the results showed that the most accurate culture time prediction has an average error of 4.461 days, and the most important cardiac gene that correlates to the culture time is Gene CASQ2 (Calsequestrin 2), whose protein plays a role in the storage and transportation of calcium ions.

The major features of this work are: 1) The existing hPSC-CM maturity quantification schemes mainly focus on deciding whether the hPSC-CMs are mature or not. For the first time, we have proposed a data-driven pipeline that is capable of quantifying the maturity level of the cell in vitro with a finer granularity. 2) The selected dominating cardiac genes can guide the design of more efficient stimulation schemes to promote the hPSC-CM maturation. 3) The data-driven approach opens the door of the hPSC-CM study to multidisciplinary researchers without a strong biology background.

The remaining of the paper is organized as follows. Section 2 introduces existing biological studies on the maturity determination methods for hPSC-CMs. Section 3 verifies the feasibility of adopting in vitro culture duration as the proxy to estimate the maturity level of hPSC-CMs. Section 4 explains the cardiac gene selection and culture duration prediction models. Section 5 discusses the results and Section 6 presents the conclusion.

2. Related Works

As compared to the human adult cardiomyocytes, the current in vitro hPSC-CMs remain largely immature. By the principle of “if you cannot measure it, you cannot improve it” [15], the slowly moving hPSC-CM maturation research is challenged by the lack of a convenient and standard method to assess hPSC-CMs’ degree of maturity quantitatively. Currently, researchers have been using various biological methods to evaluate whether hPSC-CMs are becoming more mature towards adult CMs from various aspects: morphology and structure, calcium handling, contractile function and electrophysiological properties, and gene expression. The differences between immature-CMs and mature-CMs in these four aspects are summarized in Figure 1.

2.1. Morphology and Structure

The process from human fetal CMs to adult phenotype in vivo needs 6 to 10 years [16]. However, hPSC-CMs can be generated within 15 days of in vitro differentiation [17]. To nurture more mature hPSC-CMs, researchers increased culturing time up to 120 days and found key changes in cultured hPSC-CMs [9]. The main difference can be classified into four parts: 1) shape of the cell; 2) sarcomeres: a contractile unit of the muscle fiber; 3) sarcoplasmic reticulum (SR): also a structure found within muscles, and the main function of SR is to store calcium ions Ca2+; and 4) transverse tubules (T-tubules): a cell membrane that penetrates the center of the skeletal and cardiac muscle cell.

Adult cardiomyocytes are well-aligned, rod-like, multinucleated cells. They have highly organized sarcomeres and well-developed SR and T-tubules. In immature hPSC-CMs, those features tend to be lacking. The hPSC-CMs are small, rounded, mononucleated, with disorganized and shorter sarcomeres. Moreover, hPSC-CMs have poorly developed SR and no T-tubules [10].

Data Collection and Processing Methods: During the culturing period of 120 days, the status of hPSC-CMs can be divided into two stages: the early stage (from 20 days to 40 days) and the late stage (from 80 days to 120 days). During the culture process, the hPSC-CMs are observed under the microscope. Cell images are captured and saved into a computer. The parameters of sarcomeres length, cell perimeter, cell area, percentage of multinucleation, and circularity index (ratio between the cell width and length) are obtained with appropriate image analysis software [9].

2.2. Calcium Handling

In adult CMs, T-tubules and SR are well developed to regulate Ca2+ induced Ca release (CICR) and fast excitation-contraction coupling (ECC). The sharp and uniform increase of intracellular Ca2+ concentration in adult CMs is important for the synchronized contraction in multiple sarcomeres. By contrast, in hPSC-CMs, T-tubules are absent and SR is underdeveloped with low expression for the most part of sarcoendoplasmic reticulum calcium ATPase (SERCA) and other key proteins.

Data Collection and Processing Methods: Calcium imaging technology refers to a method for monitoring calcium ion concentration in tissues using calcium ion indicators named dyes. The cultured slip is observed under the microscope after placing the coverslip cultured with hPSC-CMs in a Petri dish containing the dye solution for about 30 minutes at 37 degrees Celsius. Images are captured and quantified with appropriate image analysis software.

2.3. Contractile Function and Electrophysiological Properties

The contractile function is a fundamental status indicator of CMs. 1) The hPSC-CMs and fetal CMs display a comparable force-generation capacity (0.22 ± 0.70 mN/mm² to 11.8 ± 4.5 mN/mm² in hPSC-CMs and ∼0.4 mN/mm² in fetal CMs), while adult CMs generate much larger forces (∼51 mN/mm²) [18]. 2) The hPSC-CMs also show immaturity in their electrophysiological properties, as compared with adult CMs, including reduced electrical excitability, higher resting membrane potential ($-20$ to $-60$ mV v.s.∼$-90$ mV), low capacitance (30-50 pF v.s.∼150 pF), smaller upstroke (15–50 V/s v.s. 180–400 V/s) and conduction velocity (2.1–20 cm/s v.s. 41–84 cm/s), and presence of spontaneous beating, which is found in early fetal CMs [19].

Data Collection and Processing Methods: Contractile and electrophysiological data are monitored and collected in real-time by devices such as HEKA EPC-10 patch-clamp amplifier. The changes in contractile and electrophysiological properties are analyzed by appropriate software, e.g., Patchmaster and Igor Pro.

2.4. Gene Expression

During the culturing of hPSC-CMs, several cardiac-specific genes express a more adult heart-like expression level over time, such as CASQ2, CRYAB, MYH6, MYH7, TNNI3, and ACTC1. 1) Gene CASQ2 results in a more mature calcium handling phenotype during culturing [20]. The CASQ2 gene provides instructions for marking a protein called calsequestrin 2 found in myocytes, where it is involved in storing and transporting calcium ions. 2) CRYAB is a part of the small heart protein family and functions as a molecular chaperone that primarily binds misfolded proteins to prevent protein aggregation, inhibit apoptosis, and contribute to in-tricellular architecture. 3) In an adult heart, the genes MYH6 and MYH7 are predominantly expressed in the ventricle and atrium, respectively. As the culture time of hPSC-CMs prolongs, the gene expression level of MYH7 and MYH6 tend to mature toward the adult cardiomyocytes level. Also, MYH7 and MYH6 provide essential instructions to form part of a large protein called type II myosin that generates the mechanical force to pump blood to the rest of the body. 4) The level of TNNI3 expression gradually increases with the culture time. The TNNI3 gene provides instruction for making a protein called cardiac troponin I, which helps coordinate the heart’s contraction. 5) The ACTC1 is the major protein of the thin filament in cardiac sarcomeres. Some other important cardiac-specific genes are listed in Table 1.

Data Collection and Processing Methods: The differentiation protocol is performed over millions of cells. At different culture dates, these cells are removed for RNA sample collection, and total RNA from hPSC-CMs can be isolated using reagent by following the manufacturer’s protocol [21,22]. To obtain the gene expression profiles, 48804 genes are probed with the microarray method, which can determine whether the RNA or DNA from a particular individual contains a gene mutation. The gene expression profile of hPSC-CMs is then compared with adult CMs by conducting statistical tests.

3. The hPSC-CM Maturity Evaluation Proxy

To propose a data-driven model capable of quantifying the maturation stage of hPSC-CMs with fine granularity, we first examine the gene expression data collected during culturing.

3.1. Cardiac-specific Gene Data Collection

The miRNA expressions of hPSC-CMs (Supp_table4 in [22]) are adopted as input data to analyze the maturation process. The gene expression datasets include: (a) The mRNA expression collected for a set of in vitro culture time points in $\mathcal{T}=\{0,3,7,10,14,20,28,35,45,60,90,120\}$ days. (b) Two independent runs of the adult mRNA expression, namely adult 1 sample and adult 2 sample, where the culture time is indicated as $t^{*}>>120$, approximately two years.

On day $t\in \mathcal{T}$, three million cells were sampled for RNA collection and $\mathcal{N}=\{1,\cdots ,203\}$ cardiac genes listed in the Genomic Institute of the Novartis Research Foundation (GNF) expression atlas were probed [23]. For the i-th probed cardiac gene expression, $i\in \mathcal{N}$, the fold change value $y_i^t$ is measured, which is used to describe the degree of change from the initial gene profile of the fetal sample to the gene profile on day t.

$$y_i^t=lo{g}_2\left({\displaystyle \frac{z_i^t}{z_i^{t_0}}}\right),\forall i\in \mathcal{N},\forall t\in \mathcal{T},$$

(1)

where for the i-th gene, $z_i^t$ indicates the corresponding gene probe profile of the hPSC-CMs on day t, and $z_i^{t_0}$ is the initial profile of the fetal hPSC-CMs. The $lo{g}_2$ transformed the ratio to “fold”, i.e., “times”, where $y_i^t=1$ means $z_i^t$ is doubling $z_i^{t_0}$.

The resulting fold change data set given below will be used to verify whether culture time can serve as the proxy of the maturity level of hPSC-CMs.

(1) A set of fold change values $Y_t=\left\{y_i^t\right|i\in \mathcal{N}\}$ for hPSC-CMs on day $t\in \mathcal{T}$, where $y_i^t$ is the fold change value of the i-th gene. (2) $Y^{*}=\left\{y_i^{*}\right|i\in \mathcal{N}\}$ for adult CMs on day $t^{*}$, where $y_i^{*}$ is the fold change value of the i-th gene.

3.2. Maturation Level v.s. Culture Duration

To compare the fold changes of the heart signature genes between hPSC-CMs on day t and adult CMs on day $t^{*}$, the Pearson correlation coefficient $r_{Y_t,Y^{*}}$ has been calculated during the hPSC-CMs differentiation for two random variables $Y_t$ and $Y^{*}$.

$$r_{Y_t,Y^{*}}={\displaystyle \frac{{cov}_{Y_t,Y^{*}}}{{\delta }_{Y_t}{\delta }_{Y^{*}}}}={\displaystyle \frac{\sum _{i\in \mathcal{N}}(y_i^t-\overline{Y_t})(y_i^{*}-\overline{Y^{*}})}{\sqrt{\sum _{i\in \mathcal{N}}{(y_i^t-\overline{Y_t})}^2\sum _{i\in \mathcal{N}}{(y_i^{*}-\overline{Y^{*}})}^2}}},$$

(2)

where $cov$ is the covariance and $\sigma$ is the standard deviation. $N=203$ is the sample size, which equals to the total number of sample cardiac genes indexed by i. For hPSC-CMs on day t during in vitro differentiation, $\overline{Y_t}=\left({\sum }_{i\in \mathcal{N}}{y}_i^t\right)/N$ is the sample mean of the 203 genes. For adult CMs, $\overline{Y^{*}}=\left({\sum }_{i\in \mathcal{N}}{y}_i^{*}\right)/N$ is the corresponding sample mean.

Pearson correlation coefficient is a measure of the similarity of two random variables between -1 and 1. If $Y_t$ and $Y^{*}$ are highly correlated, i.e., $|r_{Y_t,Y^{*}}|$ is close to 1, it is reasonable to build a model such that $Y_t$ can be predicted based on the value of $Y^{*}$. An example relationship between the response $Y_t$ and the predictor $Y^{*}$ can be described as a linear function:

$$\widehat{Y_t}={\alpha }_t+{\beta }_t{Y}^{*},\forall t\in \mathcal{T},$$

(3)

where $\widehat{Y_t}$ is the predicted response of $Y_t$. The coefficients ${\alpha }_t$ and ${\beta }_t$ are the intercept and slope, respectively. ${\beta }_t$ tells how much the dependent variable $\widehat{Y_t}$ is expected to increase when the independent variable $Y^{*}$ increases by one.

For the simple linear regression model in Eq. (3), ${\alpha }_t$ and ${\beta }_t$ are selected to minimize the difference between the predicted response $\widehat{Y_t}$ and the measured response $Y_t$. In particular, to quantify how much variance remains after fitting the linear model, the sum of squared error (SSE) is adopted to measure the performance of the predictive model.

$$SSE=\sum _{i\in \mathcal{N}}{(y_i^t-\widehat{Y_t})}^2,$$

(4)

where SSE is measured by the squared differences between the predicted and actual target values. To minimize SSE, the coefficients have the following values.

$$\left\{\begin{array}{c}{\beta }_t=r_{Y_t,Y^{*}}{\displaystyle \frac{{\delta }_{Y^{*}}}{{\delta }_{Y_t}}},\\ {\alpha }_t=\overline{Y_t}-\beta \overline{Y^{*}},\end{array}\right.$$

(5)

where the above coefficients are derived by setting the partial derivative of SSE to 0. As compared with slop ${\beta }_t$, $r_{Y_t,Y^{*}}$ can be treated as the standardized slope of the simple linear regression model.

Meanwhile, the coefficient of determination, $R_{Y_t,Y^{*}}^2$, measures that “$R_{Y_t,Y^{*}}^2$ ×100 % of the variation in $Y_t$ is reduced by taking into account predictor $Y^{*}$”. $R_{Y_t,Y^{*}}^2$ tells how much of the total variance can be explained by the linear model.

$$R_{Y_t,Y^{*}}^2=1-{\displaystyle \frac{SSE}{SST}}=1-{\displaystyle \frac{SSE}{\sum _{i\in \mathcal{N}}{(y_i^t-\overline{Y_t})}^2}},$$

(6)

where the sum of squared total (SST) quantifies the total variance of the target outcome as the sum of squared distance between individual data point $y_i^t$ and the mean of the response variable $\overline{Y_t}$. Note that for the linear regression model in Eq. (3), $R_{Y_t,Y^{*}}^2={\left(r_{Y_t,Y^{*}}\right)}^2$.

As illustrated in Table 2, $R_{Y_t,Y^{*}}^2$ increases with culture time t. In particular, the squared Pearson correlation coefficient $R_{Y_{28},Y^{*}}^2=0.55$ shows a moderate positive linear relationship between the gene expression profiles of day 28 hPSC-CMs and adult CMs. Meanwhile, the value of 93.0% in Figure 2 suggests a strong linear relationship between two adult samples, where only 7% of the variation in the adult 2 sample is left to explain after taking into account the adult 1 sample in a linear way. However, knowing the hPSC-CMs are continually becoming mature is not sufficient. The quantification of the mature stage will be discussed next.

4. Data-Driven Maturity Quantification

Aiming for an effective algorithm to quantify the maturity level of hPSC-CMs, the data-driven pipeline in Figure 3 is designed to select cardiac-specific genes and then map those genes to the in vitro culture time t, which is adopted as the proxy of the maturity stage.

4.1. Data-Driven Maturity Quantification Pipeline

Data Collection: We collected two publicly available data sets with gene-specific fold change values across 12 in vitro culture time points in $\mathcal{T}$. The first data set has recorded 203 cardiac genes for 1 group of cells (Supp_table4 in [22]), and the second data set has 48804 genes recorded for 3 groups of cells (GSM873339 - GSM873374 in [14]).

Step 1: Since both data sets are collected in the same culturing environment setting, to integrate them and increase the data volume, the gene “ID” (used by data set 2) is mapped to the “ILMN_Gene” (used by data set 1) according to the GPL6884 table in [24], which can translate the “unique identifier for the probe” to the “internal gene symbol”. Since multiple IDs can map to the same ILMN_Gene, we remove the redundant entries in data set 2.

After cleaning data set 2 and taking the intersection of gene symbols in both data sets, the resulting integrated data set has 48 data records: 4 groups cells across 12 in vitro culture time points. For each record, 189 genes listed in Table A1 in Appendix are recorded. For notation simplicity, the genes that are fed into the next steps (feature selection and predictive modeling) are denoted as $\mathbf{X}=\left\{x_i\right|i\in {\mathcal{N}}_c\}$, where $x_i$ represents the fold change value of the i-th cardiac gene in ${\mathcal{N}}_c=\{1,\cdots ,189\}$.

Steps 2-3: With 48 records of $\left\{\mathbf{X},T\right\}$, the train-test split in Table 3 separates the total data set into two parts: ${\mathcal{M}}_{train}=\{1,\cdots ,M_{train}\}$ records for training and the remaining ${\mathcal{M}}_{test}=\{1,\cdots ,M_{test}\}$ for testing. Based on the training data, a subset of $\mathbf{X}$ is selected such that regression analysis can be conducted to learn a function f and predict T.

$$\widehat{T}=f\left({\mathbf{X}}^{\prime }\right),{\mathbf{X}}^{\prime }\subseteq \mathbf{X},$$

(7)

where $\widehat{T}$ is the time estimated by the predictive model f.

The purpose of feature selection is to remove unnecessary, irrelevant, and redundant genes. The remaining important cardiac genes in ${\mathbf{X}}^{\prime }$ are useful to create an accurate predictive model f, which aims to minimize the root-mean-square error (RMSE).

$$f=argminRMSE\left(T\right)=argmin\sqrt{{\displaystyle \frac{1}{M_{train}}}SSE\left(T\right)},$$

(8)

where for the j-th training record, $T_j$ is the ground truth, and $\widehat{T_j}$ is the corresponding estimated value. $SSE\left(T\right)={\sum }_{j\in {\mathcal{M}}_{train}}{(T_j-{\widehat{T}}_j)}^2$ is the sum of squared errors of predicting T, which measures the deviation between the observed value and the true value.

Note that after splitting the training data and test data, the standardization of numerical features is performed over training data, such that each gene’s fold change value has zero mean and unit variance.

$$x_{i,j}^{Normalized}={\displaystyle \frac{x_{i,j}-{\overline{x}}_i}{{\sigma }_{x_i}}},\forall j\in {\mathcal{M}}_{train}\cup {\mathcal{M}}_{test},$$

(9)

where $x_{i,j}$ is the fold change value of the i-th gene in the j-th collected record. For each cardiac gene, based on all necessary statistics of the training data (mean ${\overline{x}}_i={\sum }_{j\in {\mathcal{M}}_{train}}{x}_{i,j}/M_{train}$ and standard deviation ${\sigma }_{x_i}$), the standardization is also performed on test data. For notation simplicity, we dropped the superscript $Normalized$ in the rest of the paper.

With all of the numerical values in the dataset being standardized, no gene can dominate the objective function and make the estimator unable to learn from other features correctly as expected.

Step 4: To verify the performance of the feature selection algorithm and predictive modeling, the standardized data records in the test data set are fed into the model, and the corresponding RMSE and $R^2$ score are adopted as the performance evaluation metrics.

$$R_{T,\widehat{T}}^2=1-{\displaystyle \frac{\sum _{j\in {\mathcal{M}}_{test}}{\left(T_j-{\widehat{T}}_j\right)}^2}{\sum _{j\in {\mathcal{M}}_{test}}{\left(T_j-{\overline{T}}_{test}\right)}^2}},$$

(10)

where ${\overline{T}}_{test}=\left({\sum }_{j\in {\mathcal{M}}_{test}}{T}_j\right)/M_{test}$ is the average ground truth culture time of the testing data. $R^2$ is the fraction of the total sum of squares that is explained by the regression, and the closer $R^2$ is to 1, the better the model.

Since Steps 1 and 4 are fixed for the proposed pipeline, to improve the performance of the hPSC-CM culture time prediction, the following feature selection algorithms (Step 2) and predictive modeling (Step 3) are investigated in our study.

• Method 1 (${\mathbf{M}}_1$): Filter method and linear regression
• Method 2 (${\mathbf{M}}_2$): Wrapper method
• Method 3 (${\mathbf{M}}_3$): Embedded method
• Method 4 (${\mathbf{M}}_4$): Non-linear feature selection and non-linear regression
• Method 5 (${\mathbf{M}}_5$): Non-linear feature selection and linear regression

4.2. ${\mathbf{M}}_1$ — Filter Method (Pearson Correlation) + Linear Regression

To select the dominant cardiac genes in ${\mathcal{N}}_c$ that determine the culture duration T, the univariate feature selection technique examines the linear strength of the relationship (such as Pearson correlation coefficient) between each input feature $x_i$ and the corresponding output variable T. The ranking of input features based on their strength of relationship concerning the output variable can gain a preliminary understanding of the collected data.

The detailed steps of the correlation-based filter method include: (1) the constant and quasi-constant cardiac genes in $\mathbf{X}$, which have variance less than the threshold 0.01, i.e., ${\left({\delta }_{x_i}\right)}^2<0.01$, are removed from ${\mathcal{N}}_c$. This step can delete the cardiac genes that do not change significantly when the culture time T increases. (2) Irrelevant cardiac genes with the absolute Pearson correlation coefficient less than the threshold 0.5, i.e., $|r_{x_i,T}|<0.5$ are removed as well, as they contain little information to predict the output T. (3) When the correlation coefficient for a pair of cardiac genes $(x_{i_1},x_{i_2})$ is higher than the threshold 0.5, i.e., $|r_{x_{i_1},x_{i_2}}|>0.5$, only one gene with a higher correlation with the output variable T is kept.

As illustrated by the heatmap in Figure 4, $k^{*}=2$ cardiac genes are kept by the filter method, and each selected input feature (${\mathbf{X}}^{\prime }=\left\{x_{19},x_{175}\right\}$) has a high correlation to the output T. The correlations among the selected features are less than the predefined threshold.

After the above feature selection steps, the multiple linear regression is applied to the clean data set and the coefficient of each feature in the linear function $f_1$ is obtained as follows:

$$\begin{array}{cc}\widehat{T}=\hfill & f_1\left(x_{19},x_{175}\right)=\sum _{i\in \left\{19,175\right\}}{\beta }_i{x}_i+\alpha \hfill \\ =\hfill & 30.42x_{19}-8.32x_{175}+37.67,\hfill \end{array}$$

(11)

where the coefficients of the above linear function are selected to minimize RMSE in Eq. (8). 37.67 is the intercept and other coefficients show the change in the output variable $\widehat{T}$ for one unit of change in the input cardiac gene while holding other input genes in the model constant. Since $|{\beta }_{19}|>|{\beta }_{175}|$, $x_{19}$ has more impact on the culture time prediction as compared to $x_{175}$.

4.3. ${\mathbf{M}}_2$ — Wrapper Method (Recursive Feature Elimination + Linear Regression)

Different from ${\mathbf{M}}_1$ with sequential feature selection and regression modeling, the wrapper method uses the performance of the regression model as evaluation criteria for the feature selection scheme. In particular, the wrapper method can generate different subsets of features, and each subset is used to build a model and train the learning algorithm. The subset yielding the best performance in terms of RMSE is selected as the final features.

During the subset generation process, the recursive feature elimination (RFE) method involves multi-rounds of elimination of the input variables. In each round, RFE specifies the number of input features that should be selected to build the regression model, and then performs recursive feature elimination. The detailed procedure is given below.

Iteration k: the number of features being selected is iterated from $k=1$ to $k=min\left\{\right|{\mathcal{N}}_c|,|{\mathcal{M}}_{train}\left|\right\}=40$, where $|•|$ means the cardinality of the cardiac genes in set ${\mathcal{N}}_c$ and the cardinality of the training set ${\mathcal{M}}_{train}$. Note that the number of features being selected cannot exceed the number of records in the training dataset.

In the k-th iteration, RFE performs the following two steps.

Step ${\mathbf{M}}_{\mathbf{2}.\mathbf{1}}$: Build a full linear regression model with all of the input features in the cardiac gene set ${\mathcal{N}}_c$.

$$\widehat{T}=f_2\left(\mathbf{X}\right)=\sum _{i\in {\mathcal{N}}_c}{\beta }_i{x}_i+\alpha ,$$

(12)

where the above model aims to minimize the RMSE in Eq. (8). Note that since we are trying to build a model with $|{\mathcal{N}}_c|=189$ coefficients based on $|{\mathcal{M}}_{train}|=40$ records in the training dataset, the full model is not unique because of the small dataset: $|{\mathcal{N}}_c|>|{\mathcal{M}}_{train}|$.

Step ${\mathbf{M}}_{\mathbf{2}.\mathbf{2}}$: Rank features based on the absolute value of the coefficient. The least important feature is pruned from the current set of cardiac features. Note that since the training dataset has been standardized to unit variance, the importance of a feature increases with increasing $|{\beta }_i|$. Otherwise, the feature importance ranking has to take account of the standard deviation of each cardiac gene as well: $|{\beta }_i|{\sigma }_{x_i}$ [25].

Ending Condition for the k-th iteration: For each k-th iteration, repeat the feature elimination Steps ${\mathbf{M}}_{\mathbf{2}.\mathbf{1}}$-${\mathbf{M}}_{\mathbf{2}.\mathbf{2}}$ over the pruned set until the number of genes left is k.

After all of the iterations are completed ($k=40$), the subset size $k^{*}$ that optimizes the performance criteria (smallest RMSE) is used to select the input variables, and the corresponding optimal subset is then used to train the final model. As illustrated in Figure 5, $k^{*}=39$ selected features are ranked based on the absolute linear regression coefficient $|{\beta }_i|$ in the final model.

4.4. ${\mathbf{M}}_3$ — Embedded Method (Lasso Regularization)

Embedded methods combine the qualities of the filter and wrapper methods. It is implemented by regression algorithms with built-in feature selection schemes. Regularization is the most commonly used embedded method that introduces additional constraints to bias the regression model toward less input cardiac features. As shown below, the least absolute shrinkage and selection operator (LASSO), i.e., $\mathbf{L}\mathbf{1}$ regularization for generalized linear models, has an additional weighted penalty against the complexity of the model.

$$f_3=argmin\left\{\sum _{j\in {\mathcal{M}}_{train}}{\left(T_j-{\widehat{T}}_j\right)}^2+\lambda \sum _{i\in {\mathcal{N}}_c}\left|{\beta }_i\right|\right\},$$

(13)

where $\widehat{T}=f_3\left(\mathbf{X}\right)={\sum }_{i\in {\mathcal{N}}_c}{\beta }_i{x}_i+\alpha$. The model complexity is the sum of the absolute coefficients for all the input features. Meanwhile, $0\le \lambda \le \infty$ is adjustable and the higher the value, the more the coefficient ${\beta }_i$ is forced to shrink. With this penalty term, if a feature $x_i$ is irrelevant, LASSO penalizes its coefficient and makes it 0, and this feature is then removed from the data set.

To determine the appropriate value for $\lambda$, cross-validation (CV) is adopted and the training dataset is divided into 5 non-overlapping folds. 100 values of $\lambda$ are iterated to train Lasso models over 5 folds. As illustrated in Figure 6, the optimal ${\lambda }^{*}$ yields the best average performance in terms of the smallest fold average RMSE: ${\lambda }^{*}=0.0226$. With the chosen optimal ${\lambda }^{*}$, $k^{*}=40$ features with non-zero coefficient ${\beta }_i$ in the corresponding predictive model are ranked in Figure 7.

4.5. ${\mathbf{M}}_4$ — Non-Linear Feature Selection and Regression (XGBoost Method)

When the regression model f in Eq. (7) represents a non-linear relationship between the input cardiac genes and the output culture time, the tree-based methods, e.g., decision trees, random forest, and extreme gradient boosting (XGBoost) [26], can be applied to perform feature selection with low complexity. They can model non-linear relations well and do not require much tuning.

As illustrated in Figure 8, the decision tree can create a regression model that predicts the culture time T by evaluating a sequence of true/false questions regarding the cardiac features. Although being an effective regression model, the decision tree is very sensitive and small changes to the training set can result in significantly different tree structures.

XGBoost addresses this issue by leveraging the wisdom of crowds wherein a large number of individual trees operating as a committee will outperform any of the individual constituent trees.

$$\widehat{T}=f_4\left(\mathbf{X}\right)=\sum _{b=1}^B{f}_4^b\left(\mathbf{X}\right),f_4^b\in \mathcal{F},$$

(14)

where B is the number of trees. $f_4^b$ is a function in the functional space $\mathcal{F}$, and $\mathcal{F}$ is the set of all possible trees, which can map the data record to the green leaf node. B is a tunable parameter that can be iterated to reduce the training RMSE.

The importance of a feature in the XGBoost can be measured as the number of times a feature is used to split the data across all trees. With this feature ranking mechanism, multiple XGBoost models are trained iteratively by feeding the top k features. As shown in Figure 9, the cross-validation performance shows that the best XGBoost model only needs the top $k^{*}=27$ features. When feeding the selected top 27 cardiac genes, various parameters in the XGBoost model are tuned based on the training data, and the resulting feature importance is ranked in Figure 10.

4.6. ${\mathbf{M}}_5$ — XGBoost Method + Linear Regression

${\mathbf{M}}_5$ adopts the input features selected by ${\mathbf{M}}_4$ (XGBoost) and predicts the output with the linear regression model. This method is introduced to measure the improvement that ${\mathbf{M}}_4$ achieves by capturing nonlinearity among input features and output T. In other words, this method signifies the impact and importance of implied nonlinearity among the input features and between pairs of inputs and output. The coefficients of the selected cardiac genes in the linear regression model are shown in Figure 11.

5. Results and Discussions

The five feature selection and regression model development methods (${\mathbf{M}}_1$ to ${\mathbf{M}}_5$) are applied to select the most important cardiac-specific genes that are related to the hPSC-CM culture duration. With the selected genes, predictive algorithms are adopted to quantify the maturity level of hPSC-CMs. The resulting performance over the test data is collected.

5.1. Culture Time Prediction

The comparison between the observed culture time T and the predicted culture time $\widehat{T}$ over the test dataset is shown in Figure 12. All of the five methods have acceptable test performance in terms of a small discrepancy between T and $\widehat{T}$, as verified by the low RMSE and the $R^2$ score closer to 1 in Table 4.

For the linear regression models, ${\mathbf{M}}_{\mathbf{2}}$ (RFE with 39 features) yields the best performance, which is followed by ${\mathbf{M}}_{\mathbf{3}}$ (Lasso with 40 features), ${\mathbf{M}}_{\mathbf{1}}$ (Pearson correlation with 2 features), and ${\mathbf{M}}_{\mathbf{5}}$ (with 27 features selected by XGBoost). For these 4 linear models, the predictive performance does not strictly grow or decrease with the number of selected cardiac genes. The reason is that on one hand, the predictive performance will improve when the number of important cardiac genes being selected increases. On the other hand, the performance will drop when less important cardiac genes are included in the predictive model.

Moreover, the non-linear model ${\mathbf{M}}_{\mathbf{4}}$ (XGBoost) outperforms all of the four linear models. This is another successful application of XGBoost because it can extract useful predictive information from 27 genes while ${\mathbf{M}}_{\mathbf{2}}$ (RFE) needs 39 cardiac genes. The performance degradation from ${\mathbf{M}}_{\mathbf{4}}$ to ${\mathbf{M}}_{\mathbf{5}}$ verifies that XGBoost can capture the non-linearity correlation among the selected cardiac genes and culture time.

5.2. Cardiac Gene Selection

As illustrated in Table 4, three methods (${\mathbf{M}}_{\mathbf{1}}$, ${\mathbf{M}}_{\mathbf{3}}$, and ${\mathbf{M}}_{\mathbf{4}}$) have selected $x_{19}$ as the top gene in ${\mathbf{X}}^{\prime }$ that correlates with the culture duration. The detailed importance level of $x_{19}$ can be found in Eq. (11), Figure 7, and Figure 10. Although ${\mathbf{M}}_{\mathbf{2}}$ (RFE in Figure 5) ranks $x_{19}$ as the 8-th important gene, the training dataset is very small (only 40 data records), there is randomness in the 189 coefficients of the full model returned by Eq. (12). Consequently, the ranking of the cardiac genes returned by ${\mathbf{M}}_{\mathbf{2}}$ has randomness as well. When $x_{19}$ is not ranked as the top cardiac gene, ${\mathbf{M}}_{\mathbf{5}}$ has the worst performance.

It is thus reasonable to conclude that $x_{19}$ is the most important gene in terms of predicting the culture time of hPSC-CMs. This is a very important observation because it can guide the maturation promotion scheme design for hPSC-CMs. $x_{19}$ (Gene CASQ2: Calsequestrin 2) triggers muscle contraction and plays a critical role in excitation-contraction coupling in the heart and in regulating the rate of heartbeats.

By stimulating hPSC-CMs such that the fold change value of CASQ2 increases faster, it is likely to promote the maturation of hPSC-CMs into a more adult-like phenotype. Meanwhile, since the cultured cells are destroyed during data acquisition, the existing gene expression-based hPSC-CM maturity quantification schemes are invasive and costly, similar to the morphology and structure-based schemes and calcium handling-based schemes described in Section 2. For non-invasive maturity quantification, the contractile function and electrophysiological property-based approaches have the most potential, because lab-on-chip devices [27] can be designed and fabricated to monitor the properties that correlate with $x_{19}$ in real-time as the culture duration prolongs.

6. Conclusion

We have proposed a data-driven maturity evaluation pipeline for hPSC-CMs. Based on the public gene expression data that spans various in vitro culture time points, the feasibility of adopting culture time as a proxy for assessing hPSC-CM maturity is verified. The biology domain knowledge and machine learning algorithms have been applied to select dominating cardiac genes that correlate with the culture time, and regression algorithms can further predict the culture time with an average error of fewer than 4.5 days. This work lays the foundation for standardizing the hPSC-CM maturation quantification process and provides insights into promoting the maturation of hPSC-CMs. Since cells may vary randomly in the rate of maturation and develop the same functionalities at different paces, the correlation between the post-differentiation culture time and the maturity level might not be perfect. Our future work will be exploring a more reliable proxy of the maturation level. Moreover, the gene expression data currently available is relatively small, a novel machine learning model will be proposed to uncover the latent maturity level with a relatively small amount of data.