Radiomics is a recent technique introduced in medicine to describe characteristics of medical images quantitatively. Radiomics belongs to artificial intelligence (AI) applications, but it is based on the calculation of features using well-defined mathematical formulas applied directly to the image pixel values (or to a filtered version of the original images). The mathematical definitions of radiomic features are based on the distribution and the relationship between pixels and voxels in the images' region of interest. The concept behind this method lies in the fact that the human eye cannot appreciate all the characteristics of a medical image. Haralick et al. [20] described how the textural features, highlighting the behavior of gray levels' dependencies, can identify different areas in an image. Later, textural information was proposed as an application in medical imaging [21,22]. The improvement in hardware calculation power made these techniques able to compute a high number of medical imaging biomarkers in an acceptable span of time; those indices should help the physician during the treatment decision task, allowing a personalized care pathway for different patients. However, these biomarkers are not yet ready to be used in oncology without a robust validation or a demonstration of their reliability [23]. Among them, radiomic indices and feature signatures are increasingly present in the panorama of modern scientific literature [24,25]. The main issue and challenge up to date are to understand how to overcome the limits of this approach [26,27].

To date, in the literature, several studies have investigated the ability of radiomics features in the tumor-healthy tissue differentiation task, both for computed tomography (CT) and positron emission tomography (PET) datasets, as described by Chu et al. [28] that used feature values in a random forest classifier for diagnostic purposes. In another study [29], healthy tissues' features were used as additional information for an automatic segmentation algorithm. More recently, feature-extracted CT images were combined with BED values to predict tumor response to SBRT [30].

Despite the great work done to date, to our knowledge, there is still no characterization of the radiomic features' ability to give specific information about healthy tissue compared to the sick one when machine learning models for prognosis are involved. One of the main challenges in the field of radiomics, which makes its clinical application difficult, is the harmonization of the features to be analyzed.

Our present multicentric work aims to propose a novel concept of harmonizing the CT radiomic signal using a combination derived from both the tumor and the healthy contralateral tissue, to overcome the variability typical in each patient in different conditions (i.e., manufacturer/technical characteristics, acquisition, reconstruction protocol, and different anatomy).

#### **2. Materials and Methods**

#### *2.1. Study Design*

This study was a retrospective multicentric work. It was approved by the Area Vasta Emilia Nord (AVEN) Ethics Committee (ID: 817/2018/OSS\*/IRCCSRE). The study was also approved by the ethics committees of all the participating institutions; it was performed in accordance with the principles of Good Clinical Practice (GCP) in respect of the ICH GCP guidelines, the ethical principles contained in the Helsinki declaration and its subsequent updates. Each patient gave informed consent for joining the study.

#### 2.1.1. Patient Cohort

Patients who underwent SBRT for histologically proven diagnosis of primary earlystage NSCLC were retrospectively collected from January 2010 to December 2019. A multicenter research project named "TEXture Analysis of PET/CT in lung cancer patients treated with Stereotactic body radiation therapy (TEXAS)" was designed to involve seven Italian Centers.

Inclusion criteria were: (1) histologically proven diagnosis of NSCLC; (2) early-stage T1–T3N0M0 (TNM 7th edition); (3) patients who underwent SBRT, with treatment biological effective dose BED10 ≥ 100 Gy; and (4) age > 18 years.

Exclusion criteria were: (1) lung tumor greater than 7 cm; (2) histologically proven diagnosis of small cell lung cancer or metastasis; (3) previous thoracic irradiation; (4) presence of bone, lymph node, or visceral metastatic lesions; (5) patients with secondary pulmonary nodules from non-NSCLC or NSCLC; (6) past non-NSCLC tumors with evidence of active disease at the time of SBRT and synchronous non-NSCLC tumors (arising within six months of SBRT diagnosis of NSCLC) with the exception in both cases of non-melanomatous skin tumors.

The patient cohort was divided into training (76 patients from three centers) and external validation (41 patients from the other four centers) datasets. This strategy for the distribution of centers among datasets was made to balance the two groups according to the patients' outcomes as described in the following sections. The external validation step was a fundamental part of the study in order to confirm the performances obtained in the training phase.

#### 2.1.2. SBRT Details

Conventional computed tomography (CT) simulation scans were obtained. The radiation oncologist contoured gross tumor volume (GTV) on the CT, as part of the therapeutic pathway. A 5–10 mm isotropic margin was added to GTV to generate the planning target volume (PTV). Intensity-modulated radiation therapy (IMRT) was delivered to all patients. The dose normalization ensured that at least 95% of PTV receives 100% of the prescribed dose with a homogeneous distribution. For all patients, ipsilateral and contralateral lung, heart, chest wall, esophagus, spinal cord, and bronchial trees were contoured as organs at risk (OARs).

#### *2.2. Image Acquisition*

All patients included in the study had PET/CT images, previously acquired as part of their care pathways, and a pre-treatment CT used for planning of SBRT. The planning CT acquisition protocols and scanning devices differed among institutions, as reported in Table 1. PET image sets, corrected for attenuation, were acquired no more than three months before the start of the treatment. Patients fasted at least 6 h before the injection 18F-FDG tracer and the serum glucose level measured at the injection time was below 160 mg/mL in all patients. PET examinations were performed 60 min after the intravenous administration of the radiotracer using a specific protocol for each institution shown in Table 1.

**Table 1.** Protocol acquisition parameters for simulation CT and PET examinations stratified for centers. Whenever two scanners were used, a "|" indicated the different configurations.

#### *2.3. Image Segmentation*

Computed tomography and PET image sets were exported in DICOM format into a dedicated research computer for radiomics analysis. For the present study, gross tumor volume contouring was separately performed on the CT (manually, referred to as GTVCT) and PET (automatically, hereinafter named GTVPET) images of the pre-treatment PET/CT studies.

Two radiation oncologists with experience in lung cancer contoured each lesion on every sequential slice of the planning CT using standardized window settings for parenchyma (W = 1600 and L = −600), according to EORTC guidelines [31] for all patients. Regarding GTVPET delineation, the radiation oncologists placed a region of interest (ROI) on the area of tumor FDG uptake on PET images and an automatic contour—consisting of the region encompassed by a given fixed percent intensity level relative to the maximum registered tumor activity (40% of SUV max)—was generated. We decided to use this approach as a previous study showed that GTVPET delineation using this fixed threshold was better correlated with the gross tumor (based on pathologic examination) instead of using as basis the manually delineated GTVCT [32].

In order to perform radiomic feature harmonization, we used an ROI from the healthy tissue. This was obtained by copying the GTVCT into a healthy lung region, i.e., the contralateral lung at the same level of the GTV (named Contra\_Lung). The Contra Lung initial volume was also shifted by 0.6 and 0.3 cm in six directions for a total of 12 shifts, avoiding the inclusion of surrounding tissues of the healthy lung. These shifts had the aim of simulating the uncertainty in the positioning of the healthy ROI (for future reproducibility of the harmonization method). We chose the shifts in accordance with PET image resolution to account for a likely uncertainty in ROI positioning since PET imaging can be used to localize the GTV before the treatment.

An example picture showing the location of the ROIs mentioned above is shown in Figure 1.

**Figure 1.** Visualization of the CT ROIs in a patient. The contralateral ROI was shifted in 12 different positions (shown in red).

#### *2.4. Outcome*

In this study, PFS was considered as the primary endpoint and was converted into a binary outcome, which was set to 1 for patients who were alive and without disease progression at 24 months, 0 otherwise. PFS was defined as the time from the start of the SBRT to documented relapse or death. The use of the 2-year threshold was chosen because it could properly describe the treatment effectiveness. In fact, a preliminary analysis of the Kaplan–Meier curves of PFS after SBRT for our cohort showed that the majority of the progressions occurred in a period ranging from 2 to 3 years.

#### *2.5. Radiomics Analysis*

Our analysis followed the steps defined for our radiomic study (Figure 2), which included image preprocessing. The first phase consisted of spatial resampling to an isotropic voxel size to obtain reproducible and rotationally invariant features. Then, image range re-segmentation updated the ROI voxels according to a chosen intensity range to remove all voxels for which intensity values fall outside the selected intensity range. Finally, the images were discretized by intensity, grouping the original intensity values (256) into specific ranges (bins). The aim was to reduce image noise and computational burden. The intensity discretization process fixed the width of the re-segmentation interval and the bin width, defining a new bin for each intensity interval. Selecting the bin width allowed direct control of the absolute range represented on each bin. The image preprocessing of intensity and spatial discretization is described in Supplementary Material Table S1. Intensity discretization parameters were chosen accordingly to the guidelines proposed by Orlhac et al. [33,34].

**Figure 2.** Radiomic pipeline description of the implemented steps in our evaluation process.

#### 2.5.1. Feature Extraction

After the image preprocessing steps, feature calculation and their extraction was performed. Features (intensity-based, shape-based, and second-order) were extracted from original images and filtered images (using wavelets, Laplacian of Gaussian (LoG), and gamma modifier filters) [35]. Radiomic features were calculated using a homemade software employing the widely used pyRadiomics library in order to apply pre-determined filters to the original images and compute features from the edited results. The list of the extracted radiomic features can be found at https://pyradiomics.readthedocs.io/en/latest/ features.html, (accessed on 15 July 2022).

#### 2.5.2. Harmonization Process

The harmonization process consisted of, for our two available image modalities, calculating features for 14 different ROIs. One of them coincided with GTVCT, the other 13 with the duplicated GTVCT positioned in the healthy region and its shifts, as described in Section 2.3. This allows us to consider operators' variability in the positioning of the healthy ROI. The general idea of employing this harmonization formula was inspired by another work [36], and it is shown in Equation (1):

$$f\_{HARM}(i) = \frac{f\_{GTV}(i) - f\_{HEALTHY}(i)}{\sigma(i)} \tag{1}$$

where: *fHARM*(*i*) is the ith harmonized feature, *fHEALTHY*(*i*) is the median ith feature value calculated on the 13 healthy tissue samples for each patient and modality, and *σ*(*i*) is the difference between the 75th and the 25th percentile of the ith feature distribution. Shape features were not harmonized. We applied this harmonization only to CT data because of its intrinsic dependence on the protocol acquisition parameters. Furthermore, CT is used for a morphologic and anatomical characterization and pixel values are related to a physical characteristic of the tissue (the linear attenuation coefficient). On the other hand, PET, being a functional imaging modality, is less sensitive to low signal changes in spatial coordinates. Especially in this case, for the healthy lung, in PET pixel values, there is no useful physical information regarding a region where we do not register a signal from the radiotracer absorption.

#### 2.5.3. Feature Selection

LASSO feature selection was applied, in which the following function is minimized (Equation (2))

$$\sum\_{i=1}^{n} \left( y\_i - \beta\_0 - \sum\_{j=1}^{p} \beta\_j \mathbf{x}\_{ij} \right)^2 + \lambda \sum\_{j=1}^{p} |\beta\_j| \tag{2}$$

where: *yi* is the observed value, *βj* is the LASSO coefficients, and *λ* ∑*p j*=1 *βj* is the shrinkage penalty [37].

The parameter λ was chosen using 10-fold cross-validation (CV) computing its error. LASSO penalty brings to zero the weight coefficients (*βj*) of irrelevant features not predictive of the chosen outcome. In addition, LASSO handles sets of collinear features by increasing the weight of one of them while setting the other weights to zero. Because the considered outcome was binary, we used a binomial function for LASSO regression. In Table S2 we show the shrinkage penalties for our trained models.

#### 2.5.4. Model Building

The original and harmonized features were used to develop a supervised machine learning binary classifier. A linear support vector machine (SVM, Model 1) [38] and an Ensemble Subspace Discriminant (ESD, Model 2) [39] were trained by optimizing their performance in 10-fold cross-validation in the training dataset.

Linear SVM classifiers provide low generalization error, even with small learning sample datasets. ESD classifiers are used to decide an explicit discriminant subspace of low dimension.

The two described model types were applied to five different combinations of input features: (A) harmonized CT + PET, (B) harmonized CT, (C) original CT, (D) only original PET, and (E) harmonized CT + PET + selected clinical variables in order to assess the effect of harmonization on the performance of the predictive models. The interested reader can find more information in Text S1 in Supplementary Materials.

The clinical variables in method (E) were chosen among the available ones by using Kaplan–Meier survival curves as described in Section 2.5.5.

#### 2.5.5. Statistical Analysis

For each model, a 95% confidence interval (CI) of the AUCs was calculated for the training and external validation sets. Furthermore, accuracy (95% CI are reported), precision, and recall were calculated.

Subsequently, Kaplan–Meier survival curves were computed using the PFS to select the clinical features. A clinical feature exhibiting a *p*-value from a log-rank test less than 0.05 was considered significant and included in model E. Matlab R2021b (Mathworks, Natick, MA) and R (Vienna, Austria), available at https://www.R-project.org (accessed on 15 July 2022), were used to perform the statistical analysis.

The *p*-values related to statistical differences among the AUC values of each model were calculated using two-sided DeLong test.

#### **3. Results**

#### *3.1. Clinical Results*

