Cohort-Individual Cooperative Learning for Multimodal Cancer Survival Analysis

Survival analysis can be expressed as an estimation of the hazard function, which models the patient’s probability of death at a certain time, conditioned on personal clinical records. In an early stage, Cox’s proportional hazards regression model [31] conceptualises the hazard function as a multiplication of two components: 1) the underlying baseline hazard function, describing how the risk of event per time unit changes over time at baseline levels of covariates; and 2) the hazard ratio, measuring the impact of covariates. Since the baseline hazard function for a given cancer is constant, the hazard ratio determines whether a patient is at high risk. Based on the Cox’s regression model, subsequent survival analysis methods [3, 4, 5, 7, 6, 8] aim to predict the personalized hazard ratio using quantitative data extracted from short-term clinical indicators or long-term follow-up reports. For example, Kappen et al. [6] summarized seventeen pretreatment characteristics, such as residual tumor size, age, and thrombocytes, to predict the treatment outcome using a neural network. Ohno-Machado et al. [8] predicted the survival of AIDS patients conditioned on demographics, laboratory markers, and clinical findings, such as age, hemoglobin, and albumin.

Genomic profiles provide molecular information of tumors, which is important for cancer prognosis prediction. For example, certain genetic mutations or variations can affect tumor growth, metastasis, and response to chemotherapy or targeted therapies. Extensive studies [9, 10, 12, 11] have been conducted to predict cancer prognosis and treatment response by leveraging genomics data. Recently, Yousefi et al. [12] combined deep learning and Bayesian optimization methods to tackle high-dimensional cancer outcomes prediction tasks. Qiu et al. [11] built good predictive models with limited high-dimensional samples by using a meta-learning survival analysis framework.

Pathological images provide morphological features of tumors, which can provide valuable information about the aggressiveness of the tumor, its response to treatment, and the likelihood of disease recurrence. Some recent works [13, 14, 15, 16] constructed effective survival analysis models based on giga-pixel Whole Slide Images (WSIs). For example, Yao et al. [15] introduced the siamese MI-FCN and attention-based MIL pooling to efficiently learn imaging features from the WSI and then aggregate WSI-level information to patient-level. Zhu et al. [16] adaptively sampled and grouped hundreds of patches from each WSI into different several clusters. Then they employed an aggregation model to make patient-level predictions based on cluster-level survival prediction results.

Despite impressive performance achieved on survival analysis datasets, they have a narrow perspective and potentially overlook important aspects or correlations present in other modalities. On the contrary, multimodal models can offer improved performance, robustness, and a more comprehensive understanding of cancer prognosis.

Integrating genomics and pathology images can enhance the predictive power of survival analysis models [18, 19, 20, 32, 25, 22, 21, 24, 23] and thus attract increasing attention recently. By utilizing the cellular graph embedded in the tissue, Nakhli et al. [22] utilized a unified representation for each patient, leveraging the hierarchical organization of the tissue. Chen et al. [21] captured the interactions between features across multiple modalities by combining the unimodal feature representations using the Kronecker product. Additionally, they incorporated a gating-based attention mechanism to control the expressive power of each representation. Furthermore, Chen et al. [25] proposed an interpretable, dense co-attention mapping between WSIs and genomic features formulated in the embedding space. Xu et al. [28] introduced the optimal transport theory to match WSI patches and gene embeddings for selecting informative patches to represent gigapixel WSIs and build an interpretable co-attention module to effectively fuse multimodal data. Moreover, Zhou et al. [29] found that the generated cross-modal representation can enhance and recalibrate intra-modal representation, and thus significantly improve its discrimination for survival analysis.

Existing co-attention-based methods focus on extracting common knowledge by using multimodal interactions. Moreover, they are exposed to the risk of overfitting high-dimensional data, leading to performance degradation on unseen samples. To address these issues, multimodal knowledge is comprehensively and explicitly decomposed in our framework to facilitate effective multimodal fusion. Meanwhile, our framework leverages cohort guidance to improve the generalization ability of the decomposed knowledge components.

Genomics. Genomic profiles can identify specific genetic alterations or biomarkers associated with cancer prognosis. For example, certain genetic mutations, gene expression patterns, or alterations in DNA copy number can serve as prognostic markers, helping to predict the likelihood of survival conditions of patients. In our framework, we partition RNA sequencing (RNA-seq), Copy Number Variation (CNV), Simple Nucleotide Variation (SNV) sequences into six sub-sequences as previous methods [29, 28]. Each sub-sequence is transformed into a feature by a Self-normalizing Neural Network (SNN) [17], which is a more powerful alternative to conventional Multi-Layer Perceptron (MLP). By integrating all sub-sequence features, we obtain the genomic representation $F_{g}\in\mathbb{R}^{1\times 256}$.

