Somayyeh Azizi^*
Independent Researcher, Sendai, Japan.
DOI: 10.4236/jamp.2026.141010   PDF    HTML   XML   163 Downloads   601 Views   Citations

Abstract

This study presents a novel fractional-order mathematical model to investigate the dynamics of glioma-immune interactions under therapeutic interventions. Building upon and extending previous models by de Pillis et al. and Pillay et al., we incorporate immune checkpoint inhibition and fractional-order derivatives to capture the memory-dependent behavior of immune responses. The model distinguishes between therapy-sensitive and resistant glioma cells and includes key immune components such as natural killer cells, cytotoxic T lymphocytes, cytokines, and exhausted effectors. An optimal control problem is formulated with three control variables representing immunotherapy, virus therapy, and checkpoint blockade. The forward-backward sweep method is employed to compute optimal treatment strategies over a 50-day horizon. Numerical simulations demonstrate that the fractional-order framework significantly influences treatment outcomes, with lower fractional orders delaying immune activation and reducing therapeutic efficacy. The proposed optimal control strategy achieves superior glioma suppression and immune activation compared to monotherapies and fixed-dose combinations. These findings highlight the importance of incorporating memory effects and multi-modal control in the design of effective glioma treatment protocols.

Keywords

Glioma Modeling, Fractional-Order Differential Equations, Immunotherapy, Virus Therapy, Checkpoint Inhibition, Optimal Control, Forward-Backward Sweep Method, Immune Dynamics

Share and Cite:

1. Introduction

Gliomas are the most common and aggressive type of primary brain malignancy, originating from glial cells and accounting for most malignant glial neoplasms in adults. High-grade gliomas, such as glioblastoma multiforme (GBM), are particularly lethal due to their rapid proliferation, diffuse infiltration into surrounding brain tissue, and resistance to conventional therapies. Despite advances in surgical resection, radiotherapy, and chemotherapy, the median survival time for GBM patients remains approximately 12 to 15 months, with a five-year survival rate below 10% [1] [2]. In recent years, immunotherapy has emerged as a promising approach to glioma treatment, aiming to harness the body’s immune system to recognize and eliminate malignant glial cells. Strategies such as immune checkpoint inhibitors, adoptive T cell transfer, and oncolytic virus therapies have shown encouraging results in preclinical and early clinical studies [3] [4]. However, gliomas are known to create a highly immunosuppressive microenvironment that impairs immune cell infiltration and function. The overexpression of checkpoint molecules like PD-1 and CTLA-4, along with the accumulation of exhausted T cells and regulatory immune cells, contributes to immune evasion and limits the effectiveness of immunotherapy [5] [6]. Mathematical modeling has become an essential tool for understanding glioma-immune interactions and optimizing treatment strategies. Early models by de Pillis et al. [7] introduced validated frameworks for simulating the cell-mediated immune response to glioma growth and for exploring the effects of combined immunotherapy and chemotherapy. These models provided critical insights into how immune cells interact with glioma populations and how treatment timing and intensity influence outcomes. Building on this foundation, Pillay et al. [8] [9] developed more recent models incorporating virus therapy and fractional-order dynamics, demonstrating that memory effects in biological systems can significantly alter treatment efficacy.

Inspired by these foundational works, the present study extends the modeling framework in several important directions. First, we adopt a fractional-order differential equation system to capture the memory-dependent behavior of immune responses, which is particularly relevant in chronic and slowly evolving diseases like glioma. Second, we introduce a checkpoint inhibition mechanism to reflect the role of immune suppression pathways in glioma progression, an aspect not explicitly addressed in earlier models. Third, we formulate a multi-control optimal control problem involving immunotherapy, virus therapy, and checkpoint inhibition, allowing for the design of adaptive treatment strategies that balance therapeutic efficacy with cost and toxicity. Finally, we conduct a comparative analysis of treatment outcomes under both integer-order and fractional-order dynamics, evaluating the impact of each therapeutic component and the importance of memory effects on glioma suppression and immune activation.

The primary objective of this study is to develop a comprehensive and flexible mathematical framework for simulating glioma-immune interactions and to identify optimal, personalized treatment strategies. By integrating fractional calculus with optimal control theory, we aim to provide deeper insights into the timing, synergy, and intensity of therapeutic interventions. Our results offer valuable guidance for the design of future immunotherapeutic protocols and contribute to the broader effort of improving clinical outcomes for glioma patients.

The remainder of this paper is organized as follows. In Section 2, we present the formulation of the fractional-order glioma-immune model. Section 3 provides the equilibrium analysis of the untreated system, while Section 4 extends this analysis to the treated model under therapeutic interventions. In Section 5, we formulate and analyze the optimal control problem involving immunotherapy, virus therapy, and checkpoint inhibition. Section 6 presents numerical simulations and a comparative discussion of treatment strategies. Finally, conclusions are summarized in Section 7.

2. Model Formulation

We construct a fractional-order mathematical model to describe the complex interactions among therapy-sensitive and therapy-resistant glioma cells, healthy lymphocytes, immune effector cells (both active and exhausted), immune checkpoint molecules, cytokines, and therapeutic agents. The system is governed by a set of nonlinear fractional differential equations using Caputo derivatives of order $\alpha \in \left(0,1\right]$ , which allow the model to capture memory effects a critical feature for representing the delayed and cumulative nature of immune responses and treatment dynamics. The Caputo fractional derivative of order $\alpha \in \left(0,1\right]$ for a function $f\left(t\right)$ is defined as:

$D_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{f^{\prime }\left(s\right)}{{\left(t-s\right)}^{\alpha }}\text{d}s},$

where $\Gamma (\cdot )$ is the Gamma function. This operator incorporates the entire history of the function $f\left(t\right)$ , making it ideal for modeling biological systems with non-local temporal behavior. Let $T_s\left(t\right)$ and $T_r\left(t\right)$ denote the populations of therapy-sensitive and therapy-resistant glioma cells, respectively. The total glioma burden is given by $T\left(t\right)=T_s\left(t\right)+T_r\left(t\right)$ . The variable $N\left(t\right)$ represents healthy lymphocytes, while $C\left(t\right)$ and $E\left(t\right)$ denote active and exhausted immune effector cells. The concentrations of checkpoint molecules and cytokines are represented by $P\left(t\right)$ and $I\left(t\right)$ , respectively. Therapeutic interventions are modeled by $M\left(t\right)$ for drug therapy and $B\left(t\right)$ for oncolytic virus therapy. This formulation integrates key immunological mechanisms such as immune activation and suppression, checkpoint inhibition, immune exhaustion and recovery, cytokine-mediated feedback, and glioma resistance evolution, providing a comprehensive framework for analyzing glioma-immune dynamics under combination therapy. The full system of fractional-order differential equations is presented in Appendix A.

2.1. Untreated Model

In the absence of therapy, the system evolves according to the following equations:

$D^{\alpha }{T}_s\left(t\right)=a_s{T}_s\left(1-b_s{T}_s\right)-c_{1}N\left(t\right)T_s-\omega T_s,$

$D^{\alpha }{T}_r\left(t\right)=a_r{T}_r\left(1-b_r{T}_r\right)-c_{2}N\left(t\right)T_r+\omega T_s,$

$D^{\alpha }N\left(t\right)=w-fN+\frac{g{T}^2}{h+T^2+\varphi N}\left(1+\kappa I\right)N-pNT-\theta P\left(t\right)N,$

$D^{\alpha }C\left(t\right)=-mC-qCT+rNT-\delta C+\zeta E,$

$D^{\alpha }E\left(t\right)=\delta C-\mu E-\zeta E,$

$D^{\alpha }M\left(t\right)=0,$

$D^{\alpha }B\left(t\right)=0,$

$D^{\alpha }P\left(t\right)=\sigma T-\lambda P,$

$D^{\alpha }I\left(t\right)={\rho }_{1}N+{\rho }_{2}C-{\rho }_{3}I.$

In this setting, glioma cells grow logistically and are suppressed by lymphocytes. A fraction $\omega$ of sensitive cells mutates into the resistant phenotype. Lymphocytes are produced at a constant rate, decay naturally, and are activated by glioma presence through a hill-type function modulated by cytokines. They are suppressed by glioma-induced immunosuppression and checkpoint molecules. Effector cells are activated by lymphocyte-glioma interactions, suppressed by glioma contact, and may become exhausted or recover. Checkpoint molecules and cytokines evolve naturally without therapeutic intervention.

2.2. Treated Model with Combination Therapy

To simulate the effects of therapy, we extend the model to include drug and viral treatments, immune response delay, checkpoint inhibition, and other immunological feedback mechanisms. The treated system is described by:

$\begin{array}{l}{D}^{\alpha }{T}_s\left(t\right)=a_s{T}_s\left(1-b_s{T}_s\right)-c_{1}N\left(t-\tau \right)T_s-B\left(t\right)T_s-k_{T}M\left(t\right)T_s\\ -\rho T_s^2-\omega T_s,\end{array}$

$D^{\alpha }{T}_r\left(t\right)=a_r{T}_r\left(1-b_r{T}_r\right)-c_{2}N\left(t-\tau \right)T_r-B\left(t\right)T_r-\rho T_r^2+\omega T_s,$

$\begin{array}{l}{D}^{\alpha }N\left(t\right)=w-fN+\frac{g{T}^2}{h+T^2+\varphi N}\left(1+\kappa I\right)N-pNT-k_{N}M\left(t\right)N\\ -\theta P\left(t\right)N+ϵU_P\left(t\right)N,\end{array}$

$D^{\alpha }C\left(t\right)=-mC+\frac{j{B}^2}{k+B^2}C-qCT+rNT-k_{C}M\left(t\right)C-\delta C+\zeta E,$

$D^{\alpha }E\left(t\right)=\delta C-\mu E-\zeta E,$

$D^{\alpha }M\left(t\right)=-\gamma M+U_M\left(t\right),$

$D^{\alpha }B\left(t\right)=-\eta B+U_V\left(t\right),$

$D^{\alpha }P\left(t\right)=\sigma T-\lambda P-\xi U_P\left(t\right),$

$D^{\alpha }I\left(t\right)={\rho }_{1}N+{\rho }_{2}C-{\rho }_{3}I.$

In this extended system, glioma cells are suppressed by both immune and therapeutic agents. Lymphocytes are affected by drug toxicity and checkpoint inhibition, with $U_P\left(t\right)$ restoring immune function. Effector cells are activated and suppressed as before, but now also experience drug-induced toxicity. Checkpoint molecules are cleared by therapy, and cytokines continue to mediate immune feedback. The control functions $U_M\left(t\right)$ , $U_V\left(t\right)$ , and $U_P\left(t\right)$ represent the administration rates of drug, viral, and checkpoint inhibitor therapies, respectively. The terms $-\gamma M$ and $-\eta B$ model natural decay of drug and viruses, while the Hill-type term $\frac{j{B}^2}{k+B^2}$ captures saturation in virus-induced immune activation. The nonlinear term $\rho T^2$ reflects glioma resistance to therapy. This formulation captures the complex interplay between glial progression and immune regulation under modern combination therapies, incorporating immune suppression, checkpoint signaling, immune exhaustion and recovery, cytokine-mediated feedback, treatment delay, and biological memory.

2.3. Immune Suppression and Recovery Dynamics

To capture the multifaceted regulation of immune activity in glial malignancies, we extend the model to include checkpoint inhibition, reversible immune exhaustion, and cytokine-mediated feedback. These mechanisms reflect the dynamic balance between immune suppression and recovery, which is critical in the glioma microenvironment. We introduce the state variable $P\left(t\right)$ , representing the concentration of immune checkpoint molecules such as PD-1/PD-L1. These molecules are upregulated by glial cells and suppress lymphocyte activity through the term $-\theta P\left(t\right)N\left(t\right)$ . The checkpoint dynamics are governed by:

$D^{\alpha }P\left(t\right)=\sigma T-\lambda P-\xi U_P\left(t\right)$

Here, $\sigma$ is the production rate of checkpoint molecules by glial cells, $\lambda$ is their natural decay rate, and $\xi$ quantifies the therapeutic clearance induced by the checkpoint inhibitor control $U_P\left(t\right)$ . The term $ϵU_P\left(t\right)N\left(t\right)$ models the restoration of immune function due to checkpoint blockade. To represent immune exhaustion and its reversibility, we include the exhausted immune cell compartment $E\left(t\right)$ . Effector cells $C\left(t\right)$ transition into exhaustion at rate $\delta$ , while a fraction $\zeta$ of exhausted cells can recover and return to the active pool:

$\begin{array}{l}{D}^{\alpha }C\left(t\right)\supset +\zeta E\left(t\right)\\ D^{\alpha }E\left(t\right)\supset -\zeta E\left(t\right)\end{array}$

This reversible exhaustion mechanism reflects the potential for immune recovery, particularly under checkpoint inhibition. Additionally, we introduce the cytokine concentration $I\left(t\right)$ , representing signaling molecules such as interleukins or interferons. These are produced by lymphocytes and effector cells and decay naturally:

$D^{\alpha }I\left(t\right)={\rho }_{1}N+{\rho }_{2}C-{\rho }_{3}I$