One-hundred and seventeen early-stage NSCLC patients met the inclusion criteria. The baseline characteristics of the patients are summarized in Table 2. The median age was 78 years and there were more male (72.6%) than female patients. With a median follow-up of 29.8 months, the median PFS was 24.2 months, and 2-year PFS percentage was 51.2%. Median OS and 2-year OS percentage were 28.5 months and 64%, respectively. The clinical characteristics, including age, gender, Charlson comorbidity index (CCI), diffusing capacity of carbon monoxide (DLCO), tumor size, Eastern Cooperative Oncology Group (ECOG) performance status, and biological equivalent dose to PTV, showed no significant differences between the training and external validation cohorts (Table 2).

**Table 2.** Statical analysis of clinical variables. Abbreviations: PS: Performance status according to ECOG scale, BPCO: chronic obstructive pulmonary disease; ADK: Adenocarcinoma, SCC: squamous cell carcinoma, Fr: fraction; RT: radiotherapy VMAT: volumetric arc-therapy; IMRT: intensity modulated radiotherapy, TOMO: Tomotherapy, PTV: planning target volume. *p*-values in bold mean the statistical significance.

**Table 2.** *Cont.*

No clinical or treatment-related features were shown to be significantly related to PFS in the univariate analysis of the whole population, except for gender (*p* = 0.04 in favor of female) and lung site (right vs. left in favor of the right one, *p* = 0.006).

#### *3.2. PFS Models*

The PFS predictive performance of the linear SVM and ESD models using radiomic features and clinical features are reported in Table 3. In Figure 3, all the models are graphically compared considering their confidential intervals. Models using harmonized features and PET (A,E) achieved AUCs greater than 0.70, both in training and validation. The performances of models using CT-only harmonized features (B models) are not confirmed on the validation dataset (AUC training > 0.75, AUC validation <0.60), while adding PET features leads to better stability between the training and validation sets. Only C-type models (original CT-only features) showed a low-mean AUC (< 0.62), both in training and validation.

**Figure 3.** Performances (AUC) of the studied models. The boxplot shows the minimum, maximum, and average values of the bootstrapped 95% CIs.

**Table 3.** Models' results in terms of AUC, accuracy, precision, and recall.

\* AUCs in square brackets are their bootstrapped 95% CIs. \*\* Precision and recall are presented for class 1. \*\*\* *p*-values are calculated with respect to the conditions C1 and C2 for linear SVM and ESD models, respectively. Values in bold mean the statistical significance.

It is worth noting that both A1 and A2 models significantly outperformed C1 and C2 models, both in the training and external validation datasets (*p* = 0.0001, *p* = 0.01 and *p* = 0.02, *p* = 0.046, respectively, for linear SVM and subspace discriminant models), likewise for E and C models (E1: *p* < 0.0001, and *p* = 0.02, and E2: *p* = 0.01, and *p* = 0.02, respectively, for training and external validation datasets). C models outperformed B models, but only in the training dataset (*p* < 0.0001 and *p* = 0.01, respectively, for linear SVM and subspace

discriminant models) and D and C models had the same performance metrics. In summary, models using clinical information (E models) do not add a significant improvement to A models. This effect is also appreciable in Figure 4, where all the *p*-values among the models are represented.

**Figure 4.** *p*-values calculated using the two-sided DeLong test. Numbers in bold mean the statistical significance.

#### **4. Discussion**

In this work, a multi-centric cohort of early-stage NSCLC patients treated with SBRT was used to build and validate predictive models of PFS greater than 24 months using radiomic features from CT and PET exams and clinical information. Several existing works in the literature [40] describe acquisition protocol variabilities in multi-centric studies, which could affect the performance of radiomic models. Since radiomics computes features from the pixel values in the images, differences in acquisition protocol can lead to biased results. The rationale behind our harmonization method lies in the fact that retrospective multicenter radiomic studies are challenging but necessary, as gathering data from several centers for a centralized analysis is complex for legal, ethical, administrative, and technical reasons. Most of the time, the different centers involved do not follow standardized acquisition and reconstruction protocols; therefore, the collected data suffer from intra-, and inter-variability, making radiomic features sensitive to multicenter variability. Our novel harmonization technique aims to reduce the bias caused by the absence of standardized protocols. Generally, feature analysis is performed by calculating them inside an ROI that coincides with the lesion target. Our study aims to tackle this issue by attempting to reduce this effect using the healthy region of the patient as the baseline from which to harmonize the radiomic data computed from the lesion. From our results, the harmonization improves models' performance when it is used on CT image sets. On the other hand, we expected that for PET-only images, the harmonization method is not easily applicable due to the functional aim of this imaging modality. In such a modality, several healthy tissues (i.e., lungs) are not 18-FDG-avid, while a harmonization based on healthy radiotracer accumulation, to our knowledge, has yet to be studied. Our work investigated and evaluated the feasibility of this technique, which could be employed and better analyzed in future studies.

When using original feature values, PET features were preferred over CT during feature selection, resulting in an only PET-based model. Furthermore, in A models, two features from CT were included in the final prediction score (Log\_Sigma30mm\_GLDM\_Small-Dependence-High-Gray-Level-Emphasis (SDHGLE) and Wavelet\_LHH\_NGTDM\_Busyness), which were also selected in B models. In the same way, a subset of selected features in the A

method (Log\_sigma10mm\_GLSZM\_Size-Zone-Low-Gray-Level-Zone-Emphasis (LGLZE), Exponential\_FIRST ORDER\_Median, Square\_GLSZM\_ZoneEntropy (ZE)) is also present among the selected features in the PET-only-based method. This could mean that our approach, also given the higher performance metric of the A method, was able to merge the information hidden in the CT and PET image sets in a multi-centric cohort of patients, highlighting the importance of properly handling hybrid imaging in radiomic models.

Other harmonization techniques were previously described in the literature [41]. Among these, the ComBat harmonization method, which removes batch effects, mostly based on an empirical Bayes framework, is one of the most used. Conversely, the ComBat method has some limitations: for instance, the dimension of the homogeneous group cannot be too few (in our study, it would have not been applicable as most of the centers provided less than 15 patients). Indeed, recently, some methods have been proposed to overcome these limitations, e.g., using bootstrap and Monte Carlo technique to improve robustness in the estimation [42].

Even if Monte Carlo and bootstrap strategies aimed at overcoming the cohort size limitations, the objective of this method still focuses on removing differences in radiomic feature distributions among different labels (corresponding to different centers). Com-Bat, thus, relies on the individual distributions and changes made to feature values are dependent on a group of patients. While we know that there are data supporting the effectiveness of this method (especially in making the feature distributions uniform), we wanted to tackle the multicenter studies issue from a different angle, which is to account for the individual patients' differences (caused both by the scanner/institution protocols and their anatomy) taken directly from the lesion imaging. This renders the method easier, both computationally and for cohort eligibility reasons (which, in the ComBat method, is needed for representing the single center in terms of homogeneity being an assumption of the method). In fact, if validated further, our method can be applied even in heterogenous cohorts since it uses only the single image set of the patient.

Our approach aimed to use all the information present in the CT data, both from cancer lesions and the contralateral healthy tissue, simulating the radiologists' skill in subjectively evaluating a lesion and adding this information in quantifiable and statistical terms (through the radiomic features).

In our work, the well-known and studied concept of delta radiomics was implemented not in a temporal sense but spatially (cancer vs. healthy tissue), which is an approach that, to our knowledge, was not applied in other prognostic oncological works. Traditional radiomics uses absolute values extracted from regions of interest to predict a clinical outcome. On the other hand, delta radiomics predicts a clinical outcome through the combination of radiomic values computed from image sets acquired at different time points (i.e., radiographs to monitor follow-ups or differences between basal PET and interim PET), which is a rationale also used in clinical practice to assess lesion progression (i.e., PERCIST). In our manuscript, we decided to apply delta radiomics not between different time points but between different anatomic locations (healthy vs. tumor tissues). The assumption behind this use of delta radiomics is that each patient can have an intrinsic "baseline" value for a certain radiomic feature (caused by individual anatomy and institution protocols) that needs to be accounted for when building predictive models. Comparison between normal and tumor tissue behavior (even in terms of pixel values) is also common in clinical practice (i.e., SUV values typical of physiological metabolism or HU/density values of healthy tissue). Some authors [43–45] explain that delta radiomics—which is the use of textural indices associated with different time points or anatomical regions—is more successful than traditional radiomics. Our work aims to provide the basis for a framework where the study of simple absolute feature values can make room for the analysis of their relationship to a reference, used as a threshold or as a comparison.

There are several limitations to the current work. Our study suffered a restricted number of patients selected retrospectively. Nonetheless, the patients' number seemed reasonable at the current phase of our study. It assures the homogeneity in terms of pathologies, as including only NSCLC lesions prevented possible biases created by evaluating different diseases, even in the lung anatomical district. Our study highlighted the necessity to monitor and carefully use features related to pixel values and their relationships. In this case, we can assume that the importance of the radiomic features is not only held in their numerical value but also in the intrinsic relationship among those values. The textural indices' ability to perform more complex clinical tasks (i.e., predicting toxicity and its grade) could be further examined in the next phase of our work or in a prospective study design, which could also assess the robustness of our method. We believe that a prospective study will be able to validate these models within a cohort gathered with a better strategy.

Interestingly enough, we found that the improved performance in models employing harmonized features in the training phase was also confirmed in the validation dataset; the use of an external dataset is becoming more and more crucial to radiomics studies to assure and facilitate their introduction in clinical practice.

In our study, no clinical or treatment-related features were shown to be significantly related to PFS, except for gender and lung site. It is well known in the literature that gender is a prognostic factor for PFS [46–48]. Due to the size of our cohort, we did not find significant correlations between PFS and other studied clinical prognostic variables, such as age or histology (also due to the inclusion criteria). To our knowledge, we did not find any other study reporting a significant correlation between PFS and tumor laterality; thus, we will investigate this finding together with our model generalization power in a future prospective study. Indeed, in the literature, many studies showed a significant correlation between some clinical or treatment-related characteristics and outcomes (PFS and OS) and some created predictive models. Among the various statistical prediction models, nomograms can be accurate and feasible prognostic instruments with high utility in estimating individual patient risk and may, thus, help guide treatment decisions in clinical practice. At present, there are some nomograms, based on clinical features, developed for early-stage NSCLC treated with SBRT [49–51], but there is still need for validation of the clinical variables found in those studies and their experimental results in more robust cohorts, such as prospective ones. Therefore, a need exists for a robust recurrence-related prediction model to help select high-risk candidates who may benefit from additional systemic therapies.

In this scenario, a predictive model based on radiomics and clinical and treatmentrelated characteristics can improve the prediction of clinical outcomes, as already demonstrated by other works. We also found out that clinical variables did not improve the radiomics models, but only the proposed harmonization process statistically significantly improves the model's performance.

As previously stated, our future aim is to apply our method to a prospective multicentric cohort to further validate the framework's stability. In addition, other anatomical regions should be explored to generalize the harmonization, even when the healthy area is not so easily defined as in the lung case. Regarding the employment of this method also in PET datasets, it could be interesting to explore the feasibility of applying our harmonization to 18F-FDG-avid anatomical regions, such as the liver or the brain, which are, however, not related to a pathologic response. Such a method could be especially useful where CT-PET is the only exam included in the care pathway of the patient.

#### **5. Conclusions**

A novel strategy of CT data harmonization involving delta radiomics, considering both cancer and healthy tissue in the contralateral lung, was tested and externally validated in a multi-centric study for NSCLC patients, to initially assess its feasibility.

The radiomics models with harmonized features can predict better the selected patient outcome in our cohort, providing valuable additional information to the clinician.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/curroncol29080410/s1, Table S1: The image preprocessing of intensity and spatial discretization, Table S2: Lasso feature selection, Text S1: Models' description.

**Author Contributions:** M.B.: study, design, literature search, data analysis/interpretation, writing original draft, supervision, V.T.: study, design, literature search, data analysis/interpretation, writing original draft, A.B.: Data analysis/interpretation, data validation, reviewing draft, N.C.: literature search, data analysis/interpretation, data validation, figures, writing, M.G.: study design, data collection, S.C.: literature search, reviewing draft, P.B.: resources, data collection, S.L.M.: resources, data collection, E.P.: resources, data collection, M.A.: resources, data collection, reviewing draft, A.R.: resources, data collection, M.S.: resources, data collection, C.P.: resources, data collection, S.U.: resources, data collection, G.M.: resources, data collection, N.G.-L.: resources, data collection, L.F.: resources, data collection, C.I.: resources, data collection, M.I.: resources, data collection, P.C.: resources, data collection, supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was partially supported by the Italian Ministry of Health—Ricerca Corrente.

**Institutional Review Board Statement:** The study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of by the Area Vasta Emilia Nord (AVEN) (protocol code 817/2018/OSS\*/IRCCSRE approved on 16 July 2019).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The training weights of the models proposed in this work (namely, models E and A described in Section 2.5.4) are available in a GitHub repository at the link: https: //github.com/ausl-re/TEXAS, accessed on 15 July 2022. Further instructions on how to perform the harmonization and to use the model will be added.

