Integration of Graph Neural Network and Neural-ODEs for Tumor Dynamic
  Prediction
In anti-cancer drug development, a major scientific challenge is
disentangling the complex relationships between high-dimensional genomics data
from patient tumor samples, the corresponding tumor's organ of origin, the drug
targets associated with given treatments and the resulting treatment response.
Furthermore, to realize the aspirations of precision medicine in identifying
and adjusting treatments for patients depending on the therapeutic response,
there is a need for building tumor dynamic models that can integrate both
longitudinal tumor size as well as multimodal, high-content data. In this work,
we take a step towards enhancing personalized tumor dynamic predictions by
proposing a heterogeneous graph encoder that utilizes a bipartite Graph
Convolutional Neural network (GCN) combined with Neural Ordinary Differential
Equations (Neural-ODEs). We applied the methodology to a large collection of
patient-derived xenograft (PDX) data, spanning a wide variety of treatments (as
well as their combinations) on tumors that originated from a number of
different organs. We first show that the methodology is able to discover a
tumor dynamic model that significantly improves upon an empirical model which
is in current use. Additionally, we show that the graph encoder is able to
effectively utilize multimodal data to enhance tumor predictions. Our findings
indicate that the methodology holds significant promise and offers potential
applications in pre-clinical settings.
Omid Bazgir, Zichen Wang, Ji Won Park, Marc Hafner, James Lu
# Integration of Graph Neural Network and Neural-ODEs for Tumor Dynamic Prediction

###### Abstract

In anti-cancer drug development, a major scientific challenge is disentangling the complex relationships between high-dimensional genomics data from patient tumor samples, the corresponding tumor's organ of origin, the drug targets associated with given treatments and the resulting treatment response. Furthermore, to realize the aspirations of precision medicine in identifying and adjusting treatments for patients depending on the therapeutic response, there is a need for building tumor dynamic models that can integrate both longitudinal tumor size as well as multimodal, high-content data. In this work, we take a step towards enhancing personalized tumor dynamic predictions by proposing a heterogeneous graph encoder that utilizes a bipartite Graph Convolutional Neural network (GCN) combined with Neural Ordinary Differential Equations (Neural-ODEs). We applied the methodology to a large collection of patient-derived xenograft (PDX) data, spanning a wide variety of treatments (as well as their combinations) on tumors that originated from a number of different organs. We first show that the methodology is able to discover a tumor dynamic model that significantly improves upon an empirical model which is in current use. Additionally, we show that the graph encoder is able to effectively utilize multimodal data to enhance tumor predictions. Our findings indicate that the methodology holds significant promise and offers potential applications in pre-clinical settings.

## 1 Introduction

In the development of novel anti-cancer therapies, patient-derived xenografts (PDX) have become an important platform for addressing key questions, such as evaluating the treatment response to therapeutic agents and the combinations thereof, identifying the relevant biomarkers of response and the mechanisms of resistance development Byrne et al. (2017). Furthermore, as PDX models are obtained by surgically removing patients' tumor and implanting them in mice, co-clinical avatar trials Byrne et al. (2017) can be performed, whereby treatment of the patients occur simultaneously to the treatment of the corresponding (pre-clinical) PDX models generated from the same patients. Such avatar studies could enable for real-time clinical decision making, to identify the best treatment(s) that the patients could respond to and help deliver some of the promises of precision medicine.

Given the myriad applications of PDX models, an important computational task is the prediction of their dynamic response from the baseline -omics data and/or early tumor size data. This is a particularly challenging task due to the need to meld high-dimensional omics data measured on baseline (e.g., RNA-seq pre-treatment), with the low-dimensional but serially assessed tumor size
measurements under treatment. While empirical Zwep et al. (2021) and spline-based Forrest et al. (2020) tumor dynamic models have been proposed, there has been little progress in melding such dynamic models with high dimensional -omics data. In Zwep et al. (2021), an machine learning (ML) approach has been proposed to use a model based on least absolute shrinkage and selection operator (LASSO) to predict tumor dynamic parameters from copy-number variations (CNVs) of genes from a large PDX data set consisting of various treatments Gao et al. (2015). In Ma et al. (2021), a few-shot learning (multi-layer perceptron) approach, and in Peng et al. (2022), a heterogeneous GCN approach, have been proposed to learn drug responses from in-vitro (cell-line) data and predict in-vivo (PDX) outcomes within two response categories. While the predictions made by the models in Ma et al. (2021) and Peng et al. (2022) show promising correlations with true responses, they do not necessarily capture the tumor dynamic, which is crucial for clinical decision-making. While the proposed method shows promise, improving the predictivity of the model and including omics data remains an important topic for further research.

