AI-Driven Hybrid Ecological Model for Predicting Oncolytic Viral Therapy Dynamics

The study by Kiyokawa et al. (2021) investigated the synergistic effects of ICOVIR17, a hyaluronidase-expressing oncolytic adenovirus, on a murine glioblastoma model (005 cells) treatment in combination with anti-PD-1 immune checkpoint blockade. Key findings of the combination therapy included an increased ratio of CD8+ T cells to regulatory T cells and enhanced cytotoxicity of T cells against GBM cells ex vivo (Kiyokawa et al., 2021). The underlying mechanism was found to be that ICOVIR17 degrades hyaluronan (HA) in the extracellular matrix (ECM), thus, enhancing tumor-infiltrating immune cells, activating NF-kB signaling in macrophages, and promoting a proinflammatory tumor microenvironment. These changes were associated with increased expression of PD-L1 and an enhancement of therapeutic efficacy when combined with PD-L1 blockade.

However, some open questions from these findings include: How do the virus and anti-PD-1 therapies interact with each other? Can we identify molecular mediators that drive such interactions and complex dynamics? Can a mathematical or computational model accurately predict the contributions of tumor microenvironment (ECM) remodeling and immune activation to therapeutic efficacy, including the dynamics of immune markers such as TNF and NF-kB identified in this study?

To address these questions, we will use the data from Kiyokawa et al. (2021) to train and inform our computational modeling, particularly in exploring the interactions between oncolytic virotherapy and immune checkpoint inhibitors.

State-of-the-art OVT models in mathematical oncology and computational medicine, ranging from ordinary differential equations (ODEs) to reaction-diffusion equations (RDEs), to attractor reconstruction, have explored tumor-virus interactions, capturing complex dynamics like collective oscillations and chaos (Sherlock and Coster, 2023), (Wang et al., 2019), (Heidbuechel et al., 2020), (Ramaj and Zou, 2023), (Jenner et al., 2022). Early ODE models incorporated tumor-virus heterogeneity, capturing dormancy and recurrence, though they struggled with parameter optimization and data availability (Karev et al., 2006). RDE-based multiscale models simulated virus spread and tumor responses, revealing oscillatory dynamics but were sensitive to stochastic fluctuations (Paiva et al., 2009). Nonlinear ODE models showed complex oscillatory behaviors, including limit cycles and chaotic attractors, highlighting time delay and immune-tumor interactions, though they lacked empirical validation (Agarwal and Bhadauria, 2011). Recent ODE models with viral and tumor cell densities, optimized to experimental data ($R^{2}=0.4286$), had limitations in modeling virus replication dynamics (Jenner et al., 2018).

Wodarz and Komarova (2009) developed an ODE model emphasizing spatial constraints, but it lacked immune response and molecular insights, limiting its clinical predictive power. Mahasa et al. (2017) used Delayed-Differential Equations (DDEs) to show that partial normal-cell infection enhances viral potency, but with high parameter sensitivity and no clinical validation. Similarly, Friedman and Lai (2018) modeled combination therapy using RDEs, noting reduced efficacy with higher inhibitor doses, but lacked molecular mechanisms. Agent-based models (ABMs) have also been explored, such as the glioblastoma model by Jenner et al. (2022), which showed that stromal density impacts therapeutic efficacy, but lacked long-term clinical validation and predictive ability for adaptive dynamics.

Current OVT models, while informative, face significant challenges in clinical translatability, including limitations in parameter optimization, patient-specific dynamics, and forecasting clinically relevant biomarkers. These challenges highlight the potential of DDEs and RDEs to predict complex attractors and oscillatory dynamics but also underscore the need for enhanced clinical data-based validation and predictive robustness. In principle, the complex oscillatory dynamics of OVT ecosystems—such as limit cycles and chaotic attractors—enable adaptive timing and dosage strategies to maintain tumor-virus dynamics in a non-therapy-resistant oscillatory phase. This resembles predator-prey systems, where oscillations prevent bifurcations into resistant or eradicated states. In this context, OVT can steer malignant cells in a feedback loop between treatment-resistant and therapy-susceptible (responding) phases, adapting dynamically to therapeutic interventions.