Cytokines enhance immune activation through the amplification factor $\left(1+\kappa I\right)$ , which modulates the lymphocyte stimulation term. This feedback loop captures the self-reinforcing nature of immune signaling and its role in sustaining anti-glioma responses. Together, these components provide a comprehensive framework for modeling immune suppression and recovery, enabling the analysis of checkpoint inhibition, immune exhaustion, and cytokine signaling interact with therapeutic interventions and glioma dynamics.

2.4. Comparison with Existing Models

The proposed model significantly extends classical glioma-immune interaction frameworks by incorporating multiple layers of biological complexity relevant to glial malignancies. It is a developed extension of the mathematical model introduced by de Pillis et al. [7], incorporating additional mechanisms to more accurately represent tumor-immune dynamics under drug-based immunotherapy.

Key distinctions between our model and existing approaches include:

• Fractional-Order Dynamics: Unlike traditional models that use integer-order derivatives, our system employs Caputo fractional derivatives to capture memory effects and the non-local behavior of biological processes. This allows for a more accurate representation of delayed immune responses and long-term interactions between glial and immune cells.

• Glial Resistance Evolution: We introduce a two-compartment glial population, distinguishing between therapy-sensitive ($T_s\left(t\right)$ ) and therapy-resistant ($T_r\left(t\right)$ ) cells. This enables the modeling of mutation-driven resistance, a critical factor in glioma progression and treatment failure, which is not addressed in earlier models.

• Immune Checkpoint Inhibition: Our model explicitly incorporates the dynamics of immune checkpoint molecules ($P\left(t\right)$ ) and their suppressive effects on lymphocyte activity. The inclusion of a control function $U_P\left(t\right)$ allows for the simulation of checkpoint inhibitor therapies, a key component of modern immuno-oncology not present in the de Pillis framework.

• Reversible Immune Exhaustion: We model the transition of effector immune cells into an exhausted state ($E\left(t\right)$ ) and their potential recovery via a reactivation term ($\zeta$ ). This addition reflects recent biological findings on immune plasticity and the potential for functional restoration under immunotherapy.

• Cytokine-Mediated Feedback: The inclusion of a cytokine variable ($I\left(t\right)$ ) introduces a regulatory feedback loop that modulates immune activation. This mechanism captures the role of signaling molecules such as interleukins and interferons in enhancing immune responses, a feature absent in most existing models.

Collectively, these enhancements position our model as a comprehensive and biologically informed framework for studying glial-immune interactions under combination therapy. By integrating modern immunological insights and advanced mathematical tools, the model provides a robust platform for optimizing treatment strategies and exploring the nonlinear dynamics of glial-immune coevolution.

3. Equilibrium Analysis of the Untreated Model

In this section, we analyze the equilibrium behavior of the untreated glioma-immune interaction model. Our objective is to characterize the long-term dynamics of the system in the absence of therapeutic intervention by identifying and classifying its steady states. We begin by deriving the disease-free equilibrium (DFE), which corresponds to complete glioma eradication, using direct analytical substitution. To reinforce this result, we re-derive the DFE using a general nonlinear stability theorem, ensuring both biological clarity and mathematical rigor.

We then investigate the existence and structure of a positive glioma-persistent equilibrium (GPE), representing a biologically meaningful state in which glioma cells coexist with immune components. By reducing the model to a tractable subsystem, we derive a nonlinear scalar equation whose solutions determine the conditions under which the glioma persists. A graphical analysis of this equation provides further insight into the threshold dynamics governing glioma persistence.

Finally, we examine the local stability of both the DFE and the GPE by computing the Jacobian matrix at each equilibrium and analyzing its eigenvalues in the context of Caputo-type fractional-order systems. This comprehensive equilibrium analysis lays the foundation for understanding how therapeutic interventions may shift the system from a glioma-persistent state to a disease-free state, which will be explored in the next section.

3.1. Disease-Free Equilibrium (DFE): Direct Analytical Derivation

We now analyze the equilibrium behavior of the untreated system described in Section 2. To determine the disease-free equilibrium (DFE), we assume the complete absence of glioma cells by setting $T_s^{\text{*}}=0$ and $T_r^{\text{*}}=0$ . We then set all fractional derivatives to zero and solve for the steady-state values of the remaining variables.

From the first two equations of the model, we have:

$0=a_s{T}_s^{\text{*}}\left(1-b_s{T}_s^{\text{*}}\right)-c_1{N}^{\text{*}}{T}_s^{\text{*}}-\omega T_s^{\text{*}}\Rightarrow T_s^{\text{*}}=0,$

$0=a_r{T}_r^{\text{*}}\left(1-b_r{T}_r^{\text{*}}\right)-c_2{N}^{\text{*}}{T}_r^{\text{*}}+\omega T_s^{\text{*}}\Rightarrow T_r^{\text{*}}=0.$

Substituting $T^{\text{*}}=T_s^{\text{*}}+T_r^{\text{*}}=0$ into the lymphocyte equation yields:

$0=w-f{N}^{\text{*}}+\frac{g\cdot 0}{h+0+\varphi N^{\text{*}}}\left(1+\kappa I^{\text{*}}\right)N^{\text{*}}-p{N}^{\text{*}}\cdot 0-\theta P^{\text{*}}{N}^{\text{*}},$

$\Rightarrow 0=w-f{N}^{\text{*}}-\theta P^{\text{*}}{N}^{\text{*}}.$

Solving for $N^{\text{*}}$ , we obtain:

$N^{\text{*}}=\frac{w}{f+\theta P^{\text{*}}}.$

From the effector cell equation:

$0=-m{C}^{\text{*}}-q{C}^{\text{*}}\cdot 0+r{N}^{\text{*}}\cdot 0-\delta C^{\text{*}}+\zeta E^{\text{*}}\Rightarrow 0=-\left(m+\delta \right)C^{\text{*}}+\zeta E^{\text{*}},$

$\Rightarrow C^{\text{*}}=\frac{\zeta }{m+\delta }{E}^{\text{*}}.$

From the exhaustion equation:

$0=\delta C^{\text{*}}-\mu E^{\text{*}}-\zeta E^{\text{*}},$

$\Rightarrow \delta \cdot \frac{\zeta }{m+\delta }{E}^{\text{*}}=\left(\mu +\zeta \right)E^{\text{*}}.$

This equation is satisfied only if $E^{\text{*}}=0$ , which implies $C^{\text{*}}=0$ .

From the checkpoint equation:

$0=\sigma \cdot 0-\lambda P^{\text{*}}\Rightarrow P^{\text{*}}=0.$

From the cytokine equation:

$0={\rho }_1{N}^{\text{*}}+{\rho }_2{C}^{\text{*}}-{\rho }_3{I}^{\text{*}}\Rightarrow {\rho }_1{N}^{\text{*}}={\rho }_3{I}^{\text{*}},$

$\Rightarrow I^{\text{*}}=\frac{{\rho }_1}{{\rho }_3}{N}^{\text{*}}=\frac{{\rho }_{1}w}{{\rho }_3\left(f+\theta P^{\text{*}}\right)}.$

Substituting $P^{\text{*}}=0$ , we obtain:

$I^{\text{*}}=\frac{{\rho }_{1}w}{f{\rho }_3}.$

Since no therapy is applied in the untreated model, we also have:

$M^{\text{*}}=0,B^{\text{*}}=0.$

Therefore, the disease-free equilibrium is given by:

$\text{DFE}=\left(T_s^{\text{*}},T_r^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}},M^{\text{*}},B^{\text{*}},P^{\text{*}},I^{\text{*}}\right)=\left(0,0,\frac{w}{f},0,0,0,0,0,\frac{{\rho }_{1}w}{f{\rho }_3}\right).$

This equilibrium represents a biologically meaningful state in which both sensitive and resistant glioma cells are eradicated. The immune system maintains a stable population of lymphocytes, while effectors and exhausted immune cells are absent due to the lack of glioma-induced stimulation. Checkpoint molecules are cleared, and cytokines persist at a basal level due to ongoing lymphocyte activity.

3.2. Re-Derivation of the DFE Using a General Stability Theorem

To confirm the validity of the disease-free equilibrium and prepare for a rigorous stability analysis, we now re-derive the DFE using a general nonlinear stability theorem. We consider a nonlinear system of the form:

$\frac{\text{d}\overrightarrow{x}}{\text{d}t}=J\left(\overrightarrow{x}-\overrightarrow{p}\right)+\overrightarrow{P}\left(\overrightarrow{x}\right),$

where $\overrightarrow{x}\in {ℝ}^n$ is the state vector, $\overrightarrow{p}\in {ℝ}^n$ is an equilibrium point, $J\in {ℝ}^{n\times n}$ is a constant matrix, and $\overrightarrow{P}\left(\overrightarrow{x}\right)$ is a continuously differentiable nonlinear function. The function $\overrightarrow{P}\left(\overrightarrow{x}\right)$ satisfies the conditions:

$\overrightarrow{P}\left(\overrightarrow{p}\right)=\overrightarrow{0}\text{and}\underset{\overrightarrow{x}\to \overrightarrow{p}}{\lim }\frac{\Vert \overrightarrow{P}\left(\overrightarrow{x}\right)\Vert }{\Vert \overrightarrow{x}-\overrightarrow{p}\Vert }=0.$

This implies that $\overrightarrow{P}\left(\overrightarrow{x}\right)$ consists of at least second-order terms near $\overrightarrow{p}$ . Under these conditions, the local stability of the equilibrium point $\overrightarrow{p}$ is determined by the eigenvalues of the matrix $J$ . If all eigenvalues of $J$ have strictly negative real parts, then $\overrightarrow{p}$ is locally asymptotically stable. If at least one eigenvalue has a positive real part, then $\overrightarrow{p}$ is unstable. To apply this theorem to our model, we define the extended state vector:

$\overrightarrow{x}=\left[\begin{array}{c}{T}_s\\ T_r\\ N\\ C\\ E\\ P\\ I\end{array}\right],\overrightarrow{p}=\left[\begin{array}{c}0\\ 0\\ \frac{w}{f}\\ 0\\ 0\\ 0\\ \frac{{\rho }_{1}w}{f{\rho }_3}\end{array}\right].$

We express the system as:

$\frac{\text{d}\overrightarrow{x}}{\text{d}t}=\overrightarrow{F}\left(\overrightarrow{x}\right)=J\left(\overrightarrow{x}-\overrightarrow{p}\right)+\overrightarrow{P}\left(\overrightarrow{x}\right),$

where $\overrightarrow{F}\left(\overrightarrow{x}\right)$ represents the right-hand side of the model equations. We now verify that the nonlinear term $\overrightarrow{P}\left(\overrightarrow{x}\right)$ satisfies the required conditions. In the equations for $T_s$ and $T_r$ , the nonlinear terms include logistic growth $T_s^2$ , $T_r^2$ , and bilinear interactions with $N$ , all of which vanish at $T_s=T_r=0$ and are at least quadratic near the DFE. The lymphocyte equation includes the nonlinear term:

$\frac{g{T}^2}{h+T^2+\varphi N}\left(1+\kappa I\right)N-pNT-\theta PN,$

which vanishes at $T=0$ , $P=0$ , and is at least quadratic in $T$ , $P$ , and $I$ . Similarly, the effector cell equation includes bilinear terms $-qCT+rNT$ , which vanish at $T=0$ , and the exhaustion equation includes $\delta C$ , which vanishes at $C=0$ . The checkpoint and cytokine equations are linear in $T$ , $N$ , and $C$ , and vanish at the DFE. Therefore, all nonlinear terms in $\overrightarrow{P}\left(\overrightarrow{x}\right)$ vanish at $\overrightarrow{p}$ , and are of order two or higher near the DFE. Hence, the system satisfies the structural conditions of the general nonlinear stability theorem, and the equilibrium point:

$\text{DFE}=\left(0,0,\frac{w}{f},0,0,0,\frac{{\rho }_{1}w}{f{\rho }_3}\right)$

is again confirmed as a valid equilibrium under this general framework. In Section 3.3, we will compute the Jacobian matrix at the DFE and analyze its eigenvalues to determine the local stability of the system.

3.3. Stability of the Disease-Free Equilibrium—Untreated Model

To analyze the local stability of the disease-free equilibrium (DFE) obtained in Section 3.1, we apply the general nonlinear stability theorem introduced in Section 3.2. According to this theorem, the local behavior of a nonlinear system near an equilibrium point can be determined by examining the eigenvalues of the Jacobian matrix evaluated at that point, provided that the nonlinear terms vanish at the equilibrium and are at least of second order in a neighborhood of that point. We consider the untreated subsystem of the full model presented in Section 2, focusing on the dynamics of the therapy-sensitive glioma cells $T_s\left(t\right)$ , lymphocytes $N\left(t\right)$ , active immune effector cells $C\left(t\right)$ , and exhausted immune cells $E\left(t\right)$ . The governing equations for this subsystem are as follows:

$D^{\alpha }{T}_s\left(t\right)=a_s{T}_s\left(t\right)\left(1-b_s{T}_s\left(t\right)\right)-c_{1}N\left(t\right)T_s\left(t\right)-\omega T_s\left(t\right),$