Over the past few years, Neural-Ordinary Differential Equations (NODEs) Chen et al. (2018) has emerged as a promising deep learning (DL) methodology for making predictions from irregularly sampled temporal data. Recently, a methodology based on NODE has been developed for tumor dynamic modeling in the clinical trial setting and the generated embeddings from patients' tumor data have been demonstrated to be effective for predicting their Overall Survival (OS) Laurie and Lu (2023). While the Tumor Dynamic NODE (TDNODE) has set the mathematical foundations for modeling longitudinal tumor data Laurie and Lu (2023), a methodology for the incorporation of high-dimensional, multimodal data into such models has yet to be developed.

In this work, we propose a novel way to combine the previously developed tumor volume encoder with a heterogeneous graph encoder (see Fig. 1). The latter takes as inputs multimodal data consisting of drugs, disease and RNA-seq by incorporating graphs that encapsulate drug-gene associations, disease-gene associations and gene-gene interactions respectively. In application to the PDX data of Gao et al. (2015), we show that by leveraging the multimodal data including RNA-seq in conjunction with early tumor response data, the proposed methodology can significantly improve predictive performance of future tumor response. Finally, we summarize the findings and discuss the potential future applications of the proposed approach in enabling precision medicine in oncology.

Figure 1: **Model architecture overview.****A)** Integration of GCNs and Neural-ODEs for tumor dynamic prediction. **B) Heterogeneous Graph Encoder:** two bipartite graph attention convolution NN to extract disease-gene, and drug-gene association associations, and a graph convolution NN to extract gene-gene interactions, and integrated into emdedding space in baseline level.
Method

Within our multi-modal framework, we constructed a multi-relational network using three large datasets covering interactions between drugs, genes, and diseases. We used RNA-seq data as node features within tissue-specific knowledge graphs and further integrated it with drug targets of treatments as a heterogeneous graph Zhang et al. (2019).

This integration allowed us to learn an embedding that captures the complex relationships between genes, tumors, and treatments. The following 3 graphs were utilized to represent both biological and pharmacological knowledge:

* **gene-gene graph:** tissue-specific functional gene networks obtained from TissueNexus Lin et al. (2022), which provides gene-gene associations based on gene expression across 49 human tissues.
* **drug-gene graph:** obtained from DGIdb Cotto et al. (2018), which consolidates information on drug target and interacting genes from 30 disparate sources through expert curation and text mining.
* **disease-gene graph:** obtained from DisGeNET Pinero et al. (2016), one of the largest collections of genes and variants associated with human diseases.

For each PDX model, the pre-treatment embeddings using the above multi-modal graphs are used in conjunction with the early, observed tumor measurements to predict the future tumor dynamic profiles. The following subsections cover these two respective aspects.

### Heterogeneous Graph Encoder

We formulate the PDX representation learning task as one of graph embedding, by fusing information from a heterogeneous network that incorporates drug, disease, and gene relationships. We apply multi-layer GCN Welling and Kipf (2016) in the RNA domain to model the gene interactions (_gene-gene graph_). Additionally, we use a bipartite graph attention convolutions Wang et al. (2020); Nassar (2018) for message passing from drugs to target genes (_drug-gene graph_), as well as gene embeddings to the disease domain (_disease-gene graph_). The mathematical formulation of our proposed framework is described as follows.