We predict that multimodal therapies, like OVT combined with immune checkpoint inhibitors or photodynamic therapy with modulated light pulses, informed by predictive biomarkers identified by our model, can reinforce these oscillations by synchronizing treatment (entrainment) with the system’s natural periodicity, thus maintaining susceptibility while avoiding resistance (Agarwal and Bhadauria, 2011; Shimizu et al., 2023; Kiyokawa et al., 2021). Given this, we hypothesize that combining DDEs with the Generalized Lotka-Volterra (GLV) model, a predator-prey framework, provides an optimal basis to capture the evolutionary dynamics, temporal delays, and feedback loops critical to virus-tumor microenvironment interactions. Additionally, incorporating AI algorithms to optimize hyperparameter tuning could enhance predictive power and replicate experimental data more efficiently.

The GLV equations with time delay best-fitting the tumor-immune-virus population dynamics are given by:

• •

  $x(t)$: Population of tumor cells (prey) at time $t$.

• •

  $y(t)$: Population of virus-infected cells (predator) at time $t$.

• •

  $a$: Intrinsic growth rate of tumor cells.

• •

  $b$: Rate at which tumor cells are killed by virus-infected cells.

• •

  $c$: Intrinsic death rate of virus-infected cells.

• •

  $d$: Rate at which virus-infected cells are produced by interactions with tumor cells.

• •

  $K_{x}$: Carrying capacity of the tumor cell population.

• •

  $K_{y}$: Carrying capacity of the virus-infected cell population.

• •

  growth_mod_x: Modifier for the growth rate of tumor cells influenced by the virus-infected cell population.

• •

  growth_mod_y: Modifier for the growth rate of virus-infected cells influenced by the tumor cell population.

• •

  damping_x: Damping coefficient for the tumor cell population, representing self-limiting growth or other inhibitory effects.

• •

  damping_y: Damping coefficient for the virus-infected cell population.

• •

  delay: Time delay in the interaction between the two populations, representing the time lag between infection of tumor cells by the virus and the resultant effect on the population sizes.

The introduction of time delay terms allows the nonlinear model to account for past states when modeling the current interaction dynamics between tumor cells, infected cells, and viruses over time(Aeedian et al., 2022; Elsadany et al., 2018; Akjouj et al., 2024; Hernández-Bermejo and Fairén, 1997). Thus, the time-delay GLV model predicts self-organizing, higher-level emergence of multiscale ecological behaviors in tumor-immune-virus interactions within heterogeneous microenvironments.

Advanced computational algorithms were used to analyze OVT dynamics datasets, treating the interaction between tumor cells (prey) and virus-infected cells (predators) using the time-delay GLV equations. The GLV model is enhanced with dynamic growth rates, adaptive learning in their competitive tumor-immune environmental interactions, damping effects, and delay terms to simulate the complex interactions observed in the complex experimental data(Aeedian et al., 2022; Akjouj et al., 2024). Other stochastic and partial differential equations including reaction-diffusion systems and general time-delay differentials were first screened before settling on the GLV as the optimal model by use of statistical performance metrics (e.g., R-square, MSE fit, etc.).

The GLV model captures the predator-prey dynamics between tumor and virus-infected cells, incorporating terms for dynamic growth rates, self-limiting behavior, and interaction delays. The delay term reflects the time lag between viral infection and the subsequent effect on tumor viability. The goal is to accurately simulate the oscillatory behaviors observed in the experimental data in representing the evolutionary dynamics of OVT interactions with the tumor ecosystems.

By optimizing the GLV model parameters using Genetic Algorithms (GA) and Differential Evolution (DE), we minimized the mean squared error (MSE) between the hybrid computational model’s predictions and the experimental data. These algorithms can dynamically adjust the damping and growth rates of the model. Reinforcement Learning (RL) is then applied to further fine-tune the model, where an RL agent iteratively adjusts the parameters to improve the model fit using the R-square fitness parameter as the reward function to optimize.

