{'fraction_non_alphanumeric': 0.029250696835546812, 'fraction_numerical': 0.020987047056894574, 'mean_word_length': 4.742044812652984, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 5, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistent Homology of Tumor CT Scans Predicts Survival In Lung Cancer

Radiomics, the objective study of non-visual features in clinical imaging, has been useful in informing decisions in clinical oncology. However, radiomics currently lacks the ability to characterize the overall structure of the data. This field may benefit by incorporating persistent homology, a popular new algorithm that analyzes whole data structure. We hypothesized that persistent homology could be applied to lung tumor scans and predict clinical variables. We obtained computed tomography lung scans (n = 565) from the NSCLC-Radiomics and NSCLC-Radiogenomics datasets in The Cancer Imaging Archive. Segmentation data was used to create a cubical region centered on the primary tumor in each scan. For each scan, a cubical complex filtration based on Hounsfield units was generated. We created a feature curve that plotted the number of 0 dimensional topological features against each Hounsfield unit. The curve's first moment of the distribution was utilized as a summary statistic to predict survival in a Cox proportional hazards model. The first moment of the distribution is equivalent to the area under the curve of our topological feature curves (AUC). The Kruskal-Wallis H Test and a post-hoc Dunn's test with Bonferroni correction were used to test AUC differences among survival quartiles. After controlling for tumor image size, age, and stage, AUC, was associated with poorer survival (HR = 1.118; 95% CI = 1.026-1.218; p = 0.01). AUC was significantly higher for patients in the lowest survival quartile compared to the highest survival quartile (p < 0.001). We have shown that persistent homology can generate useful clinical correlates from tumor CT scans. Our 0-dimensional topological feature curve statistic predicts survival in lung cancer patients. This novel statistic may be used in tandem with standard radiomics variables to better inform clinical oncology decisions.


Introduction
The use of radiomics in tumor imaging has been an invaluable, noninvasive tool informing cancer diagnosis, treatment, and prognosis 1 . Geometric properties of tumors such as size, surface area, and volume are commonly studied features among more non-visual features such as texture 2 . We propose extending the utility of radiomics by studying topological properties. In contrast to the local structural focus of geometry, topology focuses more on global structure.
Since topological properties would focus on general tumor shape properties, they may be more robust to noise and capture information not found in traditional radiomic features. Intuitively, one can describe several topological differences between malignant and benign tumors. Benign tumors are well-connected and have a homogeneous shape whereas malignant tumors are more likely to have diffuse spread and necrotic cavities 3 . Quantifying the number of topological features may be useful in predicting patient survival. In fact, pathology already uses shape properties in the context of Gleason scores, which is a measure of prostate cancer severity based on prostate gland shape 4 . Persistent homology is a popular technique in the "-omics" sciences to describe the topological features of large data sets. Persistent homology has already been used in cancer biology from genomics to histology 5, 6 . A statistical measure inspired by persistent homology called the smooth Euler characteristic has been developed to predict clinical outcomes from glioblastoma tumor imaging 7 . Persistence images, an alternative representation of persistent homology, have been used in machine learning models to classify MRI images of hepatic tumors 8 . While radiomics and persistent homology have largely existed in separate worlds, they share similar challenges. In radiomics, one challenge is to find the most informative image features for analysis 9 . In persistent homology, a similar challenge is finding the most useful topological data representation for visualization, statistical comparison, and predictive modeling. We were interested in whether image features related to topology in lung CT scans were associated with survival. Since we hypothesized earlier that an increased number of topological features may correlate to a more malignant tumor, we wanted to create a topological data representation that captured the quantity of topological features of a CT tumor scan. In this project, we create a new radiomics variable using the statistical summary variables of the 0D topological feature curve, our way of representing topological feature quantity. We describe this further in section 2. To our knowledge, no work has been done in using persistent homology to characterize survival in cancer patients using such a summary statistic. We show through a discrete and continuous analysis that moment 1 of our 0D topological feature curve (i.e. area under the curve) is significantly associated with worse survival.

Persistent Homology Background
Topological data analysis (TDA) encompasses a broad set of techniques used to highlight topological patterns in data. Computing persistent homology using cubical complexes is a popular method within TDA for describing topology in imaging data 10 . Counting topological features on an image requires a binary black and white image, so we cannot count the topological features from the grayscale CT scan directly. However, we can convert the CT scan into a series of binary black and white images (i.e. filtration) from which the topological features can be counted. In the context of persistent homology, cubical complexes are essentially synonymous with these binary images.
For each CT scan, we create a series of binary images (i.e. cubical complexes), one for each Hounsfield unit filtration value. For example, if the filtration value is -900 Hounsfield units, then all pixels with units less than or equal to -900 are colored black, and all pixels above are colored white. Topological features can be counted from these binary images. A particular topological feature can be described by dimension and the range of Hounsfield unit filtration values across which the feature is found.
We show this filtration process with slice 13 of scan 1 in figure 1A. This slice is shown as it demonstrates a potential connection between tissue characteristics (calcifications) and topological features. At low Hounsfield unit thresholds, the image is mostly white. Islands of black pixels are considered connected components or 0D features. Tissue necrosis may result in a more fragmented tumor, which would increase the number of 0D features. "Lakes" of white pixels are considered holes or 1D features. Two-dimensional features are not shown in figure 1 as they only appear when considering the cubical complexes of all the slices together as a three-dimensional structure. Intuitively, they can be thought of as pockets of white pixels. Calcifications would likely manifest as 1D or 2D features late in the filtration, which is shown when the Hounsfield unit threshold is 200 in figure 1A.
Topological barcodes are one of the most popular ways to visualize persistent homology 11 . Each topological feature in a barcode diagram is given a color to represent its dimension and a horizontal bar that spans the Hounsfield unit filtration values where the feature can be found. The barcode diagram in figure 1B represents the entire cubical complex persistent homology of scan 1. Since we were interested in topological feature quantity, we developed an alternative representation of persistent homology called the topological feature curve.
The topological feature curve is a transformation of the barcode diagram that counts the number of bars (i.e. topological features) at each Hounsfield filtration value. The number of bars is then plotted against each normalized Hounsfield unit filtration value. We can construct four topological feature curves, one for each feature dimension, plus a fourth curve that counts all features regardless of dimension. The topological feature curves in Figure 1C represent all feature curves in scan 1.
We use the moments of the distribution of the resulting 0D topological feature curve as predictor variables for our survival analysis. We used raw moments. The exact mathematical calculation is described in Table 3 in the Supplement.
Similar topological summary statistics have been described in analyzing brain connectomes of individuals with ADHD using 0D Vietoris-Rips complex persistent homology features 12 . Likewise, we also focus our analysis on 0D features; however, we generate our topological features using cubical complexes. Table 1 gives a descriptive overview of the two lung scan data sets. These data sets were obtained from The Cancer Imaging Archive 13 . The NSCLC-Radiomics set (n = 421) had more patients compared to the NSCLC Radiogenomics set (n = 138) 14,15 .

Results
The NSCLC Radiogenomics had a greater proportion of patients who died during its study. There was no significant difference in age between the Radiomics cohort (68.05 ± 10.09) and the Radiogenomics cohort (69.27 ± 8.79, t = 1.27, p = 0.206). Stage had a different distribution between the two cohorts (Stage I tumors 22.1% vs 63.0%, χ 2 = 135, p < .001). There were no significant differences in sex proportion between the two data sets (χ 2 = 1.38, p = 0.239). Males were most prevalent in both sets. Moment 1 of the 0D feature curve (Z = 28592, p = 0.782) and tumor image size (Z = 30544, p = 0.36) were not significantly different between the two data sets. Both of these variables were significant in predicting survival in the univariate Cox hazard model. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2020.  Table 5 in the supplement). The moments of the 0D feature curves are presented as median [IQR] since the data was not normally distributed. Tumor image size was also not normally distributed and also presented as median [IQR]. All other data are presented as mean (SD) or as a proportion. Statistical comparisons between the two cohorts were performed using Wilcoxon rank-sum test for the non normally distributed data. Student's t-test was used for the other continuous data. Chi-squared test was used for the categorical stage and sex data.