**Bipartite Graphs Attention Convolution.** Conventional GCN assumes that all nodes belong to the same category. However, in our scenario there are heterogeneous attributes across various node types, such as genes, diseases, and drug targets. This limitation becomes evident when node attributes span across different domains. Consequently, we adopt a bipartite graph as a natural representation for modeling inter-domain interactions among distinct node types. We adapt GCN to operate within a bipartite graph, where node feature aggregation exclusively occurs over inter-domain edges. Specifically, let us denote graphs as \(\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)\), where \(\mathcal{V}\) represents the set of vertices, given by \(\left\{v_{i}\right\}_{i=1}^{N}\), and \(\mathcal{E}\) is the set of edges. We consider a bipartite graph \(\mathcal{BG}(\mathcal{U},\mathcal{V},\mathcal{E})\) defined as a graph \(\mathcal{G}(\mathcal{U}\cup\mathcal{V},\mathcal{E})\), where \(\mathcal{U}\) and \(\mathcal{V}\) represent two sets of vertices (nodes) corresponding to the two respective domains. Here, \(u_{i}\) and \(v_{j}\) denote the \(i\)-th and \(j\)-th node in \(\mathcal{U}\) and \(\mathcal{V}\), respectively, where \(i=1,2,\ldots,M\) and \(j=1,2,\ldots,N\). All edges within the bipartite graph exclusively connect nodes from \(\mathcal{U}\) and \(\mathcal{V}\) (i.e., \(\mathcal{E}=\{(u,v)|u\in\mathcal{U},v\in\mathcal{V}\}\)). The features of the two sets of nodes are denoted by \(X_{u}\) and \(X_{v}\), where \(X_{u}\in\mathbb{R}^{M\times P}\) is a feature matrix with \(\vec{x}_{u_{i}}\in\mathbb{R}^{P}\) representing the feature vector of node \(u_{i}\), and \(X_{v}\in\mathbb{R}^{N\times Q}\) is defined similarly. For the message passing \(\mathrm{MP}_{v\to u}\) from domain \(\mathcal{V}\) to \(\mathcal{U}\), we define a general bipartite graph convolution (_bg_) as:

\[bg_{\mathcal{E}}(u_{i})=\rho\left(\mathrm{agg}\left(\{W_{u_{i},v_{j}}\vec{x}_{ v_{j}}|v_{j}\in\mathcal{N}_{u_{i}}^{\mathcal{E}}\}\right)\right), \tag{1}\]

where \(\mathcal{N}_{u_{i}}^{\mathcal{E}}\) represents the neighborhood of node \(u_{i}\) connected by \(\mathcal{E}\) in \(\mathcal{BG}(\mathcal{U},\mathcal{V},\mathcal{E})\) (i.e., \(\mathcal{N}_{u_{i}}^{\mathcal{E}}\subset\mathcal{V}\)). \(W_{u_{i},v_{j}}\in\mathbb{R}^{M\times N}\) is a feature weighting kernel transforming \(N\)-dimensional features to \(M\)-dimensional features, the \(\mathrm{agg}\) is a permutation-invariant aggregation operation, and the \(\rho\) operator can be a non-linear activation function. In our work, we used element-wise mean-pooling and ReLU Nair and Hinton (2010) for \(\mathrm{agg}\) and \(\rho\) respectively.

Our bipartite graph convolution layers utilize the graph attention network Velickovic et al. (2017) as the backbone on the node features, resulting in the bipartite graph attention convolution layer
(\(bga\)). Since the attention mechanism considers features of two sets of nodes, we specifically define a learnable matrix \(W^{u}\in\mathbb{R}^{P\times S}\) (resp. \(W^{v}\in\mathbb{R}^{Q\times S}\)) for \(X_{u}\) (resp. \(X_{v}\)). The \(bga\) can be formulated as:

\[bga_{\mathcal{E}}(u_{i})=\mathrm{ReLU}\left(\sum_{v_{j}\in\mathcal{N}^{ \mathcal{E}}_{u_{i}}}\alpha_{u_{i},v_{j}}W^{v}\vec{x}_{v_{j}}\right), \tag{2}\]

The attention mechanism is a single-layer feedforward neural network, parameterized by a weight vector \(\vec{a}\) and applying the LeakyReLU non-linearity function. The attention weight coefficients can be expressed as:

\[\alpha_{u_{i},v_{j}}=\frac{\exp\left(\rho(\vec{a}^{T}[W^{u}\vec{x}_{u_{i}}\|W^ {v}\vec{x}_{v_{j}}])\right)}{\sum_{v_{k}\in\mathcal{N}^{\mathcal{E}}_{u_{i}}} \exp\left(\rho(\vec{a}^{T}[W^{u}\vec{x}_{u_{i}}\|W^{v}\vec{x}_{v_{k}}])\right)}, \tag{3}\]