$\begin{array}{l}{D}^{\alpha }N\left(t\right)=w-fN\left(t\right)+\frac{g{\left(T_s\left(t\right)+T_r\left(t\right)\right)}^2}{h+{\left(T_s\left(t\right)+T_r\left(t\right)\right)}^2+\varphi N\left(t\right)}\left(1+\kappa I\left(t\right)\right)N\left(t\right)\\ -pN\left(t\right)\left(T_s\left(t\right)+T_r\left(t\right)\right)-\theta P\left(t\right)N\left(t\right),\end{array}$

$\begin{array}{l}{D}^{\alpha }C\left(t\right)=-mC\left(t\right)-qC\left(t\right)\left(T_s\left(t\right)+T_r\left(t\right)\right)\\ +rN\left(t\right)\left(T_s\left(t\right)+T_r\left(t\right)\right)-\delta C\left(t\right)+\zeta E\left(t\right),\end{array}$

$D^{\alpha }E\left(t\right)=\delta C\left(t\right)-\mu E\left(t\right)-\zeta E\left(t\right).$

At the disease-free equilibrium, we assume that

$T_s^{\text{*}}=T_r^{\text{*}}=C^{\text{*}}=E^{\text{*}}=P^{\text{*}}=0,$

and from Section 3.1, we have

$N^{\text{*}}=\frac{w}{f},I^{\text{*}}=\frac{{\rho }_{1}w}{f{\rho }_3}.$

Let the state vector be defined as

$\overrightarrow{x}=\left[\begin{array}{c}{T}_s\\ N\\ C\\ E\end{array}\right],$

and let the vector field be denoted by

$\overrightarrow{F}\left(\overrightarrow{x}\right)=\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right],$

where

$F_1=a_s{T}_s-a_s{b}_s{T}_s^2-c_{1}N{T}_s-\omega T_s,$

$F_2=w-fN+\frac{g{T}^2}{h+T^2+\varphi N}\left(1+\kappa I\right)N-pNT-\theta PN,$

$F_3=-mC-qCT+rNT-\delta C+\zeta E,$

$F_4=\delta C-\mu E-\zeta E,$

and $T=T_s+T_r$ . At the DFE, we substitute

$T=0,N=\frac{w}{f},C=0,E=0,P=0,I=\frac{{\rho }_{1}w}{f{\rho }_3}.$

We now compute the Jacobian matrix $J=\left[\frac{\partial F_i}{\partial x_j}\right]$ evaluated at the DFE. The

partial derivatives are given by:

$\frac{\partial F_1}{\partial T_s}=a_s-2a_s{b}_s{T}_s-c_{1}N-\omega ,\frac{\partial F_1}{\partial N}=-c_1{T}_s,$

$\frac{\partial F_1}{\partial C}=0,\frac{\partial F_1}{\partial E}=0.$

Evaluating at the DFE yields:

$\frac{\partial F_1}{\partial T_s}=a_s-\frac{c_{1}w}{f}-\omega ,\frac{\partial F_1}{\partial N}=0.$

For the second component:

$\frac{\partial F_2}{\partial T_s}=0,\frac{\partial F_2}{\partial N}=-f,\frac{\partial F_2}{\partial C}=0,\frac{\partial F_2}{\partial E}=0.$

For the third component:

$\frac{\partial F_3}{\partial T_s}=rN-qC,\frac{\partial F_3}{\partial N}=r{T}_s,$

$\frac{\partial F_3}{\partial C}=-m-\delta -q{T}_s,\frac{\partial F_3}{\partial E}=\zeta .$

Evaluating at the DFE gives:

$\frac{\partial F_3}{\partial T_s}=\frac{rw}{f},\frac{\partial F_3}{\partial N}=0,\frac{\partial F_3}{\partial C}=-m-\delta ,\frac{\partial F_3}{\partial E}=\zeta .$

For the fourth component:

$\frac{\partial F_4}{\partial T_s}=0,\frac{\partial F_4}{\partial N}=0,\frac{\partial F_4}{\partial C}=\delta ,\frac{\partial F_4}{\partial E}=-\mu -\zeta .$

Substituting all partial derivatives, the Jacobian matrix at the DFE is:

$J=\left[\begin{array}{cccc}{a}_s-\frac{c_{1}w}{f}-\omega & 0& 0& 0\\ 0& -f& 0& 0\\ \frac{rw}{f}& 0& -m-\delta & \zeta \\ 0& 0& \delta & -\mu -\zeta \end{array}\right].$

Since this matrix is lower triangular, its eigenvalues are the entries on the diagonal:

${\lambda }_1=a_s-\frac{c_{1}w}{f}-\omega ,{\lambda }_2=-f,{\lambda }_3=-m-\delta ,{\lambda }_4=-\mu -\zeta .$

Given that all parameters $f$ , $m$ , $\delta$ , $\mu$ , and $\zeta$ are strictly positive, it follows that ${\lambda }_2<0$ , ${\lambda }_3<0$ , and ${\lambda }_4<0$ . The sign of ${\lambda }_1$ depends on the inequality

${\lambda }_1<0\iff \frac{c_{1}w}{f}>a_s-\omega .$

Therefore, the disease-free equilibrium is locally asymptotically stable if and only if the inequality

$\frac{c_{1}w}{f}>a_s-\omega$

is satisfied. This condition implies that the cytotoxic effect of lymphocytes on glioma cells, represented by the parameter $c_1$ , in combination with their production rate $w$ , must be sufficiently strong to overcome both the intrinsic growth rate $a_s$ of therapy-sensitive glioma cells and the mutation rate $\omega$ into resistant phenotypes. If this inequality is not satisfied, the glioma population may escape immune control, resulting in the instability of the disease-free equilibrium.

3.4. Endemic and Glioma-Persistent Equilibria—Untreated Model

When the disease-free equilibrium (DFE) becomes unstable, the system may admit one or more positive steady states in which glioma persists in coexistence with immune components. These endemic equilibria, also referred to as glioma-persistent equilibria (GPEs), represent chronic disease states where the immune system fails to fully eliminate the glioma population. In this section, we analyze the existence and structure of such equilibria in the untreated model.

3.4.1. General Endemic Equilibrium

We begin by considering the fully untreated model and derive the conditions under which an endemic equilibrium exists. This equilibrium satisfies $T_s^{\text{*}}+T_r^{\text{*}}>0$ , indicating the persistence of glioma in the system. Setting all derivatives to zero, we obtain the following algebraic system:

$0=a_s{T}_s^{\text{*}}\left(1-b_s{T}_s^{\text{*}}\right)-c_1{N}^{\text{*}}{T}_s^{\text{*}}-\omega T_s^{\text{*}},$

$0=a_r{T}_r^{\text{*}}\left(1-b_r{T}_r^{\text{*}}\right)-c_2{N}^{\text{*}}{T}_r^{\text{*}}+\omega T_s^{\text{*}},$

$0=w-f{N}^{\text{*}}+\frac{g{\left(T^{\text{*}}\right)}^2}{h+{\left(T^{\text{*}}\right)}^2+\varphi N^{\text{*}}}\left(1+\kappa I^{\text{*}}\right)N^{\text{*}}-p{N}^{\text{*}}{T}^{\text{*}}-\theta P^{\text{*}}{N}^{\text{*}},$

$0=-m{C}^{\text{*}}-q{C}^{\text{*}}{T}^{\text{*}}+r{N}^{\text{*}}{T}^{\text{*}}-\delta C^{\text{*}}+\zeta E^{\text{*}},$

$0=\delta C^{\text{*}}-\mu E^{\text{*}}-\zeta E^{\text{*}},$

$0={\rho }_1{N}^{\text{*}}+{\rho }_2{C}^{\text{*}}-{\rho }_3{I}^{\text{*}},$

$0=\sigma T^{\text{*}}-\lambda P^{\text{*}},$

where $T^{\text{*}}=T_s^{\text{*}}+T_r^{\text{*}}$ is the total glioma burden.

From the last two equations, we obtain:

$P^{\text{*}}=\frac{\sigma T^{\text{*}}}{\lambda },I^{\text{*}}=\frac{{\rho }_1{N}^{\text{*}}+{\rho }_2{C}^{\text{*}}}{{\rho }_3}.$

Solving the exhaustion equation yields:

$E^{\text{*}}=\frac{\delta }{\mu +\zeta }{C}^{\text{*}}.$

Substituting into the effector equation gives:

$C^{\text{*}}=\frac{r{N}^{\text{*}}{T}^{\text{*}}}{m+\delta +q{T}^{\text{*}}-\frac{\zeta \delta }{\mu +\zeta }}.$

Substituting this into the expression for $E^{\text{*}}$ , we obtain:

$E^{\text{*}}=\frac{\delta }{\mu +\zeta }\cdot \frac{r{N}^{\text{*}}{T}^{\text{*}}}{m+\delta +q{T}^{\text{*}}-\frac{\zeta \delta }{\mu +\zeta }}.$

Next, we return to the equation for the therapy-sensitive glioma population. Assuming $T_s^{\text{*}}>0$ , we factor out $T_s^{\text{*}}$ to obtain:

$T_s^{\text{*}}\left[a_s\left(1-b_s{T}_s^{\text{*}}\right)-c_1{N}^{\text{*}}-\omega \right]=0,$

which implies:

$N^{\text{*}}=\frac{a_s\left(1-b_s{T}_s^{\text{*}}\right)-\omega }{c_1}.$

Substituting this expression into the lymphocyte equation yields a nonlinear algebraic equation in $T^{\text{*}}$ . Using the expressions for $I^{\text{*}}$ and $P^{\text{*}}$ in terms of $N^{\text{*}}$ , $C^{\text{*}}$ , and $T^{\text{*}}$ , this equation becomes a single nonlinear equation in $T^{\text{*}}$ . Once a positive solution $T^{\text{*}}$ is found (analytically or numerically), the corresponding values of $N^{\text{*}},C^{\text{*}},E^{\text{*}},P^{\text{*}}$ , and $I^{\text{*}}$ can be computed using the expressions above. The existence of a biologically meaningful endemic equilibrium thus depends on the solvability of the nonlinear equation in $T^{\text{*}}$ and the positivity of the resulting expressions for all state variables. This equilibrium represents a persistent glioma state in which both therapy-sensitive and resistant glial cells coexist with immune components. The presence of exhausted immune cells and elevated checkpoint molecule levels reflects the immunosuppressive environment characteristic of high-grade gliomas.

3.4.2. Reduced Model and Existence of a Positive GPE

To simplify the analysis and isolate the core glial-immune dynamics, we now consider a reduced version of the model by assuming the absence of therapy-resistant glioma cells and checkpoint molecules. Specifically, we set $T_r\left(t\right)=0$ , $P\left(t\right)=0$ , and treat $I\left(t\right)=I^{\text{*}}$ as a constant. The reduced system becomes:

$D^{\alpha }T\left(t\right)=aT\left(t\right)\left(1-bT\left(t\right)\right)-c_{1}N\left(t\right)T\left(t\right),$

$D^{\alpha }N\left(t\right)=w-fN\left(t\right)+\frac{gT{\left(t\right)}^2}{h+T{\left(t\right)}^2+\varphi N\left(t\right)}\left(1+\kappa I^{\text{*}}\right)N\left(t\right)-pN\left(t\right)T\left(t\right),$

$D^{\alpha }C\left(t\right)=-mC\left(t\right)-qC\left(t\right)T\left(t\right)+rN\left(t\right)T\left(t\right)-\delta C\left(t\right),$

$D^{\alpha }E\left(t\right)=\delta C\left(t\right)-\mu E\left(t\right).$

Setting all derivatives to zero, we solve the algebraic system:

$0=a{T}^{\text{*}}\left(1-b{T}^{\text{*}}\right)-c_1{N}^{\text{*}}{T}^{\text{*}},$

$0=w-f{N}^{\text{*}}+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}+\varphi N^{\text{*}}}\left(1+\kappa I^{\text{*}}\right)N^{\text{*}}-p{N}^{\text{*}}{T}^{\text{*}},$

$0=-m{C}^{\text{*}}-q{C}^{\text{*}}{T}^{\text{*}}+r{N}^{\text{*}}{T}^{\text{*}}-\delta C^{\text{*}},$

$0=\delta C^{\text{*}}-\mu E^{\text{*}}.$

From the last equation:

$E^{\text{*}}=\frac{\delta }{\mu }{C}^{\text{*}}.$

Substituting into the third equation:

$C^{\text{*}}=\frac{r{N}^{\text{*}}{T}^{\text{*}}}{m+\delta +q{T}^{\text{*}}}.$

Then:

$E^{\text{*}}=\frac{\delta }{\mu }\cdot \frac{r{N}^{\text{*}}{T}^{\text{*}}}{m+\delta +q{T}^{\text{*}}}.$

From the first equation, assuming $T^{\text{*}}>0$ , we divide both sides by $T^{\text{*}}$ to get:

$N^{\text{*}}=\frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}.$

Substituting this into the second equation yields:

$\begin{array}{l}0=w-f\cdot \frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}+\varphi \cdot \frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}}\left(1+\kappa I^{\text{*}}\right)\cdot \frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}\\ -p\cdot \frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}\cdot T^{\text{*}}.\end{array}$

Factoring out $\frac{a\left(1-b{T}^{\text{*}}\right)}{c_1}$ , we define the function:

$\Phi \left(T^{*}\right)=\frac{a\left(1-b{T}^{*}\right)}{c_1}\left[f-\frac{g{T}^{*2}}{h+T^{*2}+\varphi \cdot \frac{a\left(1-b{T}^{*}\right)}{c_1}}\left(1+\kappa I^{*}\right)+p{T}^{*}\right],$

and seek a value $T^{\text{*}}>0$ such that $\Phi \left(T^{*}\right)=w$ . The function $\Phi \left(T\right)$ is continuous on $\left(0,1/b\right)$ , with:

$\underset{T\to 0^{+}}{\lim }\Phi \left(T\right)=\frac{af}{c_1},\underset{T\to {\left(1/b\right)}^{-}}{\lim }\Phi \left(T\right)=0.$

By the Intermediate Value Theorem, if

$0<w<\frac{af}{c_1},$

then there exists at least one $T^{\text{*}}\in \left(0,1/b\right)$ such that $\Phi \left(T^{*}\right)=w$ , confirming the

3.5. Graphical Illustration of $\Phi \left(T\right)$ —Untreated Model

To further support the analytical results derived in the previous section, we now examine the behavior of the function

$\Phi \left(T\right)=\frac{a\left(1-bT\right)}{c_1}\left[f-\frac{g{T}^2}{h+T^2+\varphi \cdot \frac{a\left(1-bT\right)}{c_1}}\left(1+\kappa I^{*}\right)+pT\right],$

which arises from the steady-state condition of the lymphocyte population in the untreated glioma-immune interaction model. This function captures the balance between glioma-induced immune stimulation and suppression, modulated by the glioma population size $T$ . The existence of a positive glioma-persistent equilibrium (GPE) depends on whether the immune recruitment rate $w$ intersects this function within the biologically feasible interval $T\in \left(0,1/b\right)$ . The function $\Phi \left(T\right)$ is continuous and strictly positive on the interval $\left(0,1/b\right)$ . At the lower bound, we observe that

$\underset{T\to 0^{+}}{\lim }\Phi \left(T\right)=\frac{af}{c_1},$

which represents the maximal immune recruitment capacity in the absence of glioma burden. As $T$ increases, the perfector $\left(1-bT\right)$ decreases linearly, reflecting the saturation of glioma growth near it carrying capacity. Simultaneously, the bracketed term in $\Phi \left(T\right)$ captures the nonlinear interplay between immune stimulation and suppression. The term $\frac{g{T}^2}{h+T^2+\varphi \cdot \frac{a\left(1-bT\right)}{c_1}}$ increases with $T$ , while the linear term $pT$ also increases. As a result, the entire bracketed expression may initially decrease due to dominant immune stimulation but eventually increases as glioma-induced suppression becomes more significant. This behavior implies that $\Phi \left(T\right)$ is unimodal: it increases near $T=0$ , reaches a peak, and then decreases toward zero as $T\to 1/b$ . The horizontal line $w$ , representing the constant rate of lymphocyte recruitment, intersects $\Phi \left(T\right)$ if and only if $w<\frac{af}{c_1}$ . In this case, the Intermediate Value Theorem guarantees the existence of at least one intersection point $T^{\text{*}}\in \left(0,1/b\right)$ such that $\Phi \left(T^{*}\right)=w$ , confirming the existence of a positive glioma-persistent equilibrium. Figure 1 illustrates this relationship. The blue curve represents the function $\Phi \left(T\right)$ , while the red dashed line corresponds to the constant immune recruitment rate $w$ . The point where the two curves intersect, marked by $T^{\text{*}}$ , indicates a steady-state glioma population that can coexist with the immune system. If such an intersection exists, the system admits a GPE, meaning the immune response is insufficient to eliminate the glioma entirely, and the glioma persists at a stable level. If no intersection occurs specifically, if $w\ge \frac{af}{c_1}$ , then the immune system is strong enough to prevent glioma persistence, and the system converges to the disease-free equilibrium. This graphical perspective provides valuable insight into the threshold dynamics of glioma persistence. It highlights the critical role of immune parameters in determining long-term outcomes and suggests that increasing the effective immune response either by enhancing cytotoxic activity or by reducing glioma-induced suppression may shift the system away from a glioma-persistent state.

Figure 1. Graph of $\Phi \left(T\right)$ and immune recruitment rate $w$ . The intersection indicates a positive glioma-persistent equilibrium.

3.6. Stability of the Glioma-Persistent Equilibrium (GPE)—Untreated Model

In this section, we investigate the local stability of the glioma-persistent equilibrium (GPE) in the untreated model. GPE corresponds to a nontrivial steady state $\left(T^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}}\right)\in {ℝ}_{>0}^4$ , where glioma cells coexist with immune components in a dynamic balance. To determine whether this equilibrium is locally asymptotically stable, we linearize the system around the GPE and analyze the eigenvalues of the Jacobian matrix evaluated at this point. Letting $\overrightarrow{x}\left(t\right)={\left[T\left(t\right),N\left(t\right),C\left(t\right),E\left(t\right)\right]}^{\top }$ , the linearized system near the GPE is given by

$D^{\alpha }\overrightarrow{x}\left(t\right)=J\left(T^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}}\right)\cdot \overrightarrow{x}\left(t\right),$

where $J$ is the Jacobian matrix evaluated at the equilibrium. The partial derivatives of the system components yield the following entries of the Jacobian matrix.

For the first equation:

$\frac{\partial F_1}{\partial T}=a\left(1-2b{T}^{\text{*}}\right)-c_1{N}^{\text{*}},\frac{\partial F_1}{\partial N}=-c_1{T}^{\text{*}}.$

For the second equation, we define the function

$\phi \left(T\right)=\frac{g{T}^2}{h+T^2},$

whose derivative with respect to $T$ is

$\frac{\text{d}\phi }{\text{d}T}=\frac{2gTh}{{\left(h+T^2\right)}^2}.$

Using this, we compute:

$\frac{\partial F_2}{\partial T}=\left(\frac{2g{T}^{\text{*}}h}{{\left(h+T^{\text{*}2}\right)}^2}-p\right)N^{\text{*}},\frac{\partial F_2}{\partial N}=-f+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}}-p{T}^{\text{*}}.$

For the third equation:

$\frac{\partial F_3}{\partial T}=-q{C}^{\text{*}}+r{N}^{\text{*}},\frac{\partial F_3}{\partial N}=r{T}^{\text{*}},\frac{\partial F_3}{\partial C}=-m-q{T}^{\text{*}}-\delta .$

For the fourth equation:

$\frac{\partial F_4}{\partial C}=\delta ,\frac{\partial F_4}{\partial E}=-\mu .$

Substituting these expressions, the Jacobian matrix evaluated at the GPE is

$J\left(T^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}}\right)=\left[\begin{array}{cccc}a\left(1-2b{T}^{\text{*}}\right)-c_1{N}^{\text{*}}& -c_1{T}^{\text{*}}& 0& 0\\ \left(\frac{2g{T}^{\text{*}}h}{{\left(h+T^{\text{*}2}\right)}^2}-p\right)N^{\text{*}}& -f+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}}-p{T}^{\text{*}}& 0& 0\\ -r{N}^{\text{*}}+q{C}^{\text{*}}& r{T}^{\text{*}}& -m-q{T}^{\text{*}}-\delta & 0\\ 0& 0& \delta & -\mu \end{array}\right].$

This matrix is lower block triangular, so its eigenvalues are the union of the eigenvalues of the two $2\times 2$ diagonal blocks. The eigenvalues of the $\left(C,E\right)$ subsystem are

${\lambda }_3=-m-q{T}^{\text{*}}-\delta ,{\lambda }_4=-\mu ,$

which are both strictly negative due to the positivity of all parameters and $T^{\text{*}}>0$ .

For the $\left(T,N\right)$ subsystem, we define the submatrix

$J_1=\left[\begin{array}{cc}a\left(1-2b{T}^{\text{*}}\right)-c_1{N}^{\text{*}}& -c_1{T}^{\text{*}}\\ \left(\frac{2g{T}^{\text{*}}h}{{\left(h+T^{\text{*}2}\right)}^2}-p\right)N^{\text{*}}& -f+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}}-p{T}^{\text{*}}\end{array}\right].$

The characteristic polynomial of $J_1$ is given by

${\lambda }^2-\text{Tr}\left(J_1\right)\lambda +\text{det}\left(J_1\right)=0,$

where

$\text{Tr}\left(J_1\right)=a\left(1-2b{T}^{\text{*}}\right)-c_1{N}^{\text{*}}-f+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}}-p{T}^{\text{*}},$

$\text{det}\left(J_1\right)=\left(a\left(1-2b{T}^{\text{*}}\right)-c_1{N}^{\text{*}}\right)\left(-f+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}}-p{T}^{\text{*}}\right)+c_1{T}^{\text{*}}\left(\frac{2g{T}^{\text{*}}h}{{\left(h+T^{\text{*}2}\right)}^2}-p\right)N^{\text{*}}.$

The eigenvalues ${\lambda }_{1,2}$ of the $\left(T,N\right)$ subsystem are

${\lambda }_{1,2}=\frac{1}{2}\left(\text{Tr}\left(J_1\right)\pm \sqrt{\text{Tr}{\left(J_1\right)}^2-4\det \left(J_1\right)}\right).$

To ensure local asymptotic stability of the GPE in the fractional-order system of Caputo type with order $\alpha \in \left(0,1\right]$ , all eigenvalues ${\lambda }_i$ must satisfy the condition

$\left|\text{arg}\left({\lambda }_i\right)\right|>\frac{\alpha \pi }{2},\forall i=1,2,3,4,$

which guarantees that the eigenvalues lie outside the sector

$\left\{\lambda \in ℂ:\left|\text{arg}\left(\lambda \right)\right|\le \frac{\alpha \pi }{2}\right\},$

the stability region for Caputo fractional-order systems. If this condition is satisfied, then the GPE is locally asymptotically stable in the sense of fractional dynamics. Otherwise, the GPE is unstable, and the system may transition to a different dynamical regime, such as glioma elimination or unbounded growth, depending on the nature of the perturbation and the global phase space structure. This completes the local stability analysis of the glioma-persistent equilibrium in the untreated model.

4. Equilibrium Analysis of the Treated Model

In this section, we analyze the equilibrium behavior of the glioma-immune interaction model under therapeutic intervention. The model incorporates both drug-induced glioma cell death and oncolytic virus therapy, as well as the dynamics of drug and virus administration, checkpoint inhibition, cytokine signaling, and immune exhaustion. A time delay $\tau$ is included to represent the lag in NK cell response to glioma presence. Our goal is to identify the disease-free equilibrium (DFE), corresponding to complete glioma eradication, and to analyze its local stability using the general stability theorem for fractional-order systems.

4.1. Disease-Free Equilibrium (DFE)—Treated Model

To determine the disease-free equilibrium, we assume that both glioma populations are completely eradicated, i.e., $T_s\left(t\right)=T_r\left(t\right)=0$ . In the absence of glial cells, the immune system is no longer stimulated, and the concentrations of drug, virus, checkpoint molecules, and cytokines decay to zero in the absence of external input. Therefore, we set

$T_s=T_r=B=M=P=0,U_M\left(t\right)=U_V\left(t\right)=U_P\left(t\right)=0.$

Substituting these into the system, we solve for the remaining variables.

From the glioma equations:

$D^{\alpha }{T}_s=0\Rightarrow T_s^{\text{*}}=0,D^{\alpha }{T}_r=0\Rightarrow T_r^{\text{*}}=0.$

From the drug and virus equations:

$D^{\alpha }M=-\gamma M+U_M=0\Rightarrow M^{\text{*}}=0,D^{\alpha }B=-\eta B+U_V=0\Rightarrow B^{\text{*}}=0.$

From the checkpoint and cytokine equations:

$D^{\alpha }P=\sigma T-\lambda P-\xi U_P=0\Rightarrow P^{\text{*}}=0,$

$D^{\alpha }I={\rho }_{1}N+{\rho }_{2}C-{\rho }_{3}I=0.$

Since $C=0$ at equilibrium (as shown below), we obtain:

$0={\rho }_{1}N-{\rho }_{3}I\Rightarrow I^{\text{*}}=\frac{{\rho }_1{N}^{\text{*}}}{{\rho }_3}.$

From the lymphocyte equation:

$D^{\alpha }N=w-fN+\frac{g{T}^2}{h+T^2+\varphi N}\left(1+\kappa I\right)N-pNT-k_{N}MN-\theta PN+ϵU_{P}N=0.$

Substituting $T=0$ , $M=0$ , $P=0$ , $U_P=0$ , we obtain:

$0=w-fN\Rightarrow N^{\text{*}}=\frac{w}{f}.$

From the effector equation:

$D^{\alpha }C=-mC+\frac{j{B}^2}{k+B^2}C-qCT+rNT-k_{C}MC-\delta C+\zeta E=0.$

Substituting $T=0$ , $B=0$ , $M=0$ , $E=0$ , we obtain:

$0=-mC-\delta C\Rightarrow C^{\text{*}}=0.$

From the exhaustion equation:

$D^{\alpha }E=\delta C-\mu E-\zeta E=0.$

Substituting $C=0$ , we obtain:

$0=-\left(\mu +\zeta \right)E\Rightarrow E^{\text{*}}=0.$

Therefore, the disease-free equilibrium of the treated model is given by