Discrete Analysis of 0D Topological Feature Curves
We first performed a discrete analysis comparing survival groups by their moments of distribution of the 0D feature curves. Patients were placed in survival groups based on their survival quartile in the combined study cohort. This discrete analysis was performed to explore possible visual differences in how our new 0D topological feature curve metric interacted with survival. Figure 6 in the supplement shows that poorer survival quartiles had taller 0D feature curves. Figure 7 in the supplement quantifies the difference in these curves showing that poorer survival quartiles had larger moments of the distribution. Fully detailed statistical comparisons among the moments of distribution values are found in table 4 in the supplement.

Cox Proportional Hazard Analysis of 0D Topological Feature Curves
The discrete analysis was performed primarily for intuition and data visualization. The discrete analysis is not as rigorous as a Cox model, which would measure the impact of the moments of the 0D topological feature curve on survival outcomes. We verified that the assumptions of proportional hazards were met through visual appraisal of the Schoenfeld residuals plot, which is shown in figure 9 in the supplement. . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020.  Table 2. Both univariate and multivariate Cox hazard models show moment 1 of the 0D feature curve is associated with poorer survival. Data is shown as hazard ratio (95% confidence interval). Bolded predictor variables were significant in both the univariate and multivariate models. Italicized predictor variables were only significant in the univariate model. The moments of distribution and tumor image size were rescaled to lie between 0 and 50 units. Linear scaling does not alter HR significance but does change the magnitude of the HR and CI to a more interpretable value. Patients with stage 0 tumors were removed from this model as the Stage 0 vs stage I variable did not converge in the model (table 5 in supplement).