Pathology. Whole Slide Images (WSIs), i.e., pathology images, describe the information about the tumor immune microenvironment and provide valuable information for cancer prognosis prediction. Considering the huge resolution of WSI, which exceeds the capacity of Convolutional Neural Networks (CNNs), we adopt a strategy of splitting tissue regions within each WSI into non-overlapping patches at 20x magnification (256 x 256 resolution for each patch). Following previous works [28, 29, 21], we utilize an ImageNet pre-trained ResNet-50 [2] model to extract a 1024-dimensional embedding for each patch, while all patch embeddings of the same WSI are collected as an embedding set. It is worth noting that such an embedding set is still high-dimensional as it typically contains tens of thousands of patches for each WSI. To further reduce information redundancy, we employ the K-means algorithm [33] to cluster all patch embeddings into $k$ groups and leverage the cluster centers as pathology features. However, due to the stochastic nature of K-means, cluster centers of different samples may be misaligned, which means that two cluster centers from the same ordinal position of two samples may exhibit completely different phenotypes. Consequently, deep models prefer to learn phenotype-independent knowledge instead of specialized knowledge tailored to each phenotype, overlooking crucial information in important patches. In our framework, we address this misalignment issue by cluster center alignment (CCA), an optimal matching between an anchor and cluster centers, as shown in Fig. 3. This involves assigning each cluster center to a feature in the anchor with maximized similarities of matched pairs. To solve this, the Hungarian algorithm [34] is employed to calculate the permutation matrix that maps the current centers to their matched ordinal positions. After that, we obtained the aligned centers by multiplying cluster centers with the permutation matrix. Aligned centers are utilized to update the anchor with a ratio of $\tau$ during training and generated the pathology representation $F_{p}\in\mathbb{R}^{1\times 256}$, ensuring accurate alignment and capturing crucial information.

Genomic profiles and pathology images provide both unique and collaborative insights into tumors. Multimodal methods that can effectively integrate these two modalities are promised to be more impressive in learning discriminative representations than unimodal ones. To achieve this goal, decomposing knowledge into modular components enables scalability and flexibility in handling high-dimensional and heterogeneous multimodal data. Previous multimodal methods [35, 36] simply decompose multimodal knowledge into common and specific components. However, synergy-new knowledge that can only be generated through multimodal collaborations-is valuable yet overlooked in these methods. To better utilize all knowledge components in multimodal data, we propose a multimodal knowledge decomposition module and cohort guidance modeling to eliminate redundancy and strengthen valuable information under the guidance of patient cohorts.

Multimodal Knowledge Decomposition. In our framework, we decompose the knowledge within genomic profiles and pathology images into four distinct components, including redundancy, synergy and uniqueness of two modalities. Concretely, the large heterogeneity gap between pathology and genomics indicates a wealth of modality-specific knowledge. Meanwhile, there still exist considerable common knowledge between them. For example, particular genetic mutations or alterations can be observed at both the molecular level and the histopathological level simultaneously. Furthermore, another essential component, synergy, can provide new insights on tumors to enrich the knowledge in multimodal data, especially for heterogeneous modalities. Emerging from multimodal interactions, it surpasses common knowledge by forming a unique realm of insights unattainable through any individual modality alone. As an example, the diagnosis of glioma typically necessitates genetic markers (IDH mutation and 1p/19q co-deletion) to categorize gliomas into subtypes and histopathological images to determine the presence of microvascular proliferation and necrosis simultaneously. Overall, a comprehensive and explicit decomposition enhances integration, representation, synergy capture, and flexibility. It enables more effective modeling and utilization of multimodal data, leading to improved performance, insights, and decision-making in cancer prognosis.

To achieve a comprehensive decomposition of knowledge, we develop four encoders within our framework. These encoders are specifically employed to capture and model different components—redundancy, synergy, and uniqueness—present within the genomic profiles and pathology images. Concretely, we employ MLP layers as modality encoders $\Phi_{p}$ and $\Phi_{g}$ to extract specific knowledge $P$ and $G$, focusing on aspects distinct from common knowledge shared by both modalities. The formula can be written as:

