A Minimal Model Framework for Robust CAR-T Cell and Oncolytic Virus Combination Therapy

TL;DR

This paper introduces a simplified mathematical model to predict the effectiveness of combining CAR-T cells and oncolytic viruses for treating glioblastoma.

Contribution

The study presents a minimal model framework using quasi-steady-state approximations to optimize combination immunotherapy.

Findings

The QSS model uses 9 parameters and achieves comparable accuracy to the full model with 11 parameters.

The QSS model is favored by the Akaike Information Criterion in most therapy conditions.

CAR-T exhaustion dynamics do not improve predictions within the 72-hour observation window.

Abstract

Glioblastoma remains one of the most lethal brain cancers. Combination therapy using CAR-T cells and oncolytic viruses shows promise, yet mechanisms underlying synergy remain poorly understood. We develop mathematical models to analyze IL-13Rα2-targeting CAR-T cells and the oncolytic virus C134 using patient-derived glioblastoma data. We present a minimal model framework for predicting combination immunotherapy outcomes. Applying timescale separation between rapid viral and slower cellular dynamics, we derive quasi-steady-state (QSS) approximations that reduce complexity while maintaining accuracy. The QSS model uses 9 parameters compared with 11 in the full model and achieves comparable fits. Model comparisons using the Akaike Information Criterion indicate that the QSS model is generally favored; it consistently yields lower AIC values for oncolytic virus monotherapy and produces lower AIC values in three of four combination therapy conditions. Models with and without CAR-T exhaustion produce identical fits, indicating that exhaustion dynamics do not improve predictions within the 72-hour observation window. Overall, our results demonstrate that simplified QSS formulations effectively capture viral dynamics and provide a practical framework for optimizing combination immunotherapies.

Figures

Equations

Full text

Introduction

1

Glioblastoma multiforme (GBM) is the most common and aggressive brain tumor in adults, characterized by rapid growth, diffuse infiltration, and high recurrence rates despite standard treatments such as surgical resection, radiotherapy, and temozolomide chemotherapy [1, 2]. Survival improvements over recent decades have been limited, reflecting the highly heterogeneous and treatment-resistant nature of GBM [3]. These challenges have motivated the development of diverse therapeutic strategies, including immunotherapies and locoregional approaches, aimed at eliminating tumor cells, stimulating anti-tumor immunity, and overcoming barriers such as the blood-brain barrier and intrinsic resistance [4].

Among these emerging approaches, immunotherapies have attracted particular attention, with Chimeric Antigen Receptor (CAR) T-cell therapy emerging as a promising strategy. CAR-T therapy involves engineering a patient’s T cells to recognize and selectively eliminate tumor cells [5]. In the context of GBM, CAR-T cells targeting IL-13Rα2 have demonstrated tumor regression [6]. However, their effectiveness in solid tumors is still limited by poor T-cell infiltration, an immunosuppressive tumor microenvironment, and antigen heterogeneity [7].

To address some of these limitations, oncolytic virus (OV) therapy has been explored as a complementary strategy. The concept of using viruses to treat cancer dates back to early 20th-century observations, when viral infections were occasionally associated with spontaneous tumor regression [8]. Modern OVs exploit this phenomenon through a dual mechanism: selective infection and lysis of malignant cells, coupled with stimulation of anti-tumor immune responses. Viral replication within tumor cells triggers cell death and the release of new virions, which spread to neighboring cancer cells. Importantly, this process induces immunogenic cell death, releasing tumor antigens and danger signals that activate the immune system and convert immunologically “cold” tumors into “hot” ones. When combined with chemotherapy or other immunotherapies, OV therapy can achieve synergistic effects, enhancing efficacy while minimizing toxicity [9–11].

Since CAR-T and OV therapies target tumors via distinct but complementary mechanisms, their combination has emerged as a promising strategy. OVs can remodel the tumor microenvironment to support CAR-T function by reducing immunosuppression, enhancing antigen presentation, and promoting chemokine-mediated recruitment of effector cells [12]. Conversely, CAR-T cells can target virus-infected tumor cells, potentially increasing viral spread and overall treatment efficacy. Understanding how these therapies interact and optimizing their combination requires quantitative frameworks that can capture the underlying dynamics.

