Having identified the significant genes and proteins within the pathways, we now sought to obtain activity scores for each pathway. The purpose was two-fold. Firstly, to identify pathway activities that correlate with event-free patient survival; and secondly, to relate the genes and proteins expression of a pathway to its activity score. The RStudio package function PROGENy was used to obtain pathway scores from the RNA sequencing data of the TCGA legacy dataset. Albeit PROGENy's composition of 14-cancer relevant pathways, specifically for this study, only the CRC pathways were analysed in depth. The only available CRC pathways in PROGENy were PI3K, MAPK, TGF-Beta, WNT, and p53.

The PROGENy activity scores for each CRC signalling pathway are shown in Figure 6. Although the TGF-Beta and PI3K appear to have an inverse relationship, most pathways exhibit a cluster of high activity scores in a small subset of patients. However, specific activity patterns are difficult to discern.

**Figure 6.** PROGENy pathway activity scores from the RNA Sequencing expression dataset. Heatmap visualising the *z*-coefficients matrix for all features in each CRC signalling pathway for all 329 patients within the RNA sequencing TCGA legacy dataset cohort. The colourmap on the left indicates the activity score values ranging from 3–−3. The Euclidean distance metric was used to pass the pairwise distance between observations between both rows and columns. The complete linkage method was used to create the hierarchical cluster tree. Data were clustered along the columns of data, then along the rows of row-clustered data.

#### Kaplan–Meier Survival Analysis on PROGENy Activity Scores

To investigate the correlation between the PROGENy pathway activity scores and patient event-free survival, a Kaplan–Meier analysis and log-rank test were used. The activities of the WNT, PI3K, TGF-Beta, and MAPK pathways were all significantly correlated with event-free survival (Table 5). In particular, we found that activity of the WNT STN exhibited the strongest association with event-free patient survival (hazard ratio and *p*-value of 1.9731 and 0.0013, respectively).

**Table 5.** The association between PROGENy pathway activity and patient event-free survival. PROGENy scores were analysed using Kaplan–Meier analysis and the log-rank test to obtain hazard ratios and *p*-values.

The event-free survival Kaplan–Meier curves for the CRC STNs MAPK, TGF-Beta, PI3K, p53, and WNT, associated with patient survival through PROGENy activity scores, are shown in Figure 7. The STNs MAPK, TGF-Beta, and WNT followed a similar trend, where a high pathway activation was associated with shorter event-free survival. From the literature, this is expected considering that both pathways, MAPK and WNT, drive cell proliferation. Activation of the WNT pathway in CRC increases the levels of β-catenin within the cytosol, causing it to further travel into the nucleus and express WNT target genes, including genes that control the cell cycle [41]. For the WNT signalling pathway, the log2 hazard ratio was 0.98 with a 95% confidence interval between 0.38 and 1.58 and a *p*-value of 1.3 × 10−3 . The high group characterised by high WNT activity consisted of 108 patients and the low group of 221 patients. Similarly, the results for TGF-Beta are in line with its immunosuppressive function. Conversely, for the other STNs, p53 and PI3K, a high pathway activity was associated with longer event-free survival. The p53 result is in line with its known functions as a tumour-suppressor. The PI3K result is counterintuitive considering that the PI3K/AKT pathway is a classical survival signalling pathway but might be a confounding factor explained by "side-effect" activation due to crosstalk.

#### *2.3. Developing Linear Regression Machine Learning Models*

Having ascertained the significant genes, proteins, and pathways, we next sought to determine how well the expression of genes and proteins of a pathway can predict its activity scores. The results should give valuable insight into which candidate genes and proteins should be prioritised in follow-up studies focusing on constructing patientspecific mechanistic models with significant value. The WNT STN was the most significant pathway (Table 5, Figure 7); thus, we focused solely on it. Because there was a limited overlap between the mRNA and proteomics data in the TCGA legacy dataset (only five patients had data for both), we focused the analysis on the CPTAC (Clinical Proteomic Tumor Analysis Consortium) dataset. Both mRNA and proteomics data were available for the entire CPTAC cohort of 79 patients.

In total, three machine learning models were developed to predict the PROGENY pathway scores using different features as inputs. The features were the genes and proteins of the WNT pathway and their associated gene or protein expression values. To minimise the risk of overfitting, 10-fold cross-validation was used to build each model. The predictive power of each model was judged based on the root mean squared error (RMSE), the correlation coefficient (R), and the *p*-value. A *p*-value less than 0.05 was taken as significant for this field of study. The RMSE is an appropriate measure of predictive power for this study because model complexity was not an issue, and RMSE measures the differences between predicted and actual PROGENy values. There is no absolute threshold; however, lower values indicate a better overall fit of how accurately the model predicts the response. An acceptable range for the RSME values is difficult to define, but Model 1, consisting of the genes of the current mechanistic model, can act as a baseline benchmark. Notably, the predicted values of Model 2 show a much better correlation with the true values and a

lower RMSE value compared to Model 1. (Figure 8B); thus, demonstrating a much better fit. We focus our analysis on three models:

Model 1 was developed based on the mechanistic WNT signalling pathway. The selected features were the eleven genes of the mechanistic models in the literature [27–32], using the mRNA expression data from the CPTAC dataset (Table 6). The gene DKK1 is typically found in the mechanistic WNT pathway; however, it was not available in the RNA CPTAC dataset; therefore, its paralogue, DKK2, replaced it [42]. Models 2 and 3 were developed using systematic feature selection based on the entire set of 119 WNT pathway genes. LASSO regression was used to identify the features based on the mRNA expression data (Model 2) or protein expression data (Model 3). To complete the model, the two most significant genes (Model 2) or proteins (Model 3) were also included as features. LASSO regression on the mRNA data of the 119 WNT network genes selected seven genes as features. Thus, Model 2 consisted of nine features: seven identified by LASSO plus the two most significant genes. LASSO regression on the protein data of the 27 proteins of the WNT network for which data were available selected five proteins as features. Thus, Model 3 consisted of seven features: five identified by LASSO plus the two most significant proteins.

How well can the mechanism-based, mRNA-based, and protein-based models predict the pathway scores? Model 1, based on the genes of the mechanistic models, serves as a benchmark and achieved a reasonable predictive power with a Spearman correlation coefficient of 0.5241, a *p*-value of 7.12 × 10−7 , and an RMSE of 1.0167. Model 2, based on LASSO identified genes from the canonical WNT pathway, resulting in a greater predictive power, compared to Model 1, with a Spearman correlation coefficient of 0.8367, a *p*-value of 7.88 × 10−22, and an RMSE of 0.6416. Despite the smaller number of input features for Model 2, compared to Model 1, Model 2 evidentially has greater predictive power, shown by the higher *p*-value, correlation coefficient, and lower RMSE. Similarly to Model 2, Model 3 used LASSO to identify proteins from the canonical WNT pathway. In total, data for only 27 proteins were available due to missing data and imputation (see materials and methods). Model 3 is the least predictive model. The associated Spearman correlation coefficient, RMSE, and *p*-value were 0.4095, 1.0509, and 1.78 × 10−4 , respectively. Most features used in these machine learning models have nonzero regression coefficients (Figure 8D–F).

**Figure 7.** Kaplan–Meier curves of event-free survival associated with CRC PROGENy pathway activity scores. Kaplan–Meier curves associated with PROGENy pathway activity scores for the signalling pathways (**A**) MAPK, (**B**) TGFß, (**C**) PI3K, (**D**) p53, and (**E**) WNT. The patients were stratified into two groups according to the expression level of each feature. The groups are represented as high (orange line) and low (blue line), where *n* indicates the total number of patients in each group. For each group, the Kaplan–Meier curve of event-free survival was tested for statistically significant differences using a log-rank test. Log2HR = log base 2 of the hazard ratio. CI95% = 95% Confidence Interval. P-val = *p*-value. RNA Sequencing data is from the TCGA legacy dataset.

**Figure 8.** Linear regression machine learning models of the mechanistic and canonical WNT signalling pathways. Scatter plot visualising the correlation between the predicted score from (**A**) Model 1, using the 11 features of the mechanistic model (*y*-axis), and the PROGENy scores (*x*-axis). (**B**) Model 2, using the nine features of the canonical WNT model and significant genes (*y*-axis), and the PROGENy scores (*x*-axis). (**C**) Model 3, using the seven features of the canonical WNT model and significant proteins (*y*-axis), and the PROGENy scores (*x*-axis). Each blue dot represents 1 of 79 patients in the cohort. The black line indicates perfect predictions. "R" denotes the Pearson correlation coefficient, and "RMSE" is the root mean square error. The second row of figures represents the corresponding linear regression correlation coefficients for Models 1–3. Presented are the regression coefficients for the (**D**) Model 1, (**E**) Model 2, and (**F**) Model 3. The features used for each machine learning model are shown on the vertical axis, with the associated correlation coefficients shown on the horizontal axis. The blue dot represents the regression correlation coefficient, and the blue line (error bar) is the associated standard error. The dashed vertical yellow line represents the midline for the linear regression coefficients at 0.

**Table 6.** Overview of the significant genes and proteins and the linear regression machine learning models. (**A**) List of the top six significant features from RNA sequencing and proteomics from TCGA Legacy datasets. (**B**) Model 1: An overview of the features used for the Mechanistic WNT pathway model. Models 2 and 3: An overview of the features used for the canonical WNT pathway, from CPTAC RNA and proteomics datasets, respectively. The number of features used and associated RMSE for each model is listed in the first column. "RMSE" = root mean squared error. Significant features are highlighted in blue.

#### *2.4. Summary*

By incorporating the results from the event-free survival analysis, the PROGENy analysis, and the development of the linear regression machine learning models together, it is apparent that each finding paves the way to a valuable tool for predicting the heterogeneity of the WNT STN activity on an individual patient level. Concretely, each finding from each analysis builds sequentially upon another. The event-free survival analysis brings forth critical and novel features for the WNT pathway that have not yet been accounted for mechanistically. Interesting, from the PROGENy analysis, we found that the WNT signalling pathway was the most significant pathway associated with event-free patient survival. This finding solidifies the idea that there is a relationship between this pathway and patient survival that should be continued to be researched mechanistically. Finally, the motivation for using linear regression machine learning was particularly because it is an appropriate tool to predict WNT pathway activity from specified features (i.e., prognostic genes and proteins) on an induvial patient level. Precisely, the machine learning models additionally identified several features that were not only significantly associated with event-free patient survival, but also important for predicting the patient-specific WNT pathway activity scores. Thus, by incorporating such features into an updated mechanistic model of the canonical WNT STN activation, one would expect to better understand the patient-specific differences in the control of pathway activity, but also to relate pathway activity to clinical outcomes on the level of individual patients.

#### **3. Discussion**

Our experimental workflow employed numerous Kaplan–Meier survival scans, a PROGENy analysis, and developed three machine learning models based on mechanistic WNT pathway models and LASSO identified genes and proteins of the extant WNT network. The models can serve as a platform for improving current mechanistic models and WNT pathway activity-based patient stratification. Conversely, all models of the WNT signalling pathway to date have focussed on a few genes and proteins of the core pathway [27–32]. Building upon this work that describes the initial steps of WNT pathway activation, our work builds upon this by identifying new features and regression models that describe the heterogeneity of the WNT pathway activity scores across patients from two CRC cohorts [27]. We have also shown that the activity scores of the WNT pathway were the most significant in terms of the association to event-free patient survival compared to four other key pathways of CRC; the MAPK, PI3K, TGF-Beta, and p53 pathways.

#### *3.1. Event-Free Survival Analysis*

This study determined novel genes and proteins significantly associated with eventfree survival in CRC across the major five CRC signalling pathways, WNT, PI3K-Akt, TP53, MAPK, and TGF-Beta. In total, 389 genes and 75 proteins were found to be significant across all signalling pathways from a sum of 872 genes and 284 proteins. A total of 59 of these features were significant in the WNT signalling pathway, and interestingly, only seven of these significant features, DVL3, DKK1, SFPR2, WNT3A, WNT3, SFPR1, and GSK3B, are found within previous mechanistic WNT models [27–32]. Thus, in order to build new and improved mechanistic models of the WNT pathway with the potential to fulfil the promise of precision medicine, future modelling should focus on these significant genes and proteins.

For the WNT signalling pathway, the number of significant genes determined from RNA sequencing was 53. The top two significant genes found were DVL3 and VANGL2. The most significant gene was DVL3, with a corresponding hazard ratio of 2.8158 and a *p*-value of 2.58 × 10−4 . WNT ligands bind to FZD receptors and the LRP5/6 co-receptor, which leads to the recruitment of cytosolic DVL3, then relays this signal to downstream signalling events that result in the translocation of ß-catenin to the nucleus and target gene expression. Thus, DVL is recruited by the receptor Frizzled and prevents the fundamental destruction of cytosolic β-catenin [43]. Previous studies demonstrated that a high expres-

sion of DVL3 in CRC acts as an unfavourable prognostic marker [44]. Additionally, Zhao et al. suggest that DVL3 is a key regulator in CRC chemoresistance and targeting it may be a potential strategy for CRC therapy [45]. Our study solidifies this further, as shown in the Kaplan–Meier curve in Figure 3, in which a high expression of DVL3 was associated with a lower event-free patient survival. VANGL2, the second most significant gene in the WNT signalling pathway, has a hazard ratio and *p*-value of 2.7692 and 1.08 × 10−3 , respectively. Studies show that VANGL2 can be activated by WNT through Frizzled receptors [46]. Furthermore, within the WNT pathway, six proteins were significantly associated with event-free patient survival. The most significant protein was PRKCA. High PRKCA protein expression correlated with longer event-free survival (Figure 5). Interestingly, the PRKCA protein was also important for explaining the WNT pathway activity scores in Model 3, with a significant negative regression coefficient (Figure 8), which is in line with the known function of PRKCA, which is to inhibit WNT signalling via several mechanisms, including the phosphorylation of β-catenin, and enhance CRC cell death, concretely fitting to the survival data presented in Figure 5 [47]. Together, these results suggest a model in which PRKCA exerts its positive effects on patient survival via repressing the WNT activity.

Finally, the results for the WNT pathway analysis identified several significant genes and proteins that were associated with event-free survival but have not yet been accounted for in current WNT mechanistic models. Most, if not all, models to date are based on a few intracellular and extracellular components of the core WNT pathway [27–32]. Focusing on describing the general mechanistic details of WNT pathway activation, these models were not intended to model patient specific differences and describe the WNT signalling heterogeneity across patients. In contrast, our analysis and machine learning models aimed at identifying genes and proteins with significant event-free survival associations and describing patient-specific differences of WNT pathway activation. In total, the remaining 52 significant features found in this study are not found in these mechanistic models. Thus, we have identified many genes that were associated with event-free survival and that constitute prime candidates to improve current WNT mechanistic models. These improved models could then be used for patient stratification and predicting response to therapy [9,33].