$\left(T_s^{\text{*}},T_r^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}},M^{\text{*}},B^{\text{*}},P^{\text{*}},I^{\text{*}}\right)=\left(0,0,\frac{w}{f},0,0,0,0,0,\frac{{\rho }_{1}w}{f{\rho }_3}\right).$

This equilibrium represents a glioma-free state in which the immune system has returned to its basal level, and all therapeutic and signaling agents have been cleared from the system in the absence of external input.

4.2. Stability of the Disease-Free Equilibrium—Treated Model

To analyze the local stability of the disease-free equilibrium (DFE), we linearize the treated glioma-immune interaction model around the equilibrium point

$\left(T_s^{\text{*}},T_r^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}},M^{\text{*}},B^{\text{*}},P^{\text{*}},I^{\text{*}}\right)=\left(0,0,\frac{w}{f},0,0,0,0,0,\frac{{\rho }_{1}w}{f{\rho }_3}\right).$

Let $\overrightarrow{x}\left(t\right)={\left[T_s\left(t\right),T_r\left(t\right),N\left(t\right),C\left(t\right),E\left(t\right),M\left(t\right),B\left(t\right),P\left(t\right),I\left(t\right)\right]}^{\top }$ , and define the vector field $\overrightarrow{F}\left(\overrightarrow{x}\right)$ as the right-hand side of the system. The Jacobian matrix $J$ is constructed by evaluating the partial derivatives of each component with respect to each variable at the DFE. We compute the nonzero partial derivatives as follows:

From the $T_s$ equation:

$\frac{\partial F_1}{\partial T_s}=a_s-\frac{c_{1}w}{f}-\rho -\omega ,\frac{\partial F_1}{\partial N}=-c_1{T}_s^{\text{*}}=0,\frac{\partial F_1}{\partial M}=-k_T{T}_s^{\text{*}}=0,\frac{\partial F_1}{\partial B}=-T_s^{\text{*}}=0.$

From the $T_r$ equation:

$\frac{\partial F_2}{\partial T_r}=a_r-\frac{c_{2}w}{f}-\rho ,\frac{\partial F_2}{\partial T_s}=\omega ,\frac{\partial F_2}{\partial N}=-c_2{T}_r^{\text{*}}=0,\frac{\partial F_2}{\partial B}=-T_r^{\text{*}}=0.$

From the $N$ equation:

$\begin{array}{l}\frac{\partial F_3}{\partial T_s}=-p\cdot \frac{w}{f},\frac{\partial F_3}{\partial T_r}=-p\cdot \frac{w}{f},\frac{\partial F_3}{\partial N}=-f-\theta P^{\text{*}}-k_N{M}^{\text{*}}+\frac{ϵU_P^{\text{*}}}{f}=-f,\\ \frac{\partial F_3}{\partial M}=-k_N\cdot \frac{w}{f},\frac{\partial F_3}{\partial P}=-\theta \cdot \frac{w}{f},\frac{\partial F_3}{\partial I}=\frac{g\kappa \cdot 0}{h+0+\varphi \cdot \frac{w}{f}}\cdot \frac{w}{f}=0.\end{array}$

From the $C$ equation:

$\frac{\partial F_4}{\partial C}=-m-\delta ,\frac{\partial F_4}{\partial E}=\zeta .$

From the $E$ equation:

$\frac{\partial F_5}{\partial C}=\delta ,\frac{\partial F_5}{\partial E}=-\mu -\zeta .$

From the $M$ , $B$ , $P$ , and $I$ equations:

$\frac{\partial F_6}{\partial M}=-\gamma ,\frac{\partial F_7}{\partial B}=-\eta ,\frac{\partial F_8}{\partial P}=-\lambda ,\frac{\partial F_9}{\partial I}=-{\rho }_3,\frac{\partial F_9}{\partial N}={\rho }_1.$

Substituting these into the Jacobian matrix, we obtain:

$J\left(\text{DFE}\right)=\left[\begin{array}{ccccccccc}{a}_s-\frac{c_{1}w}{f}-\rho -\omega & 0& 0& 0& 0& 0& 0& 0& 0\\ \omega & a_r-\frac{c_{2}w}{f}-\rho & 0& 0& 0& 0& 0& 0& 0\\ -\frac{pw}{f}& -\frac{pw}{f}& -f& 0& 0& -\frac{k_{N}w}{f}& 0& -\frac{\theta w}{f}& 0\\ 0& 0& 0& -m-\delta & \zeta & 0& 0& 0& 0\\ 0& 0& 0& \delta & -\mu -\zeta & 0& 0& 0& 0\\ 0& 0& 0& 0& 0& -\gamma & 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\eta & 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -\lambda & 0\\ 0& 0& {\rho }_1& {\rho }_2\cdot 0& 0& 0& 0& 0& -{\rho }_3\end{array}\right].$

The eigenvalues of this matrix are:

${\lambda }_1=a_s-\frac{c_{1}w}{f}-\rho -\omega ,$

${\lambda }_2=a_r-\frac{c_{2}w}{f}-\rho ,$

${\lambda }_3=-f,{\lambda }_4=-m-\delta ,{\lambda }_5=-\mu -\zeta ,$

${\lambda }_6=-\gamma ,{\lambda }_7=-\eta ,{\lambda }_8=-\lambda .$

The final eigenvalue ${\lambda }_9$ is obtained from the $2\times 2$ block involving $N$ and $I$ , which has the form:

$\left[\begin{array}{cc}-f& 0\\ {\rho }_1& -{\rho }_3\end{array}\right].$

The eigenvalues of this block are simply $-f$ and $-{\rho }_3$ , so we conclude:

${\lambda }_9=-{\rho }_3.$

We now apply the general stability theorem for fractional-order systems.

Theorem 4.1 (Stability of Linear Fractional-Order Systems).

Consider the linear autonomous system of Caputo fractional differential equations:

$D^{\alpha }\overrightarrow{x}\left(t\right)=A\overrightarrow{x}\left(t\right),$

where $A\in {ℝ}^{n\times n}$ is a constant matrix and $\alpha \in \left(0,1\right]$ . The equilibrium $\overrightarrow{x}=0$ is locally asymptotically stable if and only if all eigenvalues ${\lambda }_i\in ℂ$ of $A$ satisfy

$\left|\text{arg}\left({\lambda }_i\right)\right|>\frac{\alpha \pi }{2},\forall i=1,\cdots ,n.$

Applying this theorem to the Jacobian matrix $J\left(\text{DFE}\right)$ , we observe that all eigenvalues except ${\lambda }_1$ and ${\lambda }_2$ are strictly negative real numbers. Therefore, their arguments satisfy

$\left|\text{arg}\left({\lambda }_i\right)\right|=\pi >\frac{\alpha \pi }{2},\text{for}i=3,\cdots ,9.$

Thus, the local stability of the DFE depends entirely on the signs of ${\lambda }_1$ and ${\lambda }_2$ . For the DFE to be locally asymptotically stable, we require:

${\lambda }_1=a_s-\frac{c_{1}w}{f}-\rho -\omega <0,{\lambda }_2=a_r-\frac{c_{2}w}{f}-\rho <0.$

Rearranging these inequalities, we obtain the stability conditions:

$\frac{c_{1}w}{f}+\rho +\omega >a_s,\frac{c_{2}w}{f}+\rho >a_r.$

These conditions imply that the combined effects of immune-mediated cytotoxicity, saturation-induced glioma death, and phenotypic switching (in the case of therapy-sensitive cells) must be sufficiently strong to overcome the intrinsic proliferation rates of both glioma subpopulations. Specifically, the term $\frac{c_{1}w}{f}$ represents the effective killing rate of therapy-sensitive glioma cells by lymphocytes at equilibrium, while $\rho$ accounts for self-limiting glioma growth, and $\omega$ captures the rate at which sensitive cells convert into resistant ones. Similarly, $\frac{c_{2}w}{f}$ reflects the immune pressure on therapy-resistant cells. If both inequalities are satisfied, then all eigenvalues of the Jacobian matrix lie outside the critical sector defined by the fractional order $\alpha$ , and the disease-free equilibrium is locally asymptotically stable. In this case, the system will return to the glioma-free state following small perturbations, indicating successful therapeutic control of both glioma subtypes. Conversely, if either inequality fails, then at least one eigenvalue lies within or on the boundary of the instability sector $\left\{\lambda \in ℂ:\left|\text{arg}\left(\lambda \right)\right|\le \frac{\alpha \pi }{2}\right\}$ , and the DFE becomes unstable. This instability implies that the glioma population may escape immune and therapeutic control, leading to recurrence or persistence of the disease. This completes the local stability analysis of the disease-free equilibrium in the treated glioma-immune interaction model. In the next section, we will investigate the existence and stability of glioma-persistent equilibria under treatment and explore how therapeutic parameters influence long-term glioma dynamics.

4.3. Interpretation of Stability Conditions for the Disease-Free Equilibrium

The local stability analysis of the disease-free equilibrium (DFE) in the treated glioma-immune interaction model reveals critical thresholds that determine whether therapeutic intervention can successfully eradicate the glioma. Specifically, the DFE is locally asymptotically stable if and only if the following two inequalities are satisfied:

$\frac{c_{1}w}{f}+\rho +\omega >a_s,\frac{c_{2}w}{f}+\rho >a_r.$

These conditions provide biologically meaningful insights into the mechanisms required for glioma clearance.

The first inequality ensures that the combined effects of immune-mediated cytotoxicity ($\frac{c_{1}w}{f}$ ), saturation-induced glioma death ($\rho$ ), and phenotypic switching from therapy-sensitive to resistant cells ($\omega$ ) are sufficient to suppress the intrinsic growth rate $a_s$ of therapy-sensitive glioma cells. The second inequality imposes a similar requirement for therapy-resistant cells, which do not undergo phenotypic switching but are still subject to immune pressure and saturation effects. These thresholds highlight the importance of maintaining a sufficiently high lymphocyte recruitment rate $w$ , as it directly enhances the immune-mediated killing terms $\frac{c_{1}w}{f}$ and $\frac{c_{2}w}{f}$ . Additionally, increasing the saturation death rate $\rho$ through therapeutic strategies such as hypoxia induction or metabolic disruption can further contribute to glioma suppression. The parameter $\omega$ , representing the rate of transition from sensitive to resistant phenotypes, plays a dual role: while it reduces the sensitive population, it also feeds the resistant compartment, potentially undermining treatment if not balanced by sufficient immune pressure on $T_r$ . From a therapeutic perspective, these conditions suggest that successful eradication of glioma requires coordinated enhancement of immune function, effective drug and viral delivery, and suppression of resistance mechanisms. Failure to meet these thresholds may result in the persistence or regrowth of glioma, even in the presence of treatment.

4.4. Glioma-Persistent Equilibrium (GPE): Existence and Stability in the Treated Model

When the stability conditions for the disease-free equilibrium (DFE) are not satisfied, the system may admit one or more positive steady states in which glioma cells persist in coexistence with immune and therapeutic components. These glioma-persistent equilibria (GPEs) represent chronic disease states where the glioma is not eradicated but maintained at a nonzero level due to a dynamic balance between proliferation, immune response, and therapeutic effects. To investigate the existence of such equilibria, we consider steady-state solutions of the full treated model where $T_s^{\text{*}}>0$ and/or $T_r^{\text{*}}>0$ . At equilibrium, all time derivatives vanish, and the system reduces to a set of nonlinear algebraic equations:

$0=a_s{T}_s^{\text{*}}\left(1-b_s{T}_s^{\text{*}}\right)-c_1{N}^{\text{*}}{T}_s^{\text{*}}-B^{\text{*}}{T}_s^{\text{*}}-k_T{M}^{\text{*}}{T}_s^{\text{*}}-\rho T_s^{\text{*}2}-\omega T_s^{\text{*}},$

$0=a_r{T}_r^{\text{*}}\left(1-b_r{T}_r^{\text{*}}\right)-c_2{N}^{\text{*}}{T}_r^{\text{*}}-B^{\text{*}}{T}_r^{\text{*}}-\rho T_r^{\text{*}2}+\omega T_s^{\text{*}},$

$\begin{array}{l}0=w-f{N}^{\text{*}}+\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}+\varphi N^{\text{*}}}\left(1+\kappa I^{\text{*}}\right)N^{\text{*}}-p{N}^{\text{*}}{T}^{\text{*}}\\ -k_N{M}^{\text{*}}{N}^{\text{*}}-\theta P^{\text{*}}{N}^{\text{*}}+ϵU_P^{\text{*}}{N}^{\text{*}},\end{array}$

$0=-m{C}^{\text{*}}+\frac{j{B}^{\text{*}2}}{k+B^{\text{*}2}}{C}^{\text{*}}-q{C}^{\text{*}}{T}^{\text{*}}+r{N}^{\text{*}}{T}^{\text{*}}-k_C{M}^{\text{*}}{C}^{\text{*}}-\delta C^{\text{*}}+\zeta E^{\text{*}},$

$0=\delta C^{\text{*}}-\mu E^{\text{*}}-\zeta E^{\text{*}},$

$0=-\gamma M^{\text{*}}+U_M^{\text{*}},$

$0=-\eta B^{\text{*}}+U_V^{\text{*}},$

$0=\sigma T^{\text{*}}-\lambda P^{\text{*}}-\xi U_P^{\text{*}},$