**Acknowledgments:** A heartfelt thanks to all the people involved in the various centers in collecting and managing the data, especially Simona Marani, Maria Paola Ruggieri, Giulia Mascari, and Cinthia Aristei. A special thanks to the AIRO Lung Group, especially Vieri Scotti, Stefano Vagge, and Alessio Bruni.

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

#### **References**

### *Article* **ViSTA: A Novel Network Improving Lung Adenocarcinoma Invasiveness Prediction from Follow-Up CT Series**

**Wei Zhao 1,†, Yingli Sun 2,†, Kaiming Kuang 3, Jiancheng Yang 3,4, Ge Li 5, Bingbing Ni 4, Yingjia Jiang 1, Bo Jiang 1, Jun Liu 1,6,\* and Ming Li 2,7,\***

**Simple Summary:** Assessing follow-up computed tomography(CT) series is of great importance in clinical practice for lung nodule diagnosis. Deep learning is a thriving data mining method in medical imaging and has obtained surprising results. However, previous studies mostly focused on the analysis of single static time points instead of the entire follow-up series and required regular intervals between CT examinations. In the current study, we propose a new deep learning framework, named ViSTA, that can better evaluate tumor invasiveness using irregularly serial follow-up CT images to avoid aggressive procedures or delay diagnosis in clinical practice. ViSTA provides a new solution for irregularly sampled data. ViSTA delivers superior performance compared with other static or serial deep learning models. The proposed ViSTA framework is capable of improving performance close to the human level in the prediction of invasiveness of lung adenocarcinoma while being transferrable to other tasks analyzing serial medical data.

**Abstract:** To investigate the value of the deep learning method in predicting the invasiveness of early lung adenocarcinoma based on irregularly sampled follow-up computed tomography (CT) scans. In total, 351 nodules were enrolled in the study. A new deep learning network based on temporal attention, named Visual Simple Temporal Attention (ViSTA), was proposed to process irregularly sampled follow-up CT scans. We conducted substantial experiments to investigate the supplemental value in predicting the invasiveness using serial CTs. A test set composed of 69 lung nodules was reviewed by three radiologists. The performance of the model and radiologists were compared and analyzed. We also performed a visual investigation to explore the inherent growth pattern of the early adenocarcinomas. Among counterpart models, ViSTA showed the best performance (AUC: 86.4% vs. 60.6%, 75.9%, 66.9%, 73.9%, 76.5%, 78.3%). ViSTA also outperformed the model based on Volume Doubling Time (AUC: 60.6%). ViSTA scored higher than two junior radiologists (accuracy of 81.2% vs. 75.4% and 71.0%) and came close to the senior radiologist (85.5%). Our proposed model using irregularly sampled follow-up CT scans achieved promising accuracy in evaluating the invasiveness of the early stage lung adenocarcinoma. Its performance is comparable with senior experts and better than junior experts and traditional deep learning models. With further validation, it can potentially be applied in clinical practice.

**Citation:** Zhao, W.; Sun, Y.; Kuang, K.; Yang, J.; Li, G.; Ni, B.; Jiang, Y.; Jiang, B.; Liu, J.; Li, M. ViSTA: A Novel Network Improving Lung Adenocarcinoma Invasiveness Prediction from Follow-Up CT Series. *Cancers* **2022**, *14*, 3675. https:// doi.org/10.3390/cancers14153675

Academic Editors: Hamid Khayyam, Ali Madani, Rahele Kafieh and Ali Hekmatnia

Received: 7 June 2022 Accepted: 20 July 2022 Published: 28 July 2022

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

**Copyright:** © 2022 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/).

**Keywords:** adenocarcinoma; invasiveness; X-ray computed tomography; deep learning; temporal attention

#### **1. Introduction**

Low-dose computed tomography (LDCT) is recommended for lung cancer screening in high-risk populations based on the National Lung Cancer Screening Trial (NLST) report, which is now included in US screening guidelines [1]. Owing to LDCT, more and more early stage lung adenocarcinomas are diagnosed and treated. In clinical practice, most people require follow-up CT scans due to the indeterminate diagnosis or low probability of malignancy on baseline CT. Assessing the changes in size, CT value, and other imaging features can substantially help the diagnosis and invasiveness evaluation of early stage lung adenocarcinomas. However, the evaluation process is tedious and lacks objectivity, which means that radiologists could be overwhelmed by numerous serial CT image evaluations. Moreover, the features indicating malignancy may not be present in the early stages of lung adenocarcinoma. As we know, biological changes may precede morphological changes. Therefore, an efficient tool for objectively evaluating the changes and mining the internal patterns of lung nodules on serial CTs is of great importance.

Deep learning is a thriving data mining method in medical imaging and has obtained surprising results [2–4]. It can efficiently and automatically process medical images and has achieved promising performances on par with clinicians on various clinical tasks, including disease classification, medical image registration, and organ segmentation [5–9]. Previous studies have shown that deep learning could aid clinical decision-making for early lung cancer in disease management and invasiveness prediction [10–14]. However, most prior studies only included single-time CT scan images, while serial CT scan images were not fully investigated. Several powerful deep learning methods have been invented to process serial data, e.g., Long Short-term Memory, Gated Recurrent Unit Network, and Transformer [15–17]. Equipped with the aforementioned tools, a deep learning system can include serial images, better evaluate the biological behavior and changes, and then better predict different clinical events, such as prognosis, therapeutic effect, and subsequent growth patterns [18].

Serial deep learning models have achieved great success in serial data domains, including natural language processing, video classification, and speech recognition [15,19,20]. Nonetheless, it is important to notice that medical serial data such as electronic health records [21] or medical examinations are almost always sampled irregularly in time, separating them from the aforementioned modalities. Since the progression of the disease is strongly correlated with the time intervals between two time points, the asynchronous (sampled irregularly) nature of medical data requires special treatment. For example, by limiting sampling time intervals to 1, 3, and 6 months, deep learning methods proved effective in integrating multiple time points and improving the prediction of lung cancer treatment response [22]. However, this restriction on time intervals still limits the usage of the deep learning method in processing clinical serial data, epically for irregularly serial data.

In this article, we propose ViSTA (Visual Simple Temporal Attention), a deep learning framework capable of predicting the tumor invasiveness of pulmonary adenocarcinomas from Follow-up CT Series. The main contributions are three-fold: First, by introducing a simple temporal attention mechanism, we propose a new deep learning network, named ViSTA, to evaluate the invasiveness of early stage lung adenocarcinoma using irregularly serial CT scans images. ViSTA is able to gather information throughout the entire series and improve the prediction performance. Compared with serial analysis using traditional recurrent neural networks [22], ViSTA is not limited by different time intervals and can process completely irregularly sampled serial data. ViSTA was trained and validated on a dataset of 1121 CT scans from 282 follow-up series and evaluated on a hold-out test set of 113 CT scans from 69 follow-up series. Second, ViSTA delivers superior performances compared with other static or serial deep learning models. ViSTA also outperforms sizebased predictive methods (Volume Doubling Time [23]) by a large margin. Third, ViSTA was proven to achieve higher scores than two junior radiologists and came close to one senior radiologist in the observer study. Our results prove ViSTA's superiority in terms of processing irregularly sampled series and its great potential of being put into clinical practice in reality. Additionally, ViSTA is completely transferrable to other medical imaging tasks where analyzing serial data should yield better performances.

#### **2. Materials and Methods**

#### *2.1. Data Collection*

From January 2011 to October 2017, a search of the electronic medical records and the radiology information systems of the hospital was performed by one author (Yingli Sun). The inclusion criteria are as follows: (1) two or more available CT examinations with thin-slice (≤1.5 mm) images before resection. If there were only two CT examinations, the interval between two scans should be over 30 or more days. (2) Complete pathologic reports. The exclusion criteria for this analysis were: (1) prior treatment before surgery; (2) poor quality CT images; (3) lesions that were difficult to clearly delineate. Finally, a total of 351 nodules from 347 patients (mean age, 58.41 years ±11.79 (SD); range, 22–84 years) were enrolled in the study. Among the 351 lung nodules, 191 nodules were pathologically identified as preinvasive lesions, including 1 atypical adenomatous hyperplasia (AAH), 39 adenocarcinomas in situ (AIS), and 151 minimally invasive adenocarcinoma (MIA); whereas 160 nodules were identified as invasive adenocarcinoma (IA). In total, 1234 serials CT scans of the 351 nodules were enrolled in this study. The median interval between the first and the last CT examinations was 366 ± 500 days (range, 30–2813 days; interquartile range, 165–852 days). The 351 nodules were randomly separated into a training set (245 nodules), validation set (37 nodules), and test set (69 nodules) (see Table 1).

**Table 1.** Number of CT scans/nodules in training, validation, and test set.

#### *2.2. CT Scanning Parameters*

Preoperative chest CT in our department was performed using the following four scanners: GE Discovery CT750 HD, 64-slice LightSpeed VCT (GE Medical Systems, Chicago, IL, USA); Somatom Definition flash, Somatom Sensation-16 (Siemens Medical Solutions, Erlangen, Germany) with the following parameters: 120 kVp; 100–200 mAs; pitch, 0.75–1.5; and collimation, 1–1.5 mm, respectively. All imaging data were reconstructed using a medium sharp reconstruction algorithm with a thickness of 1–1.5 mm.

#### *2.3. Nodule Labeling, Segmentation and Imaging Preprocessing*