Diverging from $\Phi_{p}$ and $\Phi_{g}$ that receive only a single input, common and synergistic encoders $\Phi_{c}$ and $\Phi_{s}$ employed in this context necessitate two inputs to build multimodal interactions from different perspectives. In these encoders, a co-attention block is utilized to bridge the interactions across modalities and produce modality attentions to integrate knowledge from different modalities into a fused feature $C$ or $S$, as shown in Fig 4. The detailed operations can be formulated as:

where $A$ is the co-attention matrix and $fc$ indicates a fully-connected layer. Synergy $S$ is computed the same as common knowledge $C$, but with different parameters. Using the above encoders, we decompose modality representations $F_{g}$ and $F_{p}$ into genomic-specific $G$, pathology-specific $P$, common $C$, and synergistic $S$ components, respectively.

Cohort Guidance Modeling. In the above section, we employ different encoders to capture different knowledge components from multimodal data. Beside that, an effective loss function is also significant in enhancing the discrimination and generalization abilities of our model. Due to the high dimensionality of input modalities, multimodal interactions learned from each patient, as done in existing methods [25, 28, 29], may overfit task-irrelevant information, leading to reduced generalization and discrimination abilities. Therefore, we propose to utilize the correlations between patient cohorts to acquire more general knowledge components with cohort consistency, capture the heterogeneity of multimodal data, and gain a better understanding of cancer.

In our framework, we harness the cohort guidance to learn knowledge components from both knowledge and patient levels, as shown in Fig. 5. At the knowledge level, the most prominent difference between the decomposed knowledge components is their correlation to the original modality representations. Specifically, synergy is unattainable through any individual modality alone, so it is beyond the knowledge of both modalities. Redundancy is the intersection of knowledge as it is shared by both modalities. In addition, uniqueness is located in the source modality and is distinct from another modality. Therefore, we employ a set of similarity constrains between them as:

where $cos$ computes the cosine similarity between two inputs and $|\cdot|$ is the absolute operation used to constraints two inputs orthogonal to each other. We freeze the gradients of $F_{p}$ and $F_{g}$ from $l_{k}$ to prevent potential model collapse that the extracted pathology and genomic features follow similar distributions. At the patient level, task-relevant information can be obtained by distinguishing patients conditioned on their risk scores, which are valuable for extracting more discriminative representations. Therefore, we split all patients into $r$ equal groups conditioned on their ground truth survival time, while patients in the same group are considered as similar. For uncensored patient, its features are closer to the features extracted from similar patients, while distinctive to patients in other groups. For censored patients, it may be closer to the features of patients at similar or lower risk and distinctly different from those at higher risk. Such an assumption in line with the definition of contrastive learning, which has been proved to have strong feature learning capability in both supervised and unsupervised settings. Therefore, we develop a cohort contrastive learning to enhance the discrimination and generalization abilities of our model. In the training phase, we build a cohort bank to store features of historical patients in each group. Our cohort bank is actually $r$ queues with length $b$ and follows the first-in-first-out principle for dynamically feature updating. For the feature of uncensored patient, we reduce its distances to similar patients in the cohort bank, while enlarging its distances to other patients. As an example, the formula for synergy $S$ is:

where $\mathcal{S}_{+}$ is the patient group of similar risk storing in the cohort bank, while $\mathcal{S}_{-}$ is the collection of other patients. $d$ is the similarity function. Noted that $S$ can be replaced by other knowledge components $G$, $P$, or $C$ similarly. For censored patient, we can extend $\mathcal{S}_{+}$ as the collection of patients with similar or lower risk.

Combining knowledge and patient levels, the cohort guidance loss can be formulated as:

Recently, Transformer [42] has shown impressive performance on multimodal learning, and thus is employed to integrate the decomposed components for survival predictions. Specifically, we concatenate a class token $U$ with the decomposed components $G$, $P$, $C$, and $S$ as the inputs of Transformer $\Phi_{t}$. Moreover, a fully-connected layer with Sigmoid activation, denoted as $\sigma$, is appended to the output of class token for predicting the hazard function $H$:

For survival prediction, following previous works [25, 26, 28, 29], we construct an one-hot time-series vector $T$ based on $n$ equal parts of ground truth survival time. Under this setting, the original event time regression problem is simplified to a classification problem. We find out the time interval $t_{k}$ in which the ground truth event occurred as the classification label $k$ for each patient. Each patient sample is defined as a triplet $\{H,c,k\}$, where $H=\{h_{1},...h_{n}\}$ is the predicted hazard vector measuring the probability of ground truth event time located in corresponding time interval $t_{k}$. Additionally, we define the discrete survival function $f_{sur}(H,k)=\prod_{j=1}^{k}(1-h_{j})$. Following previous works [28, 29, 25], we generalizes the negative log-likelihood (NLL) with censorship to supervise the survival prediction by:

Finally, the overall loss function of our framework is:

where $\alpha$ controls the effect of our cohort guidance.

Dataset. To validate the effectiveness of the proposed methods, we evaluate our framework on five datasets from TCGA: Bladder Urothelial Carcinoma (BLCA, n=372), Breast Invasive Carcinoma (BRCA, n=956), Glioblastoma & Lower Grade Glioma (GBMLGG, n=569), Lung Adenocarcinoma (LUAD, n=453) and Uterine Corpus Endometrial Carcinoma (UCEC, n=480). These datasets contain hundreds of paired genomics, pathological images, and follow-up data collected from multiple centers, enabling center-agnostic analysis on precision oncology. We collected all diagnostic WSIs used for primary diagnosis, resulting in 2,830 WSIs with an average of 15k patches per WSI at 20x magnification (assuming 256 x 256 patches). Genomic profiles provide molecular information of individual, including RNA sequencing (RNA-seq), Copy Number Variation (CNV), Simple Nucleotide Variation (SNV), DNA methylation. Following previous works [25, 29], we use RNA-seq, CNV, and SNV sequences and further group them into six genomic sub-sequences: 1) Tumor Suppression, 2) Oncogenesis, 3) Protein Kinases, 4) Cellular Differentiation, 5) Transcription, and 6) Cytokines and Growth.

Evaluation Metrics. The concordance index (C-index) is the proportion of all comparable patient pairs for which the predicted outcome is consistent with the actual outcome. In survival analysis, the predicted outcome is considered consistent with the actual outcome if the predicted survival time is longer for individuals who indeed have a longer survival time, compared to the predicted event time for others. Moreover, Kaplan-Meier and T-test analyses are employed to assess the significance of differences in survival predictions between high- and low-risk groups.

Implementation. For each dataset, we adopted five-fold cross-validation to evaluate our model and other compared methods. Our framework was implemented using PyTorch running on a single NVIDIA GTX 3090 GPU. For optimization, we employed the SGD optimizer with a learning rate of 1e-3 to train our framework for 30 epochs.

We compare our framework with different types of survival analysis models, including genomics-based (SNN [17]), pathology-based (AttMIL [37], CLAM [38], TransMIL [39], and DTFD [40]), and multimodal (MCAT [25], M3IF [41], GPDBN [23], Porpoise [27], HFBSurv [19], SurvPath [26], MOTCat [28], and CMTA [29]) models. Here we provide brief overviews of several representative multimodal competitors:

• •

  Porpoise [27] is a multimodal fusion model that aggregates WSI patch embeddings with the generated attention scores and processes genomic profiles using a SNN network. The extracted multimodal features are simply concatenated to produce survival predictions.

• •

  SurvPath [26] learns biological pathway tokens from transcriptomics that can encode specific cellular functions and fuse two modalities using a memory-efficient multimodal Transformer to model interactions between pathway and histology patch tokens.

• •

  MOTCat [28] contains a global structure consistency, in which optimal transport (OT) is applied to match WSI patches and gene embeddings for selecting informative patches. More importantly, OT-based co-attention provides a global awareness to effectively capture structural interactions for survival prediction.

• •

  CMTA [29] explores the intrinsic cross-modal correlations and transfer potential complementary information by enhancing modality-specific representations by integrating with multimodal representations.

In addition, we build several models for a comprehensive comparison, including: 1) SNNTrans: two SNN layers followed with a transformer; 2) MaxMIL/MeanMIL: a maximum/average operation to aggregate multiple instances; and 3) DualTrans: combining SNN and TransMIL for genomic profile and pathological images, respectively, while a Transformer is employed for multimodal fusion. We unify the settings of compared models to make a fair comparison and report the re-produced results in our own environment.