Mathematical modeling has proven to be a valuable tool for understanding the dynamics and predicting outcomes of complex cancer immunotherapies [13–15]. Mechanistic ordinary differential equation (ODE) models of CAR-T cell therapy for glioma have, for instance, investigated how multiple CAR-T cells bind to individual tumor cells and how parameters such as expansion rates and antigen receptor densities influence treatment success and tumor elimination kinetics [16]. Adhikarla et al. [17] employed ODE modeling to optimize the scheduling of CAR-T therapy in combination with targeted radionuclide therapy, showing that the timing and sequencing of treatments critically affect therapeutic outcomes. Similarly, Conte et al. [18] introduced a kinetic-theory–based approach to tumor–immune modeling, deriving macroscopic dynamics from microscopic cellular interactions and highlighting how spatial structure and immune cell state transitions shape tumor-immune behavior. Mahasa et al. [27] investigated dosing strategies for CAR-T and OV combination therapy, revealing that virus-induced synergism can prevent tumor elimination and that simultaneous administration induces more sustained tumor cell reduction compared to sequential treatment. More recently, Pell [19] developed an ODE framework demonstrating that immunostimulants such as NHS-muIL12 can both boost immune activity and induce tumor PD-L1 upregulation, thereby reducing the efficacy of PD-L1 checkpoint blockade and explaining treatment failure. Building on these foundations, Conte et al. [20] developed an early data-driven combination-treatment model integrating uninfected and infected tumor cells, viral replication, and CAR-T dynamics in the context of GBM. Inspired by this approach and the foundational work of Pell et al. [19], we propose an alternative minimal model for analyzing CAR-T cell and oncolytic virus combination therapy. The model tracks four key populations: uninfected and infected tumor cells, viral particles, and CAR-T cells, including explicit dynamics for CAR-T cell exhaustion. To examine the balance between biological detail and model tractability, we also consider a reduced formulation in which exhaustion effects are represented implicitly through modified CAR-T cell death rates over a 72-hour timescale. In addition, we apply timescale separation to distinguish fast viral dynamics from the slower evolution of cellular populations, allowing for quasi-steady-state approximations.

The model is calibrated to experimental data from Conte et al. [20] across multiple effect-or-to-target ratios and multiplicities of infection. Using this framework, we analyze interactions between CAR-T and oncolytic virus therapies, evaluate the impact of model simplifications, and identify mechanisms that are essential for reproducing key experimental trends relevant to combination treatment.

Sections 2 and 3 describe the experimental data, present the full mathematical model, introduce reduced and quasi-steady-state formulations, and apply these models to combination therapy and monotherapy scenarios using data fitting. Section 4 presents quantitative model comparisons across formulations, including performance metrics and identifiability analysis. Section 5 discusses biological interpretation and implications for therapeutic optimization.

Main Model

2

Experimental Data

2.1

The experimental data used to validate our mathematical model are obtained from recent investigations by Conte et al., who systematically examined combination therapy involving oncolytic virus C134 and IL-13Rα2-targeting CAR-T cells against patient-derived glioblastoma cultures [20]. C134 is a herpes simplex virus (HSV)-based oncolytic virus that selectively infects and lyses tumor cells, utilizing cancer-specific defects in antiviral defenses. Upon infection, C134 replicates within tumor cells, triggering cell lysis and releasing new viral particles for subsequent infection cycles. The CAR-T cells were engineered to target IL-13Rα2, a receptor overexpressed in glioblastomas but with limited expression in normal brain tissue, thereby enabling selective recognition and elimination of tumor cells.

Experiments employed PBT030 primary glioblastoma cultures with high IL-13Rα2 expression. Cells were plated at 2 × 10^4^ per well with three biological replicates ( $[eqn]$ ) per condition. Cultures were initially infected with C134 at multiplicities of infection ( $[eqn]$ ), representing different viral particle-to-cell ratios. Two hours later, CAR-T cells were introduced at effector-to-target (E:T) ratios of 1:10, 1:25, or 1:50, corresponding to 2000, 800, and 400 cells per well, allowing the virus to establish infection prior to CAR-T cell engagement.

Tumor dynamics were monitored using the xCELLigence system, which continuously measures electrical impedance as cells adhere to electrode-coated surfaces. The resulting Cell Index (CI), a dimensionless proxy for viable cell counts, was recorded every 15 minutes over 72 hours, yielding high-resolution temporal profiles (> 280 data points per condition). These measurements were used to estimate model parameters for both monotherapy and combination therapy.

Model Formulation

2.2

We develop a mathematical model to describe the tumor-immune-virus dynamics in the experiments described in Section 2.1. The model tracks four populations: uninfected tumor cells $[eqn]$ , infected tumor cells $[eqn]$ , oncolytic virus particles $[eqn]$ , and CAR-T cells $[eqn]$ . We begin with a detailed formulation incorporating the main biological mechanisms, then develop reduced models motivated by timescale separation. The dynamics of each population are determined by:

Uninfected Tumor Cells: U(t).

Uninfected tumor cells follow logistic growth dynamics, expressed as $[eqn]$ , where $[eqn]$ denotes the intrinsic growth rate and $[eqn]$ represents the environmental carrying capacity accounting for space and nutrient limitations. The tumor population is reduced through two mechanisms: infection by the oncolytic virus, represented by $[eqn]$ , where $[eqn]$ denotes the infection rate coefficient and the saturation term $[eqn]$ reflects crowding effects and limited receptor availability; and elimination by CAR-T cells through antigen recognition and cytotoxic mechanisms, represented by $[eqn]$ , where $[eqn]$ denotes the killing rate coefficient. The unified saturation constant $[eqn]$ represents a threshold reflecting the finite capacity of cellular interactions at high densities, accounting for both spatial crowding and shared metabolic constraints. The dynamics of uninfected tumor cells are therefore described by:

$$[eqn]$$

Infected Tumor Cells: I(t).

Infected tumor cells are generated when uninfected tumor cells become infected by the oncolytic virus at rate $[eqn]$ . Once infected, tumor cells lose their proliferative capacity and do not undergo further division. The infected cell population is reduced through two processes: CAR-T mediated killing, represented by $[eqn]$ , where $[eqn]$ denotes the killing rate coefficient; and viral lysis at rate $[eqn]$ , which releases new viral particles into the tumor microenvironment. The dynamics of infected tumor cells are given by:

$$[eqn]$$

Oncolytic Virus: V(t).

Viral particles are produced upon lysis of infected tumor cells, with burst size $[eqn]$ contributing to viral amplification. This production is represented by $[eqn]$ . The viral population decreases through two mechanisms: consumption during infection of uninfected tumor cells at a rate $[eqn]$ , and natural clearance through degradation and non-productive binding at a rate $[eqn]$ . Combining these mechanisms, the oncolytic virus changes as:

$$[eqn]$$

CAR-T Cells: C(t).

CAR-T cells expand in response to tumor antigen stimulation, modeled within a Michaelis-Menten framework as $[eqn]$ , where $[eqn]$ denotes the proliferation rate induced by uninfected tumor cells. The CAR-T population is subject to multiple loss mechanisms: natural apoptosis at a rate $[eqn]$ ; and exhaustion due to prolonged antigen exposure, represented by $[eqn]$ , where $[eqn]$ denotes the exhaustion rate coefficient. The dynamics of CAR-T cells are represented by:

$$[eqn]$$

Collecting the above components, the complete model with CAR-T exhaustion is given by:

$$[eqn]$$

We summarize the model parameters with their biological interpretation and units in Table 1.

Parameter Estimation

2.3

We estimated model parameters by fitting predictions to experimental Cell Index measurements. The optimization employed the L-BFGS-B algorithm [21, 22] (limited-memory quasi-Newton method for bound-constrained nonlinear optimization) implemented in SciPy [23]. For each experimental condition and biological replicate, we minimized the Weighted Absolute Percentage Error (WAPE):

$$[eqn]$$

where $[eqn]$ denotes the measured Cell Index and $[eqn]$ the model prediction at time $[eqn]$ .

Model performance was evaluated using three standard metrics: the Residual Sum of Squares (RSS), the adjusted $[eqn]$ coefficient, and the Akaike Information Criterion (AIC). The RSS, defined as

$$[eqn]$$

represents the squared deviations between model predictions and experimental data. While WAPE was minimized during parameter fitting, RSS provides a standard measure of deviation, with lower values indicating better agreement between model and data.

The adjusted $[eqn]$ accounts for both the number of fitted parameters and the number of observations, and was calculated as

$$[eqn]$$

where $[eqn]$ is the number of observations, $[eqn]$ the number of fitted parameters, and $[eqn]$ the total sum of squares. A higher adjusted $[eqn]$ indicates a better model fit while accounting for model complexity.

The AIC penalizes model complexity relative to fit quality and was defined as

$$[eqn]$$

A lower AIC indicates a better trade-off between goodness of fit and model simplicity.

Model Iterations and Comparisons

3

We use computational modeling to explore CAR-T cell and oncolytic virus treatment of glioblastoma and assess therapeutic synergy under combination therapy. We first analyze combination therapy, then monotherapies, and finally quasi-steady-state approximations for each scenario.