A medical image processing and navigation software 3D Slicer (v4.8.0, Brigham and Women's Hospital, Boston, MA, USA) was used to manually delineate the volume of interest (VOI) of the included nodules at the voxel level by one radiologist (Yingli Sun, with 5 years of experience in chest CT interpretation), then the VOI was confirmed by another radiologist (Ming Li, with 12 years of experience in chest CT interpretation). Large vessels and bronchioles were excluded as much as possible from the volume of the nodule. The lung CT DICOM (Digital Imaging and Communications in Medicine) format images were imported into the software for delineation, and then the images with VOI information were extracted with NII format for next step analysis. Each segmented nodule was attributed a specific pathological label (AAH, AIS, MIA, IA), according to the detailed pathological report. Two steps were performed to preprocess CT images before path extraction. First, the whole-volume CT image was resampled to the spacing of 1 mm in all three dimensions to guarantee isotropy. Second, HU values were clipped to the range of (−1000, 400) and normalized to (0, 1) using minimum–maximum normalization. Normalization can accelerate the convergence in the training of the deep learning model and improve its generalization ability.

#### *2.4. Development of the Deep Learning Model*

We developed a deep learning model named ViSTA to classify IA/non-IA lung nodules. The overall architecture of ViSTA is presented in Figure 1. ViSTA first extracts features from CT image patches using a CNN backbone and then integrates information from time series using a lightweight attention module named SimTA [24], which is designed specifically for analyzing asynchronous time series. Details regarding the architecture of ViSTA are provided in Supplementary Section S1, and a single SimTA layer was shown in Figure S1. To avoid overoptimization, we did not heavily tune the hyperparameters of our deep learning model and simply adopted common settings. ViSTA and all its counterparts are trained end-to-end for 100 epochs using the AdamW optimizer [25]. We used a cosine decay learning schedule from 10−3 to 10−6. The batch size of each update was 32. The drop-out probability and weight decay were set at 0.2 and 0.01 to avoid overfitting.

#### *2.5. Counterpart Methods*

For comparison with ViSTA, we conducted experiments on a few of its counterparts:

#### *2.6. Evaluation and Statistical Analysis*

We evaluated the proposed ViSTA model both quantitatively and qualitatively. To evaluate each method's performance, we used a variety of metrics, including accuracy, precision, sensitivity, F1 score, and AUC. Formulas of evaluation metrics are presented in Supplementary Section S1.

To explore the visual representation and interpretability of ViSTA, we followed Simonyan, K. et al. [26] and plotted our model's saliency maps through backpropagation, and investigated the mechanism under ViSTA and where it directed its attention.

#### *2.7. Observer Study*

To further evaluate the performance of ViSTA, we conducted an observer study to compare the performance of radiologists in the same task against other models. In the observer study, all 69 CT series in the test set were evaluated by three radiologists. One is a senior radiologist with 22 years of experience, and the other two are junior radiologists with 5 and 3 years of experience, respectively. Radiologists gave the results based on the evaluation of all available serials CTs. The reviewed results were analyzed and compared with the performance of our proposed model. Radiologists' performances were evaluated using accuracy, sensitivity, precision, and F1 score.

#### **3. Results**

#### *3.1. Performance of Deep Learning Models in Predicting the Invasiveness of Early Lung Adenocarcinoma*

To validate the effectiveness of ViSTA in predicting IA/non-IA nodules, we evaluated its performance using a variety of metrics against its counterparts: VDT (cutoff value set at best Youden index or 400 days), CNN (including CNN-last, CNN-first, CNN-all-first and CNN-all-last), and CNN+LSTM.

Tables S1 and S2 show their performances on the training dataset and validation dataset. Figure 2 provide the ROC curves of all models on the test dataset. Our proposed model outperformed all deep learning models and VDT-based methods in every metric by considerable margins (best among models are highlighted with an underscore). It is worth noting that VDT is far from effective in terms of invasiveness classification. It underperformed almost all deep learning models in terms of AUC, accuracy, and F1 score. Secondly, sequential models (ViSTA and CNN+LSTM) delivered better performances

than CNN models that utilize static data points. ViSTA outperformed CNN+LSTM by considerable margins in all metrics. This performance gap can be attributed to ViSTA's suitability to analyze asynchronous time series. Unlike CNN+LSTM, which treats all time points as if they were regularly sampled, ViSTA takes time intervals into account and is better at processing follow-up series. Furthermore, we trained CNN on all time points (CNN all-first and CNN all-last) to investigate if sequential models gain superiority over larger training datasets. It turned out that ViSTA and CNN+LSTM still outperformed CNN even when it was trained on all data.

**Figure 2.** ROC curves of different models compared with performances of radiologists. The gray dotted line indicates the performance of a random classifier with no predictive ability.

#### *3.2. Performance Comparison against Radiologists*

In the observer study, we compared the performances of ViSTA and its counterparts against three radiologists (Table 2). One is a senior radiologist with 22 years of experience, and the other two are junior radiologists with 5 and 3 years of experience, respectively. All 69 follow-up series from the test set were included in the observer study. We evaluated radiologists' performances using accuracy, sensitivity and precision, and F1 score and compared them against the proposed model. In terms of metrics that require specifying threshold, we chose the threshold that delivers the best Youden Index on the validation set as the cutoff value. Figure 2 plot deep learning models' ROC curves against radiologists' metrics. In terms of accuracy and F1 score, ViSTA scored higher than the two junior radiologists (accuracy of 81.2% vs. 75.4% and 71.0%; F1 score of 81.7% vs. 73.0% and 65.5%) and came close to the senior radiologist (accuracy of 81.2% vs. 85.5%; F1 score of 81.7% vs. 84.8%).

#### *3.3. Visual Presentation Investigation*

To investigate the mechanism of ViSTA, we used a neural network visualization technique [26] to visualize the attention heatmap of the model, which was mostly attributed to the predicted results and potentially correlated to the biological behavior (Figure 3). We took the absolute value of the raw heatmap and clipped it to the range of (0, 0.01) for better visualization and interpretation. In view of the created heatmaps, we can see that the "attention" of the deep learning system was mostly focused on the nodule. Areas surrounding the nodule draw the attention of ViSTA as well, meaning that they also carry valuable information as the nodule does (Figure 3A,B). Figure 3A show a long follow-up series of 11 time points. We observed that heatmaps stay blank in the first half of the

series, during which both nodule volume and IA probability remain relatively stable. In the latter half, heatmaps begin to show along with significant increases in nodule volume and IA probability. Heatmaps are sometimes only lit up at the last time point (Figure 3B). We contribute this to the sudden increase of nodule volume between the third and the fourth time point, which provides sufficient information for the model. This argument is supported by the spike of IA probability at the fourth time point. In some rare cases, heatmaps on all time points are close to invisible (Figure 3C). We conjecture that this is because the lung nodule had almost no progression, which was proven by the fact that both nodule volume and IA probability stayed almost unchanged throughout the entire series.

**Table 2.** The performance of different models and radiologists on the test dataset. The highest among all is highlighted in bold, and the highest among models and VDT (Volume Doubling Time)-based methods is highlighted with an underscore.

**Figure 3.** Visualization investigation of ViSTA. The top row shows CT slices of each time point in the follow-up series. The middle row shows attention heatmaps extracted using the technique proposed by Simonyan, K. et al. [26]. The bottom row masks heatmaps on top of CT slices. (**A**) Attention gradually grew along with the nodule volume and IA probability as the nodule progressed to the end of the series. (**B**) The heatmap only lit up at the last time point as it is considered the one carrying valuable information. (**C**) All time points are allocated with little to no attention, which may be caused by the slow progress of the nodule.

#### **4. Discussion**

In the current study, we proposed a deep learning framework named ViSTA to predict the invasiveness of lung adenocarcinomas using serial CT images. Our results showed that models fed with serial CT images substantially and consistently outperformed models fed with single CT images. Moreover, our proposed model can effectively process asynchronous time series and outperform the traditional serial network, i.e., LSTM. Our models achieved an AUC of 86.4% and an F1 score of 81.7% in the test dataset, which were higher than those of all its counterparts. In the observer study, ViSTA achieved higher accuracy and F1 scores than two junior radiologists (accuracy of 81.2% vs. 75.4% and 71.0%, F1 score of 81.7% vs. 73.0% and 65.5%). When compared with the senior radiologist, our proposed model delivered close performance (accuracy of 81.2% vs. 85.5%, F1 score of 81.7% vs. 84.8%).

Timely and accurately assessing the biological behavior of early stage lung adenocarcinomas has been a continuous focus of attention in clinical practice. In contrast to traditional radiographic features and handcraft features, deeper and higher dimension level features mined by the deep learning method present promising advantages in many tasks, including predicting the invasiveness of the early lung adenocarcinoma. Kim et al. performed a comparison study and revealed that the predictive accuracy of the deep learning method was superior to those of the size-based logistic model [11]. We also analyzed the predictive value of VDT [27], a size-based key parameter in the differentiation of aggressive tumors from slow-growing tumors in clinical practice [24]. Not surprisingly, the performance of our proposed model substantially exceeded that of the VDT-based methods. It indirectly verified the conjecture that a deep learning system could extract and learn deeper and more valuable features, then better discover the biological behavior of the tumors and predict the invasiveness of early stage lung adenocarcinoma.

Although the deep learning method can obtain better performance, most previous studies only used single CT scan data prior to the surgery for training and extracting features, which cannot reveal and learn the internal growth pattern of the nodules. In clinical scenarios, internal growth is a vital component of Lung-RADS, a guideline to standardize image interpretation by radiologists and dictate management recommendations. Including serial CTs can facilitate medical tasks, such as differentiating benign tumors from malignant ones [28] and monitoring and predicting treatment response [22,29]. The discovery of our study supports this. By modelling serial CTs, the predictive performance of ViSTA substantially surpassed its counterparts in analyzing static data. In clinical practice, sequential medical data is generally sampled irregularly, i.e., with different follow-up time intervals. To address the irregular sampling issue, we adopted SimTA in our proposed model to process irregularly sampled time series. This lightweight module enables modeling sequential information in an efficient way. It turned out that the proposed ViSTA significantly outperformed the standard serial framework, i.e., CNN+LSTM, with considerably fewer parameters and less computation and memory footprint. ViSTA can better take advantage of the complete information of all time point CTs by modelling simple yet effective exponentially decay attention in time series. This was proved by our experiments comparing ViSTA, CNN+LSTM, and pure CNN models trained with all time point CTs (CNN-all). ViSTA's superiority over CNN-all proved that its performance gain does not come from a larger training dataset.

In the visualization analysis, we found that ViSTA can drive its attention on the nodule and the surrounding tissue and drop more attention when the probability of invasiveness increases. It can partly explain the mechanism of the deep learning system. We also found some cases where the model appeared to use features close to the nodule, such as the vasculature and parenchyma surrounding the nodule. In fact, peritumoral tissue may possess valuable information, such as tumor-infiltrating status. Features extracted from the peritumoral tissues can improve the efficiency of intramodular radiomic analysis [30,31]. However, we still cannot fully interpret whether the model incorporates other abnormalities such as background emphysema in its predictions. Further investigation using more comprehensive model attribution techniques may allow clinicians to take advantage of the

same visual features used by the model to assess the biological status of tumors. It is worth noting that some of the heatmaps in the time series are completely blue, meaning that the deep learning model allocated close to zero attention to these time points. Though this phenomenon is not completely interpretable, we argue that it can be attributed to these two facts: these time points are too far from the current one, and they lack findings informative for the deep learning model.

Even though the proposed ViSTA proved effective in processing irregularly sampled CT series in our experiments, there are several limitations left untouched. First, due to the difficulty of collecting complete lung nodule follow-up series, we only included data from a single center in this study. In clinical practice, it is preferable if the proposed method generalizes to multiple data domains. Furthermore, it is possible that a single follow-up series contains CT scans from different centers, which would be an important challenge to solve if the proposed model were to be put into clinical usage. In future studies, we will include CT series from external centers to validate the generalization performance of ViSTA. Second, the SimTA module in ViSTA models a simple temporal attention mechanism that monotonically increases weights as the time point gets closer to the current time. However, it is viable to model more complicated attention relations using deep learning models such as Transformeror Informer [15,32]. These temporal models enable capturing non-monotonic and dynamic temporal attention that could be useful in predicting invasiveness. Last but not least, even though we conducted a visual investigation on ViSTA, the interpretation of deep learning model predictions still remains a major challenge. Additionally, the final clinical decision is still up to clinicians to date. In our future research, we will further investigate the underlying mechanism of ViSTA or other similar attention mechanisms.

#### **5. Conclusions**

To summarize, we designed a deep learning model processing irregularly sampled CT series to predict the invasiveness of early stage lung adenocarcinoma from follow-up CT scans. The model achieved promising accuracy comparable with senior experts and better than junior experts and its counterparts. With further validation, the proposed model could better evaluate the invasiveness of early stage lung adenocarcinoma, avoiding aggressive procedures or delayed diagnosis and helping precise management in clinical practice.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cancers14153675/s1. References [24,33–36]. Section S1: Supplementary Methods. Figure S1: The illustration of a single SimTA layer. The colormap of the attention matrix indicates the magnitude of attention allocated in a particular timepoint. Gray represents zero attention. The darker the orange color is, the more attention the timepoint gets; Table S1: Different models' performances in predicting invasiveness on the training dataset. The highest among all is highlighted in bold; Table S2: Different models' performances in predicting invasiveness on the validation dataset. The highest among all is highlighted in bold.

**Author Contributions:** Conceptualization, K.K. and J.L.; Data curation, W.Z., Y.S. and M.L.; Formal analysis, W.Z., Y.S. and J.Y.; Funding acquisition, M.L.; Investigation, W.Z., K.K., Y.S., G.L., B.N., Y.J. and B.J.; Methodology, W.Z., K.K., Y.S., J.Y., Y.J. and B.J.; Project administration, M.L.; Resources, Y.S.; Software, K.K., B.N. and M.L.; Supervision, J.Y., J.L. and M.L.; Validation, G.L.; Visualization, B.N.; Writing—original draft, W.Z. and Y.S.; Writing—review and editing, W.Z., J.L. and M.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the National Natural Science Foundation of China 82102157 (Wei Zhao), 61976238 (Ming Li), Hunan Provincial Natural Science Foundation for Excellent Young Scholars 2022JJ20089 (Wei Zhao), Hunan Provincial Natural Science Foundation of China 2021JJ40895 (Wei Zhao), the Research Project of Postgraduate Education and Teaching Reform of Central South University 2021JGB147 (Jun Liu), 2022JGB117 (Wei Zhao), the Clinical Research Center For Medical Imaging In Hunan Province 2020SK4001 (Jun Liu), the science and technology innovation program of Hunan Province 2021RC4016 (Jun Liu), and the Clinical Medical Technology Innovation Guidance Project in Hunan Province 2020SK53423 (Wei Zhao), Science and Technology Planning Project of Shanghai Science and Technology Commission 22Y11910700 (Ming Li), Science and Technology

Planning Project of Shanghai Science and Technology Commission 20Y11902900 (Ming Li), Shanghai "Rising Stars of Medical Talent" Youth Development Program "Outstanding Youth Medical Talents" SHWJRS [2021]-99 (Ming Li).

**Institutional Review Board Statement:** This retrospective study was approved by the institutional review board of Huadong Hospital, affiliated with Fudan University (No. 2019K134).

**Informed Consent Statement:** The institutional review board waived the requirement for patients' written informed consent.

**Data Availability Statement:** Data are available for bona fide researchers who request it from the authors.

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

#### **References**

### *Article* **Development and Validation of Novel Deep-Learning Models Using Multiple Data Types for Lung Cancer Survival**

**Jason C. Hsu 1,2,3,4, Phung-Anh Nguyen 1,2,3, Phan Thanh Phuc 4, Tsai-Chih Lo 5, Min-Huei Hsu 6,7, Min-Shu Hsieh 8,9, Nguyen Quoc Khanh Le 10,11, Chi-Tsun Cheng 3, Tzu-Hao Chang 2,5,\* and Cheng-Yu Chen 11,12,\***

**Simple Summary:** Previous survival-prediction studies have had several limitations, such as a lack of comprehensive clinical data types, testing only in limited machine-learning algorithms, or a lack of a sufficient external testing set. This lung-cancer-survival-prediction model is based on multiple data types, multiple novel machine-learning algorithms, and external testing. This predicted model demonstrated a higher performance (ANN, AUC, 0.89; accuracy, 0.82; precision, 0.91) than previous similar studies.

**Abstract:** A well-established lung-cancer-survival-prediction model that relies on multiple data types, multiple novel machine-learning algorithms, and external testing is absent in the literature. This study aims to address this gap and determine the critical factors of lung cancer survival. We selected nonsmall-cell lung cancer patients from a retrospective dataset of the Taipei Medical University Clinical Research Database and Taiwan Cancer Registry between January 2008 and December 2018. All patients were monitored from the index date of cancer diagnosis until the event of death. Variables, including demographics, comorbidities, medications, laboratories, and patient gene tests, were used. Nine machine-learning algorithms with various modes were used. The performance of the algorithms was measured by the area under the receiver operating characteristic curve (AUC). In total, 3714 patients were included. The best performance of the artificial neural network (ANN) model was achieved when integrating all variables with the AUC, accuracy, precision, recall, and F1-score of 0.89, 0.82, 0.91, 0.75, and 0.65, respectively. The most important features were cancer stage, cancer size, age of diagnosis, smoking, drinking status, EGFR gene, and body mass index. Overall, the ANN model improved predictive performance when integrating different data types.

**Keywords:** lung cancer; survival; prediction models; real-world data; artificial intelligence; machine learning

**Citation:** Hsu, J.C.; Nguyen, P.-A.; Phuc, P.T.; Lo, T.-C.; Hsu, M.-H.; Hsieh, M.-S.; Le, N.Q.K.; Cheng, C.-T.; Chang, T.-H.; Chen, C.-Y. Development and Validation of Novel Deep-Learning Models Using Multiple Data Types for Lung Cancer Survival. *Cancers* **2022**, *14*, 5562. https://doi.org/10.3390/ cancers14225562

Academic Editors: Hamid Khayyam, Ali Madani, Rahele Kafieh and Ali Hekmatnia

Received: 21 September 2022 Accepted: 10 November 2022 Published: 12 November 2022

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

**Copyright:** © 2022 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/).