C-index Comparison. As shown in Tab. 1, the proposed framework achieves an average score of 72.60%, which surpasses existing methods by a significant margin. On each dataset, our model obtains 68.62%, 68.40%, 86.14%, 69.57%, and 70.26% C-index scores on BLCA, BRCA, GBMLGG, LUAD, and UCEC, which are state-of-the-art performance on most datasets, while the second best on BLCA dataset with a small margin to the best 69.10% reported by CMTA [29]. Compared to unimodal models, either genomics- or pathology-based, our framework achieves significantly improvements on all datasets. It indicates that our framework can effective integrate complementary information in genomic profiles and pathology images. Compared to multimodal models, our framework reports more robust and significant performance on survival analysis. For example, CMTA [29] performs well on the BLCA dataset but still inferior to our method on others.

Kaplan-Meier Analysis. The Kaplan-Meier (KM) analysis is a non-parametric statistic used to estimate the survival function from lifetime data. Specifically, we use the median survival time of the entire cohort to divide all patients into high- and low-risk groups, represented as red and blue curves, respectively. The p-value in KM analysis indicates the probability of observing the observed difference in survival rates between two groups, under the null hypothesis that there is no true difference between them. As shown in Fig. 6, the p-values achieved by our framework are significantly lower than 0.05 on all five datasets, which indicates a statistically significant discrimination between high- and low-risk groups. We also compare the p-values with several representative methods in Tab. 2. Our model achieves the lowest p-values on BLCA and BRCA, meanwhile is competitive to state-of-the-art methods on GBMLGG, LUAD, and UCEC. Our framework can generalize across different cancers to produce robust predictions, holding the potential to significantly advance both clinical practice and cancer research.

T-test Analysis. The T-test is a statistical analysis method used to determine if there is a significant difference between the means of two groups. In the T-test analysis, the t-value quantifies the size of the difference relative to the variability within two groups, while the p-value is the probability that the results from our results occurred by chance. As the results shown in 7, our framework can distinguish between different groups with good p-values and t-values, especially on BLCA, BRCA, and GBMLGG datasets (p-value <  0.05). However, LUAD and UCEC datasets are still challenging for our framework (p-value >  0.05).

Module Analysis. To demonstrate the effectiveness of our methods, we progressively append our modules to our baseline DualTrans, including Cluster Center Alignment (CCA), Multimodal Knowledge Discrimination (MKD), and Cohort Guidance Modeling (CGM). As shown in Tab. 3, the proposed methods effectively and consistently improve the performance of the baseline DualTrans model. CCA aligns the cluster centers after K-means, facilitating the extraction of crucial information from important patches. Moreover, MKD enables more effective multimodal fusion by comprehensively decomposing multimodal knowledge, reducing redundancy and reinforcing valuable information. Further, CGM leverages the cohort guidance to promote a general multimodal interaction learning. Collaborating these modules, our full framework achieves significant improvements over the baseline model.

Knowledge Component Visualization. In Fig. 8, we visualize the feature distributions after MKD using T-SNE [43], where the orange, green, blue and purple points indicate common, synergistic, genomic-specific, and pathology-specific knowledge, respectively. The distributions of these components differ from each other and are consistent with our intuition that common and synergistic features lie in the middle of two modality-specific features. Moreover, the features of each component are clustered as several groups, which demonstrates that patients of similar risk tend to be closer than other patients under our cohort guidance. For example, genomic-specific features are basically clustered as four groups, which is consistent with the group number used in our framework. It proves that our framework can extract discriminative representations to facilitate multimodal fusion.

Hyper-parameter Experiment. There are several hyper-parameters used in our framework, including $k$ in K-means, anchor updating ratio $\tau$, cohort bank length $b$, and the number of patient groups $r$. We provide experiments results on these hyper-parameters in Tab. 4. From these results, we conclude that $k=6$, $\tau=0.1$, $b=10$ and $r=4$ work best in our framework, and thus are employed as the main settings of our framework. Moreover, smaller values of $k$, $b$, and $r$ imply fewer features and computational loads, but the performance usually drops significantly. When exceeding adequate values, enlarging them may not offer additional improvements. For $\tau$, it controls the updating rate of anchor, which is used in cluster center alignment. After a period of training, the anchor becomes stable, so the effect of $\tau$ may not be significant.

Knowledge Attention Analysis. To further illustrate the differences between knowledge components, we visualize their attention on input data in Fig. 9. In sub-figures (a)-(c), we visualize the attention regions of redundancy, synergy, and pathology-specific knowledge on input pathology images, while the last sub-figure (d) shows attention scores on genomic profiles. For pathology images, the focusing regions vary from component to component, indicating the distinctiveness between them. For genomic profiles, the first sample focus more on synergy, whereas the second one has a stronger reliance on redundancy. It implies that two types of multimodal interactions are necessary for effective survival analysis.