For each therapeutic condition, we present the governing equations and fit our model to experimental Cell Index measurements across all dosing conditions and biological replicates. Unless stated otherwise, all model parameters are estimated via nonlinear optimization, as detailed in Section 2.3. Baseline experimental conditions employed PBT030 glioblastoma cells, with therapeutic agents administered according to the protocols described in Section 2.1.

Combination Therapy Dynamics

3.1

Full Model with CAR-T Exhaustion

3.1.1

The combination therapy case is described by model Eq. (5), which examines the effects of varying CAR-T cell dose (E:T ratios of 1:25 and 1:50) and viral load (MOI of 0.002 and 0.0008) on tumor dynamics. Tumor response is analyzed across four dose combinations, with three biological replicates per condition.

Data-fitting results for biological replicate Y1 appear in Fig. 1, with additional replicates (Y2 and Y3) in Appendix A1. Experimental Cell Index measurements (black lines) are overlaid with model predictions (red lines). The model resolves associated population dynamics, including uninfected tumor cells (blue dashed), infected tumor cells (green dashed), and CAR-T cells (yellow dashed).

Sustained tumor suppression is achieved within 24 hours across all dose combinations and maintained throughout the 72-hour observation period. The dynamics highlight nonlinear interactions between viral oncolysis, CAR-T cytotoxicity, and exhaustion, particularly in how the two therapies differentially kill infected versus uninfected tumor cells. Parameter estimates are in Table A1.

Reduced Model Without CAR-T Exhaustion

3.1.2

CAR-T exhaustion involves progressive transcriptional changes requiring sustained antigen exposure over days to weeks. Within the 72-hour experimental window, the primary loss mechanisms are direct cytotoxic killing and activation-induced cell death, both captured by the natural death term $[eqn]$ . To test whether explicit exhaustion dynamics are necessary for accurate short-term predictions, we removed the exhaustion term $[eqn]$ from the CAR-T equation in Eq. 5:

$$[eqn]$$

Model fits for biological replicate Y1 across the same four dose combinations are shown in Fig. 2, with additional replicates (Y2 and Y3) in Appendix A2. Tumor suppression within 24 hours occurs across all conditions. The reduced model captures the same tumor suppression dynamics observed in the full model, suggesting that exhaustion occurs over longer timescales than the 72-hour experimental window. Parameter estimates are in Table A2.

Monotherapy Analysis

3.2

Oncolytic Virus Monotherapy

3.2.1

For oncolytic virus monotherapy, the full model Eq. (5) explicitly tracks viral dynamics through three coupled differential equations:

$$[eqn]$$

Here $[eqn]$ denotes uninfected tumor cells, $[eqn]$ infected tumor cells, and $[eqn]$ free viral particles. The infection term incorporates Michaelis–Menten saturation kinetics through the parameter $[eqn]$ . Viral dynamics are governed by production via lysis at rate $[eqn]$ with burst size $[eqn]$ , consumption during infection at rate $[eqn]$ , and natural clearance at rate $[eqn]$ , allowing mechanistic tracking of each process.

We evaluated oncolytic virus monotherapy at three MOI levels (0.03, 0.002, 0.0008) with three biological replicates per condition. Data-fitting results for biological replicate Y1 are in Fig. 3, with additional replicates (Y2 and Y3) in Appendix A3. The model captures viral oncolysis across three orders of magnitude in viral dose. Parameter estimates are in Table A3.

CAR-T Cell Monotherapy

3.2.2

In the absence of oncolytic virus, the system in Eq. (5) reduces to a two-compartment model describing the interaction between uninfected tumor cells $[eqn]$ and CAR-T cells $[eqn]$ . The full model incorporating exhaustion dynamics is

$$[eqn]$$

To assess the role of exhaustion over the experimental timescale, we also consider a reduced formulation obtained by setting $[eqn]$ , which removes the exhaustion term:

$$[eqn]$$

Both formulations include logistic tumor growth and CAR-T-mediated killing. CAR-T cell dynamics account for antigen-dependent expansion and natural death, with the full model additionally capturing functional decline due to sustained antigen exposure. The parameter $[eqn]$ represents the half-saturation constant governing both killing efficiency and proliferative response.

Both formulations were evaluated across three effector-to-target ratios of 1:10, 1:25, and 1:50, corresponding to high, intermediate, and low CAR-T cell concentrations relative to the initial tumor burden. Model fits for biological replicate Y1 are shown in Figs. 4 and 5. Parameter estimates are reported in Tables A4 and A5. Results for additional replicates Y2 and Y3 are provided in Appendices A4 and A5.