#### **1. Introduction**

Lung cancer is the leading cause of cancer deaths worldwide [1]. Globally, there were around 2.21 million new cases of lung cancer and 1.80 million fatalities in 2020 [2]. One study reported that lung cancer incidence and mortality rates were 22.2 and 18.0 per 100,000 people in 2020, respectively [3,4]. Lung cancer can be divided clinically into two types based on histological features: non-small-cell lung cancer (NSCLC) and small-cell lung cancer (SCLC). NSCLC is the most common among them, accounting for 80–90% of lung cancers [5]. Cell deterioration and metastasis are slower in NSCLC than in SCLC. Around 70% of patients are diagnosed at an advanced stage, making surgical resection and complete treatment challenging [6,7].

Artificial intelligence (AI) has been increasingly used in medical research and clinical practice [8,9]. The accurate prediction of disease prognosis and the outcome of drug treatment, which may serve as a reference for treatment decision-making and drug selection, has become an essential topic in the clinical medicine [9,10]. Developing disease-risk and prognosis-prediction models using machine-learning or deep-learning algorithms with big data is a major area of AI-based academic research in the medical field [10,11]. Studies have used machine-learning and/or deep-learning algorithms to develop lung cancer risk and prognosis-prediction models [12–15]. Among them, Lai et al. [16] used 15 biomarkers with clinical data (including gene expression) from 614 patients to develop a deep neural network to predict the five-year overall survival of NSCLC patients.

This study aimed to develop survival-prediction models for lung cancer patients using a large number of samples, different data types, various machine-learning algorithms, and external testing. In addition to the basic clinical data (including demographic information, disease condition, comorbidity, and current medication), we examined the role of laboratory and genomic test results, which are generally not easy to obtain in predicting lung cancer survival. Moreover, we also explored the important predictors for developing prediction models.

#### **2. Methods**

#### *2.1. Study Design and Data Source*

We conducted a retrospective study in which we obtained data from the Taiwan Cancer Registry (TCR) database and the Taipei Medical University Clinical Research Database (TMUCRD). The TCR database was established in 1979 and is managed by Taiwan's Health Promotion Administration, Ministry of Health and Welfare. It covers 98% of Taiwanese cancer patients and includes diagnosis and other related information. The TMUCRD retrieved data from various electronic medical records (EHR) of three hospitals, Taipei Medical University Hospital (TMUH), Wan-Fang Hospital (WFH), and Shuang-Ho Hospital (SHH). The database contains the electronic medical record data of 3.8 million people from 1998 to 2020, including structured data (e.g., basic information of patients, medical information, test reports, diagnosis results, treatment process, surgery, and medication history) and unstructured data (e.g., progress notes, pathology reports, and medical imaging reports) [17]. This study has been approved by the Joint Institute Review Board of Taipei Medical University (TMU-JIRB), Taipei, Taiwan (approval number N202101080). All the data were anonymous before conducting analysis.

#### *2.2. Cohort Selection*

This study selected patients with lung cancer (ICD-O-3 code: C33, C34) from 2008 to 2018 in the TCR database. Exclusion criteria included individuals under 20 years old, SCLC patients, and patients who did not have any medical history in the three hospitals (TMUH, WFH, SHH). Thus, a total of 3714 patients were included in this study, including 960 patients from TMUH, 1320 from WFH, and 1434 from SHH (Figure S1 in the Supplementary Materials).

#### *2.3. Outcome Measurement*

We ascertained the study outcomes using TMUCRD EHR and vital status data from the Taiwan Death Registry (TDR) [18]. We used the diagnosis date of NSCLC as the index date, and the outcome of this study was death within two years following diagnosis. Data were censored at the date of death or loss to follow-up, insurance termination, or the study's end on 31 December 2018.

#### *2.4. Feature Selection*

Based on a literature review and consultation with clinicians, we selected features that may lead to the mortality of NSCLC patients to build prediction models. These features consisted of:

#### *2.5. Development of the Algorithms*

This study established prediction models based on four modes and different algorithms:

This study aims to predict the survival of lung cancer patients; therefore, the problem can be formulated as a classification model as it could occur in the same patients. We used possible machine-learning techniques such as logistic regression (LR), linear discriminant analysis (LDA), light gradient-boosting machine (LGBM), gradient-boosting machine (GBM), extreme gradient boosting (XGBoost), random forest (RF), AdaBoost, support vector machine (SVC), and artificial neural network (ANN). These methods are briefly introduced below.

Logistic Regression (LR): This is a discrete choice model that models the relationship between a response and multiple explanatory variables and is based on the concept of probability [19]. It is widely used and more practical in fields such as biostatistics, clinical medicine, and quantitative psychology. Its Equation (1) is:

$$y = \frac{e^{(b\_0 + b\_1 \mathcal{K})}}{1 + e^{(b\_0 + b\_1 \mathcal{K})}} \tag{1}$$

where *x* is the input value, *y* is the predicted output, *b*0 is the bias or intercept term, and *b*1 is the coefficient for input (*x*). In this study, we used the LR function with the parameter C (inverse of regularization strength) of 0.0001 to reduce the model's overfitting.

Linear Discriminant Analysis (LDA): This is generally used to classify patterns between two classes; however, it can be extended to multiple patterns. LDA assumes that all classes are linearly separable, and according to the multiple linear discrimination functions representing several hyperplanes in the feature space are created to distinguish the classes [20]. In this study, we set the parameters' *shrinkage* to '0' and the *solver* to 'lsqr' to improve estimation and classification accuracy.

Light Gradient-Boosting Machine (LGBM): This is a gradient-boosting framework that uses tree-based learning algorithms. It is designed to be distributed and efficient with the following advantages: faster training speed and higher efficiency; lower memory usage; better accuracy; support of parallel, distributed, and GPU learning; and capability to handle large-scale data [21]. The model's *class\_weight* parameter was set as 'balanced', which uses the output's value to automatically adjust weights inversely proportional to class frequencies in the input data. The *learning\_rate*, l1 regularization—*reg\_alpha*, and l2 regularization—*reg\_lambda* parameters were set as 0.05, 0.1, and 0.1, respectively.

Gradient-Boosting Machine (GBM): Gradient-boosting regression trees produce competitive, highly robust, and interpretable procedures for regression and classification. The ability of TreeBoost procedures to give a quick indication of potential predictability, coupled with their extreme robustness, makes them a useful preprocessing tool that can be applied to imperfect data [22]. The default parameters were used in this model.

Extreme Gradient Boosting (XGBoost): XGBoost, an efficient and scalable implementation of the gradient-boosting framework, is a machine-learning system for tree boosting. The scalability of XGBoost is attributed to several critical systems and algorithmic optimizations. These innovations include a novel tree-learning algorithm for handling sparse data; a theoretically justified weighted quantile sketch procedure allows the handling of instance weights in approximate tree learning [23]. The default parameters were used in this model.

Random Forest (RF): RF is an ensemble-learning method that operates by constructing many small scales of classification modules (most often decision trees) at the training time. The model outputs the class that combines the result of the individual modules based on some voting algorithms [24]. In this study, we set the parameters as follows: *n\_estimators* (the number of trees) of 500, *max\_depth* of 10, *min\_samples\_split* of 400, and *class\_weight* of 0.5 for each class.

AdaBoost: The AdaBoost algorithm is an iterative procedure that combines several weak classifiers to approximate the Bayes classifier C∗(*x*). AdaBoost builds a classifier, e.g., a classification tree that produces class labels, starting with the unweighted training sample. If a training data point is misclassified, the weight of that data point is increased (boosted). A second classifier is built using the new weights, which are no longer equal. Again, misclassified training data have their weights boosted, and the procedure is repeated [25]. The number of estimators (*n\_estimators*) used was 100.

Support Vector Machine (SVC): This is a machine-learning algorithm that can be applied to linear and nonlinear data. SVC transforms the original data to a higher dimension, from which it can use the super vectors in the training data set to find the hyperplane for categorizing the data. An SVC mainly identifies the hyperplane with the most significant margin, e.g., the maximum marginal hyperplane, to achieve higher accuracy [26]. The SVC can be represented by the following Equation (2):

$$f(\mathbf{x}) = \sum\_{i=1}^{N} (\alpha\_i^\* - \alpha\_i) \mathbf{K}(\mathbf{x}, \mathbf{x}\_i) + \mathbf{B} \tag{2}$$

where *K*(*x*, *xi*) is the kernel function, *αi*, *α*∗ *i* ≥ 0 are the Lagrange multipliers, and B is a bias term. In this study, we used a *linear* kernel for computations.

Artificial Neural Network (ANN): This is a learning algorithm vaguely inspired by biological neural networks. Computations are structured in terms of an interconnected group of artificial neurons, and these neutrons process information using a connectionist approach to computation. They are usually used to model complex relationships between inputs and outputs, find patterns in data, or capture the statistical structure [27]. The number of hidden layers with the number of neurons in each layer was set at 3 and 16, respectively. Additionally, for each layer, the *l2 regularization* of 0.01 and the 'relu' *activation* were used in the study. We set the 'softmax' activation for the output layer. We also used the 'Adam' *optimizer*, a highly performant stochastic gradient descent algorithm, and 'binary\_crossentropy' as the binary classification outcome for the *loss* function.

#### *2.6. Evaluating the Algorithms*

The training dataset contained the data of patients from TMUH and WFH. The stratified 5-fold cross-validation was applied in the training set to assess the different machinelearning models' performance and general errors. In other words, patients in the training set were divided into five groups, each used repeatedly as the internal validation set. We recruited data from SHH and used it for the external testing dataset to generalize the model.

The performance of the algorithms was measured by the area under the receiver operating characteristic curve (AUC), accuracy, sensitivity (recall), specificity, positive predictive value (PPV, precision), negative predictive value (NPV), and F1-score. We defined the best model using the highest AUC by comparing various models based on the external testing set. Furthermore, we analyzed the feature's contribution (i.e., the feature's importance) of the best model using SHAP values (SHapley Additive exPlanations) [28].

All the data processing was performed using MSSQL server 2017 (Redmond, WA, USA), and the model training and testing were performed using Python version 3.8 (Wilmington, DE, USA) with scikit-learn version 1.1 (Paris, France) [29].

#### **3. Results**

#### *3.1. Baseline Characteristics of Patients*

We identified 3714 eligible lung cancer patients diagnosed for the first time and registered at the TCR. Among those patients, 2280 patients were included in the training dataset, whereas 1434 were in the testing dataset. Demographic characteristics, comorbidities, tumor size, tumor stage, genomic tests, medication uses, and laboratory tests are presented in Table 1. The mean (standard deviation, SD) ages and BMI of cohort patients were 68 (13.7) and 23.4 (4.33), respectively. Most of the patients were male (57.5%) with late-stage lung cancer (i.e., stage IV, 54.8%), and patients were less likely to smoke (26.7%) or drink (11%). The cohort of patients had comorbidities related to hypertension (19.8%), hyperlipidemia (13.9%), COPD (16.1%), and CVD problems (11.6%). The follow-up durations for the cohort patients were a mean (SD) of 2.25 (2.47) years and a median (interquartile range (IQR)) of 1.41 [0.46–3.04] years. Detailed information is shown in Table S1 in the Supplementary Materials.

**Table 1.** Basic Characteristics of the Study Cohort.

**Table 1.** *Cont.*

**Table 1.** *Cont.*

**Note**: SD, Standard deviation; yrs., Years; IQR, Interquartile Range; BMI, Body mass index; COPD, Chronic obstructive pulmonary disease; PUD, Peptic ulcer disease; CVD, Cardiovascular; DM, Diabetes; BUN, Blood urea nitrogen; HCT, Hematocrit; HGB, Hemoglobin; K, Potassium; MCH, Mean corpuscular hemoglobin; MCHC, Mean corpuscular hemoglobin concentration; MCV, Mean corpuscular volume; Na, Sodium; PLT, Platelet; RBC, Red blood count; WBC, White blood count; a The training set included the data from Taipei Medical University and Wan-Fang hospitals; b The testing set included the data from Shuang Ho hospital.