where \(T\) and \(||\) represent the matrix transposition and concatenation operations respectively.

**Information Fusion.** We represent the heterogeneous network as an undirected graph \(\mathcal{G}(\mathcal{V},\mathcal{E})\) with the following three sets of nodes: drugs (\(\mathcal{V}^{A}\)), diseases (\(\mathcal{V}^{B}\)), and genes (\(\mathcal{V}^{C}\)). The initial features of these three sets of nodes are denoted as \(X_{\mathcal{V}^{A}}\), \(X_{\mathcal{V}^{B}}\), and \(X_{\mathcal{V}^{C}}\), respectively. The edges consist of two inter-domain sets representing drug-gene associations (\(\mathcal{E}^{AC}\)) and disease-gene associations (\(\mathcal{E}^{BC}\)), as well as an intra-domain set representing gene network connections (\(\mathcal{E}^{C\mathcal{E}}\)).

First, we apply multi-layer GCNs to the RNA domain to model the gene interactions. Importantly, the module parameters and gene node features are initialized by a pre-trained variational graph auto-encoder (VGAE). The VGAE was trained using the RNA-seq data from all the PDX models in this study, whereby each PDX model was represented as a graph with genes as the nodes of the graph. The tumor-specific graphs were obtained from TissueNexus Lin et al. (2022).

In the second step, we apply a single bipartite graph attention convolution layer to propagate the message from drugs to target genes. Conceptually, this step can be viewed as projecting information from the macro level (e.g., from the domain of drugs) to the micro level (e.g., the domain of genes). To formulate the message-passing step, we represent the hidden embeddings of node \(v^{c}_{i}\) as \(h^{(k)}_{v^{c}_{i}}\), where \(k\) is the computational step and when \(k=0\), \(h^{(0)}_{v^{c}_{i}}=x_{v^{c}_{i}}\). The \(h^{(k)}_{v^{c}_{i}}\) is computed as follows:

\[\mathrm{MP}^{(k)}_{\mathcal{V}^{A}\rightarrow\mathcal{V}^{C}}:h^{(k)}_{v^{c}_ {i}}=bga_{\mathcal{E}^{AC}}(v^{c}_{i})+h^{(k-1)}_{v^{c}_{i}}, \tag{4}\]

Similarly in the last step, we utilize the non-linear graph information captured by gene nodes to update the hidden embeddings of the disease nodes. In particular, we apply another bipartite graph attention convolution layer to project gene embeddings to the disease domain. Therefore, the third step can be viewed as an attentional pooling of the disease gene subgraph. The updated feature representations of disease nodes are computed as follows:

\[\mathrm{MP}^{(k)}_{\mathcal{V}^{C}\rightarrow\mathcal{V}^{B}}:h^{(k)}_{v^{b}_ {i}}=bga_{\mathcal{E}^{CB}}(v^{b}_{i})+h^{(k-1)}_{v^{b}_{i}}. \tag{5}\]

We concatenate the updated drug and disease embeddings into a unified representation for PDX experiments (tumor-treatment combinations) as \(\beta_{1}\), which can be used as feature inputs for downstream tasks.

### Tumor Volume Prediction

We formulate a tumor volume dynamic model using a Neural-ODE which utilizes both the embedding generated from the heterogeneous graph encoder, as well as the embedding generated from feeding early tumor volume data into an encoder. The model aims to utilize early tumor volume data together with baseline embedding (encapsulating drug, disease and RNA-seq data) in order to make personalized predictions.

**Tumor Volume Encoder.** In a similar manner to Laurie and Lu (2023), we implemented a tumor volume encoder to inform the Neural-ODE in order to make predictions individualized to a specific PDX model. This recurrent neural network (RNN) based encoder maps a short window of early observed tumor volumes (of an arbitrary length) into an embedding that we denote as \(\beta_{2}\).