Quasi-Steady-State Approximations

3.3

Viral replication and clearance occur on faster timescales than cellular processes. Oncolytic viruses exhibit rapid burst kinetics and clearance within hours, while cellular proliferation, death, and phenotypic transitions unfold over days. This separation of timescales allows a quasi-steady-state approximation, where viral populations equilibrate rapidly relative to cellular dynamics.

Combination Therapy Under QSS

3.3.1

Setting $[eqn]$ in the full model (5) and solving for viral concentration yields:

$$[eqn]$$

When viral clearance dominates reinfection $[eqn]$ , this simplifies to:

$$[eqn]$$

where $[eqn]$ represents the effective viral amplification factor. Substituting this expression into the full model reduces it to three dimensions:

$$[eqn]$$

QSS model fits for biological replicate Y1 across varying CAR-T cell doses (E:T ratios of 1:25 and 1:50) and viral loads (MOI of 0.002 and 0.0008) are shown in Fig. 6, with additional replicates (Y2 and Y3) in Appendix A6. Parameter estimates are in Table A6.

Oncolytic Virus Monotherapy Under QSS

3.3.2

For oncolytic virus therapy without CAR-T cells, we set $[eqn]$ in Eq. (16) and remove all CAR-T-related terms. The resulting susceptible-infected system describes uninfected and virus-infected tumor cell dynamics:

$$[eqn]$$

QSS model fits for biological replicate Y1 across three MOI levels (0.03, 0.002, 0.0008) are shown in Fig. 7, with additional replicates (Y2 and Y3) in Appendix A7. Parameter estimates are in Table A7.

For CAR-T cell monotherapy, the reduced model without exhaustion (Eq. 13) serves as the quasi-steady-state (QSS) representation. Since this system is identical to the original reduced model, the predicted dynamics and parameter estimates are unchanged.

CAR-T Cell Monotherapy Under QSS

3.3.3

For CAR-T cell monotherapy, the reduced model without exhaustion (Eq. 13) serves as the quasi-steady-state representation. Since no viral dynamics are present, this system is identical to the reduced model, and the predicted dynamics and parameter estimates remain unchanged.

Results

4

Performance Metrics for Model Comparison

4.1

We compared three model formulations: the Full model (F) with exhaustion dynamics Eq. (5), the model without exhaustion (W) Eq. (10), and the quasi-steady-state approximation (QSS) Eq. (16). These formulations were evaluated across all experimental conditions, including both monotherapy and combination treatment settings.

Table 2 compares two model formulations for oncolytic virus monotherapy across three viral doses. Since no CAR-T cells are present in these experiments, the Full model (F) Eq. (11) is equivalent to the model without exhaustion (W), and both include explicit viral dynamics with six parameters. The QSS approximation uses one fewer parameter Eq. (17) ( $[eqn]$ ), and produces fits that are similar to or slightly better than those of the full model Eq. (11) across all MOI levels. At the lowest viral dose (MOI = 0.0008), the QSS model gave RSS = 3.06 ± 1.45 compared to 6.12 ± 1.62 for the full model, with corresponding AIC values of −1223 ± 144 versus −1017 ± 76. At the highest dose (MOI = 0.03), both models fit the data well (adjusted $[eqn]$ ), with the full model giving RSS = 1.05 ± 0.27 and AIC = −1486 ± 75, while the QSS model showed slightly lower RSS and AIC (RSS = 0.57 ± 0.37, AIC = −1683 ± 166). These results suggest that the QSS approximation captures the main effects of viral dynamics while simplifying the model.

Table 3 compares model formulations for CAR-T monotherapy across three effector-to-target ratios. In the absence of an oncolytic virus, the model without exhaustion (W) Eq. (13) is equivalent to the QSS approximation. Both the Full model with exhaustion Eq. (12) dynamics and the W=QSS formulation Eq. (13) yielded nearly identical fits across all E:T ratios. At E:T = 1:10, both models gave RSS = 0.354 ± 0.136 and AIC values around −1800, with adjusted $[eqn]$ . This pattern was similar at E:T = 1:25 (RSS = 0.495 ± 0.296, AIC ≈ −1729) and E:T = 1:50 (RSS ≈ 0.58, AIC ≈ −1640). The exhaustion parameter showed high uncertainty in all conditions (ranging from $[eqn]$ at E:T = 1:25 to 6.555 ± 1.419 at E:T = 1:50), and including it did not change the model fit. Overall, the simpler W=QSS formulation with five parameters captures CAR-T monotherapy dynamics within the 72-hour experimental window.