4/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ; Figure 1. Methodology to generate cubical complexes. Slice 13 of CT scan 1 from the Radiomics dataset is shown in panel A. Each pixel value in a CT scan is described by a Hounsfield unit, which typically range from -1024 to 3071. We create a cubical complex by selecting a Hounsfield value as a filter. Any pixel below or equal to this filter is colored black. Any pixel above is colored white. This binary image is considered a cubical complex. Topological features can then be counted from each binary image. An example of a 0D feature is indicated by the red asterisk. An example of a 1D feature is indicated by the green asterisk. Two-dimensional features are not shown as they only appear when considering the cubical complexes of all the slices together as a three-dimensional structure. Intuitively, they can be thought of as pockets of white pixels. When the filtration reaches the maximum value, the whole image is colored black ending the filtration process. A barcode diagram (panel B) represents persistent homology by giving each topological feature a barcode. Color represents dimension, and the range represents the Hounsfield filtration unit range during which the feature existed. The Hounsfield unit horizontal axis in the barcode diagram has been normalized zero to one. The colored vertical lines represent the normalized Hounsfield units from the filtration values shown in the binary images. We create three topological feature curves (panel C) by summing the number of topological features by dimension at each Hounsfield unit filtration value. A fourth topological feature curve is generated by summing the total number of topological features regardless of dimension. The 0D topological feature curve is highlighted as we specifically use the moments of distribution of this curve for our survival analysis.

5/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ; Increasing age had a significant effect on survival (multivariate HR: 1.025; 95% CI: 1.013-1.037; z = 4.109; p < 0.001). Sex did not have a significant effect on survival (multivariate HR: 1.128; 95% CI: 0.898-1.417; z = 1.037; p = 0.30). Increasing stage also had a significant effect on survival except for stage IV, which is likely due to small sample size (n = 4). Moment 1 of the 0D topological feature curves, which also represents the AUC, had a significant effect on survival (multivariate HR: 1.118; 95% CI: 1.026-1.218; z = 2.547; p = 0.011). Moment 1 of the 0D topological feature curve significantly predicted survival even after controlling for tumor image size. Tumor image size only had a significant effect on survival in the univariate model (univariate HR: 1.015; 95% CI: 1.003-1.026; z = 2.565; p = 0.01). Figure 2 shows a forest plot of the Cox model. The nonsignificant Stage IV vs Stage I predictor variable was removed as its wide confidence interval interfered with data visualization.

Discussion
As the "-omics" sciences continue to become more present in the clinic, we need stronger analysis techniques that make meaning out of the massive amount of data. Compared to the other "-omics," radiomics may be the most useful from a clinical perspective. Obtaining repeated tumor CT and MRI scans is easier for a hospital and less invasive for a patient compared to repeated tumor sequencing (genomics) or mass spectrometry (proteomics). Therefore, it would be quite valuable if persistent homology could uncover additional meaning from imaging data. We show that even after controlling for tumor stage, age, sex, and image size, the first moment of the 0D topological feature curve is a significant predictor of survival. This effect was significant across two independent data sets indicating a resilience of our methodology to batch effect. We considered tumor image size as a confounding possibility since the number of topological features would be expected to increase with image size. However, the effect of moment 1 of the 0D feature curve on survival remained significant in the multivariate model, and about 6/18 . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ;