**Neural-ODE.** In the clinical context, an approach to model tumor dynamic data using Neural-ODE has been developed, demonstrating the ability of such a formalism to discover the right dynamical
model from longitudinal patient tumor size data across treatment arms Laurie and Lu (2023). In this work, we generalized the methodology to meld longitudinal measurements with the multi-modal data into a single deep learning framework and demonstrated in the setting of dynamic modeling of PDX data. We consider a dynamical system of the following form:

\[\frac{dy(t)}{dt}=f_{\theta}(y(t),\beta),t\in[0,T] \tag{6}\]

where \(0\) and \(T\) denote the start time of PDX experimentation and the end of the prediction time respectively, \(f_{\theta}\) is a neural network parameterized by a set of weights \(\theta\) to be learned across all PDX data, and \(\beta=[\beta_{1}||\beta_{2}]\) as the PDX embedding obtained by concatenating the outputs of the heterogeneous graph and tumor volume encoders respectively. Thus, a dynamical law represented by \(f_{\theta}\) is learned across all PDX data, with the concatenated embedding \(\beta\) serving to provide the initial condition for the Neural-ODE specific to the PDX model of interest. After simulating eqn. 6 to obtain the time evolution of state \(y(t)\), we then reconstruct the tumor volume data using a two-layer MLP reducer. The Neural-ODE approach has the benefits of handle variable-length sequences, accommodate varying observation intervals, and incorporate positional encoding for capturing temporal context. To enhance model simplicity and generalization, we keep the ODE system dimension no larger than the existing tumor growth inhibition (TGI) modelZwep et al. (2021) The TGI model captured the longitudinal tumor volume measurements, per PDX, with two empirical ODEs through estimation of three parameters: growth rate, treatment efficacy, and time-dependent resistance development (kr).

**Loss Function.** As a measure of the overall fit of the model prediction to the observed tumor volume values, we incorporated the summation of the mean squared error (MSE) and the mean absolute percentage error (MAPE) De Myttenaere et al. (2016). The MSE is used to maintain the model accuracy, and MAPE is used to balance in capturing the tumor volume growth or shrinkage trends in the training set.

## 3 PDX dataset

The primary dataset utilized in this study was obtained from a large-scale pre-clinical study conducted in PDX mice models, as detailed in Gao et al. (2015). This dataset encompassed more than 1000 PDX models, each characterized by their baseline mRNA expression levels prior to treatment. In total, the dataset covered 62 distinct treatments across six different diseases, with tumor volume measurements every 2-3 days. For our analysis, we included data from 191 unique tumors and 59 different treatments, resulting in a comprehensive dataset of 3470 PDX experiments (combinations of tumor and treatment) spanning 5 tumor types. This selection was based on the availability of RNA-seq data.

## 4 Results

To assess the effectiveness of the heterogeneous graph encoder in describing changes in tumor volume, we trained the encoder using mRECIST Therasse et al. (2000); Gao et al. (2015) response labels. Specifically, we consolidated the response categories (CR, PR, and SD) into a single response category, and PD as the second category (Table A.1 in the Appendix). To predict these response categories, we used a binary classifier using a three-layer MLP neural network, that takes the heterogeneous graph encoder embedding as the input. The dataset was split by PDX models, and standard 5-fold cross-validation experiments were performed. We evaluated the prediction performance using three metrics: balanced accuracy, area under the receiver operating characteristic (AUROC), and F1-score. Our model's performance was compared against several competing approaches, including a non-graph-based deep learning method and traditional machine learning methods. To investigate the impact of the pre-training strategy, we also implemented a variant of our model with a randomly initialized gene graph.

**GCNs outperforms baseline methods in modeling RNA-seq for treatment response prediction.** As shown in Table 1, GCN-based methods demonstrated superior performance in the responder classification task across all evaluation metrics. These models outperformed non-graph-based approaches by more than 9% in F1 score and 5% in AUROC. Additionally, we observed that
our model achieved a noticeable improvement of approximately 4% in F1 score and 3% in AUROC when utilizing gene graph pre-training. This result highlights the effectiveness of graph reconstruction pre-training in enhancing gene latent representations.

### Tumor Volume Trace Fitting