#### *3.2. A PROGENy Analysis*

The second part of this study consisted of a PROGENy pathway activity score analysis. There were several reasons why PROGENy and not classical GSEA was performed. One major reason PROGENy does not require predefined groups for comparisons, unlike GSEA, which determines whether a defined set of genes shows statistically significant, concurring differences between two groups, for example, normal and tumour tissue [14,48]. PROGENy, on the other hand, predicts specific pathway activity scores for each individual patient; thus, allowing us to build models that can explain these scores [35,49].

An advantage of the PROGENy analysis is that it reveals patient-specific differences, showing that different pathways are active in different patients (Figure 6). Interestingly, the PROGENy analysis did not reveal a striking hyperactivation of the WNT pathway (Figure 6); most values were found to be between 0 and 0.5 despite it being the main driving event in tumorigeneses within most CRC patients. This finding is explained by the fact that the pathway analysis is a differential analysis. PROGENy results are *z*-scores quantifying standardised (mean normalised and scaled) relative differences, yet all patients might have high WNT activity in absolute terms.

The results correlating the PROGENy pathway activity scores to patient event-free survival show that the WNT, PI3K, TGF-Beta, and MAPK pathways, when activated, were significantly associated with event-free patient survival. It is very interesting to note that the WNT signalling pathway was the most significant pathway associated with event-free patient survival. These results are in line with the literature, where high WNT pathway activation was associated with shorter survival [27–32]. The activation of the WNT pathway in CRC increases the levels of β-catenin within the cytosol, causing it to translocate into

the nucleus and express WNT target genes that drive cell-proliferation [41]. How can the patient-specific differences of the WNT pathway activity be explained?

#### *3.3. Development of Linear Regression Models*

The final part of this study developed three linear regression machine learning models with the aim to predict the PROGENy WNT pathway activity scores. Current mechanistic models describe the pathway components and provide deterministic insights into how the WNT pathway is activated but neglect patient-specific differences [27–32]. Thus, combining machine learning that focuses on patient-specific differences with mechanistic models has great potential for building explainable machine learning models.

Overall, Model 2 was found as the most predictive model. In comparison, the RMSE of Model 2 was 0.3751 times smaller than Model 1. Five features were identified that were both prognostic of event-free patient survival and predictive of WNT pathway activity, including DVL3, PRKCA, VANGL2, GPC4, and ROCK2. Of these, DVL3 was the only feature found in the mechanistic model, Model 1. The finding also highlights a lack of features in current mechanistic models that can explain patient-specific differences. Concretely, solidifying the need to propose such a strategy to incorporate them. In particular, our work identified the four overlapping features above (PRKCA, VANGL2, GPC4, and ROCK2) that should be the focus of mechanistic modelling in the future. Machine learning, Model 2, accounting for such features, outperformed the current mechanistic model and can accurately predict the WNT pathway activity scores.

Interestingly, the protein PRKCA was found to be both significantly associated with event-free survival and a significant predictive feature of the WNT pathway activity score (Table 6). The results in Figure 8 show that most features used in the machine learning models have nonzero regression coefficients, in particular: DVL3, FZD5, RAC1, ROCK2, GSK3B, CTB2, CBT1, and PRKCA. These features provide valuable WNT network nodes to be included in future mechanistic models of the WNT network. The presented machine learning models incorporate patient specific significant features that describe pathway activity and can thus be termed "explainable". Despite this, these prognostic features have not yet been accounted for in current WNT models. The finding that a LASSObased machine learning model (Model 2) identified new features outperforms the current mechanistic model, and can accurately predict the WNT pathway activity scores, highlights the need to improve current mechanistic models by incorporating these features.

#### **4. Conclusions**

In conclusion, we developed machine learning models of the activity of the canonical WNT signalling pathway and incorporated significant genes and proteins that are associated with event-free patient survival. This model is a valuable tool for predicting the heterogeneity of WNT pathway activity on an individual patient level. In this study, we identified several features (DVL3, FZD5, RAC1, ROCK2, GSK3B, CTB2, CBT1, and PRKCA) that were not only significantly associated with event-free patient survival but also important for predicting the patient-specific WNT pathway activity scores. We expect that integrating these features into mechanistic models of WNT pathway activation in the future will not only help to better understand the patient-specific differences in the control of pathway activity but also to relate pathway activity to clinical outcomes on the level of individual patients.

#### **5. Materials and Methods**

#### *5.1. Data Acquisition from the Literature*

All data used throughout this study was open-source data. The first datasets were from the 2016 TCGA Colorectal Adenocarcinoma, GDAC Firehose Legacy study, previously known as TCGA provisional. The datasets used for this study included Genomic from RNA sequencing (20,532 genes), proteomic from mass spectrometry (5562 proteins), and the associated clinical dataset [39,40]. The TCGA legacy datasets, due to the minimum overlap

between patients, were only used to determine prognostic genes and proteins within such patients. Pre-processing across the datasets was performed before the survival analyses were complete. Across all datasets, patients who did not have the associated censoring status or disease-free survival months available were removed. The second set of data was from the 2019 prospective CPTAC-COAD colon adenocarcinoma studies, genomic from RNA sequencing (13,482 genes), proteomic from mass spectrometry (6422 proteins), and clinical datasets were all downloaded from LinkedOmics [39,40,50,51]. Minor preprocessing across the CPTAC datasets was performed before the following analyses were complete. The pre-processing included analysing the datasets, RNA, proteomics, and clinical to ensure that the same patients were overlapping across all and, if not, removing these patients. In total, 79 patients remained across all datasets. Additionally, for the CPTAC proteomics dataset, several NAN values were present. These values were imputed for the machine learning LASSO regression analysis of the WNT STN. This was achieved in the application Perseus using a low-shifted distribution of width 0.3 and a downshift of 1.8 across the total matrix [52]. Furthermore, proteins with data points of 10% or less were removed. In total, 27/39 proteins in the WNT STN remained after imputation. All datasets analysed and pre-processed for this study can be found in the supplementary data; consequently, Table 7 and Figure 9 represent key metrics and pre-processing steps for all datasets used.

**Table 7.** Important key metrics of datasets and corresponding analyses performed for each dataset. The two studies used, TCGA Legacy and CPTAC, are located in the first column [39,40,50,51]. The "datasets" column represents the type of omic dataset used from the corresponding study in column one, with key metrics identified. The "analyses" column represents the type of analysis performed on each specified dataset. "RSEM" = RNA-sequencing by expectation-maximisation. "UQ" = upper quantile normalisation. "TMT" = tandem mass tag.

**Figure 9.** Pre-processing flow chart detailing all pre-processing steps completed prior to each analysis. Pre-processing flowchart for, (**A**) the 2016 TCGA Colorectal Adenocarcinoma, GDAC Firehose Legacy Study and (**B**) the 2019 Prospective CPTAC-COAD Colon Adenocarcinoma Study.

The five critical CRC STNs used in this study stemmed from the literature [3,11,23]. The network databases WikiPathways, KEGG, and GSEA, were used to create a list of gene-sets for each pathway including WNT, PI3K-Akt, TP53, MAPK, and TGF-Beta [21, 22,37,38,53–62]. Each gene-set consisted of every gene listed for each homo-sapiens STN pathway on WikiPathways. These gene-sets were then used as input scripts to determine which genes and proteins are prognostic and correlate to patient event-free survival. All gene sets are found in the supplementary data.

#### *5.2. Kaplan–Meier Survival Analysis*

The Kaplan–Meier survival analysis was simulated three times, the first to find the prognostic genes from RNA sequencing, the second to find the prognostic proteins from proteomics, both using the TCGA legacy datasets, and finally to determine the associations between PROGENy pathway activity scores and event-free survival using the RNA CPTAC dataset. The optimum cut-off for stratifying the patient populations into low and high groups was identified by scanning the group sizes from 10–90 to 90–10 percent splits, where 10–90 means that 10% of the patients were in the low group and 90% of the patients were in the high group and calculating the *p*-value for the overall event-free survival difference between the groups using a log-rank test with Yates' correction. The cut-offs were based on the TCGA legacy datasets: RNA sequencing (329 patients) and proteomics (73 patients), and the CPTAC RNA dataset (79 patients). The inputs for the first two simulations were based on the gene or protein sets for each CRC STN identified above. The output was a Kaplan–Meier curve for the gene or protein expressed, indicating the number of patients in the high expression or low expression group with the corresponding statistical values apparent. The inputs for the final simulation were the PROGENy CRC pathway activity scores, and the output was 5 Kaplan–Meier curves for each CRC STN. The number of patients in the high activity or low activity group was identified with the corresponding statistical values. All statistical computations and Kaplan–Meier analyses were performed in MATLAB (version R2020b Update 5 (9.9.0. 1592791), The MathWorks, Inc., Natick, MA, USA) using the statistics toolbox and the log-rank (www.mathworks.com/matlabcentral/ fileexchange/22317-logrank (accessed on 7 April 2021)) and kmplot (www.mathworks. com/matlabcentral/fileexchange/22293-kmplot (accessed on 7 April 2021)) functions from

the MATLAB (version R2020b Update 5 (9.9.0. 1592791), The MathWorks, Inc., Natick, MA, USA) file exchange [63].

#### *5.3. PROGENy Analysis*