7/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ; 30% of the variance in moment 1 of the 0D topological feature curve was not explainable by tumor image size alone (figure 10 in the supplement).
Despite these promising results, there remain limitations in our study. Though the pipeline worked well across two independent batches, the imaging parameters across the two sets were similar. It remains unknown how sensitive this technique is to imaging parameters such as CT machine model, which has been shown to affect calculation of standard radiomic features 16 . While we controlled for some important clinical variables as potential confounders, there are other relevant clinical variables such as EGFR gene status we were unable to test due to lack of data. We could only control for relevant clinical variables common to both data sets. In addition, while we know the analyzed scans were for planning radiotherapy or surgery, all of the patients' prior clinical history is a black box, which is a common limitation of public databases 17,18 .
We only explored 0D topological feature curves in our analysis. We had found similar trends using other topological features curves in predicting survival. However, adding the moments of distribution of the other three feature curves would add up to 12 additional variables to consider, which we do not believe our study is adequately powered to assess. We chose to limit our focus to 0D features as they had the strongest association with survival.
Most of these issues can be resolved by extending our pipeline to study additional CT scans from NSCLC patients with known clinical outcome across different databases. This would allow us to control for other clinical confounders and provide greater statistical power to assess moments of distribution of higher dimensional topological feature curves.
After initial therapy, patients with lung cancer obtain follow-up CT scans to survey tumor recurrence. Local and regional failure occur respectively when the primary tumor site and nearby lesions regrow after therapy. As many of the curative intent therapies include primary-tumor directed therapy (radiation or surgery), it would be worthwhile in the future to also consider local control or local progression free survival in addition to overall survival. Long term, we envision predictive modeling that incorporates relevant clinical and radiomics variables to personalize anticancer therapy for each patient. The topological feature curve summary statistics may be appropriate variables to include in such models. Since our variable only requires a CT scan and appropriate segmentation, this tool would provide universal benefit to all healthcare practice types.

Methods
A publicly available pipeline using R (v3.6.1) and Python (v3.7.6) code was built for this project [19][20][21] . All images and segmentation objects were obtained from The Cancer Imaging Archive 13 . Within TCIA, the NSCLC-Radiomics (n = 422) and NSCLC Radiogenomics (n = 211) cohorts provided CT scans and clinical data 14,15,17,18 . Of the 633 total scans, 565 also had segmentation data, which was necessary for our pipeline. Figure 4 shows the exclusion and inclusion criteria for each analysis.
Currently, R has a limitation in reading .dcm segmentation files. To work around this issue, we first converted the .dcm segmentation files into the NIfTI file format using the dcmqi open source Bash library 22 .
In Python, we used the built in glob and pathlib libraries to recursively read in the file objects produced from the R code; numpy v1.18.5 to create data structures that represented the processed tumor scans numpy objects from R 47 ; gudhi v3.0.0 to compute the persistent homology using cubical complexes 48 ; and csv v1.0 library to output the computed persistent homology as csv files to be processed again in R.
We used ITK-SNAP v3.8.0 as the DICOM viewer to visually ensure our cubical segmentation properly highlighted the tumor region of interest 49 . Figure 5 shows an overview of the data pipeline. Within R, a cubical region that encapsulated the tumor was delineated using the NIfTI segmentation files. This data object was exported to Python to compute cubical complex persistence homology (functionality in R not present at time of analysis). Cubical complex persistence homology was computed in Python and exported as .csv files to be processed in R. Within R, the raw persistence homology was transformed into topological feature curves, which plot the topological features against each Hounsfield filtration value. The first four moments of distribution were calculated for specifically 0D topological feature curves. These moments were used as predictor variables for our survival analysis.
Unpaired t-test statistics (for normally distributed data) and Wilcoxon rank sum tests (for non-normally distributed data) were used to descriptively compare patient characteristics between the Radiomics and Radiogenomics data sets. Chi-squared test was used to compare categorical variables. In the discrete analysis, all moments of distribution had a skewed exponential distribution, so non-parametric statistical tests were performed. Kruskall-Wallis H test was used to compare the median

8/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ; Figure 4. Flowchart of tumor scan inclusion and exclusion. 565 tumors had both segmentation data and CT scans allowing for computation of cubical complexes in our pipeline. *The univariate Cox model used all 565 tumor data assuming no NAs were present. However, 22 patients were missing age information, and 1 patient was missing stage information. These patients were excluded in their respective univariate Cox analysis. **The multivariate Cox model excluded 6 additional patients who were stage 0 since the Cox models did not converge (table 5 in the supplement).