$0={\rho }_1{N}^{\text{*}}+{\rho }_2{C}^{\text{*}}-{\rho }_3{I}^{\text{*}},$

where $T^{\text{*}}=T_s^{\text{*}}+T_r^{\text{*}}$ is the total glioma burden.

These equations are highly nonlinear and coupled, with feedback loops involving immune stimulation, cytokine signaling, and checkpoint inhibition. Analytical solutions are generally intractable; however, qualitative insights can be gained by examining the effective net growth rates of the glioma subpopulations:

$\begin{array}{l}{\Gamma }_s=a_s\left(1-b_s{T}_s^{*}\right)-\frac{c_{1}w}{f}-\rho T_s^{*}-\omega -\frac{B^{*}+k_T{M}^{*}}{T_s^{*}},\\ {\Gamma }_r=a_r\left(1-b_r{T}_r^{*}\right)-\frac{c_{2}w}{f}-\rho T_r^{*}-\frac{B^{*}}{T_r^{*}}+\omega \frac{T_s^{*}}{T_r^{*}}.\end{array}$

A GPE exists if ${\Gamma }_s>0$ or ${\Gamma }_r>0$ , meaning that the net proliferation of at least one glioma subpopulation exceeds the combined suppressive effects of immune and therapeutic mechanisms. The existence of such equilibria depends sensitively on the values of $w,\rho ,\omega$ , and the treatment inputs $U_M^{\text{*}},U_V^{\text{*}},U_P^{\text{*}}$ .

To assess the local stability of the GPE, we linearize the treated glioma-immune interaction model around a nontrivial steady state:

$\left(T_s^{\text{*}},T_r^{\text{*}},N^{\text{*}},C^{\text{*}},E^{\text{*}},M^{\text{*}},B^{\text{*}},P^{\text{*}},I^{\text{*}}\right)\in {ℝ}_{>0}^9,$

where at least one of the glioma populations satisfies $T_s^{\text{*}}>0$ or $T_r^{\text{*}}>0$ . Let

$\overrightarrow{x}\left(t\right)={\left[T_s\left(t\right),T_r\left(t\right),N\left(t\right),C\left(t\right),E\left(t\right),M\left(t\right),B\left(t\right),P\left(t\right),I\left(t\right)\right]}^{\top },$

and define the vector field $\overrightarrow{F}\left(\overrightarrow{x}\right)$ as the right-hand side of the system. The linearized system near the GPE is given by:

$D_t^{\alpha }\overrightarrow{x}\left(t\right)=J\left(\text{GPE}\right)\cdot \overrightarrow{x}\left(t\right),$

where $J\left(\text{GPE}\right)$ is the Jacobian matrix evaluated at the endemic equilibrium. We compute the partial derivatives of each component with respect to each variable:

- Glioma-sensitive population $T_s$ :

$\begin{array}{l}\frac{\partial F_1}{\partial T_s}=a_s\left(1-2b_s{T}_s^{\text{*}}\right)-c_1{N}^{\text{*}}-B^{\text{*}}-k_T{M}^{\text{*}}-2\rho T_s^{\text{*}}-\omega ,\\ \frac{\partial F_1}{\partial N}=-c_1{T}_s^{\text{*}},\frac{\partial F_1}{\partial M}=-k_T{T}_s^{\text{*}},\frac{\partial F_1}{\partial B}=-T_s^{\text{*}}.\end{array}$

- Glioma-resistant population $T_r$ :

$\frac{\partial F_2}{\partial T_r}=a_r\left(1-2b_r{T}_r^{\text{*}}\right)-c_2{N}^{\text{*}}-B^{\text{*}}-2\rho T_r^{\text{*}},\frac{\partial F_2}{\partial T_s}=\omega ,\frac{\partial F_2}{\partial N}=-c_2{T}_r^{\text{*}},\frac{\partial F_2}{\partial B}=-T_r^{\text{*}}.$

- Lymphocyte population $N$ : Let

$\varphi \left(T^{\text{*}},N^{\text{*}}\right)=\frac{g{T}^{\text{*}2}}{h+T^{\text{*}2}+\varphi N^{\text{*}}}\left(1+\kappa I^{\text{*}}\right),$

Then

$\frac{\partial F_3}{\partial T_s}=\left(\frac{\partial \varphi }{\partial T^{\text{*}}}\cdot \frac{\partial T^{\text{*}}}{\partial T_s}-p{N}^{\text{*}}\right)+\varphi \left(T^{\text{*}},N^{\text{*}}\right),$

$\frac{\partial F_3}{\partial T_r}=\left(\frac{\partial \varphi }{\partial T^{\text{*}}}\cdot \frac{\partial T^{\text{*}}}{\partial T_r}-p{N}^{\text{*}}\right)+\varphi \left(T^{\text{*}},N^{\text{*}}\right),$

$\frac{\partial F_3}{\partial N}=-f+\frac{\partial \varphi }{\partial N}{N}^{\text{*}}+\varphi \left(T^{\text{*}},N^{\text{*}}\right)-p{T}^{\text{*}}-k_N{M}^{\text{*}}-\theta P^{\text{*}}+ϵU_P^{\text{*}},$

$\frac{\partial F_3}{\partial M}=-k_N{N}^{\text{*}},\frac{\partial F_3}{\partial P}=-\theta N^{\text{*}},\frac{\partial F_3}{\partial I}=\frac{g{T}^{\text{*}2}\kappa }{h+T^{\text{*}2}+\varphi N^{\text{*}}}{N}^{\text{*}}.$

- CTL population $C$ :

$\begin{array}{l}\frac{\partial F_4}{\partial C}=-m+\frac{j{B}^{\text{*}2}}{k+B^{\text{*}2}}-q{T}^{\text{*}}-k_C{M}^{\text{*}}-\delta ,\\ \frac{\partial F_4}{\partial T_s}=-q{C}^{\text{*}}+r{N}^{\text{*}},\frac{\partial F_4}{\partial T_r}=-q{C}^{\text{*}}+r{N}^{\text{*}},\frac{\partial F_4}{\partial N}=r{T}^{\text{*}},\\ \frac{\partial F_4}{\partial M}=-k_C{C}^{\text{*}},\frac{\partial F_4}{\partial B}=\frac{2jk{B}^{\text{*}}{C}^{\text{*}}}{{\left(k+B^{\text{*}2}\right)}^2},\frac{\partial F_4}{\partial E}=\zeta .\end{array}$

- Exhausted CTL population $E$ :

$\frac{\partial F_5}{\partial C}=\delta ,\frac{\partial F_5}{\partial E}=-\mu -\zeta .$

- Drug and virus concentrations:

$\frac{\partial F_6}{\partial M}=-\gamma ,\frac{\partial F_7}{\partial B}=-\eta .$

- Checkpoint molecules $P$ :

$\frac{\partial F_8}{\partial T_s}=\sigma ,\frac{\partial F_8}{\partial T_r}=\sigma ,\frac{\partial F_8}{\partial P}=-\lambda .$

- Cytokines $I$ :

$\frac{\partial F_9}{\partial N}={\rho }_1,\frac{\partial F_9}{\partial C}={\rho }_2,\frac{\partial F_9}{\partial I}=-{\rho }_3.$

Let $J\left(\text{GPE}\right)$ be the Jacobian matrix composed of these partial derivatives. The local stability of the GPE is determined by the eigenvalues ${\lambda }_i$ of $J\left(\text{GPE}\right)$ . According to the fractional-order stability criterion stated in Theorem 4.1, the GPE is locally asymptotically stable if and only if all eigenvalues satisfy

$\left|\text{arg}\left({\lambda }_i\right)\right|>\frac{\alpha \pi }{2},\forall i=1,\cdots ,9.$

Due to the nonlinear and parameter-dependent structure of $J\left(\text{GPE}\right)$ , analytical expressions for the eigenvalues are generally unavailable. Therefore, numerical computation of the Jacobian and its spectrum is typically required to assess stability for specific parameter sets. If the GPE is stable, it represents a long-term outcome in which glioma persists despite therapy. If unstable, the system may transition to the DFE or exhibit complex dynamics such as oscillations, bifurcations, or biostability, depending on the global phase space structure.

4.5. Numerical Illustration of Equilibrium Dynamics

To support the analytical results derived in Sections 4.2 and 4.4, we now present numerical simulations of the treated glioma-immune interaction model under two distinct parameter regimes. These simulations illustrate the system’s long-term behavior and confirm the theoretical predictions regarding the existence and stability of the disease-free equilibrium (DFE) and the glioma-persistent equilibrium (GPE).

We consider two cases:

• Case A (DFE): The system satisfies the stability condition for the disease-free equilibrium. This corresponds to high immune efficacy and sufficient treatment strength, leading to the eradication of glioma cells.

• Case B (GPE): The system violates the DFE stability condition, resulting in the persistence of glioma at a nonzero steady state. This reflects a scenario with weaker immune response or suboptimal treatment.

To further explore the role of memory effects in the system, we simulate each case under two different fractional orders: $\alpha =1$ (classical integer-order dynamics) and $\alpha =0.5$ (strong memory effects). The fractional order captures the biological memory inherent in immune and tumor dynamics, with lower values of $\alpha$ representing slower, history-dependent responses.

Figure 2 displays the time evolution of glioma and lymphocyte populations for all four scenarios. In Case A, both $\alpha =1$ and $\alpha =0.5$ lead to glioma elimination, but the fractional-order system exhibits a delayed immune response and slower convergence to the DFE. In Case B, glioma persists in both cases, but the fractional-order model again shows a more gradual approach to the GPE, with prolonged tumor presence and delayed immune activation. These results highlight the importance of fractional-order dynamics in modeling glioma-immune interactions. They demonstrate that memory effects can significantly influence treatment outcomes, particularly in borderline or suboptimal therapeutic scenarios. Incorporating fractional behavior provides a more realistic and nuanced understanding of tumor progression and immune response over time.

4.6. Summary

In this section, we conducted a comprehensive equilibrium analysis of the treated glioma-immune interaction model incorporating immunotherapy, oncolytic virus therapy, and drug administration. The model accounts for both therapy-sensitive and therapy-resistant glioma populations, immune delay, cytokine signaling, and checkpoint inhibition, all within a fractional-order dynamical framework. We first derived the disease-free equilibrium (DFE), corresponding to complete glioma eradication, and established its local stability conditions using the fractional-order stability theorem. These conditions revealed critical thresholds involving immune efficacy, therapeutic strength, and glioma proliferation rates. When satisfied, the system converges with a glioma-free state, indicating successful treatment. Next, we investigated the existence of endemic equilibria, or glioma-persistent equilibria (GPEs), which arise when the DFE is unstable. These equilibria represent chronic disease states where glioma cells coexist with immune and therapeutic components. We derived the nonlinear algebraic system governing the GPE and discussed the biological and mathematical conditions under which such equilibria can exist.

Figure 2. Dynamics of glioma and lymphocyte populations under different fractional orders. Top: Case A (DFE) shows glioma elimination for both $\alpha =1$ and $\alpha =0.5$ , with slower convergence in the fractional case. Bottom: Case B (GPE) shows persistent glioma, with memory effects delaying immune response. Results highlight the impact of fractional dynamics on system behavior.

Finally, we analyzed the local stability of the GPE by computing the Jacobian matrix at the endemic equilibrium and applying the fractional-order stability criterion. Due to the complexity of the system, analytical expressions for the eigenvalues are generally intractable, and numerical methods are required to assess stability. The results highlight the potential for multiple long-term outcomes, including glioma persistence, depending on treatment intensity and immune dynamics. This equilibrium framework provides a foundation for designing effective therapeutic strategies. In the next section, we formulate an optimal control problem to identify time-dependent treatment protocols that minimize glioma burden while balancing treatment costs and immune preservation.

5. Optimal Control of the Treated Glioma-Immune Model

We consider the fractional-order glioma-immune interaction model with two therapeutic control inputs: the immunotherapeutic drug dosage $U_M\left(t\right)$ and the virus dosage rate $U_V\left(t\right)$ . These controls influence the dynamics of the immunotherapy agent $M\left(t\right)$ and the oncolytic virus $B\left(t\right)$ , respectively.

The goal is to determine optimal time-dependent control strategies that minimize the total glioma burden while accounting for treatment cost, drug and virus accumulation, and potential toxicity.

Let $T\left(t\right)=T_s\left(t\right)+T_r\left(t\right)$ denote the total glioma population. The control functions are assumed to be Lebesgue measurable and satisfy

$0\le U_M\left(t\right)\le U_{\text{max}},0\le U_V\left(t\right)\le V_{\text{max}},\forall t\in \left[0,t_f\right].$

The necessary conditions for optimality are derived using Pontryagin’s Maximum Principle for fractional-order systems. Let ${\lambda }_i\left(t\right)$ , for $i=1,\cdots ,9$ , be the adjoint variables corresponding to the state variables $T_s,T_r,N,C,E,M,B,P,I$ . The adjoint equations are given by

$D^{\alpha }{\lambda }_i\left(t\right)=-\frac{\partial \mathscr{H}}{\partial x_i},{\lambda }_i\left(t_f\right)=0.$

5.1. Optimal Control Problem for Immunotherapy Only