PROGENy is a machine learning based tool installed from Bioconductor as an RStudio package (RStudio Team (2020). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA URL http://www.rstudio.com/) [35,49]. A PROGENy analysis was performed on the TCGA RNA sequencing legacy dataset to find a common core of CRC pathway activities. Despite PROGENy's composition of 14-cancer relevant pathways, specifically for this study, only the CRC pathways were analysed in depth [35,49]. The available CRC pathways using PROGENy were PI3K, MAPK, TGF-Beta, WNT, and p53. The associated PROGENy CRC pathway activity scores for the entire TCGA legacy patient cohort is found in the supplementary data. PROGENy (version 1.12.0) was used in the study [35,49].

#### *5.4. Developing Linear Regression Machine Learning Models*

The three linear regression machine learning models developed in the study were performed in MATLAB using the regression learner app [64]. The three trained models were a variety of linear, stepwise linear, and robust linear regression models. Each model consisted of 79 observations, i.e., the entire CPTAC patient cohort. Cross-validation of 10-fold was used for each model. All models were developed solely from the CPTAC datasets from RNA and proteomics. The selected input features for the models were obtained through LASSO regression of both datasets. For both LASSO regression analyses, the lambda value with minimal mean squared error plus one standard deviation was applied. A threshold (*p*-value) of 0.05 was used to determine if a feature was significant. For the proteomics dataset, this resulted in 5/27 features (proteins) found to be significant, including CTBP1, CTBP2, PLCB4, PRKCA, and RAC1. Similarly, for the RNA sequencing dataset, 5/89 features (genes) were significant, DKK3, FZD5, NKD1, NOTUM, and WNT11. All associated datasets are found within the supplementary data.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/ijms22189970/s1.

**Author Contributions:** Conceptualisation, A.N. and D.F.; Data curation, A.N.; Formal analysis, A.N.; Funding acquisition, D.F.; Investigation, D.F.; Methodology, A.N.; Project administration, A.N.; Software, A.N.; Supervision, D.F.; Validation, D.F.; Visualisation, A.N.; Writing—original draft, A.N.; Writing—review and editing, A.N. and D.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the UCD Ad Astra Fellowship No 64998 and the European Union's Horizon 2020 research and innovation programme under grant agreement No 754923. The material presented and views expressed here are the responsibility of the author(s) only. The EU Commission takes no responsibility for any use made of the information set out.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The first dataset supporting this study were from the 2016 TCGA Colorectal Adenocarcinoma, GDAC Firehose Legacy study, previously known as TCGA provisional (cBioPortal for Cancer Genomics. Available Online: https://www.cbioportal.org/study/summary? id=coadread\_tcga (accessed on 1 March 2021)) [39,40]. The second and last dataset was from the 2019 prospective CPTAC-COAD colon adenocarcinoma studies. All data were downloaded from LinkedOmics. (Datasets for Colon adenocarcinoma: Prospective CPTAC-COAD. Available Online: http://linkedomics.org/cptac-colon/ (accessed on 8 April 2021)).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Review* **Treatment of Brain Metastases of Non-Small Cell Lung Carcinoma**

**Agnieszka Rybarczyk-Kasiuchnicz \*, Rodryg Ramlau and Katarzyna Stencel**

Department of Chemotherapy, Poznan University of Medical Sciences, Clinical Hospital of Lord Transfiguration, 61-848 Poznan, Poland<EMAIL_ADDRESS>(R.R<EMAIL_ADDRESS>(K.S.)

**\*** Correspondence<EMAIL_ADDRESS>

**Citation:** Rybarczyk-Kasiuchnicz, A.; Ramlau, R.; Stencel, K. Treatment of Brain Metastases of Non-Small Cell Lung Carcinoma. *Int. J. Mol. Sci.* **2021**, *22*, 593. https://doi.org/10.3390/ ijms22020593

Received: 3 December 2020 Accepted: 30 December 2020 Published: 8 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** Lung cancer is one of the most common malignant neoplasms. As a result of the disease's progression, patients may develop metastases to the central nervous system. The prognosis in this location is unfavorable; untreated metastatic lesions may lead to death within one to two months. Existing therapies—neurosurgery and radiation therapy—do not improve the prognosis for every patient. The discovery of Epidermal Growth Factor Receptor (EGFR)—activating mutations and Anaplastic Lymphoma Kinase (ALK) rearrangements in patients with non-small cell lung adenocarcinoma has allowed for the introduction of small-molecule tyrosine kinase inhibitors to the treatment of advanced-stage patients. The Epidermal Growth Factor Receptor (EGFR) is a transmembrane protein with tyrosine kinase-dependent activity. EGFR is present in membranes of all epithelial cells. In physiological conditions, it plays an important role in the process of cell growth and proliferation. Binding the ligand to the EGFR causes its dimerization and the activation of the intracellular signaling cascade. Signal transduction involves the activation of MAPK, AKT, and JNK, resulting in DNA synthesis and cell proliferation. In cancer cells, binding the ligand to the EGFR also leads to its dimerization and transduction of the signal to the cell interior. It has been demonstrated that activating mutations in the gene for EGFR-exon19 (deletion), L858R point mutation in exon 21, and mutation in exon 20 results in cancer cell proliferation. Continuous stimulation of the receptor inhibits apoptosis, stimulates invasion, intensifies angiogenesis, and facilitates the formation of distant metastases. As a consequence, the cancer progresses. These activating gene mutations for the EGFR are present in 10–20% of lung adenocarcinomas. Approximately 3–7% of patients with lung adenocarcinoma have the echinoderm microtubule-associated protein-like 4 (EML4)/ALK fusion gene. The fusion of the two genes EML4 and ALK results in a fusion gene that activates the intracellular signaling pathway, stimulates the proliferation of tumor cells, and inhibits apoptosis. A new group of drugs—small-molecule tyrosine kinase inhibitors—has been developed; the first generation includes gefitinib and erlotinib and the ALK inhibitor crizotinib. These drugs reversibly block the EGFR by stopping the signal transmission to the cell. The second-generation tyrosine kinase inhibitor (TKI) afatinib or ALK inhibitor alectinib block the receptor irreversibly. Clinical trials with TKI in patients with non-small cell lung adenocarcinoma with central nervous system (CNS) metastases have shown prolonged, progression-free survival, a high percentage of objective responses, and improved quality of life. Resistance to treatment with this group of drugs emerging during TKI therapy is the basis for the detection of resistance mutations. The T790M mutation, present in exon 20 of the EGFR gene, is detected in patients treated with first- and second-generation TKI and is overcome by Osimertinib, a third-generation TKI. The I117N resistance mutation in patients with the ALK mutation treated with alectinib is overcome by ceritinib. In this way, sequential therapy ensures the continuity of treatment. In patients with CNS metastases, attempts are made to simultaneously administer radiation therapy and tyrosine kinase inhibitors. Patients with lung adenocarcinoma with CNS metastases, without activating EGFR mutation and without ALK rearrangement, benefit from immunotherapy. This therapeutic option blocks the PD-1 receptor on the surface of T or B lymphocytes or PD-L1 located on cancer cells with an applicable antibody. Based on clinical trials, pembrolizumab and all antibodies are included in the treatment of non-small cell lung carcinoma with CNS metastases.

**Keywords:** brain metastases; treatment; non-small cell lung carcinoma; EGFR; ALK; immunotherapy

#### **1. Introduction**

Lung cancer is one of the most common malignant neoplasms and the main cause of death from malignant neoplasms in Poland. Each year, more than 22,000 new lung cancer cases are recorded [1].

Metastases of lung cancer to the brain occur in 18–61% of patients [2–4]. Improving the effectiveness of oncological treatment leads to a higher survival rate but also increases the population of patients at risk of this complication [5,6]. Metastases to the central nervous system (CNS) are different from other metastases due to the occurrence of neurological disorders and often require discontinuing systemic treatment in order to carry out palliative care [7].

#### **2. The Mechanism of Brain Metastasis Formation Is Similar to Other Organ Locations**

As a result of mutations in cancer cells, the degree of invasiveness increases [8,9] (Figure 1). The cells detach themselves from the primary tumor and penetrate the blood vessels, reaching other organs through the bloodstream. Being a very well-vascularized organ, the brain is often subjected to metastases [10]. Metastatic cells arrest at distinct sites and extravasate through vascular walls into the brain parenchyma. Cancer cells proliferate at the metastatic niche, form colonies in this parenchyma, and the subsequent proliferation of cells leads to clinically detectable metastatic lesions [8,9].

**Figure 1.** The main stages of cancer cell colonization of brain parenchyma (source: You H. et al., Front Immunol., 2019 [8]).

Untreated metastases to the CNS lead to a gradual deterioration of the patient's performance and to death within one to two months as a result of increased intracranial pressure. Avoiding or delaying these complications requires expertise in the radical and adjunctive treatment of brain metastases [11].

#### **3. Non-Small Cell Lung Carcinoma without Activating the Mutation of Epidermal Growth Factor Receptor (EGFR) or Anaplastic Lymphoma Kinase (ALK)**

A treatment method is chosen for lung carcinoma patients with CNS metastases based on their prognosis. It is determined by Radiation Therapy Oncology Group recursive partitioning analysis—RPA [12] (Table 1).

**Table 1.** Prognostic stages in patients with central nervous system metastases according to the Radiation Therapy Oncology Group (RPA). KPS: Karnofsky Performance Status.

The eligibility criteria are the patient's age, general performance, and the presence of metastases outside the CNS. The first group (class I) includes patients in good general condition, with a KI (Karnofsky Index) of 70% or more, less than 65 years old, without extracerebral metastases, and with good primary tumor control. The third class are patients in poor general condition, with a KI below 70%. The second class includes the remaining patients. In class III, radiotherapy of metastases to the CNS is not recommended due to a very bad prognosis. It is optimal to implement the best adjunctive treatment. The average period of survival is about two months [12].

In the case of a single metastasis to the brain up to four metastatic lesions to the CNS, either surgical removal or stereotactic radiotherapy (SRS) is recommended by the RPA in class I or II patients [7,13].

There is no evidence that the addition of whole-brain radiation therapy (WBRT) to stereotactic radiotherapy or surgery affects the overall survival of patients [14].

Data from a prospective study by Japanese researchers (JLGK0901) indicates that stereotaxis may be relevant in patients with more than three CNS metastases [15]. The observational trial lasted three years and included 1194 patients [15] with 1–10 newly diagnosed CNS metastases. The largest tumor volume was <10 mL and <3 cm in the longest dimension. The total cumulative volume did not exceed 15 mL, and the Karnofsky performance status was 70% or higher. All patients qualified in this way received stereotactic CNS radiotherapy. The overall survival (OS) of patients after stereotaxis was 13.9 months in 455 patients with a single metastasis, 10.8 months in 531 patients with two–four metastatic lesions, and 10.8 months in 208 patients with 5–10 metastatic lesions. An equal overall survival (OS) in patients with two–four metastatic lesions and 5–10 lesions indicates that stereotactic radiosurgery (SRS) is an important alternative to whole-brain radiotherapy in selected patients in good general condition.

Stereotaxis with or without whole-brain radiotherapy was analyzed in a third phase trial in which patients with one–four brain metastases of lung carcinoma were randomized [16]. Three hundred and sixty-four patients meeting the criteria for inclusion in the trial were analyzed. Fifty-one percent of patients received only stereotactic radiotherapy, while 49% received SRS followed by whole-brain radiation therapy (WBRT). It was shown that the age of patients significantly affects their survival. Stereotactic radiotherapy as a stand-alone treatment improves survival in patients aged 50 years or younger, with no difference in the age group over 50 years. Patients with a single metastasis experienced significantly longer survival than patients with two–four metastases. In the assessment of

cognitive disorders during treatment, patients under 50 years of age tolerated the therapy better in both arms of the trial. Patients with a single CNS metastasis, compared to patients with two–four lesions, had less severe cognitive impairment. Local disease control was better in the arm with SRS plus WBRT in both age groups.

In the third-phase QUARTZ trial [17], the role of whole-brain radiation therapy (WBRT) was assessed in patients with non-small cell lung carcinoma with inoperable CNS metastases. The patients were randomized into two groups. In one group, they received radiation therapy—WBRT 20 Gy in four fractions—and steroid therapy, and in the other group, the best adjunctive treatment without radiotherapy. The mean survival duration of patients in the radiotherapy arm was 49 days and 51 days in the optimal adjunctive treatment arm. In both groups, there were no differences in the quality of life and the use of steroids. The entire brain can be subjected to radiation therapy in a 20-Gy regimen in five fractions or 30 Gy in 10 fractions [18]. Alternative fractionation: 40 Gy in 20 fractions twice a day does not affect patients' survival times. Attempts have been made to use chemotherapy as a radiosensitizer without improving patients' survival times [19].

In patients with asymptomatic CNS metastases who have not yet received systemic treatment, the therapy sequence should be considered. In a study published in 2014 [20], patients with non-small cell lung carcinoma with asymptomatic CNS metastases (one to four lesions) received either stereotactic radiosurgery (SRS) followed by a two-drug cisplatin-based chemotherapy or chemotherapy alone. The average age of the patients was 58 years, with a mean total survival time of 14.6 months in the arm with stereotactic radiotherapy and chemotherapy; in the arm with chemotherapy alone—15.3 months. The average time to progression in the CNS was 9.4 months in the arm with SRS; in the arm with chemotherapy alone—6.6 months. The symptomatic progression of CNS lesions was more frequently observed in patients without stereotactic radiotherapy [19].

In the phase 3 trial [21], patients with CNS metastases of non-small cell lung carcinoma (NSCLC) received chemotherapy—cisplatin with vinorelbine (days 1, 8, 15, and 22) courses every 28 days, with a maximum of six courses. Whole-brain radiation therapy (WBRT)—30 Gy in 10 fractions—took place early in some patients—on days 1–12 of the first chemotherapy course and, in the second arm, after two chemotherapy courses (56 days). The objective response rate was 20% in the early radiotherapy arm and 21% in the delayed radiotherapy arm. The average survival duration in patients with delayed radiotherapy was 24 weeks and, in patients with early radiotherapy, 21 weeks. The results indicate that, during chemotherapy treatment, the implementation time of CNS palliative radiotherapy in patients with asymptomatic NSCLC brain metastases does not affect patients' survival duration [21].

In patients with symptomatic metastases of lung cancer to the CNS, the recommended dose of corticosteroids used long term in the prevention of cerebral edema is 4 mg of dexamethasone daily. Increasing the dose of the steroid to 16 mg daily does not improve the disease control but generates treatment toxicity [22].

#### **4. Non-Small Cell Lung Carcinoma with Present Epidermal Growth Factor Receptor (EGFR) Mutation and ALK Rearrangement**

The discovery of the EGFR mutation [23,24] and ALK rearrangement [24] and then the introduction of first-generation tyrosine kinase inhibitors (gefitinib—Figure 2A and erlotinib—Figure 2B) to the treatment of non-small cell lung carcinoma allowed, compared to platinum-based dual-drug chemotherapy, for longer progression-free survival (PFS), higher objective response rates (ORR), and better disease control rates (DCR) in comparison with two-drug platinum-derivative-based chemotherapy [23–25].

EGFR is one of the four members [25] of the HER family receptors, which comprise [26] EGFR/HER1/erbB1, HER2/erbB2, HER3/erbB3, and HER4/erbB4 [27]. EGFR signaling [25] is triggered by the binding of growth factors, such as Epidermal Growth Factor (EGF) [25,28], resulting in the dimerization of EGFR molecules [27]. Autophosphorylation and transphosphorylation of the receptors through their tyrosine kinase domains [25] leads to the recruitment of downstream effectors and the activation of proliferative and cell

survival signals [25]. In recent years, intensive research has been dedicated to the Epidermal Growth Factor Receptor (EGFR) [27] due to its significant role in the pathogenesis [27,29] of malignant tumors. In many types of cancers, intracellular pathways modulated by EGFR have been identified [25,28] as crucial factors influencing tumor survival and development [30]. On the other hand, EGFR has also been shown to be a promising molecular target [25–27] for potential therapeutic agents. Attempts to modify the signal transduction exerted by EGF have been made either by blocking [25] the activity of certain elements of the EGFR pathway or by direct inhibition of the EGF receptor itself [27]. Gefitinib and erlotinib target the ATP cleft [31] within the tyrosine kinase Epidermal Growth Factor Receptor (EGFR). Specific activating mutations within the tyrosine kinase domain of the EGFR molecularly correlate to the responses [23–25] to gefitinib or erlotinib (Figure 3A,B).

**Figure 2.** (**A**) Gefitinib—first-generation tyrosine kinase inhibitor; (**B**) Erlotinib—first-generation tyrosine kinase inhibitor.

In Figure 3A, the inhibitor (dark blue), representing gefitinib, occupies the ATP cleft. The locations of the two missense mutations are shown within the activating loop of the tyrosine kinase (light blue); the three in-frame deletions are all present within another loop (shown in red), which flanks the ATP cleft. Figure 3B shows a close-up view of the EGFR tyrosine kinase domain, with the critical amino acids implicated in binding the inhibitor. Specifically, 4-anilinoquinazoline compounds such as gefitinib inhibit catalysis by occupying the ATP-binding site, where they form hydrogen bonds with methionine 769 (M769) and cysteine 751 (C751) residues, whereas their anilino ring is close to the methionine 742 (M742),

lysine 721 (K721), and leucine 764 (L764) residues (all shown in green). In-frame deletions within the loop that is targeted by mutations (shown in red) are predicted to alter the positions of these amino acids relative to that of the inhibitor. Mutated residues (red) are shown within the activation loop of the tyrosine kinase (light blue).

**Figure 3.** (**A**,**B**) Clustering of mutations in the Epidermal Growth Factor Receptor (EGFR) gene (adapted from Lynch T.J. et al., The New England Journal of Medicine, 2004) [23].

Gain-of-function mutations [32] in the tyrosine kinase domain of the EGFR gene markedly increase the sensitivity to EGFR tyrosine kinase inhibitors (TKIs) [33]. It has been shown that 10–30% of all lung adenocarcinomas [34,35] contain an EGFR-activating mutation. EGFR mutations occur mostly in adenocarcinoma, younger women and girls [26], and never-smokers [23–25]. The increased prevalence [32] of EGFR mutations in the metastatic disease (early stage—14, 2% and metastatic—30, 3%) in the dataset may partially reflect referral bias [26] (Figure 4).

**Figure 4.** Distribution of the oncogenic driver mutations in non-small cell lung carcinoma (NSCLC) (adapted from Skoulidis F et al., 2019). (**A**) Early stage and (**B**) metastatic disease.

The most common oncogenic mutations are deletion in exon 19 (45–50% of all somatic EGFR mutations) and a point mutation (L858R) in exon 21 (35–45% of mutations) [25,36,37]. Ex20Ins mutations are the third-most common EGFR-activating mutations in NSCLC [38], which collectively account for approximately 4% to 10% of all EGFR mutations [35]. These mutations are predictive of the clinical activity of the EGFR TKIs [39], which yield a superior RR (response rate) [39,40] and PFS [40,41], as well as a better QoL (quality of life) [39–41] scores when compared with combination chemotherapy in the first-line setting [23]. The discovery of EGFR mutations and ALK rearrangements also contributed to the development of a new scale [37] of prognostic factors in patients with brain metastases of non-small cell lung carcinoma, taking into account the presence of EGFR mutations or ALK rearrangements. The Lung Cancer Molecular Markers Graded Prognostic Assessment

(Lung-molGPA) index facilitates making clinical decisions in this group of patients. In addition to the previous parameters [12], such as the patient's age, general performance, presence or absence of cancer outside the CNS, the number of brain metastases (one–four or >four), it also takes into account the gene status of the EGFR and ALK mutations. The higher the number of points obtained on this scale, the better the prognosis and longer survival of patients [37] (Table 2).

**Table 2.** Lung-molGPA (Lung Cancer Molecular Markers Graded Prognostic Assessment).

KPS—Karnofsky Performance Status, NA—not applicable, neg/unk—negative or unknown, pos—positive, BM—brain metastases, EGFR—Epidermal Growth Factor Receptor, and ALK—Anaplastic Lymphoma Kinase.

> First-generation tyrosine kinase inhibitors block [38,39] the EGFR receptor in a reversible manner. A better control of neoplastic disease during treatment with gefitinib or erlotinib [40,41], and the longer lives of patients, drew attention to the problem of metastatic lesions in the CNS. Lung cancer patients treated with first-generation TKI achieved a mean survival time of 33.1 months. After the diagnosis of disease progression in the CNS or in the meninges, the average survival time was 5.5 and 5.1 months. The incomplete penetration of drugs into the CNS through the blood–brain barrier causes a worse response to the first-generation TKI treatment in the brain and meninges [42]. Despite the low molecular weights of gefitinib and erlotinib, their penetration rates into the cerebrospinal fluid (1.13% and 2.77%, respectively) and the CNS concentration rates are low (3.7 ng/mL and 28.7 ng/mL, respectively) [43] (Table 3).

**Table 3.** Concentrations of the EGFR and ALK tyrosine kinase inhibitors in the cerebrospinal fluid (CSF).

Attempts have been made to increase the doses of gefitinib or erlotinib [44] or to introduce the pulsatile administration of drugs in patients with metastatic lesions in the CNS. The achieved therapeutic effects were still unsatisfactory due to the fact that higher doses of the first-generation TKI [45] increased the drug concentration index in the CNS, but the obtained effect was short-lived. A prolonged administration of high doses of erlotinib or gefitinib causes unacceptable toxicity and is not used [44–46].

Afatinib (Figure 5) is a second-generation tyrosine kinase inhibitor.

Acquired resistance occurs [47] in patients who initially benefit from EGFR-targeted therapies (first-generation tyrosine kinase inhibitors) [25,26].

**Figure 5.** Afatinib—second-generation tyrosine kinase inhibitor.

A clinical definition of acquired resistance to EGFR TKIs: acquired resistance in systemic progression (by Response Evaluation Criteria in Solid Tumors (RECIST) or World Health Organization (WHO) criteria) after a complete or partial response or >six months of stable disease after treatment with targeted therapy [48]. It irreversibly binds to the EGFR receptor and also has a higher affinity for the receptor compared to first-generation drugs. The studies LUX-Lung 3 [49] (cisplatin with pemetrexed) and LUX-Lung 6 [50] (cisplatin with gemcitabine) demonstrated the superiority of TKI over platinum-based two-drug chemotherapy with new-generation drugs. In the presented studies, patients receiving TKI compared to chemotherapy benefited from longer progression-free survival (PFS). They showed higher objective response rates (ORR) and a better disease control rate (DCR). The CNS penetration rate for afatinib is below 1%, and the CNS concentration is 0.46 ng/mL [51]. The LUX-Lung 3 and LUX-Lung 6 [49,50] studies were analyzed, taking into account asymptomatic brain metastases. The progression-free time in the LUX-Lung 3 trial [52] in patients with CNS metastases was 11.1 months in the afatinib arm and 5.4 months in the chemotherapy arm. In the LUX-Lung 6 trial, patients with CNS metastases treated with afatinib [53–55] achieved a progression-free time of 8.2 months, and in the chemotherapy arm, PFS was 4.7 months. Progression-free time in the afatinib arm compared to chemotherapy was equal in patients without brain metastases and in patients with CNS metastases [50,56,57]. The LUX-Lung 7 trial compared gefitinib with afatinib and included patients with central nervous system metastases. The mean follow-up was 27.3 months; progression-free survival for the afatinib arm was 11 months and 10.9 months for the gefitinib arm. The time to treatment failure for afatinib was 13.7 months and, for gefitinib, 11.5 months. Afatinib and gefitinib in the LUX-Lung 7 trial no difference in the overall survival (OS) [58]. Brueckl et al. (ESMO 2018 Congress, abstract 1449P) [59] presented an analysis of GIDEON, a prospective noninterventional study that was conducted in Germany to investigate the activity and tolerability of first-line afatinib in routine clinical care. Among 151 treated patients, the majority (72.8%) started treatment at an afatinib dose of ≥40 mg; 61.8% of them had dose reductions. In the group of patients starting at <40 mg, 46.2% had dose reductions, while dose increases were performed in 33.3%. The safety profile of afatinib was consistent with the known safety profile identified by the clinical trials. In spite of relatively high proportions of patients with brain metastases (approximately 30%) and uncommon *EGFR* mutations (approximately 13%), the results corroborated the clinical data for afatinib in the routine setting. The median PFS was 12.9 months, with a 12-month PFS rate of 54.6%. Seventy-three percent of patients

responded, and 90% obtained disease control. Both the ORRs and disease control rates (DCR) were independent of the type of *EGFR* mutation, the presence of baseline brain metastases, and the starting dose (Figure 6).

**Figure 6.** Overall response rates and disease control rates obtained with first-line afatinib in the noninterventional GIDEON study [59] (adapted from Brueckl et al., ESMO, 2018).

Osimertinib (Figure 7) is a third-generation tyrosine kinase inhibitor. In the AURA 3 clinical trial [60], it was compared to pemetrexed and cisplatin or carboplatin-based two-drug chemotherapy [61,62].

**Figure 7.** Osimertinib—third-generation tyrosine kinase inhibitor.

It was the second-line treatment for all patients, with the first- and second-generation EGFR TKI used in the first-line treatment. After the disease progressed, the T790M mutation determining the resistance [63] to drugs from the first- and second-generation TKI groups was determined, and patients were randomized to the Osimertinib arm or to the chemotherapy arm [64,65]. The trial also included patients with metastases to the central nervous system, without symptoms resulting from focal lesions in the CNS, who did not require treatment with steroids for at least four weeks before the start of the trial. The median treatment duration was 10.1 months for patients treated with Osimertinib (Osimertinib, *n* = 279) and 4.4 months for patients treated with chemotherapy (*n* = 140). The objective response rate (ORR) was 71% for the Osimertinib treatment and 31% for

chemotherapy-treated patients [63,64]. A subgroup analysis was performed; patients with measurable CNS lesions (one or more brain lesions) were included in the first group and patients with one or more lesions measurable and nonmeasurable in the CNS in the second group. In the first group of patients [64,65], the ORR was 70% in the Osimertinib arm and 31% in the chemotherapy arm [66]. In the second group of patients [64,65,67], the ORR was 40% in the Osimertinib arm and 17% in the chemotherapy arm. In both groups of patients, the mean response time in the CNS was 8.9 months in the treatment with Osimertinib and 5.7 months in the treatment with chemotherapy [67,68]. The mean PFS in the group of patients with measurable changes in the treatment with Osimertinib was 11.7 months, while, in chemotherapy, it was 5.6 months [69].

The EGFR T790M mutation [70] is the most common mechanism of TKI first- and second-generation resistance (detected in 50–60% of patients) [25]. It is unlikely that any erlotinib combination [70–72] will overcome this specific drug resistance mechanism.

Osimertinib, a third-generation small molecule tyrosine kinase inhibitor, is recommended in patients with the T790M resistance mutation [66,67]. It is also effective in patients with metastases to the central nervous system and the meninges [67,73,74].

In the phase 3 FLAURA [71–74] clinical trial, patients receiving Osimertinib achieved a PFS of 18.9 months and 10.2 months in the control arm (gefitinib or erlotinib). Patients with CNS metastases also benefited from treatment with Osimertinib [74,75].

Osimertinib is a third-generation tyrosine kinase inhibitor and has demonstrated high tolerability [73–75]. Some patients showed resistance to this drug, and the major mutation site is C797S on the EGFR gene (discovery of genome sequencing) [75,76]. In the future, when EGFR TKI drug resistance occurs [75,76], genetic testing could be used to select the treatment method corresponding to the resistance mechanism [74–76].

Progress in the field of molecular biology in recent years has enabled the identification of potential oncogenic pathways [77,78]. In 2007, Soda and his colleagues found an echinoderm microtubule-associated protein-like 4 (EML4) ALK fusion gene from non-smell-cell lung cancers [77]. These ALK fusion proteins can induce the constitutive activation of the ALK tyrosine kinase [77,78]. The oligomerization of domains such as the coiled-coil [77] domain of the fusion partner gives stimulation [79] ALK downstream pathways as a result [79]. The P13K-AKT-Mtor, RAS-MAPK-ERK, or JAK-STAT pathways are constitutively activated [77,79].

ALK mutations are rare and can be found in approximately 3–7% of patients with the diagnosis of NSCLC [77–79]. ALK mutations are more common in young, nonsmoking men with adenocarcinoma [78,79].

Crizotinib is an ATP-competitive, orally bioavailable ALK inhibitor [80] and was first applied for the treatment of EML4 ALK-positive NSCLC [81]. Crizotinib (Figure 8) was introduced based on the phase 3 Profile 1014 study [81] as a standard of treatment in patients with ALK-positive lung cancer.

This first-generation tyrosine kinase inhibitor has a concentration rate in the cerebrospinal fluid of 0.616 ng/mL and a penetration rate to the cerebrospinal fluid of 0.26% [51,80]. In the Profile 1014 trial [81], crizotinib achieved significantly longer PFS compared to chemotherapy (nine months vs. four months), and after 12 and 24 weeks of treatment, higher intracerebral DCR of 85% and 65% was observed in the arm with crizotinib and 45% and 25% in the arm with chemotherapy. The intracerebral control of the disease was also better in patients with metastases to the CNS—23% compared to chemotherapy. However, the isolated progression of the disease in the CNS was more frequent during treatment with crizotinib; extracerebral progression was more frequent during the treatment with chemotherapy (pemetrexed cisplatin). Patients with untreated CNS metastases or progression [82,83] of the disease were not randomized for the trial, and 20% of patients participating in the trial had CNS radiotherapy.

Nearly one-third of patients treated with crizotinib had CNS metastases in the first year of therapy. In some of these patients, it was the only location of neoplastic disease progression [83,84].

Alectinib (Figure 9) is a second-generation tyrosine kinase inhibitor used in patients with ALK-positive lung tumors and is also effective in the central nervous system [85,86].

**Figure 8.** Crizotinib—first-generation Anaplastic Lymphoma Kinase tyrosine kinase inhibitor (ALK TKI).

**Figure 9.** Alectinib—a second-generation ALK TKI.

Alectinib shows a concentration level of 2.69 nM [87,88] in the cerebrospinal fluid, and the penetration rate into the cerebrospinal fluid is 86% [51,87,88].

The ALEX phase 3 [89] clinical trial included previously untreated patients with advanced ALK-positive lung cancer. The patients received crizotinib or alectinib [89,90]. The primary endpoint was progression-free survival. The mean follow-up was 17.6 months in the crizotinib arm and 18.6 months for alectinib. The progression-free time was significantly longer in patients treated with alectinib [89,90] and was 12 months, while, in the crizotinib arm, it was 8.5 months [89,90]. In this trial, alectinib demonstrated therapeutic efficacy in patients with metastases to the central nervous system. Twelve percent of patients treated with alectinib (18 patients) and 45% (68 patients) treated with crizotinib demonstrated changes in the CNS. The one-year (12 months) CNS cumulative events (progression) level was 9.4% vs. 41.4% when comparing alectinib with crizotinib. Good intracerebral disease control coexisted with PFS—a mean average of 25.7 months for alectinib [90,91] and

10.4 months for crizotinib. Treatment toxicity was also lower in patients receiving alectinib and had a favorable safety profile [90,91].

Another ALK inhibitor was found to overcome crizotinib resistance and to better control disease in the CNS [92]. Ceritinib (Figure 10) is a second-generation ALK TKI and has a cerebrospinal fluid penetration rate of 15% [51]. Ceritinib is effective in patients with the I117N resistance mutation [92–94].

**Figure 10.** Ceritinib—a second-generation ALK TKI.

The phase 3 clinical trial ASCEND-4 [92] compared ceritinib with chemotherapy as the first-line treatment in patients with advanced lung cancer and ALK rearrangement. In patients with metastases to the CNS, the mean PFS was 10.7 months in the ceritinib arm and 6.7 months in the chemotherapy arm [92]. The overall intracranial response rate in patients with measurable CNS changes at the baseline was 72.7% for ceritinib and 27.3% for chemotherapy [92–94].

In the phase 1 ASCEND-1 trial [95], ceritinib achieved a total intracerebral ORR of 36% in ALK TKI previously treated patients and 63% in ALK TKI untreated patients (patients had baseline CNS measurable changes). In the ASCEND-2 trial the intracerebral ORR was almost 40%, and the intracerebral DCR was 85% [95–97].

Brigatinib (Figure 11) is a second-generation ALK TKI. Brigatinib was shown to be active against the G1202R mutation [98]. The G1202R mutation is resistant to first- and second-generation ALK inhibitors (crizotinib, alectinib, and ceritinib) [98,99].

It was noted that the G1202 mutation was discovered in about 50% of relapse patients following the use of brigatinib [98,99]. Brigatinib, another second-generation ALK inhibitor, demonstrated substantial activity in patients with crizotinib refractory *ALK*positive NSCLC; however, its activity in the alectinib refractory setting is unknown [98].

The phase 2 ALTA trial [100,101] evaluated the efficacy of brigatinib in patients with advanced ALK-positive non-small cell lung carcinoma previously treated with crizotinib [100,101]. Patients were randomized to two arms of the trial—in one arm, the dose was 90 mg, and, in the other arm, 180 mg for seven days, then 90 mg [101]. In patients with measurable changes in the CNS, the ORR at a higher dose of the drug was 67%, and, at a lower dose, 37%. The DCR exceeded 80% in both arms. In the case of nonmeasurable CNS metastases, the ORR and DCR were higher in patients receiving the higher dose of the drug (19% vs. 6% and 87% vs. 72%). Two-thirds of patients receiving the higher dose of the drug and having measurable lesions in the CNS had an intracerebral response lasting, on average, 16.6 months [100,101]. Brigatinib was compared with crizotinib in a phase 3 trial in patients with ALK-positive [101] lung cancer who had not been previously treated with TKI. Ninety patients had baseline CNS metastases, and 39 patients had measurable CNS

lesions with a diameter >10 mm. The intracerebral response to treatment in patients with measurable lesions was 78% in the brigatinib arm and 29% in the crizotinib arm. In the brigatinib group, 9% of patients had disease progression in the CNS, and, in the crizotinib group, 19% of patients [99,100]. Twelve-month PFS in the group of patients with metastatic lesions in the CNS at baseline was higher in the brigatinib arm—67% than in the crizotinib arm—21% [101].

**Figure 11.** Brigatinib—a second-generation ALK TKI.

Lorlatinib (Figure 12) is a third-generation ALK inhibitor with a penetration rate to the CNS of 20–30% [51]. Lorlatinib is indicated for the treatment of patients with ALK-positive metastatic non-small cell lung cancer [102] whose disease progressed on crizotinib [103] and at least one other ALK inhibitor. Lorlatinib has been shown to be active against almost all of the previously identified ALK TKI resistance mutations, including G1202R [103,104]. It is supposed to overcome the resistance of cancer cells to early-generation drugs [104].

**Figure 12.** Lorlatinib—a third-generation ALK TKI.

In a phase 1 trial [105], an intracerebral RR of 44% was achieved in the lorlatinib arm in patients with metastatic changes in the CNS for measurable and nonmeasurable lesions and 60% for measurable lesions. Approval was based on a phase 2 study [106,107] in which lorlatinib demonstrated a substantial overall and intracranial response [106,107].

#### **5. Simultaneous CNS Radiotherapy and TKI Therapy**

It was shown that lung cancer cells with the EGFR mutation are more radiosensitive [108] than those without.  Before the era of targeted lung cancer treatment, patients had either neurosurgical surgery, SRS, or whole-brain radiotherapy with the occurrence of metastases to the brain.

Two hundred and thirty patients with CNS metastases and EGFR mutations were identified and divided into two groups [110]. In one group, 116 patients received TKI (gefitinib, erlotinib, or icotinib), and, in the other group (51), TKI and simultaneous radiotherapy of the whole brain. An ORR of 52% was achieved in both groups; OS in the radiotherapy and TKI arm was 26.4 months and, for the treatment with only TKI, 21.6 months. Compared with TKIs alone, EGFR TKIs plus WBRT demonstrated intracranial progression-free survival (PFS) of 6.9 vs. 7.4 months (*p* = 0.232) and systemic PFS of 7.5 vs. 7.9 months (*p* = 0.546) [110].

In a meta-analysis of seven trials [109] involving 1086 patients with brain metastases, TKI therapy alone was compared with radiotherapy used before TKI therapy. It was shown that patients with non-small cell lung carcinoma and brain metastases who received radiotherapy prior to TKI therapy had longer intracerebral PFS and longer OS [109]. The analysis in the subgroups showed that the survival time of patients was longer in the group with one–three metastatic lesions [109], and shorter OS was obtained by patients with more metastatic lesions. The analysis confirmed that radiotherapy, by damaging the blood–brain barrier, increases the effectiveness of TKI therapy. Consequently, the combined therapy reduces relapse and improves the overall survival [108,109].

In patients with lung cancer and brain metastases, attempts were made to combine up-front CNS radiation and TKI therapy. Based [109] on the current available evidence, patients of non-small cell lung cancer with brain metastases and EGFR mutations have better OS and iPFS (intracerebral progression free survival) when they receive up-front radiotherapy and TKI than TKI alone [108–110].

The subgroup analysis [109] showed that never-smokers lived longer compared to tobacco smokers, and patients diagnosed with adenocarcinoma lived longer compared to other histopathological types. Patients with a better overall performance status (ECOG) lived longer than patients in worse general condition. In the group of patients with symptomatic brain metastases who received TKI and simultaneous whole-brain radiotherapy, significantly worse intracerebral PFS was observed compared to patients treated only with TKI [109–111].

Currently, it is not recommended to discontinue TKI therapy while radiating the whole brain. For stereotaxis (SRS), it is recommended to discontinue TKI three days before SRS and restart it three days after treatment; therefore, the interval is seven days [111].

#### **6. Immunotherapy of Lung Cancer with Brain Metastases**

Pembrolizumab is the drug of choice in the first-line treatment of patients with a PD-L1 expression in >50% of tumor cells in patients with non-small cell lung carcinoma without the EGFR or ALK mutation. In the registration trial of pembrolizumab—Keynote-024 [112], 9% of patients had CNS metastases, and, in Keynote-010 [113], 15% of patients had CNS metastases. Pembrolizumab is recommended for the first-line treatment of stage IV non-small cell lung carcinoma (including patients with stable metastatic lesions in the central nervous system) [114].

In the CheckMate 057 trial [115], nivolumab was administered to patients with nonsquamous lung cancer as a second-line treatment. Patients achieved an OS of 12.2 months; in the arm with chemotherapy, the OS was 9.4 months. Patients with stable metastatic lesions in the central nervous system were randomized for the trial [115].

In the EAP (expanded access program) 1588 trial [116], nivolumab was administered to patients with the IIIB/IV tostages of non-squamous lung cancer after progression on prior systemic therapy. Four hundred and nine patients had CNS metastases. They were neurologically stable and could receive a steroid therapy of up to 10 mg of prednisone daily. In the group of patients with metastatic lesions in the brain, the mean follow-up time was 6.1 months (0.1–21.9); the DCR was 39%, and the mean OS was 8.6 months; the CNS disease stabilized in 96 patients, 64 patients achieved a partial response, and 4 complete CNS responses during the nivolumab treatment [116]. Currently, nivolumab is recommended for the second-line treatment of stage IV non-small cell lung carcinoma [114].

Lung cancer metastases to the central nervous system pose a serious problem in oncological treatment. These lesions not only cause progression of the neoplastic disease but also manifest focal symptoms from the CNS, affecting the general condition of patients and worsening contact with them. Neurosurgery, stereotactic radiotherapy, and chemotherapy help to improve the clinical conditions of patients. Introducing new molecules into clinical practice gives a chance not only to improve the general condition of patients but also to prolong their lives.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Evaluating the Effect of Lenvatinib on Sorafenib-Resistant Hepatocellular Carcinoma Cells**

**Tingting Shi 1,\* , Hisakazu Iwama 2 , Koji Fujita 1 , Hideki Kobara 1 , Noriko Nishiyama 1 , Shintaro Fujihara 1 , Yasuhiro Goda 1 , Hirohito Yoneyama 1 , Asahiro Morishita 1 , Joji Tani 1 , Mari Yamada 1 , Mai Nakahara 1 , Kei Takuma 1 and Tsutomu Masaki 1,\***

**Abstract:** Hepatocellular carcinoma (HCC) is one of the major causes of cancer-related deaths worldwide. Sorafenib has been used as a first-line systemic treatment for over a decade. However, resistance to sorafenib limits patient response and presents a major hurdle during HCC treatment. Lenvatinib has been approved as a first-line systemic treatment for advanced HCC and is the first agent to achieve non-inferiority against sorafenib. Therefore, in the present study, we evaluated the inhibition efficacy of lenvatinib in sorafenib-resistant HCC cells. Only a few studies have been conducted on this topic. Two human HCC cell lines, Huh-7 and Hep-3B, were used to establish sorafenib resistance, and in vitro and in vivo studies were employed. Lenvatinib suppressed sorafenib-resistant HCC cell proliferation mainly by inducing G1 cell cycle arrest through ERK signaling. Hep-3B sorafenib-resistant cells showed partial cross-resistance to lenvatinib, possibly due to the contribution of poor autophagic responsiveness. Overall, the findings suggest that the underlying mechanism of lenvatinib in overcoming sorafenib resistance in HCC involves FGFR4-ERK signaling. Lenvatinib may be a suitable second-line therapy for unresectable HCC patients who have developed sorafenib resistance and express FGFR4.

**Keywords:** lenvatinib; sorafenib-resistant; hepatocellular carcinoma; FGFR4; autophagy; microRNA

#### **1. Introduction**

Hepatocellular carcinoma (HCC) is the most common primary liver cancer, with an increasing incidence over the past few decades in various populations; it is one of the major causes of cancer-related deaths worldwide [1,2]. Only a small fraction of patients is diagnosed with early stages of the disease, when curative strategies—liver resection (LR), ablative techniques, and orthotopic liver transplantation (OLT)—can be employed [3]. For patients with advanced HCC, systemic therapy is available, with prolonged overall survival (OS) rates.

Sorafenib, an oral multi-kinase inhibitor approved in 2008, targets the Raf–MEK–ERK pathway and several receptor tyrosine kinases, including vascular endothelial growth factor receptors (VEGFRs) 2 and 3, platelet-derived growth factor receptor (PDGFR), FMSrelated tyrosine kinase 3 (FLT3), Ret, and c-Kit [4]. The median survival time with sorafenib, used as a first-line systemic therapy for the past decade, was nearly 3 months longer than that with the placebo (10.7 months vs. 7.9 months; hazard ratio (HR) 0.69; *p* < 0.001) [4].

**Citation:** Shi, T.; Iwama, H.; Fujita, K.; Kobara, H.; Nishiyama, N.; Fujihara, S.; Goda, Y.; Yoneyama, H.; Morishita, A.; Tani, J.; et al. Evaluating the Effect of Lenvatinib on Sorafenib-Resistant Hepatocellular Carcinoma Cells. *Int. J. Mol. Sci.* **2021**, *22*, 13071. https://doi.org/10.3390/ ijms222313071

Academic Editor: Nam Deuk Kim

Received: 21 November 2021 Accepted: 29 November 2021 Published: 2 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

However, intrinsic and acquired resistance to sorafenib remains a huge challenge, with only approximately 30% of patients responsive to sorafenib [4,5]. The primary resistance of HCC cells to sorafenib is postulated to be associated with genetic heterogeneity—the overexpression of epidermal growth factor receptor (EGFR) or ligand may lead to sustained activation of EGFR downstream signaling and drug resistance to sorafenib [6,7]. Acquired resistance to sorafenib, which often develops within six months [5], may be associated with several factors, such as the phosphatidylinositol 3-kinase (PI3K)/Akt pathway, autophagy, epithelial-mesenchymal transition (EMT), tumor microenvironment, epigenetic regulation, microRNAs (miRNAs), and "vessel co-option"—the ability of tumors to hijack the existing vasculature in organs such as the lungs or liver, thereby limiting the need for angiogenesis [7,8].

Second- and later-line systemic treatments are needed for patients who fail to respond or are intolerant to sorafenib. The multitargeted multi-kinase inhibitors regorafenib [9] and cabozantinib [10] were approved in 2017 and 2018, respectively, as second-line drugs. Further, PD-1 immune checkpoint inhibitors, nivolumab and pembrolizumab, have been approved as second-line therapies for patients with advanced HCC [11].

Lenvatinib is a multi-kinase inhibitor that targets VEGFR 1–3, fibroblast growth factor (FGF) receptors 1–4, PDGFR-α, RET, and KIT [12]. A recent clinical trial showed that the median survival time of 13.6 months with lenvatinib (95% CI 12.1–14.9) was non-inferior to that with sorafenib (12.3 months, 10.4–13.9; HR 0.92, 95% CI 0.79–1.06) in untreated advanced HCC [13]. Thus, lenvatinib has been approved as a first-line systemic treatment for unresectable advanced HCC. Further, some evidence indicates that lenvatinib may be used as a second-line treatment for patients who are intolerant to sorafenib or following sorafenib failure [14–16].

To explore the possible advantage of lenvatinib and the underlying mechanisms in overcoming sorafenib-resistance in HCC we developed sorafenib-resistant cell lines and performed in vitro and in vivo experiments; further, we identify the expression of miRNAs associated with the effect of lenvatinib. Only a few studies have been conducted on this topic.

#### **2. Results**

#### *2.1. Lenvatinib Inhibits Sorafenib-Resistant HCC Cell Proliferation*

The IC50 of Huh-7SR cells (6.76 ± 0.48 µM) was 2.9-fold higher than that of Huh-7 cells (2.33 ± 0.22 µM), and the IC50 of Hep-3BSR cells (7.73 ± 0.27 µM) was 2.81-fold higher than that of Hep-3B cells (2.75 ± 0.44 µM) when exposed to sorafenib (Figure 1A). The lenvatinib IC50 values were not significantly different between Huh-7 (9.91 ± 0.95 µM) and Huh-7SR (10.56 ± 0.73 µM) cells; however, the lenvatinib IC50 of Hep-3BSR cells (27.49 ± 3.01 µM) was 9.85-fold that of Hep-3B cells (2.79 ± 0.19 µM) (Figure 1B). The Huh-7SR and Hep-3BSR cells were treated with 0.3, 1, 3, 10, and 30 µM lenvatinib for 96 h, and the anti-proliferative effect of lenvatinib was assessed using the cell viability assay. DMSO-treated cells were used as controls. Lenvatinib inhibited cell proliferation in the Huh-7SR and Hep-3BSR cells in a dose- and time-dependent manner. Further, 1 µM lenvatinib showed a significant effect in the Huh-7SR cells; however, 10 µM lenvatinib was required in the Hep-3BSR cells (Figure 1C). In the colony formation assay, during a long culture period (14 days), a small dose (1 µM) of lenvatinib showed an obvious antiproliferative effect in the Huh-7SR cells, but not in the Hep-3BSR cells (Figure 1D). This may be due to the partial cross-resistance of Hep-3BSR cells to lenvatinib; the anti-proliferative effect of lenvatinib in Hep-3BSR cells needs a higher dose than Huh-7SR cells. To further verify the anti-proliferative effect of lenvatinib in Huh-7SR and Hep-3BSR cells, a 3D tumor spheroid assay was performed. An in vitro dosage (10 µM) of sorafenib and lenvatinib was used in the 3D spheroid culture of both wild-type and resistant cells; this showed that lenvatinib (10 µM) has an advantageous anti-proliferative effect in all cells (Huh-7, Huh-7SR, Hep-3B, and Hep-3BSR) compared with sorafenib (10 µM). Moreover, sorafenib could partially inhibit Huh-7SR cell proliferation but not in Hep-3BSR cells (Figure 1E).

**Figure 1.** Lenvatinib inhibits the proliferation of sorafenib-resistant human hepatocellular carcinoma (HCC) cells. (**A**) Cytotoxic effect of sorafenib (IC50): Cells were treated for 48 h with increasing concentrations of sorafenib, and dimethyl sulfoxide (DMSO) was used as a negative control. IC50 of Huh-7 sorafenib-resistant (Huh-7SR) cells was 2.9-fold that of Huh-7 cells, and IC50 of Hep-3BSR cells was 2.81-fold that of Hep-3B cells. (**B**) Cytotoxic effect of lenvatinib (IC50): Cells were treated for 48 h with increasing concentrations of lenvatinib, and dimethyl sulfoxide (DMSO) was used as a negative control. Lenvatinib IC50 values of Huh-7 and Huh-7SR cells were not significantly different. Lenvatinib IC50 of Hep-3B sorafenib-resistant (Hep-3BSR) cells was 9.85-fold that of Hep-3B cells. (**C**) Anti-proliferative effect of lenvatinib measured using cell viability assay: Huh-7SR and Hep-3BSR cells were treated for 96 h with the indicated concentrations of lenvatinib or DMSO. The relative cell number was normalized with the control. (**D**) After 1 µM lenvatinib or DMSO treatment, the colonies formed were fixed and stained with crystal violet. (**E**) Three-dimensional tumor spheroid assay was performed to evaluate the effect of lenvatinib and sorafenib on the wild-type and sorafenib-resistant HCC cell proliferative ability. The representative images of spheroids are shown (scale bar: 5th day, 71 µm; 8th and 15th days, 107 µm). Spheroid diameters were measured on 15th day. Data are represented as means ± standard error of mean (SEM) or means ± standard deviation (SD) of at least three independent experiments. \* *p* < 0.05, \*\* *p* < 0.01, \*\*\* *p* < 0.001.

#### *2.2. Lenvatinib Induces Apoptosis and Cell Cycle G1 Phase Arrest, Decreases Invasion and Migration Ability, and Regulates the Expression of Angiogenesis-Related Proteins in Huh-7SR Cells*

To determine whether lenvatinib affected the apoptosis and cell cycle in Huh-7SR cells, these cells were treated with lenvatinib, and FCM was performed. The Huh-7SR cells were treated with 10 µM lenvatinib or DMSO for 24 h, and the DMSO-treated cells were used as controls. The proportion of early apoptotic cells (lower right quadrant) among lenvatinib-treated cells was significantly higher than that among untreated cells (Figure 2A). Additionally, the expression of caspase-cleaved cytokeratin 18, cCK-18, a marker of apoptosis and necrosis, was upregulated in the lenvatinib-treated cells (Figure 2B). Furthermore, in terms of cell cycle progression, the cell population in the G0/G1 phase significantly increased, whereas cells in the S and G2/M phases decreased, suggesting that lenvatinibtreated cells were arrested in the G1 phase (Figure 2C). Invasion and wound healing assays showed that lenvatinib (10 µM) decreased invasion and migration ability after incubation for 24 h (Figure 2D,E). Further, in the Huh-7SR cells, lenvatinib affected the regulation of the expression of angiogenesis-related proteins: dipeptidyl peptidase-4 (DPPIV/CD26) [17–19] and pigment epithelium-derived factor (PEDF) [20] were upregulated, and plasminogen activator inhibitor-1 (PAI-1) [21] was downregulated, after incubation with 10 µM lenvatinib for 24 h (Figure 2F).

#### *2.3. Lenvatinib Induces Cell Cycle G1 Phase Arrest, Decreases Migration Ability, and Regulates the Expression of Angiogenesis-Related Proteins in Hep-3BSR Cells*

The Hep-3BSR cells were treated with 10 µM lenvatinib or DMSO for 24 h, and the DMSO-treated cells were used as controls. No significant differences were found between the average proportion of early apoptotic cells in the lenvatinib-treated cells and the control (Figure 3A). Moreover, no changes were observed in the expression of cCK-18 between the lenvatinib-treated and untreated cells (Figure 3B). In cell cycle progression, the cell population in the G0/G1 phase significantly increased, whereas the cells in the S and G2/M phases decreased, suggesting that the lenvatinib-treated cells were arrested in the G1 phase (Figure 3C). The invasion assay suggested no difference between the control and treatment groups (Figure 3D). The wound healing assay indicated that 10 µM lenvatinib decreased migration ability after incubation for 24 h (Figure 3E). Additionally, the expression of angiogenesis-related proteins such as endostatin [22], thrombospondin-1 (TSP-1) [23,24], interleukin-8 (IL-8) [25], along with PAI-1, were downregulated after incubation with 10 µM lenvatinib for 24 h (Figure 3F). The expression of PAI-1 was downregulated in both the Huh-7SR and Hep-3BSR cells after lenvatinib treatment.

#### *2.4. Lenvatinib Effects against Sorafenib-Resistance in HCC Cells May Be through the FGFR4-ERK Signaling Pathway*

The expression levels of EGFR, p-Akt, and p-ERK were found to be upregulated in the sorafenib-resistant cells compared to the wild-type cells (Figure 4A). These results were similar to those reported in previous studies [6–8], showing that the overexpression of EGFR activated the Akt and ERK signaling pathway to increase HCC cell survival and proliferation, which, in turn, induced sorafenib resistance. The kinase inhibition profile for lenvatinib and sorafenib indicated that lenvatinib has an advantage in FGFR4 inhibition compared with sorafenib (Table S1). Previous studies have shown that Huh-7 and Hep-3B have high FGFR4 expression [26,27]. In the Huh-7SR and Hep-3BSR cells, 10 µM sorafenib or lenvatinib decreased the expression of p-Akt, Akt, FGFR4, p-ERK, and cyclinD1, whereas the levels of p-mTOR and mTOR remained unchanged in the Hep-3BSR cells. Meanwhile, lenvatinib increased the expression of EGFR (Figure 4B). The effects of lenvatinib in the two resistant cells were mainly through FGFR4-ERK signaling; however, the difference between 10 µM sorafenib and lenvatinib treatment was not significant. After decreasing the dosage, lenvatinib showed a superior inhibition of ERK compared with the same dosage of sorafenib (Figure 4C). In addition, the low expression of FGFR4 was associated with a longer overall survival probability based on the TCGA database analysis (Figure 4D).

**Figure 2.** In Huh-7SR cells, lenvatinib induces apoptosis and cell cycle G1 phase arrest, reduces invasion and migration, and modulates the expression of proteins associated with angiogenesis. (**A**) Huh-7SR cells were treated with 10 µM lenvatinib or dimethyl sulfoxide (DMSO) for 24 h, and the level of apoptosis was measured by staining with annexin V and propidium iodide (PI) using flow cytometry. Lenvatinib treatment increased the proportion of early apoptotic cells in the Huh-7SR population. Lower right square represents early apoptosis. (**B**) The expression of caspase-cleaved cytokeratin 18 (cCK-18) was determined using Enzyme Linked Immunosorbent Assay (ELISA) after 24 h of treatment with 10 µM lenvatinib. (**C**) Huh-7SR cells treated with 10 µM lenvatinib or DMSO were analyzed using flow cytometry to determine the number of cells in each phase of the cell cycle (left panel). Representative cell cycle histograms are presented (right panel). Lenvatinib blocked the cell cycle at the G1 phase. (**D**) Invasion ability of lenvatinib-treated Huh-7SR cells was decreased. (**E**) Wound-healing assay comparing the motility of Huh-7SR cells treated with lenvatinib or DMSO. The wound-healing area was analyzed using the ImageJ software. (**F**) Representative expression of angiogenesis-related proteins: plasminogen activator inhibitor-1 (PAI-1), dipeptidyl peptidase-4 (DPPIV/CD26), and pigment epithelium-derived factor (PEDF), in Huh-7SR cells incubated with lenvatinib or DMSO for 24 h. Data are presented from three independent experiments. \* *p* < 0.05, \*\* *p* < 0.01.

**Figure 3.** In Hep-3BSR cells, lenvatinib triggers cell cycle G1 phase arrest, reduces migration, and modulates the expression of proteins associated with angiogenesis. (**A**) Hep-3BSR cells were treated with 10 µM lenvatinib or dimethyl sulfoxide (DMSO) for 24 h, and the level of apoptosis was measured by staining with annexin V and propidium iodide (PI) using flow cytometry. Lenvatinib treatment did not significantly change the average proportion of early apoptotic cells in the Hep-3BSR population. Lower right square represents early apoptosis. (**B**) The expression of caspase-cleaved cytokeratin 18 (cCK-18) was determined using Enzyme Linked Immunosorbent Assay (ELISA). After 24 h of treatment with 10 µM lenvatinib, cCK-18 levels remained changed. (**C**) Hep-3BSR cells treated with 10 µM lenvatinib or DMSO were analyzed using flow cytometry to determine the number of cells in each phase of the cell cycle (left panel). Representative cell cycle histograms are presented (right panel). Lenvatinib blocked the cell cycle at the G0/G1 phase. (**D**) Invasion ability of lenvatinib-treated Hep-3BSR cells was not significantly changed. (**E**) Wound-healing assay comparing the motility of Huh-7SR cells treated with lenvatinib or DMSO. The wound-healing area was analyzed using ImageJ software. (**F**) Representative expression of angiogenesis-related proteins, endostatin, thrombospondin-1 (TSP-1), interleukin-8 (IL-8), and PAI-1, in Hep-3BSR cells incubated with lenvatinib or DMSO for 24 h. Data are presented from three independent experiments. \* *p* < 0.05, \*\* *p* < 0.01.

**Figure 4.** The underlying mechanisms of lenvatinib in sorafenib-resistant HCC cells may be through the FGFR4-ERK signaling pathway. (**A**) Proteins expressed in wild-type and sorafenib-resistant HCC cells. (**B**) Huh-7SR and Hep-3BSR cells were treated with 10 µM sorafenib or lenvatinib for 6 h. Protein expression was analyzed using Western blot. (**C**) Huh-7SR and Hep-3BSR cells were treated with 1 and 5 µM sorafenib or lenvatinib for 6 h, and the levels of FGFR4, p-ERK, and ERK were analyzed using Western blot. C.: control; Sor: sorafenib; Len: lenvatinib. (**D**) A lower FGFR4 level is associated with longer overall survival of HCC patients based on TCGA database analysis (http://www.oncolnc.org, accessed on 5 October 2021).

#### *2.5. Different Autophagic Responsiveness between Huh-7SR and Hep-3BSR Cells*

Huh-7 and Hep-3B wild-type cells are both sensitive to lenvatinib (Figure S1). The IC50 data indicated that the Hep-3BSR cells showed partial cross-resistance to lenvatinib compared with the Huh-7SR cells. Moreover, lenvatinib inhibited Hep-3BSR cell proliferation mainly through cell cycle arrest, whereas in Huh-7SR, cell proliferation was inhibited through both apoptosis and cell cycle arrest. In addition, although mTOR is a key molecule in the regulation of autophagy, its activation remained changed after sorafenib or lenvatinib treatment in the Hep-3BSR cells. Autophagy plays a dual role in cancer development and is also associated with multidrug resistance in cancer cells [7,28]; thus, it was important to investigate if the autophagic responsiveness to sorafenib and lenvatinib is different between the two cell lines. In the Huh-7SR cells, 10 µM sorafenib or lenvatinib significantly increased the expression of microtubule-associated protein light chain 3-II (LC3-II); however, in the Hep-3BSR cells, the expression barely changed. Furthermore, in the Huh-7SR cells, while LC3-II levels increased, the expression of p62 decreased, indicating that 10 µM sorafenib or lenvatinib altered the autophagy in the Huh-7SR cells. Meanwhile, the levels of caspase-7 and caspase-3 were decreased, and the levels of PARP and cleaved-PARP were increased in the Huh-7SR cells, but not in the Hep-3BSR cells (Figure 5A). The Huh-7SR cells showed higher autophagic responsiveness to sorafenib and lenvatinib than the Hep-3BSR cells did. Lenvatinib induced high autophagic responsiveness that may induce autophagic cell death, and this may be one of the reasons that the Hep-3BSR cells showed partial cross-resistance to lenvatinib. Moreover, a lower p62 level is associated with the longer overall survival of HCC patients based on TCGA database analysis (Figure 5B).

**Figure 5.** Autophagic response of sorafenib-resistant cells to sorafenib or lenvatinib. (**A**) Huh-7SR and Hep-3BSR cells were treated with 10 µM sorafenib or lenvatinib for 6 h. LC3, p62, caspase-7, caspase-3, and PARP expression were analyzed using Western blot. C.: control; Sor: sorafenib; Len: lenvatinib. (**B**) A lower p62 level is associated with longer overall survival of HCC patients based on TCGA database analysis (http://www.oncolnc.org, accessed on 5 October 2021).

#### *2.6. Lenvatinib Affects microRNA Expression in Sorafenib-Resistant Cells*

A customized microarray platform was used to analyze the expression of 2555 miR-NAs in the lenvatinib-treated or control Huh-7SR cells. Treatment with 10 µM lenvatinib for 24 h upregulated the expression (fold change > 1.5) of 43 miRNAs, such as miR-575, miR-663a, miR-491-5p, miR-4465, miR-371b-5p, and miR-718, and suppressed the expression (fold change < 0.67) of 23 miRNAs, such as miR-4448, miR-106b-3p, miR-197-5p, and miR-130b-3p (Table 1). Unsupervised hierarchical clustering analysis was conducted by calculating Pearson's centered correlation coefficient, and the results indicated that the lenvatinib-treated Huh-7SR cells clustered together (Figure 6). Additionally, filtration

(FDR < 0.001) indicated that lenvatinib upregulated the expression of 16 miRNAs and suppressed the expression of six miRNAs (Figure S2).

**Table 1.** Statistical results and chromosomal locations of miRNAs that exhibited a fold change (FC) > 1.5, FC < 0.67, or *p* < 0.001 in HCC Huh-7SR cells treated with lenvatinib when compared with untreated cells.

**Table 1.** *Cont.*

**Figure 6.** Lenvatinib affects miRNA expression. Hierarchical clustering of differentially expressed miRNAs from Huh-7SR cells incubated with 10 µM lenvatinib or DMSO for 24 h. Fold Change >1.5 or <0.67, *p* < 0.001.

#### *2.7. Altered miRNA in HCC and Normal Tissues and the Relationship with Overall Survival of HCC Patients Based on TCGA Database Analysis*

Based on the altered miRNA expressions, 372 HCC tissues and 50 normal tissues from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) database were analyzed. miR-130b-3p and miR-1292-5p are highly expressed, while the levels of miR-491-5p and miR-1247-3p are decreased in HCC tissues compared with normal tissues. No differences were found in the levels of miR-1228-5p and miR-431-3p between the HCC tissues and the normal tissues (Figure 7A). After lenvatinib treatment, miR-491-5p and miR-1247-3p were upregulated, while miR-130b-3p and miR-1292-5p were downregulated in the Huh-7SR cells. In addition, lower levels of miR-130b, miR-106b, and miR-874, and higher levels of miR-487, are associated with the longer overall survival of HCC patients (Figure 7B).

**Figure 7.** Altered microRNA expression and the relationship with overall survival in HCC. (**A**) Expression distribution of altered miRNAs after lenvatinib treatment in HCC and normal tissues based on TCGA database analysis, where different colors represent different groups. \*\* *p* < 0.01, \*\*\*\* *p* < 0.001, ns, no significant. (**B**) Lower levels of miR-130b, miR-106b, and miR-874, and higher levels of miR-487 are associated with longer overall survival of HCC patients (http://kmplot.com/analysis/, accessed on 5 October 2021).

#### *2.8. Lenvatinib Inhibits Huh-7 Sorafenib-Resistant Cell Proliferation In Vivo*

Next, we examined the effect of lenvatinib in a nude mice xenograft model by injecting Huh-7SR cells. Tumor growth was significantly inhibited in the group treated with lenvatinib 20 mg/kg/day (5 days/week) compared with the sorafenib 30 mg/kg/day (5 days/week) and control groups (Figure 8A). Moreover, there was no difference in body weight among these three groups. Hematoxylin and eosin (H&E) staining and immunohistochemical staining of the ki-67, cyclin D1, and CD31 proteins in the subcutaneous xenograft model suggested that lenvatinib treatment decreased the levels of ki-67, cyclin D1, and the staining area of CD31, and suppressed cell proliferation and angiogenesis (Figure 8B).

**Figure 8.** Lenvatinib inhibits Huh-7 sorafenib-resistant cell proliferation in vivo. (**A**) Lenvatinib 20 mg/kg/day (5 days/week) significantly inhibits tumor growth compared with control and sorafenib 30 mg/kg/day (5 days/week) (n = 8 per group). (**B**) Hematoxylin and eosin (H&E) staining and immunohistochemical staining of ki-67, cyclin D1 and CD31 proteins in the subcutaneous xenograft model. Ki-67-positive cells and cyclin D1-positive cells in the lenvatinib-treated groups were reduced in numbers compared with that in the control group. CD31 staining area in the lenvatinib-treated groups was decreased compared with that in the control group. Data are presented as mean ± standard deviation (SD). \* *p* < 0.05, \*\* *p* < 0.01, ns, no significant.

#### **3. Discussion**

The multi-kinase inhibitor, sorafenib, has been used as first-line therapy for patients with progressive unresectable HCC for a decade. However, resistance to sorafenib limits patient response and presents a major hurdle during HCC treatment. Additionally, data from the REFLECT trial indicated that lenvatinib was the first agent to achieve non-inferiority against sorafenib [13]. Therefore, in the present study, we evaluated the inhibition efficacy of lenvatinib in sorafenib-resistant HCC cells (Huh-7SR and Hep-3BSR cells). Key points include the following: (i) lenvatinib suppressed sorafenib-resistant HCC cell proliferation, mainly by inducing G1 cell cycle arrest; (ii) the underlying advantage of lenvatinib in overcoming sorafenib resistance may occur through the FGFR4-ERK signaling pathway; (iii) along with HBV DNA, poor autophagic responsiveness may be a contributing factor toward partial cross-resistance; and (iv) miRNA alterations may contribute to the inhibition of sorafenib-resistant HCC cell growth and angiogenesis (Figure 9).

**Figure 9.** Proposed model for the underlying mechanism of lenvatinib in overcoming sorafenib resistance. Overexpression of EGFR leads to the activation of ERK and Akt signaling to induce sorafenib resistance and promote sorafenib-resistant cell proliferation. The underlying advantage of lenvatinib is the inhibition of FGFR4 compared with sorafenib. Huh-7SR and Hep-3BSR cells have high FGFR4 expression and the FGFR4-ERK signaling pathway is the major pathway to overcome increased EGFR-induced ERK activation. Compared with Huh-7SR cells, Hep-3BSR cells showed poor autophagic responsiveness to lenvatinib, which may contribute to their partial cross-resistance to lenvatinib.

The ERK signaling pathway plays a key role in anti-tumor effects and multi-kinase inhibitor resistance [28]. In both primary sorafenib-resistance and acquired resistance, the activation of ERK is major process that promotes cell proliferation. In addition, previous studies have revealed that HBx also activates the ERK signaling pathway in HCC [29]. The kinase inhibition profiles of lenvatinib and sorafenib indicated that lenvatinib is more effective at inhibiting FGFR4 (Table S1) [30,31]. Additionally, Huh-7SR and Hep-3BSR cells both expressed FGFR4 and EGFR. Further experiments showed that lenvatinib presented a better performance in ERK signaling inhibition compared with sorafenib. Thus, our findings suggest that lenvatinib overcomes sorafenib resistance mainly through the inhibition of the FGFR4-ERK signaling pathway. Previous studies also indicated that lenvatinib could strongly inhibit the activation of ERK, the downstream signaling molecules of FGFR4, compared with sorafenib and regorafenib [32], and high FGFR4 levels (positive immunohistochemistry >10% of tumor cells) were an independent predictor of a response to lenvatinib [33]. Moreover, lenvatinib enhanced the antitumor immune response of anti-programmed cell death-1 (PD-1) in HCC by blocking FGFR4 [34]. However,

lenvatinib alone could inhibit HCC cancer stem-like cells through FGFR1-3 signaling, but not FGFR4 signaling [35]. Interestingly, lenvatinib could also increase the expression of EGFR, as previously reported [36,37]. The activation of EGFR may contribute to sorafenib or lenvatinib resistance. The combination therapy of lenvatinib and gefitinib (an EGFR inhibitor) may be an option for the approximately 50% of advanced HCC patients with high EGFR expression [37]. Combination therapy or new kinase inhibitors (such as ERK inhibitors) may represent a promising strategy for sorafenib-resistant HCC patients.

Moreover, we found that Huh-7SR cells showed higher autophagic responsiveness to sorafenib and lenvatinib than Hep-3BSR cells, and this may contribute to Hep-3BSR's partial cross-resistance to lenvatinib. Autophagy plays neutral, tumor-suppressive, or tumor-promoting roles in cancer development, and is also associated with apoptosis and multidrug-resistance [38]. Previous studies demonstrated that the autophagic responsiveness to sorafenib is distinct between Hep3B and Huh7 wild-type cells and the sensitivity to sorafenib is also different [39], where sorafenib-induced autophagy improves the death rate of HCC cells [40]. Similarly, in sorafenib-resistant cells, the different autophagic responsiveness to lenvatinib may also be associated with altered sensitivity to lenvatinib.

miRNAs are a class of endogenous, small, noncoding-RNA molecules that regulate aspects of the post-transcriptional modulation of gene expression, such as cell proliferation, differentiation, metabolism, and cell death [41,42]. Additionally, several reports have shown that miRNAs can regulate the sensitivity of HCC cells to multi-kinase inhibitor drugs, such as sorafenib, by modifying diverse molecular processes [43,44]. In this study, we reported altered cell miRNA expression in Huh-7SR cells after treatment with lenvatinib. A total of 43 miRNAs were upregulated, while 23 miRNAs were downregulated. Further data set analysis from TCGA indicated that miR-130b-3p and miR-1292-5p were highly expressed, while the levels of miR-491-5p and miR-1247-3p were decreased in HCC tissues compared with normal tissues. In addition, miRNA expression was reversed after lenvatinib treatment. In many cases, miRNAs can have oncogenic effects or act as suppressors in different cancers, and some may also have a dual function. Liao et al. indicated that miR-130b-3p was upregulated in HCC and was correlated with a poor prognosis. Overexpressed miR-130b-3p was found to enhance the angiogenesis capacity of HCC cells [45]. Decreased miR-491-5p and highly pyruvate kinase M2 (PKM2) expression were associated with unfavorable clinical features and the poor prognosis of HCC patients [46]. In addition, overexpressed miR-491-5p inhibited HCC cell proliferation and migration by targeting SEC61 translocon alpha 1 subunit1 (SEC61A) [47]. The inhibition effect of lenvatinib in sorafenib-resistant HCC cells occurred partially through miRNA regulation. However, due to the scarcity of information on some of the miRNAs, further studies are needed on their function.

Only a few clinical studies have been conducted on second-line or further-line treatment of lenvatinib. Jefremow et al. suggested that in seven patients with a later line lenvatinib treatment, partial remission (PR) was shown in four of the seven patients, stable disease (SD) in two of the seven, and mixed response with overall tolerable safety in one of the seven [48]. Chen et al. enrolled 40 patients who received lenvatinib after sorafenib. The median overall survival (OS) was 9.8 months, and the objective response rate was 27.5%; moreover, the clinical outcomes of lenvatinib treatment in later lines were similar [15]. However, a report by Tomonari et al. indicated that the objective response rate was 33.3% in the second line, and 20.0% in the third line. Additionally, sorafenib-resistant HCC cells show partial cross-resistance to lenvatinib by the decreased response to FGFR signaling pathways compared with wild-type cells [16] (Table 2). Although the number of patients was one of the limitations of the study, an effective tendency was shown.

**Table 2.** Lenvatinib as second-line or further-line treatment after sorafenib.

This study has the following several limitations: (i) We only used two sorafenibresistant HCC cell lines; (ii) expanded clinical studies are needed, including in sorafenibresistant HCC patients with HBV infection, sorafenib-resistant HCC patients with high or low FGFR4 expression, and an investigation of a promising strategy, such as combination therapy with lenvatinib. Hence, we plan to perform further experiments to investigate the above limitations in our future research on this topic.

In conclusion, our study shows that (i) the advantage of lenvatinib in overcoming sorafenib-resistance may be through the FGFR4-ERK signaling pathway; (ii) HBV DNA and poor autophagic responsiveness may be the reasons for partial cross-resistance; (iii) miRNA alterations may contribute to inhibiting sorafenib-resistant HCC cell growth and angiogenesis. The present study and previous clinical data (Table 2) provide evidence that lenvatinib may be a suitable second-line therapy for unresectable HCC patients who express FGFR4 and are sorafenib resistant. Combination therapy could be a promising way to expand the efficiency of first-line multi-kinase inhibitors and our study is one of the initial investigations of this topic. Drug resistance in HCC patients is a huge barrier to overcome. Attenuating first-line treatment resistance and expanding the efficiency of second-line therapy could potentially prolong the survival time of patients. Thus, further studies and clinical trials are needed.

#### **4. Materials and Methods**

#### *4.1. Chemicals*

Sorafenib (Nexavar®) was obtained from Bayer Pharmaceutical Corporation (West Haven, CT, USA) and lenvatinib from Chem Scene (Monmouth Junction, NJ, USA). Matrigel® Matrix (356230) was purchased from Corning Inc. (Lowell, MA, USA).

#### *4.2. Cell Lines and Cell Culture*

Human HCC cell lines, Huh-7 and Hep-3B, were obtained from the Japanese Research Resources Bank (Tokyo, Japan). Huh-7 cells were maintained in Dulbecco's Modified Eagle's Medium (DMEM; Gibco-Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (FBS; Wako, Tokyo, Japan) and penicillin/streptomycin (100 mg/L; Invitrogen, Tokyo, Japan). Hep-3B cells were cultured in Modified Eagle's Medium (MEM; Gibco-Invitrogen) supplemented with 10% FBS and penicillin/streptomycin. To establish sorafenib-resistant HCC cell lines, Huh-7 and Hep-3B cells were cultured with increasing doses of sorafenib, from 1 to 6 µM, for six months. Cells were grown in a humidified incubator at 5% CO2 and 37 ◦C. Huh-7 and Hep-3B cell lines as well as Huh-7 sorafenibresistant (Huh-7SR) and Hep-3BSR sublines were authenticated using short tandem repeat (STR) profiling (BEX Co., Ltd., Tokyo, Japan).

#### *4.3. Cell Viability Assay*

Cell viability assays were performed using the Cell Counting Kit-8 (Dojindo Laboratories, Kumamoto, Japan) according to the manufacturer's instructions. Briefly, the cells

were seeded in 96-well plates at a concentration of 5000 cells/100 µL/well and incubated in a normal growth medium for 24 h. Subsequently, the cells were grown for an additional 24, 48, 72, or 96 h with sorafenib (0.3, 1, 3, 10, and 30 µM), lenvatinib (0.3, 1, 3, 10, and 30 µM), or dimethyl sulfoxide (DMSO). The medium was replaced with 100 µL of fresh medium containing 10% CCK-8 reagent, and the cells were incubated at 37 ◦C for 3 h. The absorbance was measured at 450 nm using a multi-grating microplate reader SH-9000Lab (CORONA Electric Co., Ltd., Ibaraki, Japan). The experiments were repeated three times.

#### *4.4. Three-Dimensional (3D) Tumor Spheroid Assay (3D Culture)*

Multicellular spheroids were generated using the liquid overlay technique (Figure S3) [49].

#### *4.5. Flow Cytometry Analysis of Cell Cycle*

Cell cycle progression was evaluated using a Cell Cycle Phase Determination Kit (Cayman Chemical Company, Ann Arbor, MI, USA). Cells (1.0 × 106 cells/100-mm dish) were treated with 10 µM lenvatinib or DMSO for 24 h. Cells were trypsinized and resuspended in phosphate-buffered saline (PBS) at a density of 106 cells/mL. Approximately 1.0 × 106 cells were stained in 100 µL of PBS with 10 µL of RNase A (250 µg/mL) and 10 µL of propidium iodide (PI) stain (100 µg/mL) and incubated at room temperature in the dark for 30 min. Flow cytometry (FCM) was performed to compare the proportion of lenvatinib-treated and control cells in each phase of the cell cycle. FCM was performed using a Cytomics FC 500 flow cytometer (Beckman Coulter, Brea, CA, USA) equipped with an argon laser (488 nm), and the percentage of cells was analyzed using the Kaluza software version v2.1 (Beckman Coulter). The experiments were repeated three times.

#### *4.6. Apoptosis Analysis*

Lenvatinib-mediated apoptosis was analyzed using an FCM and Annexin V-FITC Early Apoptosis Detection kit (Cell Signaling Technology, Beverly, MA, USA). Cells (1.0 × 106 cells/ 100-mm dish) were treated with 10 µM lenvatinib or DMSO for 24 h. Apoptotic and necrotic cells were analyzed by double staining with FITC-conjugated annexin V and PI according to the manufacturer's instructions. FCM was conducted using a Cytomics FC 500 flow cytometer equipped with an argon laser (488 nm) to compare the proportion of apoptotic cells in the lenvatinib-treated and control groups, and data were analyzed using the Kaluza software version v2.1. The experiments were repeated three times.

#### *4.7. Apoptosis Analysis by Enzyme Linked Immunosorbent Assay (ELISA)*

ELISA was performed to analyze the levels of caspase-cleaved cytokeratin 18 (cCK18) using the M30-Apoptosense ELISA kit (Peviva Ab, Bromma, Sweden). Cells (5000 cells/well) were seeded in 96-well plates and treated with 10 µM lenvatinib or DMSO for 24 h. Subsequently, the cells were lysed in polyoxyethylene octyl phenyl ether (Wako) and analyzed according to the manufacturer's instructions. The experiments were repeated three times.

#### *4.8. Antibody Arrays to Analyze Angiogenesis-Related Proteins*

Cells (1.0 × 106 cells/100-mm dish) were treated with 10 µM lenvatinib or DMSO for 24 h at 37 ◦C and lysed with a protease inhibitor cocktail, PRO-PREP complete protease inhibitor mixture (iNtRON Biotechnology, Seongnam, Korea). The Human Angiogenesis Array Kit (R&D Systems) was used to analyze the angiogenesis-related proteins in the lenvatinib-treated and control cells according to the manufacturer's protocol. Each array was repeated three times to validate the results.

#### *4.9. Invasion Assay*

For invasion assays, 1 × 104 cells were cultured in 100 µL of serum-free DMEM or with 10 µM lenvatinib, placed into the upper membrane chamber (Cell Biolabs, San Diego, CA, USA.) and precoated with extracellular matrix (ECM) proteins. The lower compartment was filled with 150 µL of DMEM containing 10% FBS. After incubation for 24 h in an

atmosphere of 5% CO2 at 37 ◦C, the cells in the upper compartment were removed, and the invested cells in the lower compartment were lysed. Fluorescence was read using an SH-9000Lab plate reader (CORONA) at 480 nm/520 nm.

#### *4.10. Wound Healing Assay*

Cells were seeded in 6-well plates and incubated in DMEM containing 10% FBS until they reached subconfluence. Scratches were introduced in the cell monolayer using a plastic pipette tip. After washing with PBS, 1.5% FBS media containing 10 µM lenvatinib or DMSO was added. The scratched area was photographed using a light microscope after 24 h.

#### *4.11. miRNA Microarray*

Total RNA of cells was extracted using the miRNeasy Mini Kit (QIAGEN) according to the manufacturer's instructions. After confirming the purity and quantity of each sample using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and an RNA 6000 Nano kit (Agilent Technologies), respectively, the samples were labeled using a miRCURY Hy3 Power Labeling kit (Exiqon A/S, Vedbaek, Denmark) and hybridized to a human miRNA Oligo Chip (v.21; Toray Industries, Inc., Tokyo, Japan). Chips were scanned using a 3D-Gene Scanner 3000 (Toray Industries). The 3D-Gene extraction software version 1.2 (Toray Industries) was used to calculate the raw signal intensity of the images. The raw data were analyzed using GeneSpring GX 10.0 software (Agilent Technologies) to assess the differences in miRNA expression between the samples. Global normalization was performed on raw data obtained above the background level. Differentially expressed miRNAs were determined using Welch's *t*-test.

#### *4.12. Colony Formation Assay*

Cells were trypsinized for 3 min and resuspended at a density of 1 × 103/mL. Five hundred microliters were seeded into 6-well plates, and 1.5 mL of DMEM containing 10% FBS and 1 µM lenvatinib or DMSO was added to each well. The plates were incubated at 37 ◦C in 5% CO2, and the medium was changed every three days until conspicuous colonies were observed. Colonies were fixed and stained with 0.1% crystal violet at room temperature for 5 min, and positive colony formation (>50 cells/colony) was evaluated by counting the number of colonies.

#### *4.13. Data Set Analysis*

All patients of LIHC were retrieved from the TCGA data portal (https://portal.gdc. cancer.gov/, accessed on 3 September 2021). The full clinical dataset was downloaded (up to 3 September 2021) and this study meets the publication guidelines provided by TCGA (http://cancergenome.nih.gov/publications/publicationguidelines, accessed on 3 September 2021) [50,51]. OncoLnc (http://www.oncolnc.org, accessed on 5 October 2021) [52] and KM plotter (http://kmplot.com/analysis, accessed on 5 October 2021) [53] were used for survival analysis.

#### *4.14. Western Blot*

The cells were lysed with a PRO-PREP complete protease inhibitor mixture (iNtRON Biotechnology, Korea) and collected supernatants. Protein concentration was measured using a NanoDrop 2000 spectrofluorometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA). Briefly, protein aliquots (10 µg) were separated on precast protein gels (4–20% Mini-PROTEAN TGX Gels; Bio-Rad, Hercules, CA, USA) and transferred onto nitrocellulose membranes. The membranes were blocked with 2% skimmed milk (GE Healthcare) in TBST with 0.1% Tween 20 (cat. No. T9142, Takara Bio Inc., Kusatsu, Shiga, Japan) for 30 min and were then incubated overnight at 4 ◦C with the following primary antibodies: anti-ß-Actin (#66009-1-lg, Proteintech, dilution 1:5000); anti-LC3A/B (#12741, CST, dilution 1:1000); SQSTM1/p62 (#5114, CST, dilution 1:1000); Caspase-3 (#14220, CST, dilu-

tion 1:1000); PARP (#9542, CST, dilution 1:1000); Caspase-7 (#12827, CST, dilution 1:1000); EGFR (#PAI-1110, Thermo Fisher Scientific, dilution 1:1000); p-Akt (ser473) (#4060, CST, dilution 1:2000); Akt (#4685, CST, dilution 1:1000); p-mTOR (SER2448) (#5536, CST, dilution 1:1000); mTOR (#2983, CST, dilution 1:1000); cyclinD1 (#2978, CST, dilution 1:1000); p-ERK1/2 (#4370, CST, dilution 1:1000); ERK1/2 (#4695, CST, dilution 1:1000); p-MEK1/2 (#9154, CST, dilution 1:1000); MEK1/2 (#8727, CST, dilution 1:1000); anti-ß-Actin (#4967, CST, dilution 1:1000); FGFR4 (#8562, CST, dilution 1:1000) in 5% serum (cat. No. 9048-46-8, FUJIFILM Wako, Osaka, Japan). After washing with TBST, the membranes were incubated for 1 h with corresponding horseradish peroxidase (HRP)-conjugated secondary antibodies: anti-mouse (#7076, CST, dilution 1:2000) or anti-rabbit (#7074, CST, dilution 1:2000). The signal was visualized using a chemiluminescent (ECL) kit (cat. No. 45-000-999; Cytiva) and imaged using ImageQuant LAS 4010 (GE Healthcare).

#### *4.15. Xenograft Model Analysis*

The animal study was approved (approval No. 20627, 20640 and 21681) by, and was conducted in accordance with the guidelines set by, the Committee on Experimental Animals of the Kagawa University. Female athymic mice (BALB/c-nu/nu; 5 weeks old; 15–17 g) were purchased from Japan SLC (Shizuoka, Japan). The mice were maintained under specific pathogen-free conditions using a laminar airflow rack and had continuous free access to sterilized (γ-irradiated) food (CL-2; CLEA Japan, Inc., Tokyo, Japan) and autoclaved water. Mice were subcutaneously inoculated with 5 × 106 Huh-7SR cells in the right flank. When the xenografts were palpable with an approximate diameter of 5 mm, we randomly assigned the animals to three groups of eight mice each. These groups were orally administered 20 mg/kg/day lenvatinib, 30 mg/kg/day sorafenib, or vehicle (DMSO and saline) for five days per week. Tumor volume (mm3 ) was calculated as tumor length (mm) × tumor width (mm)2/2. All animals were sacrificed on day 10 of treatment.

#### *4.16. Immunohistochemistry*

Specimens were fixed overnight in 10% formalin, embedded in paraffin, and cut into sections 5 µm thick. Briefly, slides were deparaffinized in xylene, rehydrated through a graded alcohol series, and rinsed in PBS. Depending on the protein target to be revealed, antigen retrieval was achieved by boiling the sections in either 10 mM sodium citrate buffer (pH 6.0) or 1 mM ethylenediaminetetraacetic acid (EDTA; pH 8.5) buffer for 10 min, followed by a 20-min cool down at room temperature. After a blocking step using 5% goat serum and an Avidin-Biotin Blocking Kit (Vector Laboratories, Burlingame, CA, USA), the slides were incubated with specific primary antibodies (cyclin D1, ki-67, CD31) overnight at 4 ◦C. To quench the endogenous peroxidase activity, slides were incubated for 10 min with 3% hydrogen peroxide, and subsequently, the biotin-conjugated secondary antibody was applied at a 1:500 dilution for 30 min at room temperature. Immunoreactivity was visualized using the Vectastain Elite ABC kit (Vector Laboratories) and 3,3-diaminobenzidine or Vector NovaRed (Vector Laboratories) as the chromogen. Finally, the slides were counterstained with hematoxylin.

#### *4.17. Statistical Analysis*

GraphPad Prism software version 6.0 (GraphPad Software, San Diego, CA, USA) was used for all analyses. Unpaired Student's *t*-test and Wilcox test were used to determine statistical significance between different groups. Two-way analysis of variance (ANOVA) or mixed ANOVA was performed to test the comparisons and corrected using Tukey's post hoc test. One-way ANOVA was performed before the Tukey's post hoc test to test the comparisons. Values were considered statistically significant at a *p*-value < 0.05.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/ijms222313071/s1.

**Author Contributions:** Data curation, T.S. and H.I.; investigation, T.S., K.F., H.K., Y.G. and N.N.; methodology, T.S., M.Y., M.N. and K.T.; project administration, T.S., A.M. and T.M.; supervision, S.F., H.Y. and J.T.; writing—original draft preparation, T.S.; writing—review and editing, T.S. and T.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of the Department of Laboratory Animal Science of Kagawa University, Kida, Japan (approval No. 20627, 20640 and 21681).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data supporting the conclusions of the present study have been documented in this article.

**Acknowledgments:** We thank Kayo Hirose, Keiko Fujikawa, Megumi Okamura, and Fuyuko Kokado for their assistance.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **HDAC6-Selective Inhibitor Overcomes Bortezomib Resistance in Multiple Myeloma**

**Sang Wu Lee † , Soo-Keun Yeon † , Go Woon Kim, Dong Hoon Lee, Yu Hyun Jeon, Jung Yoo, So Yeon Kim and So Hee Kwon \***

> College of Pharmacy, Yonsei Institute of Pharmaceutical Sciences, Yonsei University, Incheon 21983, Korea<EMAIL_ADDRESS>(S.W.L<EMAIL_ADDRESS>(S.-K.Y<EMAIL_ADDRESS>(G.W.K<EMAIL_ADDRESS>(D.H.L<EMAIL_ADDRESS>(Y.H.J<EMAIL_ADDRESS>(J.Y<EMAIL_ADDRESS>(S.Y.K.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +82-32-749-4513

† These authors contributed equally to this work.

**Abstract:** Although multiple myeloma (MM) patients benefit from standard bortezomib (BTZ) chemotherapy, they develop drug resistance, resulting in relapse. We investigated whether histone deacetylase 6 (HDAC6) inhibitor A452 overcomes bortezomib resistance in MM. We show that HDAC6-selective inhibitor A452 significantly decreases the activation of BTZ-resistant markers, such as extracellular signal-regulated kinases (ERK) and nuclear factor kappa B (NF-κB), in acquired BTZ-resistant MM cells. Combination treatment of A452 and BTZ or carfilzomib (CFZ) synergistically reduces BTZ-resistant markers. Additionally, A452 synergizes with BTZ or CFZ to inhibit the activation of NF-κB and signal transducer and activator of transcription 3 (STAT3), resulting in decreased expressions of low-molecular-mass polypeptide 2 (LMP2) and LMP7. Furthermore, combining A452 with BTZ or CFZ leads to synergistic cancer cell growth inhibition, viability decreases, and apoptosis induction in the BTZ-resistant MM cells. Overall, the synergistic effect of A452 with CFZ is more potent than that of A452 with BTZ in BTZ-resistant U266 cells. Thus, our findings reveal the HDAC6-selective inhibitor as a promising therapy for BTZ-chemoresistant MM.

**Keywords:** HDAC6; bortezomib-resistance; HDAC6-selective inhibitor; bortezomib; carfilzomib; multiple myeloma; LMP2; combination therapy

#### **1. Introduction**

Multiple myeloma (MM) is characterized by abnormally proliferating plasma cells derived from B cells [1] and ranks second as the cause of death from hematological malignancy [2]. Typical symptoms of MM include anemia, bone destruction, hypercalcemia, infection, and renal failure due to impaired immune function [3]. Recently, the survival rate of MM patients has improved due to the development of autologous stem cell transplantation and novel therapeutic agents that include proteasome inhibitors (PIs) and immunomodulatory drugs [4]. Despite recent advances in MM treatment, most patients fall into cyclic relapse and eventually develop refractory disease due to residual MM [5]. Many researchers have reported the recurrence and drug resistance in MM patients who have previously received chemotherapy, which remains a major obstacle for curing MM. Thus, because there are currently no effective therapies to treat chemotherapy-resistant MM, novel effective treatments for MM need to be identified.

Bortezomib (BTZ, Velcade) is the first therapeutic PI that was approved by the US Food and Drug (FDA) Administration in 2003 for the treatment of MM [6]. BTZ binds to the catalytic subunit of the proteasome with high affinity and inhibits the ubiquitinproteasome system (UPS) [7,8]. Consequently, BTZ indirectly regulates the apoptotic pathway and cell cycle progression, affecting the downstream signaling of UPS [9]. The nuclear factor kappa B (NF-κB) pathway is a typical cell signaling regulated by the UPS. The NF-κB inhibitor (I-κB) masks the nuclear localization signal of NF-κB in the cytoplasm.

**Citation:** Lee, S.W.; Yeon, S.-K.; Kim, G.W.; Lee, D.H.; Jeon, Y.H.; Yoo, J.; Kim, S.Y.; Kwon, S.H. HDAC6-Selective Inhibitor Overcomes Bortezomib Resistance in Multiple Myeloma. *Int. J. Mol. Sci.* **2021**, *22*, 1341. https://doi.org/ 10.3390/ijms22031341

Academic Editors: Daniela Grimm and Valentina De Falco Received: 17 November 2020 Accepted: 26 January 2021 Published: 29 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

When I-κB is degraded by the UPS, NF-κB translocates into the nucleus and transcribes its target genes, such as anti-apoptotic genes. However, BTZ inhibits activation of the NF-κB pathway by preventing the degradation of I-κB [10,11]. BTZ also activates both intrinsic and extrinsic apoptotic pathways [12]. In addition, BTZ causes an accumulation of the ubiquitin-conjugated proteins that enhances endoplasmic reticulum stress with cytotoxicity, resulting in apoptosis induction [13].