Table 4 compares the Full (F) Eq. (5) and quasi-steady-state (QSS) Eq. (16) model formulations across four combination therapy conditions varying in effector-to-target ratio and viral dose. Both models provided good fits to the data (adjusted $[eqn]$ ). The QSS model employs fewer parameters ( $[eqn]$ versus 11 for F) and generally yields lower RSS and AIC values. At E:T = 1:25 with MOI = 0.002, the QSS model had RSS = 0.486±0.258 and AIC = −1697±133. At the same E:T ratio but lower viral dose (MOI = 0.0008), both models performed similarly. At E:T = 1:50 with MOI = 0.002, the QSS model gave RSS = 0.387 ± 0.150 and AIC = −1832 ± 105. At E:T = 1:50 with MOI = 0.0008, both models performed similarly, with the QSS model slightly better. Including the model without exhaustion (W) in Table A8 shows that W is competitive at E:T = 1:25 with MOI = 0.0008 (RSS = 0.517±0.219, AIC = −1629±116), but overall differences between the three formulations are small. The QSS model provides a good balance of fit and parameter simplicity.

These results show two main trends across the experimental conditions. First, viral dynamics occur on a faster timescale than tumor growth. In OV monotherapy, the QSS formulation Eq. (17) generally fits the data as well as or better than the full model Eq. (17), reducing RSS by roughly two-thirds at intermediate MOI levels while using one fewer parameter. This indicates that the explicit differential equations governing viral production and clearance can be replaced by a quasi-steady-state approximation, assuming rapid equilibration of viral dynamics, without reducing predictive accuracy. Second, exhaustion dynamics do not noticeably affect therapeutic outcomes in this system. In CAR-T monotherapy, models with Eq. (12) and without exhaustion Eq. (13) gave nearly identical fits across all E:T ratios, and including the exhaustion parameter did not improve the model. In combination therapy, these patterns remain: the QSS formulation Eq. (16) performs similarly to the full model Eq. (5) using 9 parameters instead of 11, and the model without exhaustion Eq. (10) gives comparable results (Table A8). Overall, the QSS model reduces model complexity and focuses on the slower dynamics that dominate the experimental timescale while keeping the system interpretable.

Parameter Estimation Across Model Formulations

4.2

Oncolytic Virus Monotherapy

4.2.1

We compared the Full model Eq. (5), which includes explicit viral dynamics with six parameters, to the quasi-steady-state (QSS) Eq. (17) approximation with five parameters for OV monotherapy (Table A9). In the absence of CAR-T cells, the Full model is equivalent to the model without exhaustion, as the exhaustion term applies only to CAR-T dynamics; these formulations are denoted F=W in Table A9.

Tumor growth rates were consistent across models and viral doses, ranging from $[eqn]$ to 0.060 day^−1^. The saturation threshold $[eqn]$ converged to approximately 0.10 in both models at low and intermediate MOI, with minimal variability.

The viral infection rate $[eqn]$ showed model-dependent behavior. In the Full model Eq. (5),, $[eqn]$ ranged from 1.03 to 1.54 and increased with viral dose, consistent with mass-action kinetics. In the QSS model Eq. (17), estimates reached the upper bound of 1.50 at low and intermediate MOI, indicating identifiability issues when viral dynamics are not modeled explicitly. The infected cell lysis rate $[eqn]$ differed between formulations at low and intermediate MOI: the QSS model produced values two- to four-fold higher than the Full model (0.142–0.168 vs. 0.035–0.078 day^−1^). At high MOI, estimates converged (QSS: 0.085, Full: 0.074 day^−1^), reflecting that the QSS formulation combines viral production and clearance into an effective lysis rate, whereas the Full model treats these processes separately.

Parameters unique to the Full model Eq. (5), also exhibited identifiability challenges. The burst size $[eqn]$ varied from 13 to 24 virions per cell without a clear dose-dependent trend. The viral clearance rate $[eqn]$ reached its upper bound (~2.0 day^−1^), corresponding to a half-life of ~8 hours, indicating that decay occurs faster than the measurement interval can capture.

In contrast, the QSS parameter $[eqn]$ provides a simpler alternative. It increased with viral dose, from 0.178±0.033 at MOI = 0.0008 to 0.746±0.523 at MOI = 0.03, effectively capturing the net oncolytic effect without requiring explicit estimation of fast viral kinetics.

CAR-T Cell Monotherapy

4.2.2