In the absence of virus therapy, we set $U_V\left(t\right)=0$ . The control variable is the immunotherapy dosage $U_M\left(t\right)$ . The objective is to minimize the total glioma population and drug accumulation while penalizing treatment intensity. The performance index is defined as

$J\left(U_M\right)={\displaystyle {\int }_0^{t_f}\left(A_1{\left(T_s\left(t\right)+T_r\left(t\right)\right)}^2+A_{2}M{\left(t\right)}^2+\frac{1}{2}{B}_1{U}_M{\left(t\right)}^2\right)\text{d}t}.$

The Hamiltonian is given by

$\begin{array}{l}\mathscr{H}=A_1{\left(T_s+T_r\right)}^2+A_2{M}^2+\frac{1}{2}{B}_1{U}_M^2\\ +{\lambda }_1\left[a_s{T}_s\left(1-b_s{T}_s\right)-c_{1}N\left(t-\tau \right)T_s-B{T}_s-k_{T}M{T}_s-\rho T_s^2-\omega T_s\right]\\ +{\lambda }_2\left[a_r{T}_r\left(1-b_r{T}_r\right)-c_{2}N\left(t-\tau \right)T_r-B{T}_r-\rho T_r^2+\omega T_s\right]\\ +{\lambda }_3\left[w-fN+\frac{g{\left(T_s+T_r\right)}^2}{h+{\left(T_s+T_r\right)}^2+\varphi N}\left(1+\kappa I\right)N-pN\left(T_s+T_r\right)-k_{N}MN-\theta PN\right]\\ +{\lambda }_4\left[-mC+\frac{j{B}^2}{k+B^2}C-qC\left(T_s+T_r\right)+rN\left(T_s+T_r\right)-k_{C}MC-\delta C+\zeta E\right]\\ +{\lambda }_5\left[\delta C-\mu E-\zeta E\right]+{\lambda }_6\left(-\gamma M+U_M\right)+{\lambda }_7\left(-\eta B\right)\\ +{\lambda }_8\left[\sigma \left(T_s+T_r\right)-\lambda P\right]+{\lambda }_9\left[{\rho }_{1}N+{\rho }_{2}C-{\rho }_{3}I\right].\end{array}$

The adjoint equations are defined by

$D^{\alpha }{\lambda }_i\left(t\right)=-\frac{\partial \mathscr{H}}{\partial x_i},{\lambda }_i\left(t_f\right)=0,\text{for}i=1,\cdots ,9.$

The optimality condition is obtained by minimizing $\mathscr{H}$ with respect to $U_M$ :

$\frac{\partial \mathscr{H}}{\partial U_M}=B_1{U}_M+{\lambda }_6=0,$

$U_M^{\text{*}}\left(t\right)=\text{min}\left(\text{max}\left(0,-\frac{{\lambda }_6\left(t\right)}{B_1}\right),U_{\text{max}}\right).$

The corresponding adjoint system is derived by differentiating the Hamiltonian with respect to each state variable. The full expressions for the adjoint equations are provided in Appendix B.1.

5.2. Optimal Control Problem for Virus Therapy Only

In the absence of immunotherapy, we set $U_M\left(t\right)=0$ . The control variable is the virus dosage $U_V\left(t\right)$ . The objective is to minimize the glioma population and virus accumulation while penalizing treatment intensity. The performance index is defined as

$J\left(U_V\right)={\displaystyle {\int }_0^{t_f}\left(A_1{\left(T_s\left(t\right)+T_r\left(t\right)\right)}^2+A_{3}B{\left(t\right)}^2+\frac{1}{2}{B}_2{U}_V{\left(t\right)}^2\right)\text{d}t}.$

The Hamiltonian $\mathscr{H}$ has the same structure as in Section 5.1, with $U_M\left(t\right)=0$ and $U_V\left(t\right)$ active. The only control-dependent term becomes

${\lambda }_7\left(-\eta B+U_V\right).$

The adjoint system is defined by

$D^{\alpha }{\lambda }_i\left(t\right)=-\frac{\partial \mathscr{H}}{\partial x_i},{\lambda }_i\left(t_f\right)=0,\text{for}i=1,\cdots ,9.$

The optimality condition is obtained by minimizing $\mathscr{H}$ with respect to $U_V$ :

$\frac{\partial \mathscr{H}}{\partial U_V}=B_2{U}_V+{\lambda }_7=0,$

$U_V^{\text{*}}\left(t\right)=\text{min}\left(\text{max}\left(0,-\frac{{\lambda }_7\left(t\right)}{B_2}\right),V_{\text{max}}\right).$

The full expressions for the adjoint equations are provided in Appendix B.2.

5.3. Optimal Control Problem for Combined Therapy

In this case, both control variables are active: the immunotherapy dosage $U_M\left(t\right)$ and the virus dosage $U_V\left(t\right)$ . The objective is to minimize the total glioma population, drug and virus accumulation, and treatment intensity over the fixed time interval $\left[0,t_f\right]$ . The performance index is defined as

$\begin{array}{l}J\left(U_M,U_V\right)={\displaystyle {\int }_0^{t_f}(A_1{\left(T_s\left(t\right)+T_r\left(t\right)\right)}^2+A_{2}M{\left(t\right)}^2+A_{3}B{\left(t\right)}^2\begin{array}{c}\\ \end{array}}\\ +\frac{1}{2}{B}_1{U}_M{\left(t\right)}^2+\frac{1}{2}{B}_2{U}_V{\left(t\right)}^2)\text{d}t.\end{array}$

The Hamiltonian $\mathscr{H}$ follows the same structure as in Section 5.1, with both control terms active. The control-dependent components of the Hamiltonian are

${\lambda }_6\left(-\gamma M+U_M\right),{\lambda }_7\left(-\eta B+U_V\right).$

The adjoint system is defined by

$D^{\alpha }{\lambda }_i\left(t\right)=-\frac{\partial \mathscr{H}}{\partial x_i},{\lambda }_i\left(t_f\right)=0,\text{for}i=1,\cdots ,9.$

The optimality conditions are obtained by minimizing $\mathscr{H}$ with respect to both control variables:

$\frac{\partial \mathscr{H}}{\partial U_M}=B_1{U}_M+{\lambda }_6=0\Rightarrow U_M^{\text{*}}\left(t\right)=\text{min}\left(\text{max}\left(0,-\frac{{\lambda }_6\left(t\right)}{B_1}\right),U_{\text{max}}\right),$

$\frac{\partial \mathscr{H}}{\partial U_V}=B_2{U}_V+{\lambda }_7=0\Rightarrow U_V^{\text{*}}\left(t\right)=\text{min}\left(\text{max}\left(0,-\frac{{\lambda }_7\left(t\right)}{B_2}\right),V_{\text{max}}\right).$

The full expressions for the adjoint equations corresponding to the combined therapy case are provided in Appendix B.3.

5.4. Numerical Implementation

To solve the optimal control problems formulated in Sections 5.1-5.3, we employ a forward-backward sweep method adapted for fractional-order systems. The optimality system consists of the state equations (governed by Caputo fractional derivatives), the adjoint equations derived via Pontryagin’s Maximum Principle, and the control characterizations.

Let $\alpha \in \left(0,1\right]$ represents the fractional order of the system. The Caputo fractional derivative of a sufficiently smooth function $x\left(t\right)$ is defined as

$D^{\alpha }x\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t\frac{x^{\prime }\left(s\right)}{{\left(t-s\right)}^{\alpha }}\text{d}s},$

where $\Gamma (\cdot )$ is the Gamma function. This non-local operator captures memory effects and long-range temporal interactions, which are essential for modeling biological processes such as immune response and glial disease progression.

To numerically approximate the Caputo derivative, we employ the L1 finite difference scheme on a uniform time grid ${\left\{t_n=nh\right\}}_{n=0}^N$ , where $h=t_f/N$ is the time steps. The L1 approximation of $D^{\alpha }x\left(t_n\right)$ is given by

$\begin{array}{l}{D}^{\alpha }x\left(t_n\right)\approx \frac{1}{h^{\alpha }\Gamma \left(2-\alpha \right)}\left[x\left(t_n\right)-x\left(t_{n-1}\right)\right]\\ +{\displaystyle {\sum }_{k=1}^{n-1}\frac{x\left(t_k\right)-x\left(t_{k-1}\right)}{h^{\alpha }}}\left[{\left(n-k+1\right)}^{1-\alpha }-2{\left(n-k\right)}^{1-\alpha }+{\left(n-k-1\right)}^{1-\alpha }\right].\end{array}$

This approximation is applied to both the state and adjoint systems. The resulting optimality system is solved iteratively using a forward-backward sweep algorithm, as described in the remainder of this section.

5.5. Model Parameters

The parameters employed in the treated glioma-immune model and the associated optimal control simulations are selected based on relevant biological literature, estimated from experimental data, or chosen to ensure biologically realistic and numerically stable dynamics. These parameters include glial cell proliferation rates, immune interaction coefficients, pharmacokinetic constants, and weights in the objective functional. All numerical experiments presented in Section 6 are conducted using these baseline values unless otherwise stated. The complete list of parameter definitions and values is provided in Appendix C.

6. Numerical Results and Discussion

In In This section presents the numerical results obtained from the optimal control problems formulated in Section 5 and solved using the forward-backward sweep method described in Section 5.4. All simulations are based on the treated model introduced in Section 2.2, which incorporates fractional-order dynamics, immune response delays, and therapeutic interventions including immunotherapy, virus therapy, and checkpoint inhibition. Parameter values used in all simulations are listed in Appendix C.

The objective is to evaluate the therapeutic impact of different treatment strategies on glioma suppression and immune activation. Simulations are conducted over a fixed time horizon $\left[0,50\right]$ , with the fractional order $\alpha \in \left(0,1\right]$ capturing memory-dependent biological behavior. The computed optimal control functions $U_M\left(t\right)$ , $U_V\left(t\right)$ , and $U_P\left(t\right)$ are used to generate the corresponding state trajectories.

6.1. Untreated Dynamics

Figure 3 shows the untreated glial-immune dynamics for $\alpha =1$ and $\alpha =0.5$ in subplots (a) and (b), respectively. The model tracks therapy-sensitive glial cells $T_s\left(t\right)$ , resistant glial cells $T_r\left(t\right)$ , natural killer cells $N\left(t\right)$ , cytotoxic T cells $C\left(t\right)$ , and checkpoint molecules $P\left(t\right)$ . In the integer-order case, the immune system responds more promptly, partially containing glioma growth. In contrast,

Figure 3. (a) and (b) show the untreated glial-immune dynamics for $\alpha =1$ and $\alpha =0.5$ , respectively. The model tracks therapy-sensitive glial cells $T_s\left(t\right)$ , resistant glial cells $T_r\left(t\right)$ , natural killer cells $N\left(t\right)$ , cytotoxic T cells $C\left(t\right)$ , and checkpoint molecules $P\left(t\right)$ . Lower $\alpha$ values delay immune activation, leading to prolonged glial persistence and higher checkpoint accumulation.

the fractional-order model exhibits delayed immune activation, leading to prolonged glial persistence and increased checkpoint accumulation. This highlights the importance of memory effects in glioma-immune interactions.

6.2. Immunotherapy Only

Figure 4 illustrates the system dynamics under optimal immunotherapy alone for $\alpha =1$ and $\alpha =0.5$ . The control $U_M\left(t\right)$ enhances immune activation, increasing cytotoxic T cells $C\left(t\right)$ and cytokine levels $I\left(t\right)$ . However, therapy-resistant glial cells $T_r\left(t\right)$ and checkpoint molecules $P\left(t\right)$ remain elevated. The fractional-order

Figure 4. (a) and (b) show glial-immune dynamics under immunotherapy alone for $\alpha =1$ and $\alpha =0.5$ , respectively. Immunotherapy reduces therapy-sensitive glial cells $T_s\left(t\right)$ and increases cytokine levels $I\left(t\right)$ , but resistant cells $T_r\left(t\right)$ and checkpoint molecules $P\left(t\right)$ persist. Lower $\alpha$ slows immune activation, limiting treatment effectiveness.

case shows a slower immune response, resulting in limited glioma suppression. These results suggest that immunotherapy alone may not be sufficient for complete control of glial malignancies, especially under memory-affected dynamics.

6.3. Combined Therapy

Figure 5 presents the system’s behavior under combined immunotherapy and virus therapy. This approach leads to a more substantial reduction in both $T_s\left(t\right)$ and $T_r\left(t\right)$ , along with stronger activation of immune components. The synergy between immune stimulation and viral cytotoxicity enhances treatment efficacy in both integer and fractional-order cases, although the latter still shows a slightly delayed response.

Figure 5. (a) and (b) show the effects of triple therapy drug, virus, and checkpoint inhibitor on glial-immune dynamics for $\alpha =1$ and $\alpha =0.5$ , respectively. Both glial populations $T_s\left(t\right)$ and $T_r\left(t\right)$ are strongly suppressed. Checkpoint inhibition reduces $P\left(t\right)$ , enabling immune recovery and enhanced cytokine activation $I\left(t\right)$ . Lower $\alpha$ slows immune response, but the combined therapy remains effective.

6.4. Optimal Control Profiles