The rationale for RL among other machine learning algorithms was its ability as an open-ended adaptive optimization tool for dynamic treatment regimens and to fine-tune parameters of the GLV model for accurately predicting OVT dynamicsYang et al. (2023). The RL’s data-driven iterative learning in a delayed system adjusts parameters in response to changes in the GLV model behavior, which mimics adaptive responses in complex ecosystem interactions. Specifically, for the GLV equations, RL dynamically learns the predictive features of predator-prey interactions, accounting for time delays and non-linear feedback loops crucial in capturing oscillatory behaviors and dynamic patterns inherent in the ecosystem-based modeling of tumor-immune-virus interactions. In evidence to this rationale, Yauney and Shah (2018) demonstrated that RL agents optimize temozolomide dosing frequency by personalizing strategies to reduce tumor size and minimize toxicity, showing potential for improving oncolytic virotherapy by dynamically adjusting parameters based on patient-specific responses.

Experimental data were collected from cell viability assays and growth assays from our previous findings in (Kiyokawa et al., 2021). The cell viability data represented the population of virus-infected cells (predators), and the growth assay data represented the population of tumor cells (prey). These datasets were preprocessed by normalizing each dataset to the range $[0,1]$ using min-max scaling.

• •

  In Vitro Cell Killing Data:

  • –

    Viability Assay: Measures 005 (murine glioblastoma) cell death over 4 days, providing temporal dynamics of tumor cell viability under ICOVIR17 and combination therapies.

  • –

    Growth Assay: Captures growth inhibition at 24, 72, and 120 hours, providing early and long-term response kinetics with metrics such as mean, standard deviation, and percentiles. Thus, early responses (24–72 hours) and long-term effects (up to 11 days) were modeled to predict tumor clearance and survival outcomes.

• •

  Nanostring Immune Panel Data: Gene expression profiles of glioblastoma samples treated with ICOVIR17 and anti-PD-1 sourced from 080119 PanCancer and Myeloid_NormalizedData Kiyokawa.xlsx . Comparative analyses (ICOVIR17 vs. Control, Combo vs. ICOVIR17, Combo vs. Anti-PD-1) reveal immune regulatory networks and key mediators.

The integrated GLV hybrid model combines tumor cell dynamics and immune regulatory insights, testing synergistic, additive, or antagonistic interactions against molecular expressions from animal survival data (Nanostring panel) to predict optimal therapeutic outcomes.

Post-hyperparameter optimization, a feature importance analysis was conducted on the Nanostring gene expression data, using Random Forests (RF). This analysis identifies key predictive genes influencing the oscillatory dynamics of tumor-virus interactions, providing insights into potential therapeutic vulnerabilities or biomarkers for precision medicine.

g:Profiler was used for gene ontology and gene set enrichment analysis, identifying functional terms, pathways, and regulatory elements associated with the top 30 feature importance gene list. It uses a hypergeometric test to calculate the statistical significance of enrichment for each term. By default, g:Profiler applies the g:SCS multiple testing correction algorithm, which accounts for the hierarchical structure of gene ontologies, ensuring robustness and low false discovery.

The hybrid model effectively learned the optimal parameters to fit the dynamics of the system. The optimized parameters in Table 1 capture the predator-prey dynamics. The negative values of $a$ and $d$ represent natural declines in prey (tumor cells) and predator (virus-infected cells) populations, while negative $b$ and $c$ indicate antagonistic effects where predators infect prey, and prey depletion incurs a cost to predators. A delay of 0.3 reflects a time lag in predator-prey interactions, reflecting biological processes like immune responses or virus propagation.

The positive $K_{x}$ suggests a carrying capacity for prey, while negative $K_{y}$ implies the predator population is constrained or unstable without prey. Damping parameters ($damping\_x$, $damping\_y$) indicate how quickly oscillations in prey and predator populations stabilize, while growth modifiers ($growth\_mod\_x$, $growth\_mod\_y$) reflect the net growth trends under specific conditions. The parameters collectively minimize prediction error, modeling tumor (prey) and virus-infected cell (predator) interactions accurately. Table 2 indicates that the PPO RL agent achieves stable training with a low policy gradient loss (1.21e-10) and effective exploration-exploitation balance (entropy loss: -5.68) while maintaining efficiency at 73 fps over 8192 timesteps.