We compared the Full model Eq. (12), which incorporates exhaustion dynamics with six parameters, to the model without exhaustion Eq. (13), which uses five parameters, for CAR-T monotherapy (Table A10). Since no oncolytic virus is present, the quasi-steady-state approximation reduces to the model without exhaustion, and only two formulations were compared for this therapeutic modality.

Tumor growth rates were consistent between formulations, ranging from $[eqn]$ to 0.254 day^−1^, corresponding to doubling times of 2.7–3.1 days. CAR-T-related parameters displayed dose-dependent trends. The killing rate of uninfected tumor cells, $[eqn]$ , decreased from 10.0 day^−1^ at E:T = 1:10 to 8.0 day^−1^ at E:T = 1:50, likely reflecting reduced per-cell cytotoxic efficiency when effectors engage multiple targets sequentially. This accounts for the time required for immunological synapse formation, granule release, and detachment before the next engagement. The saturation threshold $[eqn]$ showed a similar pattern, decreasing from 45.2 to 13.8 across the same E:T range, suggesting killing efficiency saturates at lower tumor densities when fewer effectors are available. Conversely, the natural death rate $[eqn]$ increased from 0.016 day^−1^ at E:T = 1:10 to 0.117 day^−1^ at E:T = 1:50, possibly reflecting enhanced CAR-T survival at higher effector densities due to abundant antigen stimulation and autocrine survival signals. These parameters were identical between model formulations.

The proliferation rate $[eqn]$ differed substantially between models. In the Full model Eq. (12), $[eqn]$ day^−1^, whereas in the model without exhaustion, $[eqn]$ day^−1^. This discrepancy arises because the exhaustion term reduces CAR-T function, forcing the optimizer to increase proliferation to match observed tumor killing. Proliferation rates from the model without exhaustion align with known T cell division times of 12–24 hours, whereas Full model Eq. (12) estimates exceed physiological rates. The exhaustion parameter, $[eqn]$ , showed high uncertainty and did not improve model fit, confirming it cannot be reliably estimated from tumor growth data over this timescale.

Combination Therapy

4.2.3

We compared all three formulations for combination therapy: Full model with exhaustion Eq. (5) using eleven parameters, model without exhaustion Eq. (10) using ten parameters, and quasi-steady-state approximation using nine parameters Eq. (16). Parameter estimates are shown in Table A11.

The tumor growth rate of $[eqn]$ day^−1^ and infected cell lysis rate of $[eqn]$ day^−1^ were consistent across all three formulations, confirming these parameters can be reliably estimated. The CAR-T proliferation rate of $[eqn]$ day ^−1^ showed better consistency than in monotherapy, with all estimates in the plausible range for activated T cells.

Other parameters varied considerably across formulations. The CAR-T killing rate of uninfected cells $[eqn]$ differed up to 21-fold within identical conditions. At E:T = 1:50 with MOI = 0.002, the Full model estimated $[eqn]$ day^−1^, the model without exhaustion returned $[eqn]$ day^−1^, and the quasi-steady-state model predicted $[eqn]$ day^−1^. This variability suggests that CAR-T killing effects, including direct cytotoxicity and their interaction with viral dynamics, cannot be separated based on tumor-level data alone. The killing rate of infected cells $[eqn]$ and saturation threshold $[eqn]$ showed similar patterns.

The CAR-T natural death rate $[eqn]$ varied across formulations. The quasi-steady-state model consistently gave higher estimates of 0.876–0.976 day^−1^ than the Full model of 0.143–0.497 day^−1^ or the model without exhaustion of 0.137–1.029 day^−1^. The viral infection rate of $[eqn]$ also varied without clear dose-dependent trends.

Viral dynamics parameters showed the same identifiability issues seen in monotherapy. The burst size of $[eqn]$ virions/cell and clearance rate of $[eqn]$ day^−1^ varied widely. The exhaustion parameter of $[eqn]$ was lower than in CAR-T monotherapy but remained poorly constrained.

The quasi-steady-state parameter $[eqn]$ varied across conditions—stable in some with E:T = 1:50, MOI = 0.002: $[eqn]$ but variable in others with E:T = 1:25, MOI = 0.002: $[eqn]$ —suggesting that interactions between viral and CAR-T effects are absorbed differently depending on the therapeutic balance.

Taken together, these results show that the quasi-steady-state formulation achieves comparable or more stable estimates for key parameters of $[eqn]$ while avoiding poorly constrained terms of $[eqn]$ , supporting the use of simpler models for reliable biological inference.

Discussion and Conclusion