Figure 6 displays the optimal control inputs over a 50-day treatment period. Drug therapy $U_M\left(t\right)$ (red) ramps up and maintains moderate intensity, virus therapy $U_V\left(t\right)$ (blue) is administered in bursts, and checkpoint inhibition $U_P\left(t\right)$ (green) peaks early and gradually declines. Markers highlight key transitions in each control. This adaptive strategy balances therapeutic effectiveness with treatment cost and potential toxicity.

Figure 6. Optimized control profiles for drug $U_M\left(t\right)$ , virus $U_V\left(t\right)$ , and checkpoint inhibitor $U_P\left(t\right)$ . Drug therapy (red) ramps up and sustains moderate intensity, virus therapy (blue) is applied in bursts, and checkpoint inhibition (green) is strongest early and tapers off. Markers highlight key points in each control strategy.

6.5. System Dynamics Under Optimal Control

Figure 7 shows the full system dynamics under the optimal control strategy for $\alpha =1$ . The triple therapy leads to rapid and sustained suppression of both glial populations $T_s\left(t\right)$ and $T_r\left(t\right)$ , a sharp reduction in checkpoint molecules $P\left(t\right)$ , and robust immune activation through increased $N\left(t\right)$ , $C\left(t\right)$ , and $I\left(t\right)$ . The system stabilizes toward a glioma-free state, demonstrating the effectiveness of the optimized treatment plan.

Figure 7. System dynamics under optimal control for $\alpha =1$ . Triple therapy leads to strong suppression of both glial populations $T_s\left(t\right)$ , $T_r\left(t\right)$ , rapid checkpoint reduction $P\left(t\right)$ , and sustained immune activation via $N\left(t\right)$ ,$C\left(t\right)$ , and cytokines $I\left(t\right)$ .

6.6. Comparative Glioma Dynamics

Figure 8 compares the dynamics of $T_s\left(t\right)$ and $T_r\left(t\right)$ under four treatment strategies: no treatment, immunotherapy only, triple therapy, and optimal control. The untreated case results in unchecked glial proliferation. Immunotherapy and triple therapy slow disease progression but fail to fully suppress the glioma. Only the optimal control strategy achieves rapid and sustained reduction of both glial cell types, confirming the superiority of coordinated, multi-modal treatment.

6.7. Comparative Immune Dynamics

Figure 9 presents the immune response under the same four treatment strategies. Subplots (a)(d) show the dynamics of $N\left(t\right)$ , $C\left(t\right)$ , $I\left(t\right)$ , and $E\left(t\right)$ , respectively. Optimal control (green) consistently produces the strongest immune activation and lowest exhaustion levels. In contrast, untreated and immunotherapy-only cases show weaker responses. These results emphasize the importance of combining therapies to fully engage the immune system and achieve effective glioma suppression.

Figure 8. Comparison of therapy-sensitive $T_s\left(t\right)$ and resistant $T_r\left(t\right)$ glial cell dynamics under different treatment strategies. Optimal control achieves the most effective suppression of both cell types, while untreated and immunotherapy-only scenarios allow continued growth.

Figure 9. Comparative immune dynamics under four treatment strategies: untreated (gray), immunotherapy only (red), triple therapy (blue), and optimal control (green). Each subplot displays a key immune component over time: (a) $N\left(t\right)$ natural killer (NK) cells, (b) $C\left(t\right)$ cytotoxic T lymphocytes (CTLs), (c) $I\left(t\right)$ cytokines, and (d) $E\left(t\right)$ exhausted CTLs. Optimal control leads to the strongest immune activation and lowest exhaustion, while untreated and immunotherapy-only scenarios show weaker responses.

7. Conclusions

In this study, we developed and analyzed a fractional-order mathematical model of glioma-immune interactions incorporating immune response delays and multiple therapeutic strategies, including immunotherapy, virus therapy, and checkpoint inhibition. By formulating and solving an optimal control problem, we evaluated the effectiveness of different treatment approaches in suppressing glioma growth and enhancing immune activation.

Our numerical simulations revealed several key insights:

• Memory effects matter: Lower fractional orders ($\alpha <1$ ) significantly delay immune responses, leading to prolonged glioma persistence and reduced treatment efficacy. This underscores the importance of incorporating fractional dynamics when modeling biological systems with memory.

• Immunotherapy alone is limited: While immunotherapy enhances immune activity, it fails to fully eliminate glioma cells, particularly in the presence of resistant populations and immune exhaustion.

• Combination therapy is more effective: The joint application of immunotherapy and virus therapy produces synergistic effects, resulting in stronger immune responses and more substantial glioma suppression.

• Optimal control enhances outcomes: The optimized treatment strategy dynamically adjusts drug, virus, and checkpoint inhibitor administration over a 50-day period. This approach achieves the most effective glioma reduction while maintaining immune activation and minimizing treatment intensity.

• Immune dynamics are treatment-sensitive: Comparative analysis of immune components under different strategies highlights the critical role of coordinated therapy in achieving robust and sustained immune engagement.

Overall, our results demonstrate that integrating fractional-order modeling with optimal control provides a powerful framework for designing and evaluating personalized glioma treatment strategies. Future work may extend this model to include spatial dynamics, patient-specific parameters, and real clinical data to further enhance its predictive and translational potential.

Appendices

Appendix A. Full Model Equations

The fractional-order model describing the glioma-immune dynamics under therapeutic interventions is given by the following system of Caputo fractional differential equations:

$D_t^{\alpha }{T}_s\left(t\right)=f_1\left(T_s,T_r,N,C,I,E,M,B,P\right),$

$D_t^{\alpha }{T}_r\left(t\right)=f_2\left(T_s,T_r,N,C,I,E,M,B,P\right),$

$D_t^{\alpha }N\left(t\right)=f_3\left(T_s,T_r,N,C,I,E,M,B,P\right),$

$D_t^{\alpha }C\left(t\right)=f_4\left(T_s,T_r,N,C,I,E,M,B,P\right),$

$D_t^{\alpha }E\left(t\right)=f_5\left(T_s,T_r,N,C,I,E,M,B,P\right),$

$D_t^{\alpha }M\left(t\right)=f_6\left(U_M\left(t\right)\right),$

$D_t^{\alpha }B\left(t\right)=f_7\left(U_V\left(t\right)\right),$

$D_t^{\alpha }P\left(t\right)=f_8\left(T_s,T_r,C,P,U_P\left(t\right)\right),$

$D_t^{\alpha }I\left(t\right)=f_9\left(C,E\right).$

Here, $D_t^{\alpha }$ denotes the Caputo fractional derivative of order $\alpha \in \left(0,1\right]$ , and the functions $f_i$ represent nonlinear interactions as defined in the main text.

Appendix B. Adjoint Systems for the Treated Model

The adjoint system is derived from the Hamiltonian $H$ using Pontryagin’s Maximum Principle. For the state vector:

$x\left(t\right)={\left[T_s\left(t\right),T_r\left(t\right),N\left(t\right),C\left(t\right),E\left(t\right),M\left(t\right),B\left(t\right),P\left(t\right),I\left(t\right)\right]}^{\top },$

the corresponding adjoint vector is:

$\lambda \left(t\right)=\left[{\lambda }_1\left(t\right),{\lambda }_2\left(t\right),{\lambda }_3\left(t\right),{\lambda }_4\left(t\right),{\lambda }_5\left(t\right),{\lambda }_6\left(t\right),{\lambda }_7\left(t\right),{\lambda }_8\left(t\right),{\lambda }_9\left(t\right)\right]{.}^{\top }$

Each adjoint equation satisfies:

$D_t^{\alpha }{\lambda }_i\left(t\right)=-\frac{\partial H}{\partial x_i},{\lambda }_i\left(t_f\right)=0,\text{for}i=1,\cdots ,9.$

We present the adjoint systems for each of the three treatment scenarios considered in Sections 5.1-5.3.

B.1. Adjoint System for Immunotherapy Only

Only the immunotherapy control $U_M\left(t\right)$ is active. Virus therapy $U_V\left(t\right)$ and checkpoint inhibition $U_P\left(t\right)$ are set to zero. The adjoint system simplifies by removing terms involving $B\left(t\right)$ , $P\left(t\right)$ , and their controls.

B.2. Adjoint System for Virus Therapy Only

Only virus therapy control $U_V\left(t\right)$ is active. Immunotherapy $U_M\left(t\right)$ and checkpoint inhibition $U_P\left(t\right)$ are inactive. The adjoint system excludes terms involving $M\left(t\right)$ , $P\left(t\right)$ , and their controls.

B.3. Adjoint System for Combined Therapy

This scenario represents the fully treated model, where all three controls immunotherapy $U_M\left(t\right)$ , virus therapy $U_V\left(t\right)$ , and checkpoint inhibition $U_P\left(t\right)$ are active. The complete adjoint system is given by:

$\begin{array}{c}{D}_t^{\alpha }{\lambda }_1=-[2A_{1}T+{\lambda }_1\left(a_s\left(1-2b_s{T}_s\right)-c_{1}N\left(t-\tau \right)-B-k_{T}M-2\rho T_s-\omega \right)\\ +{\lambda }_2\omega +{\lambda }_3\left(\frac{2gT\left(1+\kappa I\right)N}{h+T^2+\varphi N}-\frac{2g{T}^3\left(1+\kappa I\right)N}{{\left(h+T^2+\varphi N\right)}^2}-pN\right)\\ +{\lambda }_4\left(-qC+rN\right)+{\lambda }_8\sigma ],\end{array}$

$\begin{array}{c}{D}_t^{\alpha }{\lambda }_2=-[2A_{1}T+{\lambda }_2\left(a_r\left(1-2b_r{T}_r\right)-c_{2}N\left(t-\tau \right)-B-2\rho T_r\right)\\ +{\lambda }_3\left(\frac{2gT\left(1+\kappa I\right)N}{h+T^2+\varphi N}-\frac{2g{T}^3\left(1+\kappa I\right)N}{{\left(h+T^2+\varphi N\right)}^2}-pN\right)\\ +{\lambda }_4\left(-qC+rN\right)+{\lambda }_8\sigma ],\end{array}$

$\begin{array}{c}{D}_t^{\alpha }{\lambda }_3=-[{\lambda }_3(-f+\frac{g{T}^2\left(1+\kappa I\right)\left(h+T^2+\varphi N\right)-g{T}^2\varphi \left(1+\kappa I\right)N}{{\left(h+T^2+\varphi N\right)}^2}\\ \begin{array}{c}\stackrel{}{}\\ \end{array}-pT-k_{N}M-\theta P)+{\lambda }_{4}rT+{\lambda }_9{\rho }_1],\end{array}$

$D_t^{\alpha }{\lambda }_4=-\left[{\lambda }_4\left(-m+\frac{j{B}^2}{k+B^2}-qT-k_{C}M-\delta \right)+{\lambda }_5\delta +{\lambda }_9{\rho }_2\right],$

$D_t^{\alpha }{\lambda }_5=-\left[{\lambda }_4\zeta -{\lambda }_5\left(\mu +\zeta \right)\right],$

$D_t^{\alpha }{\lambda }_6=-\left[2A_{2}M-{\lambda }_1{k}_T{T}_s-{\lambda }_3{k}_{N}N-{\lambda }_4{k}_{C}C+{\lambda }_6\gamma \right],$

$D_t^{\alpha }{\lambda }_7=-\left[2A_{3}B-{\lambda }_1{T}_s-{\lambda }_2{T}_r+{\lambda }_4\left(\frac{2jkBC}{{\left(k+B^2\right)}^2}\right)+{\lambda }_7\eta \right],$

$D_t^{\alpha }{\lambda }_8=-\left[-1-{\lambda }_3\theta N-{\lambda }_8\lambda \right],$

$D_t^{\alpha }{\lambda }_9=-\left[{\lambda }_3\left(\frac{g{T}^2\kappa N}{h+T^2+\varphi N}\right)-{\lambda }_9{\rho }_3\right].$

Appendix C. Parameter Values and Initial Conditions

The table below lists the parameters used in the simulations, along with their descriptions, values, and units.

Parameter	Description	Value	Unit
$\alpha$	Fractional order of the system	0.5 - 1.0	-
$T_s\left(0\right)$	Initial therapy-sensitive glioma cells	0.5	dimensionless
$T_r\left(0\right)$	Initial therapy-resistant glioma cells	0.2	dimensionless
$N\left(0\right)$	Initial NK cell population	1.0	dimensionless
$C\left(0\right)$	Initial CTL population	0.1	dimensionless
$I\left(0\right)$	Initial cytokine level	0.05	dimensionless
$E\left(0\right)$	Initial exhausted CTLs	0.01	dimensionless
$M\left(0\right)$	Initial drug concentration	0.0	dimensionless
$B\left(0\right)$	Initial virus concentration	0.0	dimensionless
$P\left(0\right)$	Initial checkpoint molecule level	0.1	dimensionless
${\beta }_1$	Glioma proliferation rate	0.3	${\text{day}}^{-1}$
${\delta }_1$	Immune-induced glioma death rate	0.1	${\text{day}}^{-1}$
$\gamma$	Cytokine production rate	0.05	${\text{day}}^{-1}$
$\mu$	Natural decay rate of immune cells	0.01	${\text{day}}^{-1}$
$\eta$	Checkpoint inhibition efficacy	0.2	dimensionless
...	...	...	...

Note: Additional parameters and their sources are provided in the main text or cited literature.