#### *3.2. The Performances of Different Prediction Models*

The performances of different prediction models are shown in Table 2. In Mode 1, the highest AUC of 0.88 was observed for the ANN model (i.e., accuracy, 0.82; precision, 0.90; recall, 0.75; and F1-score, 0.64), followed by the GBM and RF models with an AUC of 0.83 and 0.82, respectively. In Mode 3, the best performance was found with an AUC of 0.89 for the ANN model (i.e., accuracy, 0.83; precision, 0.89; recall, 0.81; and F1-score, 0.64). The following AUCs were observed 0.85 for LGBM, GBM, and 0.84 for RF models. Moreover, when considering all features in Mode 4, we found that the best model was the ANN model with an AUC of 0.89 (i.e., accuracy, 0.82; precision, 0.91; recall, 0.75; and F1-score, 0.65). Figures 1 and 2 show the ROC curves of different prediction models in four modes. Detailed information on the various models' measurements (i.e., sensitivity, specificity, PPV, NPV, accuracy, and F1-score) is shown in Table S2 in the Supplementary Materials.

**Table 2.** Performance of various Prediction Models by Modes.

**Note**: LR, Logistic Regression; LDA, Linear Discriminant Analysis; LGBM, Light Gradient Boosting Machine; GBM, Gradient Boosting Machine; XGBoost, Extreme Gradient Boosting; RF, Random Forest; SVC, Support Vector Machine; ANN, Artificial Neural Network; \*, Best model based on AUC values.

Figure 3 shows the top 20 important features of the ANN model in Mode 4. The most important features were the cancer stage, size, age of diagnosis, smoking, and EGFR gene. In other words, patients with advanced cancer stage, large cancer size, older age, and smoking behavior had a higher risk of death within two years. The SHAP value presented the important features of the GBM model in Mode 4 and was consistent with the ANN model, such as cancer stage, age at diagnosis, cancer size, and smoking status (Figure S2 in the Supplementary Materials).

**Figure 1.** The Performance of the Prediction Models in the Testing dataset by different Modes. **Note**: (**A**), Mode 1; (**B**), Mode 2; (**C**), Mode 3; (**D**), Mode 4.

**Figure 2.** *Cont.*

**Figure 2.** The Performance of the ANN Prediction Models in the Testing dataset by different Modes. **Note**: (**A**), Mode 1; (**B**), Mode 2; (**C**), Mode 3; (**D**), Mode 4.

**Figure 3.** Feature Importance of the ANN Prediction Model in Mode 4. **Note**: BMI, Body mass index; EGFR, Epidermal growth factor receptor; WBC, White blood cell; PD-L1, Programmed death-ligand 1; COPD, Chronic obstructive pulmonary disease; CCI, Charlson comorbidity index.

#### **4. Discussion**

In recent years, the prediction of cancer patients' survival has attracted the medical community's attention in various countries because it can facilitate medical decision making, strengthen the relationship between doctors and patients, and improve the quality of medical care. Rapid progress in the development of AI based on machine learning has led to more diversified applications of AI in the field of precision medicine. Based on previously published studies on machine-learning algorithms to build prediction models for the survival of lung cancer patients [12,14–16], this study further compared the performance of various novel machine-learning algorithms. In addition, we also analyzed the relationship between the diversity of features and the accuracy of prediction results and determined the most important features affecting lung cancer survival.

Studies using multiple data types and multiple novel machine-learning algorithms simultaneously are limited. In previous studies on lung cancer prediction, most of them used a single machine-learning (e.g., RF [30]) or deep-learning (e.g., NN [14–16]) algorithm or a few basic machine-learning algorithms (e.g., LR, SVM, decision tree, RF, GBM [12,31]) to develop prediction models. Our results showed that the ANN model had the highest AUC value (it was the most suitable tool for survival prediction). In contrast, the AUC value of the traditional LR algorithm exhibited the lowest performance (it had the lowest predictive ability). Lai Y.H. et al. [16] presented a deep neural network to predict the overall survival of NSCLC patients. They obtained a good predictive performance (AUC = 0.82, accuracy = 75.4%) by integrating microarray and clinical data. While only using basic clinical data (demographics, comorbidities, and medications), our predicted model demonstrated a higher performance (ANN, AUC, 0.88; accuracy, 0.82; precision, 0.90, recall, 0.75, and F1-score, 0.64). Furthermore, when combining other variables, such as laboratory and genomic tests, the AUC values of the predicted model were better (based on the external testing, the AUCs of the ANN model in Mode 1 and Mode 4 were 0.88 and 0.89, respectively; the AUCs of LGBM model in Mode 1 and Mode 4 were 0.81 and 0.86, respectively; the AUCs of the RF model in Mode 1 and Mode 4 were 0.82 and 0.85, respectively).

In this study, we explored the variables that might affect the predictive performance of the survival model. As expected, these variables were highly correlated to the mortality of lung cancer patients, such as advanced cancer stage, tumor size, age at diagnosis, and smoking and drinking status [32]. Our findings also showed that lymphocytes, platelets, and neutrophils tests were associated with the likelihood of lung cancer survival [33]. Thus, lymphocytes play an essential role in producing cytokines, inhibiting the proliferation of cancer cells, and provoking cytotoxic cell death [34]. In words, a decrease in lymphocyte count may predict worse survival in cancer patients. Neutrophils are recruited with cytokines released by the tumor microenvironment, enhancing carcinogenesis and cancer progression [35]. Platelets modulate the tumor microenvironment by releasing factors contributing to tumor growth, invasion, and angiogenesis [36]. Another study by Wang J. et al. [37] reported that lung cancer patients with a higher BMI have prolonged survival compared to those with a lower BMI. The same was true for our study's results, which may be due to the poor nutrition and weight loss caused by respiratory diseases [38], such as COPD.

There are limitations to this study. First, although the study used data from various clinical settings (e.g., TMUH and WFH for establishing the prediction model and SHH for conducting an external test) located in the north of Taiwan, the results may not directly apply to lung cancer patients in other regions. Future studies may need to consider validating the model using data from other areas. Second, this study used retrospective data for development and validation. Further experiments with a prospective study design in clinical settings are needed. Third, to obtain a highly accurate prediction, we developed the machine-learning algorithms with binary outcomes (i.e., survival and death) rather than expected continuous outcomes (i.e., length of survival) for the NSCLC patients. Further studies should be conducted with larger sample sizes to deal with continuous outcomes for lung cancer survival.

#### **5. Conclusions**

In summary, to observe the expected survival of NSCLC patients during a two-year period, we designed an artificial neural network model with high AUC, precision, and recall. Moreover, integrating different data types (especially laboratory and genomic data) led to better predictive performance. Further research is necessary to determine the feasibility of applying the algorithm in the clinical setting and explore whether this tool could improve care and outcomes.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cancers14225562/s1, Figure S1: Cohort Selection Process; Figure S2: Feature Importance of the GBM Prediction Model of Mode 4; Table S1: Detailed Demographic Characteristics of Cohort Patients; Table S2: Detailed Performance of various Prediction Models by Modes.

**Author Contributions:** T.-H.C., P.-A.N. and J.C.H. conceptualized and designed the study. P.-A.N., P.T.P. and T.-C.L. collected the data, performed the analysis, and drafted the manuscript. C.-Y.C. and T.-H.C. provided suggestions for the research design and article content. M.-H.H., M.-S.H., N.Q.K.L., C.-T.C. and J.C.H. reviewed all data and revised the manuscript critically for intellectual content. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by Taiwan Ministry of Science and Technology grants (grant numbers: MOST109-2321-B-038-004; MOST110-2321-B-038-004). The funders had no role in the study design, data collection and analysis, publication decision, or manuscript preparation.

**Institutional Review Board Statement:** This study has been approved by the TMU-Joint Institutional Review Board (Project number: TMU-JIRB N202101080).

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The authors obtained data from the Taiwan Cancer Registry (TCR) database and the Taipei Medical University Clinical Research Database (TMUCRD).

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

#### **Abbreviations**

#### **References**

### *Article* **Prognostication in Advanced Cancer by Combining Actigraphy-Derived Rest-Activity and Sleep Parameters with Routine Clinical Data: An Exploratory Machine Learning Study**

**Shuchita Dhwiren Patel 1,\*, Andrew Davies 2, Emma Laing 3, Huihai Wu 3, Jeewaka Mendis 4 and Derk-Jan Dijk 5,6**

**Simple Summary:** Survival prediction is an important aspect of oncology and palliative care. Measures of night-time relative to daytime activity, derived from a motion sensor, have shown promise in patients receiving chemotherapy. Measuring rest-activity and sleep may, therefore, result in improved prognostication in advanced cancer patients. Fifty adult outpatients with advanced cancer were recruited, and rest-activity, sleep, and routine clinical variables were collected just over a one week period, and used in machine learning models. Our findings confirmed the importance of some well-established survival predictors and identified new ones. We found that sleep-wake parameters may be useful in prognostication in advanced cancer patients when combined with routinely collected data.

**Abstract:** Survival prediction is integral to oncology and palliative care, yet robust prognostic models remain elusive. We assessed the feasibility of combining actigraphy, sleep diary data, and routine clinical parameters to prognosticate. Fifty adult outpatients with advanced cancer and estimated prognosis of <1 year were recruited. Patients were required to wear an Actiwatch® (wrist actigraph) for 8 days, and complete a sleep diary. Univariate and regularised multivariate regression methods were used to identify predictors from 66 variables and construct predictive models of survival. A total of 49 patients completed the study, and 34 patients died within 1 year. Forty-two patients had disrupted rest-activity rhythms (dichotomy index (I < O ≤ 97.5%) but I < O did not have prognostic value in univariate analyses. The Lasso regularised derived algorithm was optimal and able to differentiate participants with shorter/longer survival (log rank *p* < 0.0001). Predictors associated with increased survival time were: time of awakening sleep efficiency, subjective sleep quality, clinician's estimate of survival and global health status score, and haemoglobin. A shorter survival time was associated with self-reported sleep disturbance, neutrophil count, serum urea, creatinine, and C-reactive protein. Applying machine learning to actigraphy and sleep data combined with routine clinical data is a promising approach for the development of prognostic tools.

**Keywords:** biomarkers; circadian; machine learning; palliative care; prognosis; survival

### **1. Introduction**

Prognostication (i.e., estimation of survival) is an important aspect of the management of patients with cancer. It is of particular importance in advanced cancer where it has

**Citation:** Patel, S.D.; Davies, A.; Laing, E.; Wu, H.; Mendis, J.; Dijk, D.-J. Prognostication in Advanced Cancer by Combining Actigraphy-Derived Rest-Activity and Sleep Parameters with Routine Clinical Data: An Exploratory Machine Learning Study. *Cancers* **2023**, *15*, 503. https://doi.org/ 10.3390/cancers15020503

Academic Editors: Hamid Khayyam, Ali Madani, Rahele Kafieh and Ali Hekmatnia

Received: 10 August 2022 Revised: 23 December 2022 Accepted: 6 January 2023 Published: 13 January 2023

**Copyright:** © 2023 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/).

immediate implications for clinicians' decisions about the treatment of the cancer, treatment of co-morbidities, so-called "ceilings of care", and referral to palliative care services [1,2]. Furthermore, it has implications for patients (and families) in terms of current decisionmaking, advance care planning, and "getting one's affairs in order".

Healthcare professionals are inaccurate prognosticators, often overestimating survival [3], and the accuracy of estimates is inversely related to survival [2]. Healthcare professionals are relatively good at predicting if patients will die within a couple of days, but not so good at predicting if patients will live for a couple of months or longer.

Various prognostic tools/algorithms have been developed to improve prognostication in patients with cancer [2,4]: these tools vary in their content (e.g., objective items only; subjective items only; objective and subjective items). However, none of these tools have been shown to be consistently better than clinicians' predictions of survival [2]. Current prognostication tools often include measures such as performance status, symptoms, venous blood sample data, and clinician-predicted survival [2,5]. The integration of other physiological and behavioural parameters, such as rest-activity rhythms ("diurnal or circadian") and sleep parameters are yet to be considered in prognostic models. (The term 'circadian' is meant to refer to rhythms that persist in constant conditions. Rhythms assessed in the presence of environmental rhythms, as in the present study, are referred to as diurnal or 24 h rhythms, although increasingly these rhythms are also referred to as 'circadian')

Sleep-wake cycles and circadian rhythms have a key role in sustaining normal body function and homeostasis [6]. Deterioration of rest-activity rhythmicity (loss of rhythmicity) and fragmentation of the sleep-wake cycle may be a marker of deterioration of health and, indeed, a predictor of illness including cancer, as well as cancer survival [7–9].