5

Our results demonstrate that minimal model formulations perform comparably to mechanistically detailed alternatives for predicting combination immunotherapy outcomes. The quasi-steady-state (QSS) model Eq. (16) uses nine parameters in combination therapy compared to eleven in the Full model Eq. (5), yet achieves comparable fits. This makes it better suited for applications relying on parameter estimates from tumor data, including treatment optimization, dose scheduling, and outcome prediction.

Incorporating distinct CAR-T exhaustion states did not improve model fit within the 72-hour observation period, reflecting the distinct temporal scales governing T cell dysfunction. Phenotypic exhaustion develops over 2–4 weeks through extensive epigenetic and transcriptional reprogramming [24]. Within a 72-hour window, CAR-T cells have not yet undergone this commitment; instead, functional decline is driven by activation-induced cell death (AICD), which occurs within hours to days. The minimal formulation with effective death rates captured the dominant processes observable within this temporal window without requiring additional exhaustion parameters.

Similar identifiability issues arise for viral dynamics. The viral clearance rate consistently reached parameter bounds, indicating that decay occurs faster than can be resolved by daily measurements. Burst size $[eqn]$ estimates varied widely and lacked consistent trends across conditions. The QSS approximation Eq. (17) accounts for this timescale separation by replacing poorly constrained parameters with a single term representing the net oncolytic effect. The monotonic increase of this parameter with viral dose confirms that meaningful therapeutic information is preserved despite the reduced representation.

Predictive performance depends not on mechanistic detail alone, but on whether parameters can be constrained by available data. Phan et al. [25] showed that complex models require more data for parameter identification, while Gevertz and Kareva [26] demonstrated that identifiability analysis can guide both experimental design and model selection. Under data-limited conditions, representing the leading mechanisms supported by observations is more informative than introducing biologically plausible but unresolvable structure. For extended treatment courses beyond 2–4 weeks, models could incorporate reversible and terminal exhaustion states; however, for acute assays, minimal formulations are sufficient.

Comparing model formulations allows separation of robustly identifiable parameters from those sensitive to model structure. The tumor growth rate was the most consistent across all contexts, showing minimal variation, reflecting its clear contribution to observed tumor dynamics during control and early treatment phases. The infected cell lysis rate in combination therapy was also stable across formulations, indicating that virus-mediated cell death can be accurately captured. CAR-T proliferation showed improved consistency in combination therapy relative to monotherapy, suggesting that interacting mechanisms provide additional information to constrain this parameter.

The impact of model complexity on parameter behavior is evident when examining estimates across formulations. The model without exhaustion Eq. (10) produced CAR-T proliferation rates ( $[eqn]$ day^−1^), whereas the Full model Eq. (5) yielded inflated values ( $[eqn]$ day^−1^) that would be difficult to sustain over the 72-hour observation period. This discrepancy reflects non-identifiability between proliferation and exhaustion terms: when exhaustion suppresses CAR-T activity, the optimizer compensates by inflating proliferation rates. The exhaustion parameter itself showed high uncertainty and did not improve fit, confirming it cannot be resolved from tumor data at this timescale.

In contrast to robustly estimated parameters, CAR-T killing rates for both uninfected $[eqn]$ and infected cells $[eqn]$ varied substantially between models, preventing reliable interpretation of direct cytotoxic effects from tumor data alone. Viral burst size $[eqn]$ and clearance $[eqn]$ remained non-identifiable wherever explicit kinetics were included. Despite this variability, several parameters exhibited consistent dose-dependent trends: the CAR-T killing rate decreased with increasing tumor burden across all formulations, consistent with limits on consecutive killing, while the QSS parameter $[eqn]$ increased with viral dose, matching the expected relationship between viral load and therapeutic impact.

The fact that multiple formulations fit the data equally well highlights a fundamental limitation of inference from aggregate tumor measurements. When several mechanisms contribute to a common outcome, individual effects may be non-identifiable even if the net response is well captured. Processes that proved non-identifiable, including detailed viral kinetics and CAR-T exhaustion, would require complementary data: direct viral titer measurements, flow cytometry of CAR-T phenotypes and exhaustion markers, or single-cell functional assays. Parameter estimates from detailed models should not be over-interpreted as true biological rates when simpler formulations achieve similar accuracy with improved reliability.

Our results highlight that useful models balance mechanistic detail with what the data can reliably resolve. The most informative model is not necessarily the most detailed, but rather one whose structure supports robust parameter estimation and meaningful biological interpretation.

Supplementary Material

1