9/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted December 7, 2020. ; https://doi.org/10.1101/2020.12.06.20244863 doi: medRxiv preprint Figure 5. Overview of data pipeline. Using the tumor segmentation data, a cubical region that encapsulated the tumor was isolated from each scan. This data object was created in R and passed on to Python to compute persistent homology. Persistent homology was computed using cubical complexes. The pink arrow in the third panel points to an example of a 1D feature. The full persistent homology can be described in a bar code diagram where the horizontal axis represents the Hounsfield filtration value. Each bar in the diagram represents a single topological feature, and its range represents persistence. We transform this object into topological feature curves, which plot topological features against the normalized Hounsfield unit filtration values. Finally, we calculated the moments of distributions for 0D feature curves to be predictor variables in our survival analysis. moments values of all survival groups, and a post hoc Dunn's test with Bonferroni correction was used to make pairwise comparisons.
For the continuous analysis, univariate and multivariate Cox proportional hazard modeling were performed for the moments, tumor image size, cancer stage, age, and sex. We included tumor image size as a covariate in the multivariable Cox Hazard model to control for potential confounding with our 0D topological feature curve. Increasing the number pixels of an image may increase the topological features without necessarily changing the overall topological structure of the tumor image. Tumor image size was calculated as the product of number of pixels per slice and number of slices. All survival times were directly provided in the NSCLC-Radiomics data set. This survival time was not directly provided in the NSCLC-Radiogenomics data set, but the date of imaging and date of death or last follow-up were provided. Survival times in the NSCLC-Radiogenomics were calculated by subtracting the date of death/last follow-up from date of imaging. The stage 0 predictor variable did not converge since none of the 6 patients with stage 0 cancer died during the study, so these patients were excluded (table 5 in the  supplement).
Log-rank test was used to compare all four moment 1 quartiles in the Kaplan Meier curves, and a post hoc log-rank test was used to make pairwise comparisons. All post hoc statistical tests p-values received Bonferroni correction. Summary statistics are stated as mean ± SD or median [IQR]. All code to reproduce results and figure is publicly available on https://github.com/eashwarsoma/TDA-Lung-Phom-Reproducible with detailed instructions. Figure 6. Median 0D feature curves are larger for patients in poorer survival quartiles. Each survival group had approximately one quarter of the combined study cohorts (n = 145, n = 138, n = 141, n = 141, from worst to best survival). Within each survival group, the median 0D feature count value for each Hounsfield units was collected into one vector and plotted. The Hounsfield units were normalized to 0 to 1.

Moments of Distribution Formulas
Moment Formula . CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

15/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2020. ; https://doi.org/10.1101/2020.12.06.20244863 doi: medRxiv preprint Figure 8. Groups were divided into quartiles based on survival. Each 0D feature curve was analyzed for the first four moments of distribution. The range of the moment values for each group is provided here. Color depth represents magnitude of moment value.

16/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2020. ; https://doi.org/10.1101/2020.12.06.20244863 doi: medRxiv preprint Figure 9. Schoenfeld residuals were visually appraised to ensure assumptions of proportionality were met. While moment 1's Schoenfeld residuals in the figure are significant, we believe it is a chance result as visually there does not seem to be a trend in the residuals over time. The global Schoenfeld residuals were also nonsignificant Figure 10. Correlation plot comparing tumor image size and moment 1 of 0D topological feature curves. Increasing the image size is expected to increase the topological features. Tumor image size was calculated by multiplying the number of pixels per CT scan slice by the number of slices. Both variables were linearly scaled to 0 to 50. We generated a correlation plot to check how much variance provided by moment 1 of the 0D feature curve was unexplained by tumor image size. About 30% of the variance in moment 1 is unexplained by tumor image size. Despite the large overlap, moment 1 remained significant even after controlling for tumor image size in our Cox model.

17/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2020. ; https://doi.org/10.1101/2020.12.06.20244863 doi: medRxiv preprint Author contributions statement E.S. conceived the project. E.S., A.L., and R.W. designed project plan. E.S. and A.L. wrote code for the project. E.S. wrote manuscript. J.S. supervised all previous steps. All authors reviewed the manuscript.

18/18
. CC-BY-NC 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted December 7, 2020. ; https://doi.org/10.1101/2020.12.06.20244863 doi: medRxiv preprint