As shown in Table 3, the model demonstrates robustness by effectively capturing complex oscillatory dynamics, through empirical validation, by a low MSE (mean square error) and high R-square values for the fitness between the actual data and model predictions. Further, the use of innovative optimization algorithms, including reinforcement learning, genetic algorithms, and differential evolution, mimic real-world biological processes to optimize the timing (i.e., time-delays in tumor-virus interactions), dosage, and parameter tuning in OVT dynamics. The predictive power of these algorithms captured the oscillatory patterns in the combination therapy dynamics, as revealed in Figure 1. These methods allow adaptive adjustments to clinically relevant therapy parameters, minimizing error and improving predictive accuracy for precision medicine.

As shown in Table 4, the discovery of biomarkers, particularly related to immune response pathways like T cell activation and TNF/NF-kB, aligns with molecular signatures identified in prior studies, supporting functional validation and improving therapeutic targeting. Despite using only 10 importance features out of over 1000 genes from the Nanostring dataset, the model’s agnostic approach identified statistically significant pathways matching prior results.

As shown in Table 5 and Figure 2, the RF-predicted salience map highlighted the most enriched features and their corresponding rankings. The model predicts that photodynamic therapy synergizes with immune checkpoint inhibitors and OVT by activating comparable immune and signaling mechanisms, such as T-cell-mediated responses and TNF/NF-KB pathways. Interestingly, our findings suggest that combination therapy of oncolytic virotherapy (ICOVIR17), a hyaluronidase-expressing adenovirus, combined with anti-PD-1 immune checkpoint blockade, elicits immune responses similar to those observed with photodynamic therapy. As seen in the gProfiler gene set enrichment analysis:

Photodynamic Therapy AP-1 Pathway: The identification of ”Photodynamic therapy induced AP-1 survival signaling” (WP:WP3611) with genes such as BID, PDGFRA, and TRAF2 highlights the model’s relevance and predictive power to the therapy context. This pathway is critical for cell survival and apoptosis regulation in response to therapy progression.

Immune Responses: The top hits, including ”positive regulation of lymphocyte-mediated immunity” and ”positive regulation of leukocyte-mediated immunity,” highlight pathways involving T cell regulation and leukocytes. Genes such as TFRC (iron regulation in T cells), CD81 (immune modulation), and TRAF2 (involved in NF-kB signaling and cell survival) are essential for immune response activation and regulation.

Functional Validity: The statistical significance (adjusted p-values ranging from 0.0333 to 0.0013) confirms the robustness of these pathways in the combination therapy model. The results demonstrate how the model can extract meaningful biological insights while remaining agnostic to input data.

Other signature genes such as CD81, IL7R, IL18, PRF1, and TFRC were identified as significant predictors of T cell regulation by G: profiler, particularly influencing T cell activation, proliferation, and cytokine production. These genes likely indicate activation of cytotoxic CD8+ T cells, given their roles in immune regulation and effector function.

Positive regulation of the NF-kappaB pathway was associated with genes such as TFRC, TRAF2, PRKCI, and IL18, emphasizing their role in promoting inflammatory and survival signals in T cells. Gene signatures like BID, TFRC, TRAF2, and IFIT2 predicted involvement in apoptotic processes, providing a balance between growth and decay in these predator-prey feedback loops. Furthermore, IL18 and NF-kappaB are key to neuro-immune and brain-gut-immune axes, with IL18 driving interferon-gamma signaling critical for OVT. IFN-gamma secretion could enhance antigen presentation, and drive Th1-mediated immunity in targeted OVT. A long-short term memory network (LSTM)-based feature importance analysis on the parameter-optimized GLV model further validated these findings, highlighting contributions from TNF signaling, specifically TNFR1, interleukins like IL19 and IL17, NF-kB signaling, extracellular matrix remodeling via MMP9 metalloproteinases, and the mTOR pathway (See Supplementary Information, Table-8).

This interdisciplinary study integrates AI techniques such as RL, genetic algorithms GA, DE, and feature importance analysis via RF to forecast, optimize, and interpret complex OVT dynamics. By utilizing data-driven, explainable AI approaches, including salience maps (feature importance), the hybrid ecological model enhances predictions of OVT dynamics, for advancing precision oncology and systems medicine. These models, paired with experimental and patient-specific data, offer the ability to predict therapy responses. As shown in Figure 1, the model forecasts the complex oscillatory behaviors of multi-agent interactions and identifies critical biomarkers for personalized treatment. The model provides a good estimate of the ecological (population) dynamics, closely matching the amplitude and frequency of oscillations.