Several studies in cancer patients have incorporated actigraphy to objectively assess daytime activity, 24 h variation in rest-activity, as well as nocturnal and daytime sleep [7]. A number of actigraphy-derived parameters have been used to quantify rest-activity rhythms in this population including acrophase (time of peak activity), amplitude (peak to nadir difference, i.e., height of activity rhythm peak), mesor (average activity over a 24 h period), and the "dichotomy index" (I < O). Of these parameters, the I < O is one of the most commonly studied rest-activity measures in cancer studies. The I < O has been identified as an independent prognostic biomarker for overall survival, particularly in patients with metastatic colorectal cancer [10,11]. The I < O is defined as the percentage of the activity counts measured when the patient is in bed that are inferior to the median of the activity counts measured when the patient is out of bed [12]. An I < O of ≤97.5% is indicative of a disrupted rest-activity circadian rhythm (i.e., increased fragmented sleep and reduced daytime activity patterns) [7]. However, the I < O has not been used to prognosticate per se, either alone or in combination with other items. Furthermore, few studies have explored the potential of actigraphy-derived sleep parameters as prognostic markers in advanced cancer patients [13].

The first aim of this study was to investigate the feasibility of using I < O and other actigraphy-derived parameters as stand-alone items, to prognosticate in patients with advanced cancer. The second aim of the study was to determine whether the I < O and other actigraphy and sleep parameters should be combined with established prognostic indicators, e.g., Eastern Cooperative Oncology Group performance status (ECOG-PS), modified version of the Glasgow Prognostic Score (mGPS), Prognosis in Palliative Care Study (PiPS) –B, as well as putative prognostic variables from routine clinical data derived from blood samples, to improve prognostic accuracy. To achieve this second aim we deployed regularised regression, a supervised machine learning approach which overcomes some of the limitations of classical multiple regression, to identify effective prognostic indicators and develop more robust prognostic algorithms [14].

#### **2. Materials and Methods**

#### *2.1. Study Design and Setting*

The study was a prospective observational study conducted in a medium-sized district general hospital/cancer centre in the United Kingdom. The study was sponsored by the Royal Surrey County Hospital and received ethical approval from the London–Bromley REC (reference number—16/LO/0243). The study was registered on the CancerTrials.gov registry (reference number—NCT03283683). The study was funded by the Palliative Care Research Fund (Prof. Davies—Royal Surrey County Hospital), including an unrestricted donation from the family of Mr. John Spencer.

#### *2.2. Study Participants*

Participants were recruited from outpatients at the study site. All patients that met the criteria for the study were eligible for entry into the study (convenience sampling, consecutive recruitment). The inclusion criteria were: (a) age ≥ 18 years; (b) diagnosis of locally advanced/metastatic cancer; (c) clinician estimated prognosis of more than 2 weeks but less than 1 year; and (d) known to a specialist palliative care team. The exclusion criteria were: (a) cognitive impairment; (b) physical disability that affected general activity; and (c) physical disability that affected non-dominant arm movement.

Patients were diagnosed with locally advanced/metastatic cancer according to NHS guidelines, which consider TNM staging. All patients who met the inclusion criteria were deemed eligible for entry into the study. Potentially eligible patients were identified by the clinical team and approached by a member of the research team and invited to participate in the study. Any patient referred to the specialist palliative care team was expected to die within the next twelve months (as per the General Medical Council definition for end-of-life care [14]).

#### *2.3. Routine Data Collection*

Written informed consent was obtained from participants prior to entry into the study. The initial review (day 0) involved a collection of routine clinical data: patient demographics, information about cancer diagnosis/treatment, information about co-morbidities/ medication, assessment of Eastern Cooperative Oncology Group performance status (ECOG-PS) (by clinician and patient) [15], and completion of the Abbreviated Mental Test Score [16], the Memorial Symptom Assessment Scale—Short Form (MSAS-SF) [17], and the Global Health Status question from the PiPS-B algorithm [18]. The participant's pulse was measured (as part of the PiPS), and a venous blood sample was taken to measure haemoglobin, white blood cell count (WBC), neutrophil count, lymphocyte count, platelet count, sodium, potassium, urea, creatinine, albumin, alanine aminotransferase (ALT), alkaline phosphatase (ALP), and C-reactive protein (CRP). The final review (day 8) involved further assessment of ECOG-PS (by clinician and patient), completion of the MSAS-SF, the Pittsburgh Sleep Quality Index (PSQI) [19], and a patient acceptability questionnaire. The blood test results were used to complete the PiPS-B scoring algorithm, and serum CRP and albumin were used to calculate the mGPS [20].

#### *2.4. Wrist Actigraphy and Consensus Sleep Diary*

Wrist actigraphy was used to measure physical activity and standard sleep measures. Participants were fitted with the Actiwatch Spectrum Plus® (Philips Respironics, Bend, OR, USA) on the non-dominant arm after the initial review (day 0) and were instructed to wear the device for eight consecutive 24 h periods. The Actiwatch Spectrum Plus® is a CE-marked device with an accelerometer (i.e., motion sensor) that samples movement at 32 Hz [21] with a sensitivity of 0.025 G (at 2 count level). Participants were also given a Consensus Sleep Dairy in order to provide confirmatory information about specific sleep parameters (e.g., number of awakenings, time of final awakening) [22]: the "diary" was completed for eight consecutive sleep periods. The Actiwatches were configured and data were retrieved using device-specific software (Actiware version 6.0.9: Philips Respironics, Bend, OR, USA). The Actiwatches were adjusted to provide an epoch length (sampling interval) of one minute, which is the most common epoch length used in studies of cancer patients [23]. The Consensus Sleep Diary was used in conjunction with the Actiwatch to assist in actigraphy data interpretation (i.e., determine the major sleep/wake periods) [24].

The data from the Actiwatches was downloaded into an Excel spreadsheet, and the following rest-activity parameters were calculated using a study specific SAS programme (SAS® Version 9.4 Statistical Analysis Software, SAS Institute, Cary, NC, USA): I < O, r24 (an autocorrelation coefficient at 24 h, that is "a measure of the regularity and reproducibility of the activity pattern over a 24 h period from one day to the next") [25], mean daily activity (MDA), and mean activity during daytime wakefulness. MDA was calculated as the average number of wrist movements per minute throughout the recording time [25], and the mean duration of activity during wakefulness was calculated as the mean activity score (counts/minute) during the time period between two major sleep period intervals [26]. In addition, the following sleep parameters were calculated both automatically from the Actiwatches (using the Actiware sleep scoring algorithm) and manually from the sleep diary [27]: bedtime (BT), get-up time (GUT), time in bed (TIB), sleep onset latency (SOL), total sleep time (TST), sleep efficiency (SE), wake after sleep onset (WASO), and number of awake episodes (NA). The sleep parameters derived manually solely from the sleep diary were: time tried to sleep, time of final awakening and terminal awakening (TWAK) [22]. See Table 1 for definitions of the sleep parameters.

**Table 1.** Definitions of actigraphy-derived sleep/consensus sleep diary parameters [22,26,28].

#### *2.5. Follow-Up*

During the study period (from time of first patient recruited to six months after last patient recruited), participants' survival status (and date of death, if applicable) was determined every three months by reviewing the hospital clinical records, and/or contacting the general practitioner.

#### **3. Statistical Analyses**

The sample size for the study (*n* = 50) was derived from guidance on sample sizes for feasibility studies (and represents the upper range) [29]. Statistical support was provided by statisticians, within the Research Design Service South-East (based in the Clinical Trials Unit at the University of Surrey). Descriptive statistics were used to explain much of the data (e.g., mean and standard error; median and range). The Intraclass Correlation Coefficient (ICC) was used to assess the robustness of I < O as a marker of the restactivity rhythm, and its stability throughout the actigraphy recording. The Spearman's Rank correlation coefficient was used to measure the association between I < O and other actigraphy-derived parameters. The Spearman's rank correlation '*r*' values were defined as follows: 0 ≤ *r* < 0.3 indicated a negligible correlation, 0.3 ≤ *r* < 0.5 a low correlation, 0.5 ≤ *r* < 0.7 a moderate correlation, 0.7 ≤ *r* < 0.9 a high correlation, and 0.9 ≤ *r* ≤ 1 a very high correlation [30]. Kaplan–Meier plots, a non-parametric statistical method, were used to estimate the probability of survival past a given time point along with the log rank test to compare the survival distribution of two groups. Statistical significance was evaluated at 5%.

The "per protocol set" refers to participants that wore the Actiwatch for the eight consecutive 24 h periods with the corresponding sleep diary, whilst the "full analysis set" refers to participants that wore the Actiwatch for at least three consecutive 24 h periods (i.e., 72 h) and completed the corresponding sleep diary for the actigraphy rest-activity and sleep analysis, or for at least three consecutive or non-consecutive nights in the sleep diary for the subjective sleep analysis (i.e., calculation of the sleep diary parameters).

#### **4. Machine Learning Methods and Data Analysis**

Cox regression has been the standard approach to survival analysis in oncology. However, Cox regression has a number of limitations. In particular, it is not an adequate approach for situations in which the number of predictors is high relative to the number of observations, as is the case in this feasibility study. We therefore opted to use simple alternative methods that can (1) adequately deal with situations in which the number of predictors is large relative to the number of observations and (2) yield models that are interpretable, i.e., are not 'black box models'. Penalised (Regularised) regression models represent such an approach.

A supervised machine learning algorithm was used to develop a predictive model, where the collated subjective and objective parameters (i.e., routine clinical data and actigraphy-derived rest-activity and sleep parameters) were individual predictor variables and survival was the 'response' variable [31]. Sixty-six predictor variables were tested for potential predictive value (Appendix A, see Table A1 for descriptive statistics of the numerical predictor variables). Overall survival was defined as the time from initial review (day 0) to death or until 14 May 2020 for patients that remained alive until the end of the study.

#### *4.1. Machine Learning Dataset*

All patients recruited into the study (*n* = 50) were used for the machine learning analysis. The predictor variables were classified into the relevant variable type (e.g., binary, categorical\_nominal, etc.) and entered into a .csvfile in Excel. Binary variables, such as 'use of opioid analgesia' were transformed into dummy variables (0 or 1). Categorical\_ordinal variables with a numerical ranking, such as ECOG-PS were labelled using the 'LabelEncoder' approach, where the output integer value from the LabelEncoder function was

used to reflect the ordering of the original integer. Categorical\_ordinal variables with non-numerical values, such as PSQI sleep disturbance, were assigned a numerical ranking. Numerical\_continuous variables involving sleep/wake times were entered in the 24 h format. Missing data values were imputed with the average of the group or with the corresponding subjective/objective data from the same participant. Missing data accounted for <4% of the dataset.

#### *4.2. Regularised Regression Methods*

Regularised regression was used to reduce "overfitting" and aid the generalisability of the model. 'Regularisation' corresponds to a penalty that limits the overall weight that can be assigned across all predictor variables in the model, which reduces model complexity (compared to traditional multivariate regression). For some regularised regression approaches, the penalty can drive the weight of a variable to zero, effectively selecting the optimal combination of predictor variables that can be used to predict the given outcome.

Here, three regularised multivariate regression methods were applied and compared: ridge regression, least absolute shrinkage and selection operator (Lasso) and elastic net. The ridge regression algorithm includes all the predictor variables, shrinking the coefficients towards (but not set at) zero in a continuous manner [32]. The Lasso-derived algorithm combines the method of shrinkage with the sub-selection of predictor variables, using a penalty '*L*1 norm' [32,33], creating a 'sparse' model (i.e., selecting only a few variables from the dataset) [32]. The elastic net algorithm is broadly a combination of the ridge and Lasso [34]. This method simultaneously performs continuous shrinkage and feature selection, selecting groups of correlated variables, using a penalty of '*L*1 norm' and '*L*2 norm' [34]. Highly correlated predictor variables are averaged and entered into the model to remove any deviances caused by extreme correlations [35]. Since survival data are censored, i.e., at the end of the observation period some participants may still be alive, we applied regularised Cox regression using the glmnet package in R.

#### *4.3. Model Development*

The models were validated using a *k*-fold (10 folds used) cross-validation approach [32]. For each of the 50 individuals, the predicted survival was based on a model which was constructed on '*k* − 1 subjects, i.e., the model was blind to the participant and the participant did not contribute to the estimation of the prediction. All analyses were carried out within the statistical computing environment R (version 3.6.2). For machine learning, ridge, Lasso and elastic net (alpha = 0.5) regression the package glmnet (version 2.0) was used. Here, an exhaustive search for lambda able to produce the minimum Mean Cross-Validated Error (CVM) was performed. All subjects were used as the training set to build a final model, then *k*-fold cross-validation for performances (CVM) was performed. Analyses were performed with different settings of elastic net mixing parameter (alpha), which were elastic net (alpha = 0.5), Lasso (alpha = 0.99) and ridge (alpha = 0.01). The models generated a predicted hazard, which was compared to the actual survival in days using Pearson's correlation coefficient. To estimate the intra-variable variation in their contribution to the predictor, we computed the mean cross-validated error of the weights of each of the variables that were consistently identified in all 50 participants.

#### **5. Results**

