Bayesian Multi-view Graph Convolutional Network (BMGCN) for Integrative Multi-Omics Analysis with Survival Outcomes and Zero-Inflated Data

Bayesian Multi-view Graph Convolutional Network (BMGCN) for Integrative Multi-Omics Analysis with Survival Outcomes and Zero-Inflated Data Akhilesh Kaushal Department of Pediatric Hematology, University of Arkansas for Medical Sciences [email protected] October 15, 2025 Abstract Modern biomedical research generates vast, multi-modal datasets (multi-omics) from the same patient cohorts, offering an unprecedented opportunity to understand complex diseases. However, integrating these heterogeneous data views to predict clinical outcomes like patient survival presents significant statistical challenges. These challenges include data heterogeneity, high dimensionality, inherent zero-inflation due to technical dropouts or biological absence, and the need to incorporate prior biological knowledge. We propose the Bayesian Multi-view Graph Convolutional Network (BMGCN), a deep generative framework designed to address these challenges. BMGCN factorizes the data into shared and view-specific latent representations, enabling both data integration and the identification of view-specific signals. It employs graph-convolutional encoders to integrate prior biological network knowledge, a zero-inflated likelihood to accurately model sparse omics data, and a spike-and-slab prior for Bayesian view selection to identify modalities most relevant to the outcome. Finally, a semi-parametric Cox proportional hazards module allows the model to handle right-censored survival data directly. We detail the full generative model, derive the variational inference objective, and outline a comprehensive validation strategy. BMGCN provides a powerful, interpretable, and flexible framework for integrative multi-omics analysis. 1 Introduction The ability to measure multiple molecular layers—such as gene expression (RNA-seq), DNA methylation, microRNA, and proteomics—from the same subjects has revolutionized systems biology and precision medicine. Each of these “omics” views provides a unique window into the biological state of a system. The central hypothesis of multi-omics integration is that a unified analysis of these views can reveal emergent biological insights and yield more accurate predictive models than any single view alone [Ritchie et al., 2015]. However, the integrative analysis of such data is fraught with challenges. First, the data are heterogeneous, with different views having unique statistical properties, dimensionalities, and levels of noise. Second, many omics datasets, particularly those from single-cell and sequencing-based assays, are characterized by a high degree of sparsity and zero-inflation, where an excess of zero values complicates modeling [Lambert, 1992]. These zeros can represent true biological absence or technical artifacts (e.g., “dropouts”), and distinguishing between them is crucial. Third, predicting clinical endpoints, such as patient survival, requires specialized statistical models that can properly handle right-censored data, where the event of interest is not observed for all subjects [Cox, 1972]. Finally, ignoring the rich web of prior biological knowledge, such as protein-protein interaction networks or gene regulatory pathways, is a missed opportunity to guide the model toward more robust and interpretable solutions. To address these multifaceted challenges, we introduce the Bayesian Multi-view Graph Convolutional Network (BMGCN). BMGCN is a deep generative framework that provides a principled solution for integrating multi-omics data with censored survival outcomes. The design of BMGCN is guided by the following core objectives: 1. Latent Abstraction: To learn a low-dimensional representation of the data by factorizing the signal into latent variables that are shared across all omics views (Zi) and variables that are specific to each view (Z(k) i). This captures both common and view-specific biological processes. 1 2. View Selection: To automatically determine the relevance of each omics view for predicting the survival outcome. We achieve this through a Bayesian spike-and-slab prior on the view-specific latents (γk), which can effectively “turn off” uninformative views, enhancing interpretability [George and McCulloch, 1993]. 3. Principled Sparsity Modeling: To explicitly model the excess zeros common in omics data using a zero-inflated mixture likelihood [Lambert, 1992]. This distinguishes between true biological absence (a point mass at zero) and low-level biological signal (a continuous component). 4. Censored Outcome Modeling: To directly model time-to-event data through an integrated Cox proportional hazards module that operates on the learned latent representations [Cox, 1972]. 5. Integration of Biological Priors: To leverage prior knowledge in the form of biological networks (A(k)) by using Graph Convolutional Network (GCN) encoders, which regularize the model and encourage biologically meaningful latent representations [Kipf and Welling, 2017]. This paper details the statistical foundation of BMGCN. We first present the notation and model overview, followed by a detailed description of the full generative process. We then derive the variational inference framework and the Evidence Lower Bound (ELBO) objective function. Finally, we describe the architecture of the graph-based encoders and outline a rigorous strategy for model validation and diagnostics. 2 Methods 2.1 Notation and Core Entities Let our dataset consist of Nsubjects and Komics views. The core entities of BMGCN are defined as follows: Observed Data:  X(k) i∈Rpk: The feature vector for subject iin view k, where pkis the number of features in that view.  A(k)∈ {0,1}pk×pk: An optional adjacency matrix for view k, representing prior biological knowledge (e.g., a protein-protein interaction network).  (Ti,∆i): The observed time for subject i, where Ti∈R+is the survival or censoring time and ∆i∈ {0,1}is the event indicator (∆i= 1 if the event was observed, ∆i= 0 if censored). Latent Variables:  Zi∈RL: A shared latent vector for subject i, capturing variation common to all views.  Z(k) i∈RLk: A view-specific latent vector for subject iin view k.  γk∈ {0,1}: A binary variable indicating whether view kis active in predicting survival.  πk∈(0,1): The prior inclusion probability for view k. 2.2 The BMGCN Generative Model As illustrated in Figure 1, BMGCN employs a dual-encoder architecture that simultaneously learns a shared latent representation Z—capturing cross-omics biological signals—and view-specific latent representations Z(k)for each omics modality k. These latent variables are then intelligently integrated: the view-specific components are dynamically gated by Bayesian selection indicators γk, which modulate their contribution based on relevance to the survival outcome. The resulting composite representation feeds into two parallel output modules: a survival head for risk prediction and a decoder for multi-omics reconstruction. Formally, BMGCN defines a full generative process for the observed data, which factorizes into the following key steps for each subject i= 1, . . . , N: 2 1. Sample View Inclusion Probabilities and Indicators: For each view k= 1, . . . , K, we first draw a prior inclusion probability from a Beta distribution and then sample the binary indicator variable: πk∼Beta(α0, β0) (1) γk|πk∼Bernoulli(πk) (2) 2. Sample Latent Variables: A shared latent vector is drawn from a standard normal prior. The view-specific latent vectors are drawn from a spike-and-slab prior, conditioned on the inclusion variable γk. If γk= 0, the latent vector is fixed at zero (the “spike”); otherwise, it is drawn from a standard normal (the “slab”). Zi∼ N(0,IL) (3) p(Z(k) i|γk) = (1 −γk)δ0(Z(k) i)+γkN(Z(k) i|0,ILk) (4) where δ0is the Dirac delta function at zero. 3. Generate Omics Features via Decoders: The observed features for each view are generated from the combined shared and view-specific latent variables. A decoder network f(k) ψmaps the concatenated latent vector [Zi,Z(k) i] to the parameters of a zero-inflated Gaussian distribution. (µ(k) i,ρ(k) i)=f(k) ψ([Zi,Z(k) i]) (5) p(X(k) ij |...) = ρ(k) ij δ0(X(k) ij ) + (1 −ρ(k) ij )N(X(k) ij |µ(k) ij , σ2 k) (6) Here, µ(k) iis the mean vector, ρ(k) iis the vector of zero-inflation probabilities (passed through a sigmoid function to ensure they are in (0,1)), and σ2 kis a view-specific variance parameter. 4. Generate Survival Outcome: A latent risk score ηifor each subject is computed by a survival head network gθ, which takes the shared latent and the gated view-specific latents as input. The observed survival data (Ti,∆i) are assumed to follow a Cox proportional hazards model conditioned on this risk score. ηi=gθ([Zi, γ1Z(1) i, . . . , γKZ(K) i]) (7) (Ti,∆i)∼Cox(ηi) (8) The likelihood for the Cox model is given by the partial log-likelihood: l(η) = N X i=1 ∆i ηi−log X j∈Ri eηj (9) where Ri={j:Tj≥Ti}is the set of subjects at risk at time Ti. This semi-parametric form is powerful because it does not require specifying a baseline hazard function h0(t). When multiple events occur at the same time (ties), approximations such as Breslow’s or Efron’s can be used. The full joint distribution over all observed and latent variables is given by: p(X,T,∆,Z,{Z(k)},γ,π) = "N Y i=1 p(Ti,∆i|ηi) K Y k=1 p(X(k) i|Zi,Z(k) i)p(Z(k) i|γk)p(Zi)#× "K Y k=1 p(γk|πk)p(πk)#(10) 2.3 Variational Inference and the ELBO The posterior distribution over the latent variables, p(Z,{Z(k)},γ|X,T,∆), is intractable due to the non-linearities in the model and the Cox likelihood. We therefore resort to amortized variational inference 3 (VI) [Kingma and Welling, 2013, Rezende et al., 2014]. We introduce a tractable variational distribution q(·) to approximate the true posterior and maximize the Evidence Lower Bound (ELBO), L, with respect to the variational parameters. We use a mean-field variational family that factorizes across subjects and latent variables: q(Z,{Z(k)},γ) = N Y i=1 q(Zi) K Y k=1 q(Z(k) i)q(γk) (11) where each factor is parameterized by an encoder network: q(Zi|Xi) = N(µi(Xi),diag(σ2 i(Xi))) (12) q(Z(k) i|X(k) i,A(k))=N(µ(k) i(X(k) i,A(k)),diag((σ(k) i)2(X(k) i,A(k)))) (13) q(γk) = Bernoulli(τk) (14) where τkis the posterior inclusion probability for view k. The ELBO is defined as: L=Eq[log p(X,T,∆,Z,{Z(k)},γ)] −Eq[log q(Z,{Z(k)},γ)] (15) This objective can be decomposed into several intuitive terms: L=Lrecon +Lsurv − LKL (16) where: 1. Reconstruction Term (Lrecon): The expected log-likelihood of the omics data. Lrecon = N X i=1 K X k=1 Eq[log p(X(k) i|Zi,Z(k) i)] (17) The piecewise nature of the zero-inflated log-likelihood makes this term directly computable without explicit sampling of discrete zero-inflation masks. 2. Survival Term (Lsurv): The expected Cox partial log-likelihood. Lsurv =Eq  N X i=1 ∆i ηi−log X j∈Ri eηj  (18) 3. KL Divergence Regularizer (LKL): A sum of Kullback-Leibler (KL) divergences that penalize deviations of the variational posterior from the prior. LKL = N X i=1 KL(q(Zi)∥p(Zi)) (19) + N X i=1 K X k=1 Eq(γk)[KL(q(Z(k) i)∥p(Z(k) i|γk))] (20) + K X k=1 KL(q(γk)∥p(γk)) (21) The KL terms for the Gaussian latents have a closed-form solution. The term involving the discrete γkis approximated as τk·KL(q(Z(k) i)∥ N(0,I)), effectively weighting the KL penalty by the posterior probability that the view is included. 2.4 Graph-Convolutional Encoders To incorporate prior biological network information (A(k)), the encoder networks that parameterize q(Z(k) i) are implemented as Graph Convolutional Networks (GCNs) [Kipf and Welling, 2017]. This 4 allows the model to learn feature representations that are regularized by known biological interactions, improving robustness and interpretability. For each view k, we first normalize the adjacency matrix to ensure stable learning: ˆ A(k)= ( ˜ D(k))−1/2˜ A(k)(˜ D(k))−1/2(22) where ˜ A(k)=A(k)+I(adding self-loops) and ˜ D(k)is the diagonal degree matrix of ˜ A(k). The GCN takes the full data matrix for view k,X(k)∈RN×pk, as input features for the nodes (e.g., genes) of the graph. The GCN then learns feature embeddings by propagating and transforming information between nodes in the graph. The propagation rule for layer lis: H(l+1) = ReLU( ˆ A(k)H(l)W(l)) (23) where H(0) = (X(k))⊤∈Rpk×N,W(l)is a learnable weight matrix, and H(l)∈Rpk×dlis the feature embedding matrix at layer l. After Llayers, we obtain a final feature embedding matrix H(L)∈Rpk×dL. To obtain sample-specific parameters for q(Z(k) i), we apply a readout function. Specifically, the final feature embeddings are used to parameterize a sample-specific network that maps the input data to the latent space parameters: h(k) i= FFNenc(X(k) i,H(L)) (24) This sample-specific embedding h(k) iis then mapped to the variational parameters: µ(k) i=Wµh(k) i(25) log σ(k) i=Wσh(k) i(26) 2.5 Optimization and Implementation The ELBO is optimized with respect to the parameters of the encoders, decoders, and survival head using stochastic gradient ascent. We employ the reparameterization trick for Gaussian latents to obtain low-variance gradient estimates [Kingma and Welling, 2013]. For the discrete view-selection variables γk, we use a hard-concrete relaxation with a straight-through gradient estimator during training [Jang et al., 2017, Maddison et al., 2017]. Training is performed in mini-batches. For the Cox survival term, which is defined over the full cohort, we use mini-batch approximations, such as constructing risk sets only from subjects within the current batch or using a memory bank. While this introduces a small bias, it is essential for scalability. Key implementation details for ensuring stable training include:  Using the log-sum-exp trick for stable computation of the Cox likelihood denominator.  Applying KL annealing to gradually introduce the KL regularization term, preventing posterior collapse [Bowman et al., 2015].  Using modern optimizers like AdamW and employing learning rate schedulers [Kingma and Ba, 2014, Loshchilov and Hutter, 2017].  Initializing the view-selection logits (τk) to favor inclusion initially to promote stable learning. 5 BMGCN Experimental Pipeline Algorithm 1 BMGCN End-to-End Training and Validation Procedure Phase 1: Global Setup 1: Initialize: Create results directory and configure logging 2: Set random seeds for NumPy, PyTorch, etc., to ensure reproducibility 3: Detect and set compute device (e.g., CUDA if available, else CPU) Phase 2: Data Simulation 4: Simulate Data: 5: X, A ←Simulate Npatients with Komics views {X(k)}K k=1 and adjacency matrices {A(k)}K k=1 6: T, ∆←Simulate survival outcomes with pre-defined ground-truth informative views 7: Create F-fold stratified cross-validation splits based on event indicator ∆ Phase 3: Model Definition 8: Define Model Architecture: 9: Encoders: Define qΦ(Z|X) and qΦ(Z(k)|X(k), A(k)) as GCNs or MLPs 10: Decoders: Define pΨ(X(k)|Z, Z(k)) mapping latents to Zero-Inflated Gaussian parameters (µ(k), ρ(k)) 11: Survival Head: Define gΘ(·) mapping latents [Z, {γk·Z(k)}] to a risk score η 12: View Selection: Define learnable logits γlogits ∈RKto produce probabilities τk=σ(γlogits,k) Phase 4: Cross-Validation and Training 13: for each fold f= 1 to Fdo 14: Load training and validation data splits for the current fold 15: Initialize a new instance of the BMGCN model and transfer to the compute device 16: Initialize optimizer (e.g., AdamW) and a learning rate scheduler 17: for epoch e= 1 to Edo 18: β←min(1.0, e/anneal epochs) ▷Anneal KL regularization weight 19: Set model to training mode 20: for each batch {Xb, Tb,∆b}in the training data loader do 21: ▷— Forward Pass — 22: Infer latents Zb,{Z(k) b}, selectors γb, reconstructions {µ(k) b, ρ(k) b}, and risk scores ηb 23: ▷— Compute ELBO Loss — 24: Lrecon ←PkZeroInflatedGaussianNLL(X(k) b, µ(k) b, ρ(k) b) 25: Lsurv ←CoxPartialLikelihood(ηb, Tb,∆b) 26: LKL ←KL(q(Zb)∥p(Z)) + PkEq(γk)[KL(q(Z(k) b)∥p(Z(k)|γk))] 27: LELBO ← Lrecon +Lsurv −β· LKL 28: ▷— Backward Pass — 29: if LELBO is finite then 30: optimizer.zero grad() 31: LELBO.backward() 32: torch.nn.utils.clip grad norm (model.parameters(), max norm=1.0) 33: optimizer.step() 34: if eis a validation epoch then 35: Evaluate Concordance Index on the validation set 36: Save model checkpoint if it has the best validation performance so far 37: Load the best performing model checkpoint for the current fold 38: Record the final C-index and learned view inclusion probabilities {τk}for the fold Phase 5: Diagnostics and Reporting 39: Aggregate and save all cross-validation metrics to a summary file (e.g., Excel) 40: Generate Plots: 41: Plot cross-validation C-index (train vs. validation) across folds 42: Plot learned view inclusion probabilities {τk}across folds 43: Plot comparison of final learned {τk}vs. ground-truth view relevance 44: Generate UMAP visualizations of the shared latent space Z, colored by survival time and event status 45: Generate posterior predictive check plots comparing original vs. reconstructed omics data distributions 46: Output: Final trained models, performance metrics tables, and diagnostic figures 6 BMGCN: Multi-Omics Integration with Bayesian View Selection for Survival Prediction Inputs Omics View 1 X(1) (RNA-seq) Omics View 2 X(2) (Methylation) Omics View K X(K) (Proteomics) Encoders Shared GCN/MLP Encoder qϕ View-Specific GCN Encoder qϕ(1) View-Specific GCN Encoder qϕ(2) View-Specific GCN Encoder qϕ(K) Z(1) Z(2) Z(K) Z View Selection γ1 γ2 γK Spike&Slab gates select relevant views Interpretable Latent Space γ1Z(1) γ2Z(2) γKZ(K) Outputs Survival Prediction Head gθ Decoder pψ Predicted Risk Score (η) Survival Outcomes Reconstructed Omics Data (b X(k)) Zero-Inflated Gaussian Likelihood Feature Importance & Sparsity Unsupervised Learning Figure 1: BMGCN overview: multi-omics inputs, dual encoders, Bayesian view selection, and multi-task outputs. 3 Validation and Diagnostics Strategy To ensure BMGCN is robust, generalizable, and interpretable, we propose a multi-faceted validation strategy.  Simulation Studies: We will generate synthetic multi-view data with known ground-truth latent structures, pre-defined informative vs. uninformative views, and known survival dependencies. This will allow us to quantitatively assess the model’s ability to recover the true latent variables (cosine similarity), correctly perform view selection (accuracy of τkvs. ground truth), and accurately predict survival (Concordance Index).  Cross-Validation: On real-world datasets, we will use a k-fold cross-validation scheme, stratified by event status, to evaluate predictive performance on held-out data. The primary metric will be the Concordance Index (C-index), which measures the fraction of concordant pairs of subjects. A secondary metric will be the integrated Brier score, which assesses the calibration of survival probability predictions over time.  Posterior Predictive Checks (PPCs): To assess goodness-of-fit, we will sample from the posterior predictive distribution. We will compare the distributions of simulated data (e.g., zero frequencies, means, variances) against the observed data to detect model misspecification.  Cox Diagnostics: We will examine Schoenfeld residuals to test the proportional hazards assumption of the Cox model, a critical assumption for the validity of the survival module.  Ablation Studies: To quantify the contribution of each model component, we will perform ablation studies by systematically removing key features: (1) the graph structure (setting A(k)=I), (2) the zero-inflation module (ρ=0), and (3) the spike-and-slab selectors (forcing γk= 1 for all views).  Latent Space Analysis: We will visualize the learned shared latent space Ziusing techniques like UMAP to check if it clusters subjects by known phenotypes, treatment groups, or survival outcomes, providing a qualitative check on its biological relevance [McInnes et al., 2018]. 4 Empirical Results This section presents the empirical evaluation of the Bayesian Multi-view Graph Convolutional Network (BMGCN) on a simulated multi-omics dataset with known ground truth. Following initial analyses that indicated suboptimal performance, key architectural enhancements were implemented, including a true GCN encoder, a decoder with learnable variance, and an attention-based survival head. The results below demonstrate the impact of these improvements on survival prediction accuracy and view selection capability. 7 4.1 Survival Prediction Performance The model exhibits strong and consistent performance across all cross-validation folds. The mean validation Concordance Index (C-Index) falls in a high-performance range, indicating reliable survival prediction. Key metrics are summarized as follows:  Mean Validation C-Index: 0.7788 ±0.0385  Another Result: Some value. As shown in Table 1, the validation C-Index exceeds 0.75 in most folds, with peak performance of 0.824. The low standard deviation confirms model robustness. Critically, the learned view-inclusion probabilities (τk) converge decisively to 0 or 1 across all folds, demonstrating training stability. Table 1: Per-Fold Cross-Validation Results of the Final Model Metric Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Train C-Index 0.901 0.910 0.893 0.904 0.905 Validation C-Index 0.812 0.722 0.782 0.754 0.824 τ0(True) 1.000 1.000 1.000 1.000 1.000 τ1(False) 0.000 0.000 0.000 0.000 0.000 τ2(True) 1.000 1.000 1.000 1.000 1.000 τ3(False) 0.000 0.000 0.000 0.000 0.000 τ4(True) 1.000 1.000 1.000 1.000 1.000 τ5(False) 0.000 0.000 0.000 0.000 0.000 τ6(True) 1.000 1.000 1.000 1.000 1.000 Figure 2: Cross-validation performance (left) and view inclusion probabilities per fold (right). The model shows high predictive accuracy and perfect separation of view relevance. 4.2 View Selection Accuracy BMGCN achieves 100% view selection accuracy, correctly identifying all four informative views (0, 2, 4, 6) and excluding the three non-informative ones. As shown in Table 2 and Figure 3, the learned inclusion probabilities τkare exactly 1.000 for relevant views and 0.000 for irrelevant ones, reflecting absolute model confidence. 4.3 Model Diagnostics Latent Space Visualization. The UMAP embedding of the shared latent space Z(Figure 4) reveals a structured representation strongly correlated with survival time and event status, confirming that the latent space captures clinically meaningful information. 8 Table 2: Final View Selection Performance View Index Ground Truth Learned τkCorrectly Identified 0 Relevant 1.000 Yes 1 Not Relevant 0.000 Yes 2 Relevant 1.000 Yes 3 Not Relevant 0.000 Yes 4 Relevant 1.000 Yes 5 Not Relevant 0.000 Yes 6 Relevant 1.000 Yes Figure 3: Comparison of learned view-inclusion probabilities (τk) against ground truth, demonstrating flawless identification of relevant and irrelevant views. Figure 4: UMAP of the shared latent space, colored by survival time (left) and event status (right). The space shows clear, clinically relevant structure. Data Reconstruction. Posterior predictive checks (Figure 5) confirm that the enhanced decoder with learnable variance accurately reconstructs the zero-inflated and sparse characteristics of the original omics data, resolving earlier reconstruction issues. 9 A.6 KL Divergence for the Spike-and-Slab Prior For the view-specific latent Z(k) i, the prior is a spike-and-slab: p(Z(k) i|γk) = (1−γk)δ0(Z(k) i)+γkN(Z(k) i| 0,I) The variational posterior is a Gaussian: q(Z(k) i)=N(µ(k) q,diag((σ(k) q)2)) The KL divergence KL(q(Z(k) i)∥p(Z(k) i|γk)) is not standard because the prior is a mixture. We handle this by taking the expectation over the discrete variable γk. Eq(γk)hKL(q(Z(k) i)∥p(Z(k) i|γk))i=Eq(γk)"Zq(Z(k) i) log q(Z(k) i) p(Z(k) i|γk)!dZ(k) i# =q(γk= 0) ·KL(q(Z(k) i)∥δ0)+q(γk= 1) ·KL(q(Z(k) i)∥ N(0,I)) (56) The term KL(q(Z(k) i)∥δ0) is problematic because the delta function is not a density with respect to the Lebesgue measure. In practice, as noted in the main text, when γk= 0, the prior forces Z(k) i=0. Therefore, the model approximates this by only applying the KL penalty against the Gaussian slab when the view is included. This leads to the practical approximation: Eq(γk)hKL(q(Z(k) i)∥p(Z(k) i|γk))i≈τk·KL(q(Z(k) i)∥ N(0,I)) (57) where τk=q(γk= 1). This approximation is used in Equation (20). A.7 Derivation of the Variational Objective for the Spike-and-Slab Prior The exact KL divergence for the view-specific latent, marginalized over γk, is: KL q(Z(k) i)∥p(Z(k) i)=Zq(Z(k) i) log q(Z(k) i) p(Z(k) i)dZ(k) i =Zq(Z(k) i) log q(Z(k) i)dZ(k) i−Zq(Z(k) i) log Zp(Z(k) i|γk)p(γk)dγkdZ(k) i =−H hq(Z(k) i)i−Eq(Z(k) i)hlog (1 −πk)δ0(Z(k) i)+πkN(Z(k) i|0,I)i (58) where Hhq(Z(k) i)iis the entropy of the variational distribution. The intractability arises from the log of a sum inside the expectation. The common variational approach is to introduce an auxiliary variational distribution r(γk) for the discrete variable and derive a lower bound. However, in the context of the overall ELBO, this can become very complex. The approximation used in the paper, Eq(γk)hKL(q(Z(k) i)∥p(Z(k) i|γk))i, is itself a variational approximation that treats γkas independent. This is known as the “fully-factorized” or “mean-field” assumption across the discrete and continuous latents. The approximation in Equation (57) (using τk·KL(q∥ N)) is a further simplification that ignores the KL divergence to the spike (delta function) because it is infinite for any non-degenerate q. This is justified by the model’s structure: when γk= 0, the survival head gθreceives 0, so the model is incentivized to set q(Z(k) i) to δ0to minimize the reconstruction loss, making the infinite KL penalty moot. The τkweighting then acts as a regularizer only when the view is likely to be included. A.8 Derivation of the Optimal Variational Distribution for πk The appendix mentions a variational distribution q(πk) but does not derive its optimal form. For a conjugate model, we can find a closed-form update. From the joint distribution in Equation (33), the conditional posterior for πk(before variational approximation) is proportional to: p(πk|γk)∝p(γk|πk)p(πk) ∝πγk k(1 −πk)1−γk·πα0−1 k(1 −πk)β0−1 ∝πα0+γk−1 k(1 −πk)β0+1−γk−1(59) 16 This shows that the conditional posterior for πkis Beta(α0+γk, β0+1−γk). In variational inference, the optimal q⋆(πk) that maximizes the ELBO is the one that minimizes the KL divergence to the true conditional posterior. Under the mean-field assumption, the optimal q⋆(πk) is: q⋆(πk)∝exp Eq(γk)[log p(πk|γk)] ∝exp Eq(γk)[log Beta(πk|α0+γk, β0+1−γk)](60) Since the expectation of γkunder qis τk, we can substitute to find: q⋆(πk)∝πα0+τk−1 k(1 −πk)β0+1−τk−1 = Beta(πk|α0+τk, β0+1−τk) (61) Therefore, the optimal variational distribution for πkis q(πk) = Beta(α0+τk, β0+1−τk). The corresponding KL divergence term KL(q(πk)∥p(πk)) in Equation (39) can then be computed in closed form using the KL divergence between two Beta distributions. A.9 Graph Convolutional Layer: Signal Processing Interpretation The propagation rule H(l+1) = ReLU( ˆ A(k)H(l)W(l)) can be analyzed from a signal processing perspective on graphs. The normalized adjacency matrix ˆ A(k)can be decomposed via its eigenvectors: ˆ A(k)=UΛU⊤, where Λis a diagonal matrix of eigenvalues and Uis the matrix of eigenvectors. This is the graph Fourier transform. A graph convolution in the spatial domain (like the GCN layer) is equivalent to a multiplication in the spectral domain. The operation ˆ A(k)H(l)is a form of low-pass filter. It smooths the feature representation of each node by averaging it with the features of its neighbors. The eigenvalues in Λ are in the range [−1,1], and the largest eigenvalues correspond to the lowest frequencies (most global, smoothest signals on the graph). The weight matrix W(l)then acts as a learnable filter that can amplify or attenuate different frequency components. The ReLU non-linearity introduces the capacity to learn more complex, non-linear feature interactions. This mathematical interpretation justifies why GCNs are effective for biological data: they promote smoothness of the latent representations with respect to the prior biological network, meaning that genes/proteins that are known to interact will have similar representations in the latent space, leading to more biologically plausible and robust models. A.10 Hard-Concrete Relaxation for γk The main text mentions using a “hard-concrete relaxation with a straight-through gradient estimator” for γk. Here, we provide the mathematical formulation of this trick. The binary variable γk∈ {0,1}is replaced during training with a continuous relaxation ˜γk∈(0,1) sampled from a Hard-Concrete distribution. The sampling process is: u∼Uniform(0,1) s=σ1 β(log u−log(1 −u) + log α) ˜s=s·(γmax −γmin)+γmin ˜γk= min(1,max(0,˜s)) (62) where αis the location parameter (often derived from the logit of τk), βis a temperature parameter, and [γmin, γmax] is a small interval (e.g., [−0.1,1.1]) to stretch the output. During the forward pass, ˜γkis used. For the backward pass (gradient computation), the “straightthrough” estimator is used: the gradient is computed as if the operation was an identity function. That is, if y= hard concrete(x), then ∂y/∂x = 1. This allows gradients to flow through the discrete sampling process, enabling end-to-end training. As training progresses, the temperature βis annealed to zero, causing the Hard-Concrete distribution to collapse to a Bernoulli distribution, effectively discretizing ˜γkback to γk. 17 A.11 Derivation of Gradients for Model Parameters The appendix derives the ELBO with respect to the variational parameters. It is also crucial to derive the gradients with respect to the model parameters ψ(decoder) and θ(survival head). The ELBO is: L=Eq[log p(X,T,∆,Z,{Z(k)},γ)] −Eq[log q(Z,{Z(k)},γ)] (63) The second term (the entropy of q) does not depend on the model parameters ψ, θ. Therefore, the gradient with respect to ψand θcomes only from the first term, the expected joint log-likelihood. ∇ψ,θL=∇ψ,θEq[log p(X,T,∆,Z,{Z(k)},γ)] (64) This expectation can be broken down as shown in Equation (34). The gradients for the reconstruction term Lrecon with respect to ψare derived in Section A.3. The gradients for the survival term Lsurv with respect to θare derived in Section A.4. The key point is that these gradients can be estimated using Monte Carlo sampling from the variational distribution q, combined with the reparameterization trick for the Gaussian latents and the straight-through estimator for the discrete latents, making the entire model trainable via standard stochastic gradient descent. A.12 Reparameterization Trick for Gaussian Latents To obtain low-variance, unbiased gradient estimates for the Gaussian latent variables, we use the reparameterization trick. For a latent variable z∼ N(µ,Σ), instead of sampling zdirectly, we sample an auxiliary noise variable ϵ∼ N(0,I) and compute z=µ+Σ1/2ϵ, where Σ1/2is the matrix square root. For a diagonal covariance Σ= diag(σ2 1, . . . , σ2 d), this simplifies to element-wise multiplication: z=µ+σ⊙ϵ, where ⊙denotes the Hadamard (element-wise) product. The gradient of an expectation Eq(z)[f(z)] with respect to µand σcan then be computed as: ∂ ∂µ Eq(z)[f(z)] = Eϵ∂f(z) ∂z·∂z ∂µ=Eϵ∂f(z) ∂z(65) ∂ ∂σ Eq(z)[f(z)] = Eϵ∂f(z) ∂z·∂z ∂σ=Eϵ∂f(z) ∂z·ϵ(66) This allows gradients to flow through the sampling process, enabling efficient optimization via stochastic gradient descent. 18