As shown in Table 4, the GLV model predicted photodynamic therapy signatures as its top 10 importance features using random forests for the ICOVIR17 combination model’s experimental data fit. The overlap of photodynamic therapy-induced AP-1 survival signaling pathways with the mechanisms identified for the combination therapy suggests shared molecular response mechanisms underlying these approaches to therapy. This highlights the role of critical dynamics between decay and growth processes, since signals of apoptosis regulation (BID, TRAF2) and cell survival or growth pathways (PDGFRA), are observed in the shared context. Such overlap indicates that targeting these pathways could optimize therapeutic efficacy across multiple treatment modalities, and fine-tune the oscillations of the predator-prey ecosystem. Further, it emphasizes their importance as predictive biomarkers and potential intervention points.

Our hybrid GLV model predictions strongly align with the findings by Kiyokawa et al. (2021), particularly the identification of T cell activation markers, TNF and NF-kB as predictive mediators. The model independently highlighted these features as key regulators of the tumor-virus-immune interactions. This agreement validates the computational model, supporting its ability to accurately identify critical regulators of the therapeutic processes and their complex dynamics in a context-agnostic manner, bridging in vitro cell dynamics and in vivo immune responses. As seen in Table 4, TRAF2 emerged as one of the top signals in salience mapping (Feature importance). TRAF2 plays a key role in TNF signaling, linking it to NF-kB activation, a central mediator of inflammation, immune responses, and cell survival (Table 5) (Borghi et al., 2016). BID (apoptosis regulator) and PDGFRA (growth signaling) complement these processes by balancing cell death and proliferation, highlighting their relevance in OVT’s oscillatory dynamics. Thus, these gene expression signatures themselves support the oscillatory dynamics of the model, wherein the OVT-immune-cancer interactions are in a predator-prey feedback loop, causing the tumor to adaptively oscillate between growth and decay, toward a ’stable attractor’ (e.g., a limit cycle) during the therapy progression. Thus, an intuitive and explainable model like the hybrid GLV by capturing predator-prey dynamics, reflects the interplay of immune activation and suppression among viral and tumor interactions, as shown by its alignment with the Nanostring markers, enhancing biological interpretability.

None of the current state-of-the-art models have achieved the integration of predictive biomarkers as features with a strong statistical fit, demonstrated by the presented model’s residual variance between the predicted and actual values (i.e., low MSE of 0.02) and an R² value of 0.83 and 0.87 for the prey and predator, respectively (Table 3). This highlights the hybrid GLV model’s ability to accurately capture key biological dynamics while maintaining predictive robustness as an explainable AI system, marking a significant advancement in modeling therapy responses in computational medicine.

By integrating multimodal and real-time (longitudinal) data into these predictive models, we enabled dynamic fine-tuning of OVT dosing regimens, potentially leading to optimized therapeutic efficacy, reduced side effects, and improved patient care. The multimodal features can be further enhanced in the prospective model scaling with longitudinal tumor-immune markers (such as from multiomics, and liquid biopsies or blood profiles), to adaptively adjust therapy parameters. By integrating continuous feedback mechanisms, clinicians can use model outputs to customize treatment plans dynamically, accounting for each patient’s unique tumor microenvironment. Furthermore, the revealed molecular patterns from the gene expression studies can be extended to time-resolved, longitudinal RNA-seq and multiomics in future studies. The identified patterns and gene signatures provide potential targets for therapeutic vulnerabilities and biomarkers of OVT progression.

Given the suitability of the GLV framework for longitudinal data forecasting, it could benefit from time-delay embedding attractor reconstruction and recurrent neural networks (RNNs) such as long short-term memory networks (LSTMs) to predict patterns and temporal dependencies(Huang et al., 2020). These approaches could refine predictions and enable the identification of subtle patterns and multimodal features in vivo or patient-specific responses to OVT, with more time-resolved data. In evidence, we used LSTM as a feature importance mapping method akin to RF and found that it validated our RF-based salient features, intersecting with markers regulating T-helper cell differentiation, extracellular matrix remodeling factors, and most signals being part of the TNF and NF-kB pathways, as shown in the supplementary information in in Tables 7 and 8.