A total of 50 patients were recruited to the study, and 49 participants completed the study (Figure 1): the full analysis set consisted of 44 participants, whilst the per protocol set consisted of 37 participants. See Table 2 for characteristics of the participants. A total of 46 participants were followed up for 12 months (40 in the full analysis set, 33 in the per protocol set), and 34 died within this time period (28 in the full analysis set, 22 in the per protocol set). Unless otherwise stated, the following results relate to the full analysis set.

**Figure 1.** Study flow chart.

**Table 2.** Participant characteristics.

Note: Percentages may not sum to 100 due to rounding.

#### *5.1. Acceptability of Actigraphy and Sleep Diary Acceptability*

Actigraphy data were missing from one participant due to a technical problem. Fortytwo (84%) participants reported that the Actiwatch was "comfortable to wear", and only four (8%) reported that the Actiwatch interfered with their normal activities. No adverse

effects were reported from using the Actiwatch. Fourteen (28%) participants reported that the Consensus Sleep Diary was difficult to complete, and two (4%) subjects reported that the diary interfered with their normal activities.

#### *5.2. Univariate Analyses of Actigraphy Parameters*

5.2.1. Characteristics of the Dichotomy Index (I < O) and Correlation with Other Actigraphy and Sleep Parameters

Table 3 shows the results for the I < O. Forty-two (95%) participants had anI<O of ≤97.5%, indicating a disrupted rest-activity circadian rhythm [7]. The I < O can be considered a stable variable since the intraclass correlation coefficient for values obtained over eight days using the per protocol set, was 0.93 (95% CI: 0.88–1.00; *p* < 0.0005), which is considered an "excellent" correlation [36]. In fact, there was a "high" positive correlation between the I < O for the first three days (72 h) and for the full eight days (Spearman's correlation: *r* = 0.82; *p* < 0.0005) [31]. Moreover, there was a "high" positive correlation between the I < O on weekdays and on the weekend (Spearman's correlation: *r* = 0.76; *p* < 0.0005). Additionally, there was a "very high" positive correlation between the I < O calculated using 24 h of data, and theI<O calculated using 20 h of data, i.e., excluding the one-hour periods before/after going to bed, and the one-hour periods before/after getting out of bed (Spearman's correlation: *r* = 0.98; *p* < 0.0005).

**Table 3.** Dichotomy Index (I < O) data.

There was a "moderate" positive correlation between the I < O and the r24 (Spearman's correlation: *r* = 0.66; *p* < 0.0005), and the mean activity during wakefulness (Spearman's correlation: *r* = 0.51; *p* < 0.0005). However, there was only a "low" positive correlation between the I < O and the mean daily activity (Spearman's correlation: *r* = 0.43; *p* = 0.003). Other standard actigraphy parameters correlated with the I < O were SE, i.e., number of minutes of sleep divided by total number of minutes in bed (Spearman's correlation: *r* = 0.47, "low" correlation; *p* = 0.001), and WASO, i.e., number of minutes awake after sleep onset during sleep period (Spearman's correlation: *r* = −0.51, "moderate" correlation; *p* < 0.0005).

#### 5.2.2. I < O: Predictor of Survival and Correlation with ECOG-PS

Amongst participants that completed one year of follow-up (*n* = 40), there was no significant difference in overall survival between those separated into two groups (based on the median I < O; log rank test, *p* = 0.917), or four groups (based on the quartiles of the I<O; log rank test, *p* = 0.838). However, I < O had a "moderate" negative correlation with the physician assessed ECOG-PS (Spearman rank correlation: *r* = −0.63; *p* < 0.0005). The ECOG-PS was an independent prognostic indicator in this cohort of patients (log rank test, *p* < 0.0005). The median survival for participants with an ECOG-PS of 1 (end of study) was 141 days, ECOG-PS of 2 was 135 days, ECOG-PS of 3 was 57 days, and ECOG-PS of 4 was 17 days.

#### 5.2.3. Autocorrelation Coefficient at 24 h (r24)

The median r24 was 0.16 (range 0.04–0.37). Amongst participants that completed one year of follow-up (*n* = 40), there was no significant difference in overall survival between those separated into two groups (based on the median r24; log rank test, *p* = 0.318), or four groups (based on the quartiles of the r24; log rank test, *p* = 0.800).

#### 5.2.4. Other Actigraphy Parameters

None of the other actigraphy-derived sleep parameters were associated with a decreased overall survival: (a) TIB (log rank, *p* = 0.574: based on group median of 9 h 29 min); (b) TST (log rank, *p* = 0.147: based on normative cut-off value of ≥6.5 h [28]; (c) SOL (log rank, *p* = 0.283: based on normative cut-off value of ≤30 min [28]; (d) SE log rank, *p* = 0.224: based on normative cut-off value of ≥85% [28]; (e) WASO (log rank, *p* = 0.549: based on normative cut-off value of >30 min [28]; and (f) NA (log rank, *p* = 0.972: based on group median of 23 episodes).

#### *5.3. Multivariate Predictors of Survival: Machine Learning Results*

In the machine learning dataset, 46 participants had died within the specified time period of follow-up (i.e., by 14 May 2020). The Lasso model selected 22 predictor variables, with 14 variables consistently selected in all 50 participants during the process of validation (Figure 2). These involved eight predictor variables associated with greater survival time and six predictor variables, associated with a reduced survival time. The predictor variables associated with increased survival time, i.e., smaller hazard (in order of the coefficient associated with the predictor variable) were: later sleep diary time of final awakening, later actigraphy get up time, longer PiPS-B clinician's estimate of survival, better PSQI subjective sleep quality, greater PiPS-B global health status score (indicating better health), better actigraphy sleep efficiency, and higher haemoglobin values. The variables associated with reduced survival time were more frequent PSQI sleep disturbance wake middle of the night/early morning, higher neutrophil count, higher serum urea, serum creatinine, and serum C-reactive protein. On the contrary, a larger MSAS-SF total symptom distress was associated with a lower risk of death and a higher I < O was associated with a worse prognosis. The predicted median hazard was 0.00052, and the model was able to successfully differentiate between participants with a shorter/longer overall survival (log rank *p* < 0.0001) (Figure 3). Figure A1 shows the correlation between the actual survival and predicted hazard (Pearson's correlation coefficient *r* = −0.5; *p* = 0.0002).

The ridge model consistently identified 28 predictor variables in all 50 participants (Figure 4). During the process of validation, the top 10 variables consistently selected involved seven predictors associated with longer survival time and three predictors associated with shorter survival time. The seven predictor variables associated with longer survival time (in order of the coefficient associated with the predictor variable) were: actigraphy get-up time, sleep diary time of final awakening, sleep diary get-up time and PSQI usual get-up time; PiPS-B clinician's estimate of survival, PSQI subjective sleep quality and PiPS-B global health status score. The 3 predictor variables associated with shorter survival time were: use of opioid analgesia, modified Glasgow Prognostic Score and physicianassessed ECOG-PS day 8. The predicted median hazard was 0.44; however, there was no significant difference in overall survival when a median split was applied (log rank, *p* = 0.0914) (Figure 5). Figure A2 shows the correlation between the actual survival and predicted hazard (Pearson's correlation coefficient *r* = −0.5; *p* = 0.0002).

**Figure 2.** The mean cross—validated error (CVM) of predictor variables for hazard selected by the Lasso model.

**Figure 3.** Kaplan–Meier curve comparing survival probability predicted by the Lasso-derived algorithm (log rank, *p* < 0.0001).

**Figure 4.** The mean cross—validated error (CVM) of predictor variables for hazard selected by the ridge model.

**Figure 5.** Kaplan–Meier curve comparing survival probability predicted by the ridge-derived algorithm (log rank, *p* = 0.0914).

The elastic net model selected 10 predictor variables, with 6 variables being consistently selected during the process of validation: the two consistently selected predictor variables associated with longer survival time were (in order of the coefficient associated with the predictor variable): later actigraphy get-up time and greater PiPS-B global health status score; the 4 consistently selected predictor variables associated with shorter survival time were: higher serum urea, neutrophil count, serum C-reactive protein, and serum creatinine (Figure A3). The predicted median hazard was 0.408, but there was no significant difference in overall survival (log rank, *p* = 0.9877) (Figure A4). Figure A5 shows

the correlation between the actual survival and predicted hazard (Pearson's correlation coefficient *r* = −0.08; *p* = 0.5808).

#### **6. Discussion**

The results of this study show that univariate approaches to survival prediction, based on, for example, the I < O, are not very powerful; whereas, multivariate approaches appear to hold promise. To the best of our knowledge, this is the first study describing the application of supervised machine learning methods, involving a combination of actigraphy-derived rest-activity and sleep parameters, and data collected in routine clinical practice (i.e., simple questionnaires such as the MSAS-SF, ECOG-PS, PSQI, venous blood sampling) to prognosticate patients with advanced cancer, receiving supportive and palliative care [37]. Our study confirmed certain established novel predictors and identified some for survival in this group of patients and points to the importance of sleep characteristics for prognostication. The results of the study also confirm that clinicians are inaccurate prognosticators [3], since 11 (24%) participants were still alive at 1 year (despite the inclusion criteria of clinician estimated prognosis of more than 2 weeks but less than 1 year).

The literature had suggested that actigraphy-derived parameters, and the I < O index in particular, could be used as predictors because a low I < O is associated with increased morbidity (worse symptoms, worse quality of life), and with decreased survival [7]. At the outset of this study, we therefore focused on the I < O and other parameters describing the robustness of the rest-activity. We indeed observed a very high prevalence (i.e., 95%) of disrupted rest-activity rhythms in these advanced cancer patients, which is much higher than the reported prevalence of 19.1–54.9% [7]. This disparity undoubtedly reflects different populations, with our population having more advanced disease (and worse performance status) than previous studies [11,38]. However, in the univariate analyses of the data in our study there was no direct association between I < O and survival. Furthermore, other actigraphy-derived parameters, when used in isolation, are also not very accurate in the population.

However, the results of the study suggest that novel models developed through machine learning can facilitate improvements in prognostication. Penalised regression methods implement a feature selection strategy, providing a combination of subjective and objective predictor variables of survival that are ranked based on their contribution to the model. The models manage collinearity within the dataset, which is particularly useful in datasets involving terminal cancer patients, where often the number of features exceeds the relative sample size. The best performing method was Lasso regression which reduces the coefficients of variables with a minor contribution to zero and thereby creates a simple 'model' with only a few variables. Sleep parameters were amongst the most important variables, not only in the Lasso model but also in the more complex elastic net and ridge models. These measures primarily represented positive predictors of survival. Sleep diary final awakening (lasso and ridge) and actigraphy-derived GUT (all models) were found to have particular prognostic relevance in our study, suggesting that a later sleep diary determined 'time of final awakening' and a later actigraphy-derived 'get-up time' are associated with a lower risk of death and improved survival. Furthermore, actigraphyderived SE, which may be considered an objective measure of sleep quality, was selected as a positive predictor of survival in the lasso model (i.e., greater sleep efficiency was associated with enhanced survival) for our population. Whilst actigraphy-derived sleep quality, as opposed to sleep quantity, has been reported to have prognostic significance in advanced breast cancer patients [13], we identified quantitative sleep measures as important contributors to survival prediction.

Studies have reported actigraphy-derived circadian disruption [10,12,39] and fragmented sleep [13] to have prognostic implications in cancer patients, yet little is known about the prognostic impact of subjective sleep measures. A recent study identified the PSQI sleep duration component as a prognostic indicator in a cohort of advanced hepatobiliary/pancreatic cancer patients [40], yet a novel finding in our study was the selection of other sleep parameters from the PSQI: (1) usual get up time and (2) subjective sleep quality, where a later get up time and very good sleep quality are associated with longer overall survival, and (3) PSQI sleep disturbance components—pain, cannot breathe comfortably and wake up in the middle of the night or early morning—were associated with poorer survival. Furthermore, subjective sleep parameters, as opposed to actigraphy-derived sleep parameters, were more commonly identified in all participants in the ridge model.

Venous blood sample measurements were also significant contributors to predicting survival in our study. Previous studies have reported moderate evidence for the prognostic significance of an elevated C-reactive protein (CRP) and leucocytosis being associated with a shorter survival [1,5,18,41]. Whilst our study was able to echo these findings, we were able to further identify novel biomarkers, such as an elevated urea and serum creatinine, that may also be associated with a poorer survival, and raised haemoglobin that may be associated with a lower risk of death. Blood sampling is generally deemed 'inappropriate' when patients are in their last days/weeks of life [42], regardless only one of the 94 patients screened for our study, declined participation. Our findings endorse further evaluation of biological parameters from venous blood sample data, as they may be beneficial to improving prognostication in these patients.

Although our multivariate findings controversially imply that a higher I < O is associated with a shorter predicted survival time, all participants in our population had poor health, i.e., an I < O of <99%, which has recently been identified as an optimal cut-off for distinguishing between healthy controls and patients with advanced cancer [43]. Further inspection of our data identified that all our participants, whether they had shorter or longer survival had disrupted rest-activity rhythms, equally both groups had moderate symptom distress as measured by TMSAS, inevitably expected in an advanced cancer population. Therefore, whilst it may be a simple way of quantifying rest-activity rhythms, I<O may be a more meaningful prognostic indicator during the earlier trajectories of cancer, as opposed to the progressive stages.