We aimed to evaluate capability of the Neural-ODE in capturing the longitudinal tumor volume data trend in comparison with the state-of-the-art TGI model for PDX proposed by Zweg et al. (2021). We utilized all available tumor volume data for each PDX experiment in this experiment. Our approach involved encoding the longitudinal tumor volume data into a latent space using an RNN-based encoder. Subsequently, we used this latent space as an part of the initial condition to solve an ODE system using the Neural-ODE model to reconstruct the dynamic tumor volume data. The population level results are shown in Panel (A) of Figure 2. Notably, the Neural-ODE model outperformed the TGI model with R2 of 0.96 compared to 0.71, and Spearman correlation of 0.96 compared to 0.86.

### Tumor Volume Trace Prediction

We evaluated the ability of our model (presented in the Figure 1) to predict of future tumor volume dynamic based on a limited observed longitudinal tumor data. We selected the observation windows

\begin{table}
\begin{tabular}{l c c c} \hline Method & Balanced Accuracy & AUROC & F1 \\ \hline Pretrained GCN & **0.681** & **0.741** & **0.623** \\ GCN & 0.652 & 0.713 & 0.582 \\ MLP & 0.621 & 0.688 & 0.525 \\ SVM/RF & 0.618 & 0.676 & 0.532 \\ \hline \end{tabular}
\end{table}
Table 1: The summary of model performance on treatment response prediction using RNA-seq data. The mean over 5-fold cross validation is reported.

Figure 2: **Tumor dynamic prediction: A)** A comparison on tumor volume trace fitting between tumor growth inhibition (TGI) model proposed by Zweg et al. (2021) and our proposed Neural-ODE system. **B)** Tumor volume trace prediction using the proposed model for 3 different PDXs on 3 different row of curves. The blue curves represent the ground-truth, the green curves represent our proposed model prediction, and the red curves represent the TGI model predictions. Each column shows, one observation window, and the highlighted blocks are the observations windows where our proposed model reaches to correct classification of mRECIST response category (same as the ground-truth).
of 7, 14, 21, and 28 days, to simulate real-world scenarios where early observations are used to forecast the (future unseen) tumor volume trajectory.

We assessed the predictive performance of our model in two ways. Firstly, we employed R2 to quantify the accuracy of our model in predicting unseen tumor volumes. The results in Table 2 indicate the following: 1) The embedding learned from the heterogeneous graph encoder enhances the predictive performance of our proposed model, and 2) as the observation window size increases, our proposed model captures the unseen tumor dynamic more accurately. Additionally, as it is demonstrated in Panel (B) of Figure 2, the model effectively captures the tumor dynamic trend. However the due to noise inherent in the tumor volume measurements and the clinical significance of mRECIST response category prediction, we also assessed our proposed model's predictive performance as a classifier. The mRECIST categories are derived from the predicted tumor volume time series by applying response criteria. This evaluation measures the model's performance in correctly classifying the treatment responses based on the predicted tumor volume dynamics. Figure 3 summarizes the classification results, revealing an observable trend in which incorporating the heterogeneous graph encoder embedding improves the prediction of response categories across all observation windows.

## 5 Conclusion

In summary, we proposed a novel approach for tumor dynamic prediction that integrates RNA-seq, treatment, disease and longitudinal tumor volume data in an Neural-ODE system in a pre-clinical, PDX setting. We demonstrated that the use of Neural-ODE vastly improved the ability of the model to capture PDX tumor data than a previously proposed TGI model, as well as the benefit of adding the graph encoder to enrich the longitudinal data. As an area for further work, disentangling how the model predictions arise from the multimodal data using explainability techniques and/or attention weights is an important topic to advance our scientific understanding of the complex interplay between gene expression profiles, tumor location and drug targets. This methodology holds significant promise and warrants further evaluations, including in the clinical setting.

\begin{table}
\begin{tabular}{c c c} \hline Observation window & w/o graph encoder & w graph encoder \\ \hline
7 days & 0.233 & **0.302** \\
14 days & 0.456 & **0.479** \\
21 days & 0.586 & **0.608** \\
28 days & 0.652 & **0.659** \\ \hline \end{tabular}
\end{table}
Table 2: Predictive performance of our proposed model using different observation windows quantified with R2. The mean over 5-fold cross validation is reported.

Figure 3: Predictive performance of our proposed model as a classifier for mRECIST categories, with and without the heterogeneous graph encoder and considering different lengths of observation windows.