This study does not fully explore how the observed oscillatory patterns and complex dynamics could directly influence tumor-specific outcomes in optimizing combination OVT. A deeper analysis of how these dynamics correlate with clinical metrics, such as in vivo model tumor progression, immune response, and viral persistence, would enhance the translational relevance of the preclinical model in paving dynamic treatment regimens.

However, the model incorporates temporal features adaptable for personalized therapies, such as (patient-specific) time delays, and carrying capacities, which could further strengthen the predictive power for precision oncology. Various dynamical systems approaches can be exploited using these temporal features to enhance the model (McGehee et al., 2008; Crutchfield et al., 1982). For instance, attractor reconstruction uses delay-embedded trajectories to analyze the long-term behavior and stability of systems. By embedding time-delay dynamics into phase space, it becomes possible to visualize attractor geometries steering the ecological dynamics, and identify limit cycles or their critical transitions to strange attractors(Shin and Cho, 2023). Thereby, attractor reconstruction can help detect the stability and bifurcations of these oscillatory systems(Evers et al., 2024; Chisholm and Filotas, 2009; McGehee et al., 2008; Crutchfield et al., 1982). This enables enhanced therapeutic predictions and real-time monitoring. For instance, the time-delay parameters and growth rates predicted by the model provide explainable features in the temporal dynamics of tumor-virus interactions, such as the lag between viral infection and tumor response(Shin and Cho, 2023). These parameters can guide the timing and frequency of viral dosing or combination therapy, ensuring interventions align with the predicted peak effectiveness. By tailoring these inputs to the patient-specific tumor growth-decay patterns, clinicians can potentially optimize therapeutic efficacy while minimizing side effects.

This integrated hybrid model combines GLV equations with AI-driven optimization techniques, including RL, GA, DE, and feature importance analysis via random forest (RF), to model OVT dynamics, fine-tune parameters, and identify critical genes predicting complex tumor-immune interactions. The predator-prey oscillatory dynamics captured by the model could help inform adaptive treatment schedules and personalized strategies by adjusting therapy timing (frequency) with tumor behavior, immune response, and viral interactions, while patient-specific data can further optimize therapy. In specific, our hybrid GLV model forecasts T cell activation markers and cytokine signals, such as TNF and NF-kB pathways, identified by Kiyokawa et al. (2021), as key mediators of OVT-immune-tumor dynamics.

Biosignatures, including CD81, TRAF2, IL18, and BID as predicted by RF salience maps, and TNF, IL4, MMP9, MMP12, among other tumor microenvironment remodeling pathways identified by LSTM feature extraction (See Supplementary Information; Table 7), can guide personalized cancer immunotherapies by targeting immune checkpoint regulators, modulating T-cell differentiation, and influencing apoptosis. These markers also suggest targeting brain-immune axis inflammatory signals, supporting the development of immunomodulatory drugs. Thus, our proof-of-concept model provides a promising platform for real-time, longitudinal monitoring, adaptive control, and optimization of combination OVT in precision medicine, with predictive biomarkers informing novel targeted therapies to improve quality patient care and outcomes.

The Long Short-Term Memory (LSTM) network extracts feature importance by its ability to learn temporal dependencies in the data. Gene expression data from the Nanostring panel was preprocessed by transposing the dataset to organize genes as columns and aggregating it to match the length of the target model output with an 80:20 training-testing split. The LSTM model had 2 layers (one LSTM layer with 50 units and one Dense layer) and used the Adam optimizer (learning rate = 0.001), with early stopping (patience = 10) and a batch size of 32. The LSTM was trained on the data to predict features of the optimized GLV model. The absolute mean of the learned weights in the LSTM layer was calculated as a proxy for feature importance.

The LSTM salience map is a potentially robust approach for feature importance due to its ability to capture long-range patterns in time-series data. However, given our small sample size, Random Forest (RF) feature predictions were prioritized for the main results, while LSTM insights still demonstrated predictive power in identifying functionally validated immune-inflammatory regulatory signals and ECM remodeling pathways. We propose their scalability and translatability to preclinical in vivo studies with larger datasets.