EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6804 ISSN 1307-5543 – ejpam.com Published by New York Business Global Fractional Discrete-Time Modeling and Analysis of Oncolytic Adenovirus Therapy with Tumor-Specific Immune Response Amal T. Alshammari1,2, Normah Maan1,∗, Mahmoud A. M. Abdelaziz3 1 Department of Mathematical Sciences, Faculty of Science, Universiti Teknologi Malaysia 2 Department of Mathematics, Faculty of Science, University of Hafr Al Batin, Saudi Arabia 3 Department of Mathematics, Faculty of Arts and Sciences, Najran University, Najran, Saudi Arabia Abstract. Oncolytic viruses (OVs) are garnering increasing attention for their ability to directly target malignant cells while simultaneously stimulating the immune response against cancer. This study presents a novel discrete-time fractional-order mathematical framework to investigate the dynamics of oncolytic adenovirus therapy in conjunction with tumor-specific immune responses. The model captures the intricate interactions between viral infection processes and the immune system’s role in modulating tumor progression. To assess the effectiveness of oncolytic viral ther- apy, local stability and bifurcation analyses are conducted at the model’s equilibrium points. A set of local bifurcations is examined, and the necessary and sufficient conditions for detecting these bifurcations are derived using an algebraic criterion method. Numerical simulations support the theoretical results, indicating that increasing the viral infection rate and carefully managing time steps with immune response can achieve stable and tumor-suppressive outcomes. Given the limi- tations of achieving complete tumor eradication through genetically modified adenovirus therapy alone, this study explores the application of chaos control strategies to maintain the stability of the system dynamics. 2020 Mathematics Subject Classifications: 92D25, 92C50, 34A08, 37N25 Key Words and Phrases: Oncolytic virotherapy, fractional-order model, discrete-time dynamics, bifurcation, chaos control 1. Introduction Oncolytic virotherapy is rapidly emerging as a promising treatment for cancer. This innovative approach harnesses the power of oncolytic viruses to combat malignant cells. One defining feature of oncolytic viruses is their ability to selectively infect and replicate ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6804 Email addresses: theyab@graduate.utm.my (A. T. Alshammari), normahmaan@utm.my (N. Maan), maabdelaziz@nu.edu.sa (M. A. M. Abdelaziz) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 2 of 24 within tumor cells, either due to their natural properties or through genetic modifications. This targeted replication damages tumor cells while sparing healthy tissue [1, 2]. Clin- ical and experimental studies demonstrate significant progress in developing genetically engineered cancer-targeting viruses [3, 4]. Currently, a considerable number of oncolytic viruses derived from over ten types of viral vectors are undergoing clinical trials at various stages [5]. In addition to their direct tumor-destructive capabilities, oncolytic viruses can also induce cancer cell death by activating immune pathways and angiogenesis. The main challenge associated with this therapeutic strategy lies in the potential neutralization of viruses by pre-existing antibodies or antiviral immune responses. The interaction between oncolytic virotherapy and the immune system remains incom- pletely understood and is a focus of ongoing research. Studies on oncolytic viruses examine two types of immune responses: virus-specific, which blocks infection, and tumor-specific, which reflects the body’s reaction to tumor presence [6, 7]. Oncolytic viruses can be ge- netically engineered to selectively infect cancer cells, replicating until cell rupture (lysis) occurs and releasing new viral particles to infect neighboring cells. However, because of the antiviral immune response, the circulation time of the virus in the bloodstream is limited, challenging the sustainability of the lysis process. To overcome this, it is essential to de- sign viruses that can evade immune attack or bypass tumor immune-evasion mechanisms, thereby leveraging both the immune system’s benefits and targeted viral modifications for effective therapy. Adenovirus (Ad) is one of the most widely studied oncolytic viruses due to its po- tent ability to lyse cancer cells and stimulate immune responses. Its strong immuno- genicity enables it to rejuvenate antitumor immunity in cancer patients and disrupt the immunosuppressive tumor microenvironment [8]. In infected cells, oncolytic adenoviruses induce immunogenic cell death, which is essential for initiating adaptive antitumor re- sponses and establishing immune memory [9]. Surface modification of Ad with polymers extends circulation time and enhances tumor targeting. Thavasyappan et al. [10] classified polymer-modified oncolytic adenoviruses (OAds), highlighting their potential for systemic delivery and sustained tumor-specific lysis. For example, masking Ad capsids with non- immunogenic polymers, such as PEGylation, reduces innate immune responses and lowers plasma IL-6 levels by 95% within 6 hours after intravenous injection [11]. Clinical tri- als have shown that adenoviral therapy is safe, though often insufficient as a stand-alone treatment. Chemical conjugation of polymers with OAds enhances therapeutic stability, protects against immune recognition, and prevents antibody neutralization, thereby op- timizing therapeutic efficacy [12]. These properties allow adenoviruses to revitalize the antitumor immune response, making them powerful candidates for integration into can- cer therapy. Mathematical modeling provides an effective tool to guide the translation of these advances into clinical practice. The interactions among oncolytic viruses, immune responses, and the tumor microen- vironment are complex. Mathematical models are powerful instruments for elucidating A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 3 of 24 these interactions and for designing effective treatment strategies. Different modeling ap- proaches allow the identification of novel dynamical behaviors, leading to more realistic representations of tumor–virus–immune dynamics. Over the past two decades, several mathematical models have been developed to study oncolytic virotherapy. Some models have been formulated using systems of ordinary differential equations (ODEs) [13–16]. For instance, [13] examined three simple ODE models to explore tumor–immune interac- tions. Later, Ashyani et al. [17] extended one of these models by incorporating virus- and tumor-induced immune responses into a single variable. Notably, the models in [13, 17] did not include the free virus population, potentially providing an incomplete description of virotherapy dynamics. Phan and Tian [18] addressed this by introducing a state vari- able for the free virus population, enabling the study of how innate immune responses affect infected cancer cells and viral populations. In 2020, Al-Tuwairqi et al. [19] further advanced this model by adding parameters for immune activation and tumor eradication mediated by cytokine release from natural killer cells. In 2021, Nono et al. [20] modeled the immune system as a single effector population representing the host immune response. In both cases [19, 20], the free virus remained exposed to immune activity, influencing treatment outcomes. Parallel work has also explored systems of partial differential equa- tions (PDEs), incorporating spatiotemporal tumor distribution. Friedman et al. [21] proposed a PDE model of virotherapy with host immunity, focusing exclusively on innate immune responses. The model in [19] was later extended with diffusion terms for viral density to identify optimal treatment strategies [22]. More recently, Aljahdaly et al. [23] reformulated this PDE framework to investigate the interplay between naive and activated immune components. In this paper, we develop a new system of ODEs to model the interactions between oncolytic adenoviruses, cancer cells, and tumor-specific immune responses. We then de- rive a fractional-order discrete-time version of this model to capture additional biological features. Our analysis focuses on three aspects: (i) the impact of adenovirus and im- mune responses on cancer cell elimination; (ii) the role of biological memory, incorporated through fractional-order dynamics; and (iii) the discrete-time nature of the system, which allows realistic evaluation of treatment intervals and captures complex dynamical patterns. Because biological data collection is often discontinuous, discrete-time models provide a practical and accurate framework. By integrating modified viruses with adaptive immu- nity, our model offers insights into combining virotherapy with immunotherapy. In [24], a mathematical model of oncolytic virotherapy was proposed that focused on viral infection dynamics and tumor reduction without immune system involvement. That study demonstrated that a slight increase in the basic reproduction number R0 could induce chaotic tumor growth near the virus-free equilibrium. While viral therapy alone can reduce tumor burden, complete eradication requires specific conditions. Incorporat- ing immune responses changes tumor dynamics, influencing both treatment stability and viral spread. Motivated by this, we extend the model in [24] by including immune re- sponses and considering modified adenoviruses designed for controlled viral spread and A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 4 of 24 tumor lysis before immune clearance. Our model examines adenovirus dynamics under adaptive immune attack, with equilibrium analysis aligned with clinical findings in [10, 12]. The remainder of this paper is organized as follows: Section (2) introduces the mathe- matical model. Section (3) derives the fractional-order discrete-time form of the oncolytic virus–tumor–immune system. Section (4) presents equilibrium points and stability analy- sis. Section (5) discusses bifurcation analysis. Section (6) provides numerical simulations supporting the theoretical findings. Section (7) develops control strategies to regulate chaotic dynamics. Finally, Section (8) offers the conclusions. 2. Mathematical Model Figure (1) illustrates the mathematical model describing the interactions among tumor cells, adaptive immune responses, and the early phase of genetically modified adenovirus propagation within tumor populations. The model consists of ordinary differential equa- tions that capture the selective targeting of tumor cells by viral particles together with the tumor-specific immune response. Tumor regression occurs through two mechanisms: (i) di- rect tumor cell death due to viral replication and (ii) immune response stimulation through immunogenic cell death [25]. Oncolytic virotherapy initiates the antitumor immune re- sponse by presenting tumor-associated antigens and promoting immune cell infiltration. The model is designed to predict: (1) the optimal tumor specificity of oncolytically mod- ified adenovirus for maximum tumor reduction, (2) the impact of viral propagation on adaptive immune responses, and (3) the tumor’s overall response to adenovirus infection. As noted earlier, oncolytic viruses are generally associated with two distinct types of immune responses: virus-specific and tumor-specific. Virus-specific immunity can hin- der the effectiveness of viral therapy by blocking infection. To overcome this challenge, adenoviruses can be modified to evade immune recognition by masking their surface pro- teins, while simultaneously enhancing their ability to specifically target tumors [26, 27]. This modification allows the virus to avoid immune clearance while preserving its binding capacity to cellular receptors. A distinctive feature of the proposed model is that the therapeutic adenovirus works synergistically with adaptive immune cells to target and kill cancer cells, without immune interference impeding viral spread within the tumor microenvironment. 2.1. Model Assumptions The model’s biological assumptions, derived from the preceding discussion and the scientific literature, are as follows: • The system is divided into four populations: uninfected tumor cells U(t), infected tumor cells I(t), free virus particles V (t), and responsive antitumor immune cells M(t). A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 5 of 24 Figure 1: Interaction between immune cells and oncolytically modified adenovirus within tumor cells. • The term rU(t) (1− (U(t) + I(t))/k) represents the logistic growth rate of the unin- fected cancer cell population U(t), reflecting a biologically realistic scenario in which tumor growth slows as tumor burden increases. Here, k is the carrying capacity [28, 29]. • In the viral treatment process, uninfected tumor cells are assumed to proliferate more rapidly than infected tumor cells due to the shorter lifespan of the latter. Therefore, logistic growth is applied only to uninfected cells [30]. • Antitumor immune cells are assumed to consist primarily of CD8+ T cells, which recognize and eliminate both infected and uninfected cancer cells by detecting tumor- associated antigens [30]. • Virus-specific immunity is assumed to be absent during oncolytic virotherapy, as modified viral capsid or envelope proteins evade immune recognition while main- taining receptor binding. Thus, only active tumors stimulate an immune response. • Antitumor immunity is represented by a positive, nonlinear, increasing, and concave growth term of the form ρM(t)U(t) ω+U(t) , where ρ is the immune recruitment rate and ω is the immune threshold parameter inversely proportional to the steepness of the immune response curve [31]. • The immune response term models controlled immune cell proliferation, preventing uncontrolled population growth [31]. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 6 of 24 • Immune cells die naturally at a rate µ. Their natural turnover is generally higher than the loss due to interactions, as immune cells are continuously generated and eliminated [32]. • Finally, upon successful lysis of an infected tumor cell, a burst of newly produced virus particles is released, which can infect neighboring uninfected cells. 2.2. Equations of the Model Based on these assumptions, we propose the following nonlinear dynamical model of oncolytic virotherapy with tumor-specific immune response: dU(t) dt = rU(t) ( 1− U(t) + I(t) k ) − βU(t)V (t)− ηM(t)U(t), dI(t) dt = βU(t)V (t)− δI(t)− ηM(t)I(t), dV (t) dt = bδI(t)− γV (t), (1) dM(t) dt = ρM(t)U(t) ω + U(t) − µM(t). Here, U(t), I(t), V (t), and M(t) denote the concentrations of uninfected tumor cells, infected tumor cells, free adenovirus particles, and tumor-specific immune cells, respec- tively. In the first equation, r is the tumor growth rate, and k is the carrying capacity of tumor cells. Tumor cells are infected by free virus particles V (t) at an infection rate β. In the second equation, the term βU(t)V (t) denotes the infection of tumor cells, while δI(t) represents virus-induced lysis of infected cells. The parameter η denotes the rate at which immune cells eradicate both uninfected and infected tumor cells. In the third equation, δ is the death rate of infected cells, and γ is the clearance rate of free virus particles due to non-specific binding or defective particle formation. The parameter b represents the virus burst size, i.e., the number of new virus particles released per lysed cancer cell. Finally, the fourth equation describes the adaptive antitumor immune response. The Michaelis–Menten term models saturation in immune cell proliferation, with ρ as the recruitment rate and ω as the immune threshold parameter, while µ is the natural death rate of immune cells [33, 34]. 3. Fractional-Order Oncolytic Virus–Tumor–Immune Model in Discrete-Time In oncology, the challenge of cancer continues to drive advancements in conventional therapies as well as the development of novel approaches. This aggressive disease, known for its ability to metastasize to distant organs, requires deeper understanding for more effective treatments. Fractional-order dynamical systems provide a powerful framework A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 7 of 24 for capturing biological effects that are often missed in classical integer-order models. While integer-order derivatives can describe systems with predictable behavior, they are less effective in cases involving uncertainty, long-term memory, or nonlocal interactions, which are common in real biological processes. In such cases, nonlocal operators that account for memory and power-law effects are more appropriate. Given the complexities of modeling oncolytic virotherapy and the tumor-specific im- mune response, our proposed fractional-order model, based on Caputo’s definition [35], offers a suitable framework for capturing the dynamics of adenovirus-based therapy: DαU(t) = rU(t) ( 1− U(t) + I(t) k ) − βU(t)V (t)− ηM(t)U(t), DαI(t) = βU(t)V (t)− δI(t)− ηM(t)I(t), DαV (t) = bδI(t)− γV (t), (2) DαM(t) = ρM(t)U(t) ω + U(t) − µM(t). Here, Dα = dα dtα represents the Caputo fractional derivative of order α, where 0 < α ≤ 1 and t > 0. For large populations, discrete-time models often provide a more practical and realistic framework than continuous ones [36, 37]. This is especially relevant in cancer therapy, where treatment interventions and new tumor growth occur at distinct intervals, as in viral oncology therapy. To extend our analysis, we discretize the fractional-order model (2), thereby enhancing its suitability for numerical simulations. Several discretization methods exist, including Euler, Runge–Kutta, predictor–corrector, and nonstandard finite difference techniques. Elsayed et al. [38] proposed the piece- wise constant arguments approximation, which generalizes the Euler approach for nonlin- ear discrete-time models. Following this methodology, we discretize system (2), defining Un = U(n), In = I(n), Vn = V (n), and Mn = M(n) for n ≥ 0. The resulting discrete-time fractional-order model is: Un+1(t) = Un + s α Γ(1 + α) [ rUn ( 1− Un + In k ) − βUnVn − ηMnUn ] , In+1(t) = In + s α Γ(1 + α) [βUnVn − δIn − ηMnIn] , Vn+1(t) = Vn + s α Γ(1 + α) [bδIn − γVn] , (3) Mn+1(t) = Mn + s α Γ(1 + α) [ ρMnUn ω + Un − µMn ] . Here, s > 0 represents the time step size. The initial conditions are U0 > 0, I0 > 0, V0 > 0, and M0 > 0. The discrete fractional-order model (3) introduces two additional A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 8 of 24 parameters not present in the original ODE system: the fractional-order parameter α and the time step size s. These new parameters can lead to richer and more complex dynamical behaviors that are not captured by the classical model. Notably, as α → 1 in (3), the Euler discretization of system (1) is recovered. 4. Equilibria and Stability This section investigates the existence and stability of equilibrium points in the dis- cretized fractional-order model (3). Consider an equilibrium point (U∗, I∗, V ∗,M∗) of (3), obtained by setting the right-hand sides to zero: rU∗ ( 1− U∗ + I∗ k ) − βU∗V ∗ − ηM∗U∗ = 0, βU∗V ∗ − δI∗ − ηM∗I∗ = 0, bδI∗ − γV ∗ = 0, (4) ρM∗U∗ ω + U∗ − µM∗ = 0. Solving (4) yields five equilibria: E0, E1, E2, E3, and E4. The trivial equilibrium E0 = (0, 0, 0, 0) and the virus-free equilibrium E1 = (k, 0, 0, 0) (without immune response) always exist. The immune-present, virus-free equilibrium is E2 = ( µω ρ− µ , 0, 0, M∗ 2 ) , M∗ 2 = r [k(ρ− µ)− µω] ηk(ρ− µ) , which exists only if M∗ 2 > 0. The immune-free equilibrium is E3 = ( γ bβ , I∗3 , bβ γ I∗3 , 0 ) , I∗3 = rγδ(bkβ − γ) bβ2(rγ + bkβδ) , and exists if I∗3 > 0. The coexistence equilibrium is E4 = ( µω ρ− µ , I∗4 , bδ γ I∗4 , M ∗ 4 ) , which exists when I∗4 > 0 and M∗ 4 > 0, where I∗4 = γ [k(ρ− µ)(r + δ)− rµω]− bkβδµω (ρ− µ)(rγ + bkβδ) , M∗ 4 = δ [bβµω − γ(ρ− µ)] γη(ρ− µ) . Equilibria are biologically admissible only when all components are positive. While E0 and E1 always exist, the existence of E2, E3, and E4 depends on basic reproduction numbers, as outlined below. The basic reproduction number is a threshold quantity that measures the average number of secondary infections produced by one infected individual in a fully susceptible A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 9 of 24 population [39]. In oncolytic virotherapy, it corresponds to the expected number of newly infected cancer cells generated by a single infected cell. If the basic reproduction number is less than one, the virus-free equilibrium is stable; if greater than one, infection persists and the virus-free equilibrium is unstable. In our model, two virus-free equilibria arise, E1 and E2, with associated reproduction numbers R0 (no immune response) and R1 (with immune activity), respectively. Using the next-generation matrix method [40] we obtain R0 = kbβ γ , R1 = kbβµωδ γ [k(ρ− µ)(r + δ)− rµω] . Define R2 = k(ρ−µ) µω > 0, so that M∗ 2 = r η ( 1− 1 R2 ) . Thus, R2 > 1 if M∗ 2 > 0, i.e., E2 exists if R2 > 1. When R2 > 1, R1 > 0 as well. Moreover, R2 > 1 corresponds to U∗ 1 > U∗ 2 , reflecting the tumor-reducing effect of the immune response. Also, I∗3 = rδγ2(R0 − 1) bβ2(rγ + bkβδ) , so R0 > 1 is necessary and sufficient for the existence of E3. Straightforward calculations show that I∗4 > 0 when R1 < 1, and M∗ 4 > 0 when R0 > R2. Therefore, E4 exists if R1 < 1 and R0 > R2. The stability conditions for E0 and E1, following [41], are summarized below. Theorem 1. For system (3), the following statements hold. (i) E0 is unstable. (ii) E1 is asymptotically stable if and only if s < min  α √ 2Γ(1 + α) r , α √√√√ 4Γ(1 + α)( ρk ω+k − µ ) , α √ 4Γ(1 + α) δ + γ ± √ (δ − γ)2 + 4bδβk  . Proof. (i) At E0, J (E0) =  1 + sα Γ(1 + α) r 0 0 0 0 1− sα Γ(1 + α) δ 0 0 0 sα Γ(1 + α) bδ 1− sα Γ(1 + α) γ 0 0 0 0 1− sα Γ(1 + α) µ  . The eigenvalues are λU = 1 + sα Γ(1 + α) r, λI = 1− sα Γ(1 + α) δ, λV = 1− sα Γ(1 + α) γ, λM = 1− sα Γ(1 + α) µ. Since r > 0, s > 0, and 0 < α ≤ 1, we have λU > 1, hence E0 is unstable. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 10 of 24 (ii) At E1, J (E1) =  1− sα Γ(1 + α) r − sα Γ(1 + α) r − sα Γ(1 + α) βk − sα Γ(1 + α) ηk 0 1− sα Γ(1 + α) δ sα Γ(1 + α) βk 0 0 sα Γ(1 + α) bδ 1− sα Γ(1 + α) γ 0 0 0 0 1− sα Γ(1 + α) ( ρk ω + k − µ )  . Its eigenvalues are λ1 = 1− sα Γ(1 + α) r, λ2 = 1− sα Γ(1 + α) ( ρk ω + k − µ ) , λ3,4 = 1− sα 2Γ(1 + α) ( δ + γ ± √ (δ − γ)2 + 4bδβk ) . Thus, |λ1| < 1 if s < α √ 2Γ(1+α) r , |λ2| < 1 if s < α √ 4(ω+k)Γ(1+α) k(ρ+µ)−µω , and |λ3,4| < 1 if s < α √ 4Γ(1+α) δ+γ± √ (δ−γ)2+4bδβk . Otherwise, E1 is unstable. For the remaining equilibria E∗ i , i = 2, 3, 4, stability is determined using the Schur–Cohn criterion [42]. For n ≥ 3, define the determinants ∆i(µ, x) = ∣∣∣∣∣∣∣∣∣∣∣  1 a1 a2 · · · ai−1 0 1 a1 · · · ai−2 0 0 1 · · · ai−3 ... ... ... . . . ... 0 0 0 · · · 1 ±  an−i+1 an−i+2 · · · an−1 an an−i+2 an−i+3 · · · an 0 . .. . .. . . . . .. . .. an−1 an · · · 0 0 an 0 · · · 0 0  ∣∣∣∣∣∣∣∣∣∣∣ , i = 1, . . . , n. Let the characteristic polynomial of the Jacobian at x0 be Fµ(λ) = a0λ n + a1λ n−1 + · · ·+ an−1λ+ an = 0, (5) with a0 = 1 and ai = ai(µ), i = 1, . . . , n. The equilibrium x0 is asymptotically stable if all eigenvalues of J(µ0, x0) lie inside the unit circle. The general n-dimensional nonlinear discrete-time system is xi+1 = fµ(xi), (6) where xi+1, xi ∈ Rn, i is the iteration index, fµ is the nonlinear vector field, and µ ∈ Rm is the parameter vector. If any eigenvalue lies outside the unit circle, x0 is unstable. We use the following form of the Schur–Cohn test. Theorem 2. [42] The polynomial F (λ) has all roots in the open unit disk if and only if (a) F (1) > 0 and (−1)nF (−1) > 0. (b) ∆± 1 > 0, ∆± 3 > 0, . . . , ∆± n−3 > 0, ∆± n−1 > 0 when n is even, or ∆± 2 > 0, ∆± 4 > 0, . . . , ∆± n−3 > 0, ∆± n−1 > 0 when n is odd. Proposition 1. For any equilibrium E∗ i , i = 2, 3, 4, of system (3), let F (λ) = λ4 + ai1λ 3 + ai2λ 2 + ai3λ+ ai4 be the characteristic polynomial of the Jacobian. The Jacobian at E∗ i is J(E∗ i ) =  1− sα Γ(1+α) ( r k (2U∗ i + I∗i − k) + βV ∗ i + ηM∗ i ) − sα Γ(1+α) ( rU∗ i k ) − sα Γ(1+α) (βU∗ i ) − sα Γ(1 + α) (ηU∗ i ) sα Γ(1+α) (βV ∗ i ) 1− sα Γ(1+α) (δ + ηM∗ i ) sα Γ(1+α) (βU∗ i ) − sα Γ(1+α) (ηI∗i ) 0 sα Γ(1+α) (bδ) 1− sα Γ(1+α) γ 0 sα Γ(1+α) ( ρM∗ i ω+U∗ i )( ω ω + U∗ i ) 0 0 1− sα Γ(1+α) ( µ− ρU∗ i ω+U∗ i )  . A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 11 of 24 Applying Theorem 2 with ∆± 1 = |1| ± |ai4|, ∆± 3 = ∣∣∣∣∣∣ 1 ai1 ai2 0 1 ai1 0 0 1 ∣∣∣∣∣∣± ∣∣∣∣∣∣ ai2 ai3 ai4 ai3 ai4 0 ai4 0 0 ∣∣∣∣∣∣ , the equilibrium E∗ i is asymptotically stable if 1 + ai1 + ai2 + ai3 + ai4 > 0, 1− ai1 + ai2 − ai3 + ai4 > 0, 1± ai4 > 0, ±ai4[ai1(ai1 ± ai3)− (1± ai4)(ai2 ± ai4)] + (1± ai2)(1± ai4)∓ ai3(ai1 ± ai3) > 0. Otherwise, E∗ i is unstable. 5. Analysis of Bifurcation In this section, we analyze bifurcations of model (3). Bifurcation diagrams visually illustrate how the system dynamics evolve in response to parameter variations. In biological systems, bifurcations often signal critical transi- tions, where small changes in parameters, such as reproduction or infection rates, can cause population collapse or the emergence of new stable states. For oncolytic virotherapy with tumor-specific immunity, several key parameters govern the global dynamics. By examining how system behavior changes under parameter variation, we obtain insights into potential treatment outcomes. We therefore focus on codimension-1 and codimension-2 bifurcations in (3). 5.1. Codimension-1 Bifurcations We use algebraic criteria to establish existence conditions for codimension-1 Neimark–Sacker, flip, and fold bifurcations of (3). 5.1.1. Neimark–Sacker Bifurcation The Neimark–Sacker bifurcation (NSB) in discrete-time systems is the analogue of the Hopf bifurcation in continuous- time systems and is crucial for detecting quasiperiodic orbits. In a supercritical NSB, a stable focus loses stability as a parameter varies, giving rise to quasiperiodic behavior; in a subcritical NSB, a stable focus surrounded by an unstable invariant closed curve destabilizes and the curve disappears. Mathematically, an NSB occurs when a complex-conjugate pair of roots of (5) lies on the unit circle, all other roots lie strictly inside, and the critical pair crosses the unit circle with nonzero speed. This is formalized below. Theorem 3. [43] For system (6), a NSB occurs at µ = µ0 if the following hold. (C11) Eigenvalue assignment: ∆− n−1(µ0) = 0, (−1)nFµ0 (−1) > 0, Fµ0 (1) > 0, ∆+ n−1(µ0) > 0, ∆± j (µ0) > 0, for j = n− 3, n− 5, . . . , 1 (or 2) when n is even (or odd, respectively). (C12) Transversality: d∆− n−1(µ) dµ ∣∣∣ µ=µ0 ̸= 0. (C13) Nonresonance: cos ( 2π m ) ̸= 1− Fµ0 (1)∆ − n−3(µ0) 2∆+ n−2(µ0) for m = 3, 4, 5, . . . Applying Theorem 3, model (3) undergoes a NSB with respect to the time step size parameter s if −ai4[ai1(ai1 − ai3)− (1− ai4)(ai2 − ai4)] + (1− ai2)(1− ai4) + ai3(ai1 ± ai3) = 0, 1− ai1 + ai2 − ai3 + ai4 > 0, 1 + ai1 + ai2 + ai3 + ai4 > 0, ai4[ai1(ai1 + ai3)− (1 + ai4)(ai2 + ai4)] + (1 + ai2)(1 + ai4)− ai3(ai1 + ai3) > 0, 1± ai4 > 0, ∂ ∂s {−ai4[ai1(ai1 − ai3)− (1− ai4)(ai2 − ai4)] + (1− ai2)(1− ai4) + ai3(ai1 ± ai3)} ̸= 0. (7) Thus, a NSB occurs at values of s satisfying (7). A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 12 of 24 5.1.2. Flip Bifurcation A flip (period-doubling) bifurcation (FPB) occurs when a real eigenvalue crosses −1, creating a cycle of period two from a period-one orbit, with all other eigenvalues remaining inside the unit circle and the crossing being transversal. Theorem 4. [44] For system (6), a FPB occurs at µ = µ0 if (C21) Eigenvalue assignment: Fµ0 (−1) = 0, Fµ0 (1) > 0, ∆± n−1(µ0) > 0, and ∆± j (µ0) > 0 for j = n − 3, n − 5, . . . , 1 (or 2) when n is even (or odd, respectively). (C22) Transversality: ∑n i=1 a ′ i(−1)n−i∑n j=1(n− j + 1)(−1)n−jaj−1 ̸= 0, where a′i = dai(µ) dµ ∣∣ µ=µ0 . Applying Theorem 4, (3) undergoes a FPB at the time step size s if 1− ai1 + ai2 − ai3 + ai4 = 0, 1 + ai1 + ai2 + ai3 + ai4 > 0, ±ai4[ai1(ai1 ± ai3)− (1± ai4)(ai2 ± ai4)] + (1± ai2)(1± ai4)∓ ai3(ai1 ± ai3) > 0, 1± ai4 > 0, ∂ ∂s (−ai1 + ai2 − ai3 + ai4) ̸= 0. 5.1.3. Fold Bifurcation A fold (saddle–node) bifurcation (FDB) occurs when a real eigenvalue crosses +1, typically corresponding to the collision/creation of fixed points and qualitative changes in dynamics. Theorem 5. [42] For system (6), a FDB occurs at µ = µ0 if (C31) Eigenvalue assignment: Fµ0 (1) = 0, (−1)nFµ0 (−1) > 0, ∆± n−1(µ0) > 0, and ∆± j (µ0) > 0 for j = n− 3, n− 5, . . . , 1 (or 2) when n is even (or odd, respectively). (C32) Transversality: ∑n i=1 a ′ i(−1)n−i∑n j=1(n− j + 1)(−1)n−jaj−1 ̸= 0, where a′i = dai(µ) dµ ∣∣ µ=µ0 . By Theorem 5, (3) undergoes a FDB in s if  1 + ai1 + ai2 + ai3 + ai4 = 0, 1− ai1 + ai2 − ai3 + ai4 > 0, ai4[ai1(ai1 + ai3)− (1 + ai4)(ai2 + ai4)] + (1 + ai2)(1 + ai4)− ai3(ai1 + ai3) > 0, 1± ai4 > 0, (−ai1 + ai2 − ai3 + ai4) ′ ̸= 0. 5.2. Codimension-2 Bifurcations We next derive algebraic conditions for codimension-2 flip–Neimark–Sacker and fold–Neimark–Sacker bifurca- tions of (3). 5.2.1. Flip–Neimark–Sacker Bifurcation The following theorem provides conditions for a flip–Neimark–Sacker (flip–NS) bifurcation. Theorem 6. [45] For system (6), a flip–NS bifurcation occurs at µ = µ0 if (C41) Eigenvalue assignment: F (−1) = 0, ∆− n−2(µ0, x0) = 0, F (1) > 0, ∆+ n−2(µ0) > 0, ∆± l (µ0, x0) > 0, (−1)n−1 n∑ k=1 ( (−1)n−k k∑ ι=1 (−1)k−ιaι−1 ) > 0, A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 13 of 24 with l = n− 4, n− 6, . . . , 1 (or 2) when n is odd (or even). (C42) Transversality: ∂∆− n−2(µ, x) ∂µj ∣∣∣∣∣ µ=µ0 ̸= 0, n∑ i=1 a′ij(−1)n−i ̸= 0, for j = 1, 2, where a′ij = ∂ai/∂µj at µ = µ0. (C43) Nonresonance: cos ( 2π m ) ̸= 1− Fµ0 (1)∆ − n−4(µ0) 4∆+ n−3(µ0) for m = 3, 4, 5, . . ., with ∆± k (µ) = 1 if k ≤ 0. Applying Theorem 6, (3) admits a flip–NS bifurcation in the two parameters (β, s) if 1− ai1 + ai2 − ai3 + ai4 = 0, 1− ai3 + ai4(ai1 − ai4) = 0, 1 + ai1 + ai2 + ai3 + ai4 > 0, 1 + ai3 − ai4(ai1 + ai4) > 0, 1− ai1 + ai2 − ai3 > 0, ∂ ∂s [1− ai3 + ai4(ai1 − ai4)] ̸= 0, ∂ ∂β [1− ai3 + ai4(ai1 − ai4)] ̸= 0, ∂ ∂s (−ai1 + ai2 − ai3 + ai4) ̸= 0, ∂ ∂β (−ai1 + ai2 − ai3 + ai4) ̸= 0. (8) 5.2.2. Fold–Neimark–Sacker Bifurcation A fold–Neimark–Sacker (fold–NS) bifurcation occurs when a fold (saddle–node) and a Neimark–Sacker bifurcation occur simultaneously. Theorem 7. [45] For system (6), a fold–NS bifurcation occurs at µ = µ0 if (C51) Eigenvalue assignment: F (1) = 0, (−1)nF (−1) > 0, ∆+ n−2(µ0, x0) = 0, ∆− n−2(µ0, x0) > 0, ∆± i (µ0, x0) > 0, i = n− 4, n− 6, . . . , 1 (or 2) for n odd (or even), n∑ i=0 (n− i)ai > 0. (C52) Transversality: ∂∆+ n−2(µ, x) ∂µj ∣∣∣∣∣ µ=µ0 ̸= 0 (j = 1, . . . ,m), n∑ i=0 a′ij(−1)n−i ̸= 0 for j = 1, 2, where a′ij = ∂ai/∂µj at µ = µ0. (C53) Nonresonance: cos ( 2π m ) ̸= 1 + (∑n−1 l=0 ∑l i=0 ai ) ∆+ n−4(µ0, x0) 2∆− n−3(µ0) for m = 3, 4, 5, . . ., with ∆± k (µ, x) = 1 if k ≤ 0. Applying (C51)–(C53) in Theorem 7, (3) exhibits a fold–NS bifurcation in (β, s) if 1 + ai1 + ai2 + ai3 + ai4 = 0, 1 + ai3 − ai4(ai1 + ai4) = 0, 1− ai1 + ai2 − ai3 + ai4 > 0, 1− ai3 + ai4(ai1 − ai4) > 0, 4 + 3ai1 + 2ai2 + ai3 > 0, ∂ ∂s [1 + ai1 + ai2 + ai3 + ai4] ̸= 0, ∂ ∂β [1 + ai1 + ai2 + ai3 + ai4] ̸= 0, ∂ ∂s (ai1 + ai2 + ai3 + ai4) ̸= 0, ∂ ∂β (ai1 + ai2 + ai3 + ai4) ̸= 0. (9) Remark 1. In the stability and bifurcation analyses, we denote by ai1, ai2, ai3, ai4 the coefficients of the charac- teristic polynomial at the equilibrium E∗ i (i = 2, 3, 4). We avoid listing explicit formulae due to their length and complexity. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 14 of 24 6. Numerical Simulations In this section, we perform numerical simulations to validate the theoretical results presented above. The com- putations were carried out using Maplesoft (2023 release) and MATLAB R2024a. We present five simulation cases illustrating fold (FDB), flip (FPB), Neimark–Sacker (NSB), flip–NS, and fold–NS bifurcations. The parameters s and β are chosen as primary bifurcation parameters due to their prominent influence on stability and global dy- namics: s is the time–step size, directly affecting temporal resolution and the effective frequency of interactions in the model, while β is the infection rate, a key determinant of the efficacy of oncolytic virotherapy. Varying these parameters reveals stability shifts and transitions among bifurcation regimes. Case 1. Consider r = 3.09, k = 0.9, η = 1, δ = 0.5, b = 0.9, ρ = 0.5, µ = 0.3, ω = 0.6, γ = 0.729, s = 0.78, and α = 0.99. Let β ∈ (0, 1.5). When β = 1.1, an equilibrium E3 = (0.9, 0.11, 0.09, 0) is established, with R0 = R1 = R2 = 1. As β reaches the critical value 1.1, a pitchfork-type FDB emerges. Figure 2 displays the pitchfork bifurcation diagram at β = 1.1: as β increases beyond this threshold, a stable disease-free state develops and cancer cells disappear. The plotted points represent numerically computed steady states of U (uninfected tumor cells), I (infected tumor cells), and V (virus concentration) as β varies, obtained via long-time iteration/root-finding. When β = 0.6, the presence of two vertically aligned red points for U indicates two coexisting equilibria with distinct uninfected tumor cell densities for the same β (multistability): the ultimate outcome depends on initial conditions. The variable U shows a more prominent bifurcation near β ≈ 1.05, whereas I and V change more gradually. Biologically, once β exceeds a threshold, increased infectivity drives a sudden collapse in U and a corresponding rise in I and V , but the latter often grow from near-zero baselines, making their transitions visually less abrupt. Mathematically, this reflects the dominant nonlinear term βUV near the bifurcation point. Overall, when β > 1.1, trajectories converge to a tumor-free state, indicating that sufficiently high infection rates can eradicate cancer cells. Phase portraits in Figure3(a–c) confirm convergence to tumor elimination for β > 1.1. Figure 2: Pitchfork bifurcation diagram of model (3) at E3. (a) β = 1.0 (b) β = 1.1 (c) β = 1.2 Figure 3: Phase portraits corresponding to Figure 2. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 15 of 24 Case 2. Consider r = 3.09, k = 1.04, β = 0.9, η = 1, δ = 0.5, b = 0.9, ρ = 0.5, µ = 0.3, ω = 0.6, γ = 0.729, and α = 0.8. Let s ∈ (0, 0.85). When s = 0.7, an equilibrium E3 = (0.9, 0.11, 0.07, 0) is established with R0 = 1.15, R1 = 0.19, and R2 = 1.15. A flip (period-doubling) bifurcation occurs as s reaches 0.7. Figure 4 shows the FPB diagram for s ∈ (0, 0.85): the cancer burden tends to increase as s grows. Biologically, larger s (less frequent effective interventions) can degrade control, shifting the system from stable regulation to oscillations and then more irregular dynamics. Phase portraits in Figure 5(a–d) illustrate the progression from controlled oscillations to rising cancer cell populations as s increases, emphasizing the need to optimize s to sustain suppression (smaller s favors frequent intervention and improved control). Figure 4: FPB diagram of model (3) at E3. (a) s = 0.6 (b) s = 0.7 (c) s = 0.8 (d) s = 0.85 Figure 5: Phase portraits corresponding to Figure 4. Case 3. Consider r = 2, k = 2.1, β = 2, η = 1, δ = 0.5, b = 1.9, ρ = 0.5, µ = 0.3, ω = 0.2, γ = 0.9, and α = 0.99. Let s ∈ (0, 0.5). When s = 0.07, an equilibrium E4 = (0.3, 0.5, 0.5, 0.13) is established with R0 = 8.8, R1 = 0.28, and R2 = 7. As s increases to the critical value 0.2, E4 loses stability via an NSB. Figure 6 shows A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 16 of 24 the NSB over s ∈ (0, 0.5). Larger s makes regulation more difficult; Figure 7(a–c) shows intensifying instability with increasing s, and chaotic attractors in Figure 7(d) further highlight loss of control. Biologically, increasing s induces quasiperiodic/chaotic fluctuations in tumor levels; thus, more frequent interventions (smaller s) are required to maintain stability. Figure 6: NSB diagram of model (3) at E4. (a) s = 0.1 (b) s = 0.2 (c) s = 0.3 (d) s = 0.5 Figure 7: Phase portraits corresponding to Figure 6. Case 4. Consider r = 2, k = 1.5, η = 1, δ = 0.5, b = 0.9, ρ = 0.5, µ = 0.3, ω = 0.6, γ = 0.729, and α = 0.99. Solving the semi-system (8) yields the critical flip–NS point (s∗, β∗) = (1.6, 1.66). At (s∗, β∗), the equilibrium E3 = (0.488, 0.572, 0.353, 0) satisfies R0 = 3.1, R1 = 0.71, and R2 = 1.7. Figures 8(a,b) display flip–NS diagrams with respect to s and β, respectively. Figure 8(a) shows E3 is stable for s < 1.6 and loses stability as s increases near (1.6, 1.7). In contrast, in Figure 8(b), β produces instability on (0.7, 0.82), then stability for β > 0.82 up to the critical pair (s∗, β∗), where stability is again lost. Thus, increasing infection may eliminate the tumor, but enlarging s can reintroduce instability. Joint tuning of (s, β) is therefore necessary; values in the band (β, s) ∈ (0.82, 1.6) A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 17 of 24 ensure tumor disappearance in this case. Phase portraits in Figure 9 reflect these observations: (a) stable focus at (0.82, 1.6), (b) an unstable invariant circle as s → s∗, and (c) chaotic attractors as (s, β) ∈ (1.6, 1.7). (a) Variation in s (b) Variation in β Figure 8: Flip–NS bifurcation diagrams of model (3) at E3. (a) s = 1.55 (b) s = 1.601 (c) s = 1.68 Figure 9: Phase portraits corresponding to Figure 8. Case 5. Consider r = 2, k = 2, η = 1, δ = 0.7, b = 0.9, ρ = 0.5, µ = 0.3, ω = 0.6, γ = 0.6, and α = 0.99. Solving the first two equations in (9) yields critical values β∗ = 1.9 and s∗ = 0.4, which satisfy (9). For β ∈ (0.8, 2.5) and s ∈ (0, 0.5), the equilibrium E4 = (0.3, 0.5, 0.5, 0.13) arises at (s∗, β∗) with R0 = 5.7, R1 = 0.99, and R2 = 2.2. As s and β pass (s∗, β∗), E4 loses stability through a fold–NS bifurcation. The diagrams in Figure 10 on the β–(U, I, V,M) and s–(U, I, V,M) planes reveal the emergence of an unstable invariant circle at criticality. As s and β increase further, chaotic attractors appear and persist, accompanied by growth in cancer cell counts; see Figure 11(a–c), with chaotic dynamics highlighted in Figure 11(d). These results suggest that maintaining s and β below their fold–NS thresholds is essential to preserve stability and enable the immune response and virotherapy to suppress tumor growth effectively. The intersection of fold and NS mechanisms underscores the need for precise parameter tuning. 7. Chaos Control In our simulations, chaotic behavior of cancer cell populations is undesirable because it complicates effective treatment planning. Chaos control seeks to stabilize the dynamics, keeping trajectories within predictable bounds A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 18 of 24 (a) Variation in β (b) Variation in s Figure 10: Fold–NS bifurcation diagrams of model (3) at E4. (a) β = 1.4 (b) β = 1.95 (c) β = 2.5 (d) β = 3.99 Figure 11: Phase portraits corresponding to Figure 10. to support reliable cancer management. This section applies two approaches, state-feedback and a hybrid control strategy, to regulate chaotic dynamics. 7.1. State Feedback Strategy State-feedback control provides an effective means to regulate chaotic systems [46]. The idea is to transform the chaotic map into a (piecewise) linearized, optimally regulated system via a feedback controller that minimizes an upper bound on the state variables; control is applied under specified conditions to restore stability. With A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 19 of 24 appropriate feedback, the controlled version of model (3) is Un+1(t) = Un + sα Γ(1 + α) [ rUn ( 1− Un + In k ) − βUnVn − ηMnUn ] − j1(Un − U∗ 4 ) , In+1(t) = In + sα Γ(1 + α) [βUnVn − δIn − ηMnIn]− j2(In − I∗4 ) , Vn+1(t) = Vn + sα Γ(1 + α) [bδIn − γVn]− j3(Vn − V ∗ 4 ) , (10) Mn+1(t) = Mn + sα Γ(1 + α) [ ρMnUn ω + Un − µMn ] − j4(Mn −M∗ 4 ) . Here j1(Un−U∗ 4 ), j2(In−I∗4 ), j3(Vn−V ∗ 4 ), and j4(Mn−M∗ 4 ) are the feedback control inputs with gains j1, j2, j3, j4. The Jacobian of (10) yields the characteristic polynomial L(λ) = λ4 + k1λ 3 + k2λ 2 + k3λ+ k4 = 0, (11) where k1 = −(f44 + f33 + f22 + f11) , k2 = f11f22 + f11f33 + f11f44 − f12f21 − f14f41 + f22f33 + f22f44 − f23f32 + f33f44, k3 = −f11f22f33 − f11f22f44 + f11f23f32 − f11f33f44 − f12f21f33 + f12f21f44 − f12f24f41 − f13f21f32 + f14f22f41 + f14f33f41 − f22f33f44 + f23f32f44, k4 = f11f22f33f44 − f11f23f32f44 − f12f21f33f44 + f12f24f33f41 + f13f21f32f44 − f13f24f32f41 − f14f22f33f41 + f14f23f32f41, and f11 = 1− sα Γ(1 + α) ( ηM∗ 4 + βV ∗ 4 − r + r(I∗4 + 2U∗ 4 ) k ) − j1, f12 = − sαrU∗ 4 Γ(1 + α) k , f13 = − sαβU∗ 4 Γ(1 + α) , f14 = − sαηU∗ 4 Γ(1 + α) , f21 = sαβV ∗ 4 Γ(1 + α) , f22 = 1− sα Γ(1 + α) (δ + ηM∗ 4 )− j2, f23 = sαβU∗ 4 Γ(1 + α) , f24 = − sαηI∗4 Γ(1 + α) , f32 = sαbδ Γ(1 + α) , f33 = 1− sαγ Γ(1 + α) − j3, f41 = sαρωM∗ 4 Γ(1 + α) (ω + U∗ 4 ) 2 , f44 = 1− sα Γ(1 + α) ( µ− ρU∗ 4 ω + U∗ 4 ) − j4. By the Schur–Cohn theorem [42], all roots of (11) lie inside the unit disk iff  L(1) = 1 + k1 + k2 + k3 + k4 > 0, L(−1) = 1− k1 + k2 − k3 + k4 > 0, ∆+ 2 = −k34 − (k2 + 1)k24 + (k21 + k1k3 + 1)k4 − k1k3 − k23 + k2 + 1 > 0, ∆− 2 = k34 − (k2 + 1)k24 − (k21 − k1k3 − 2k2 + 1)k4 + k1k3 − k23 − k2 + 1 > 0. (12) (Using Theorem 2), E4 is asymptotically stable if there exist feedback gains (j1, j2, j3, j4) satisfying (12); otherwise, E4 is unstable. Using the parameter values for the NSB at E4 in Case 3, the stability region bounded by the marginal surfaces L(1), L(−1), ∆+ 2 , and ∆− 2 is shown in Figure 12. For illustration, we select (j1, j2, j3, j4) = (0, 0.7, 0.5, 0); the controlled trajectories are stable (see Figure 13). A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 20 of 24 Figure 12: Stability region for state-feedback control. Figure 13: Phase portraits of the controlled system (10). 7.2. Hybrid Control Strategy We also employ a hybrid control to manage chaos arising from Neimark–Sacker bifurcations [47]. As shown in the numerical simulations for Case 3, system (3) undergoes an NSB at E4. Introducing a convex combination of the uncontrolled and “one-step-updated” states yields Un+1(t) = ν ( Un + sα Γ(1 + α) [ rUn ( 1− Un + In k ) − βUnVn − ηMnUn ]) + (1− ν)Un, In+1(t) = ν ( In + sα Γ(1 + α) [βUnVn − δIn − ηMnIn] ) + (1− ν)In, Vn+1(t) = ν ( Vn + sα Γ(1 + α) [bδIn − γVn] ) + (1− ν)Vn, (13) Mn+1(t) = ν ( Mn + sα Γ(1 + α) [ ρMnUn ω + Un − µMn ]) + (1− ν)Mn, with 0 < ν < 1. This strategy merges parameter perturbation and feedback, and suitable ν can shift, delay, or suppress the NSB at E4. The Jacobian of (13) at the positive equilibrium is Jh =  1− h11 sα Γ(1 + α) −h12 sα Γ(1 + α) −h13 sα Γ(1 + α) −h14 sα Γ(1 + α) h21 sα Γ(1 + α) 1− h22 sα Γ(1 + α) h23 sα Γ(1 + α) −h24 sα Γ(1 + α) 0 νbδ sα Γ(1 + α) 1− γν sα Γ(1 + α) 0 h41 sα Γ(1 + α) 0 0 1− h44 sα Γ(1 + α)  , A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 21 of 24 where h11 = ( ηM∗ 4 + βV ∗ 4 − r + r(I∗4 + 2U∗ 4 ) k ) ν, h12 = νrU∗ 4 k , h13 = νβU∗ 4 , h14 = νηU∗ 4 , h21 = νβV ∗ 4 , h22 = (ηM∗ 4 + δ)ν, h23 = νbηU∗ 4 , h24 = νηI∗4 , h41 = νρωM∗ 4 (ω + U∗ 4 ) 2 , h44 = ν ( µ− ρU∗ 4 ω + U∗ 4 ) . Let S = sα Γ(1 + α) . The characteristic polynomial is H(λ) = λ4 + (−4 + ζ1S)λ 3 + ( 6− 3ζ1S + ζ2S 2 ) λ2 + ( −4 + 3ζ1S − 2ζ2S 2 + ζ3S 3 ) λ +1− ζ1S + ζ2S 2 − ζ3S 3 + ζ4S 4 = 0, (14) with coefficients ζ1 = γν + h11 + h22 + h44, ζ2 = (−bδ h23 + γh11 + γh22 + γh44) ν + (h44 + h22)h11 + h14h41 + h22h44 + h12h21, ζ3 = (( (h11 + h22)γ − bδ h23 ) h44 + (h11h22 + h12h21 + h14h41)γ − bδ (h11h23 − h13h21) ) ν + h44(h11h22 + h12h21)− (h12h24 − h14h22)h41, ζ4 = −ν(b ((h13h24 + h14h23)h41 + h44(h11h23 − h13h21)) δ + ((h12h24 − h14h22)h41 − h44(h11h22 + h12h21)) γ) . By the Schur–Cohn theorem [42], stability of (13) is ensured if ζ4S4 > 0, S4ζ4 − 2S3ζ3 + 4S2ζ2 − 8Sζ1 + 16 > 0, − ( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 )3 + ( −S2ζ2 + 3Sζ1 − 7 )( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 )2 + ( (Sζ1 − 4)2 + (Sζ1 − 4)(S3ζ3 − 2S2ζ2 + 3Sζ1 − 4) + 1 )( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 ) −(Sζ1 − 4) (S3ζ3 − 2S2ζ2 + 3Sζ1 − 4)− (S3ζ3 − 2S2ζ2 + 3Sζ1 − 4)2 + S2ζ2 − 3Sζ1 + 7 > 0,( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 )3 + ( −S2ζ2 + 3Sζ1 − 7 )( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 )2 + ( −(Sζ1 − 4)2 + (Sζ1 − 4)(S3ζ3 − 2S2ζ2 + 3Sζ1 − 4) + 2S2ζ2 − 6Sζ1 + 11 )( S4ζ4 − S3ζ3 + S2ζ2 − Sζ1 + 1 ) +(Sζ1 − 4)(S3ζ3 − 2S2ζ2 + 3Sζ1 − 4)− (S3ζ3 − 2S2ζ2 + 3Sζ1 − 4)2 − S2ζ2 + 3Sζ1 − 5 > 0. (15) Therefore, E4 is asymptotically stable if there exists ν ∈ (0, 1) satisfying (15); otherwise, E4 is unstable. Using the parameter values at the NSB of E4 from Case 3, the bifurcation diagram of the controlled system (13) versus ν is shown in Figure 14. Choosing ν = 0.2 stabilizes the trajectories; see Figure 15. Figure 14: Bifurcation diagram of the controlled system (13) with respect to ν. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 22 of 24 Figure 15: Phase portraits of the controlled system (13). 8. Conclusion Host immunity can both hinder and help oncolytic virotherapy. While antiviral responses may neutralize virions and reduce intratumoral spread, tumor-specific immunity can assist in clearing malignant cells while sparing normal tissue. Recent evidence also indicates that virus-mediated tumor lysis can prime strong antitumor immunity, improving outcomes for appropriately engineered vectors. Adenovirus (AdV) remains a leading oncolytic platform because of its safety profile, genetic tractability, and capacity to stimulate immunogenic cell death. We developed a discrete-time, fractional-order model of AdV therapy coupled to a tumor-specific immune response. The fractional term captures memory effects and aligns with discretely sampled biological data. Analyti- cally, we derived conditions for biologically admissible equilibria and established stability criteria using Schur–Cohn and Neimark–Sacker tests. Numerically, we mapped codimension-1 (fold/FDB, flip/FPB, NSB) and codimension-2 (flip–NS, fold–NS) transitions that separate clinically distinct regimes, ranging from stable tumor control to oscil- latory or chaotic progression. The model suggests clear conditions under which therapy can succeed. Coexistence with a controlled tumor burden and, in favorable windows, elimination is achievable when the reproduction numbers satisfy R1 < 1 and R0 > R2, ensuring that the immune-augmented virus can propagate within the tumor while remaining effectively checked by clearance mechanisms. The infection rate β should exceed an efficacy threshold but need not be arbitrarily large: in our simulations, tumor elimination emerged once β ≳ 1.1 at s = 0.78, whereas excessively large β combined with an unfavorable dosing cadence (large s) drove the system through fold–NS interactions into quasiperiodicity and chaos. Thus there is a practical window for β rather than a monotone “more is better.” Treatment cadence matters as much as potency. The discrete time-step s plays the role of an effective dosing or monitoring interval. Smaller s, corresponding to more frequent intervention, stabilizes dynamics and can prevent the NS and flip routes to chaos observed at larger s. Joint tuning of (β, s) is therefore essential: flip–NS and fold–NS curves demarcate narrow safe corridors in the parameter plane, and staying within these corridors maintains stability of the relevant equilibria. These findings translate into practical guidance. Therapy design strategies that raise effective intratumoral infectivity (for example, polymer masking or receptor retargeting) increase β and b/γ, helping achieve R0 > R2 while remaining inside the NS/flip stability window. Scheduling that favors more frequent administrations and monitoring (smaller s) prevents bifurcation cascades and sustains tumor suppression. Because viral therapy alone may be insufficient for complete eradication in typical regimes, combining AdV with immune-modulating agents (to increase ρ or reduce µ), cytotoxic pulses, or targeted radiotherapy can enlarge the stable region while reducing total conventional dosing. There are limitations. Parameters were explored in nondimensional form around fixed baselines; patient-specific calibration, stochastic variability, and spatial heterogeneity warrant future study. The fractional framework is a natural vehicle for memory-aware control: data-assimilated updates of (α, β, s) and real-time bifurcation tracking could guide adaptive dosing that remains inside the stability corridors identified here. In summary, maintaining R1 < 1 and R0 > R2, selecting β above its efficacy threshold but within the fold/NS boundaries, and enforcing a sufficiently small time step s yields stable, tumor-suppressive dynamics in this model. The fractional discrete-time formulation offers a practical, biologically informed tool for designing such regimens and for integrating oncolytic virotherapy synergistically with immune and conventional treatments. Acknowledgements The authors are thankful to Universiti Teknologi Malaysia for providing the facilities in this research. All authors have read and agreed to the published version of the manuscript. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 23 of 24 Funding This research received the Research Management Center (UTM) for financial support through research grants of vote Q.J130000.2554.21H19. Conflict of Interest: The authors declare that they have no conflict of interest. References [1] Nasser Hashemi Goradel, Alexander T Baker, Arash Arashkia, Nasim Ebrahimi, Sajjad Ghorghanlu, and Babak Negahdari. Oncolytic virotherapy: Challenges and solutions. Current Problems in Cancer, 45(1):100639, 2021. [2] Shyambabu Chaurasiya, Nanhai G Chen, and Yuman Fong. Oncolytic viruses and immunity. Current opinion in immunology, 51:83–90, 2018. [3] Jian Zhang, Weijie Lai, Qiang Li, Yang Yu, Jin Jin, Wan Guo, Xiumei Zhou, Xinyuan Liu, and Yigang Wang. A novel oncolytic adenovirus targeting wnt signaling effectively inhibits cancer-stem like cell growth via metastasis, apoptosis and autophagy in hcc models. Biochemical and Biophysical Research Communications, 491(2):469–477, 2017. [4] Jeong Heo, Ja-Der Liang, Chang Won Kim, Hyun Young Woo, I-Lun Shih, Tung-Hung Su, Zhong-Zhe Lin, So Young Yoo, Stanley Chang, Yasuo Urata, et al. Safety and dose escalation of the targeted oncolytic adenovirus obp-301 for refractory advanced liver cancer: Phase i clinical trial. Molecular Therapy, 2023. [5] Jing Cai and Guangmei Yan. The identification and development of a novel oncolytic virus: alphavirus m1. Human gene therapy, 32(3-4):138–149, 2021. [6] Maria Eugenia Davola and Karen Louise Mossman. Oncolytic viruses: how “lytic” must they be for therapeutic efficacy? Oncoimmunology, 8(6):e1581528, 2019. [7] Ulrich M Lauer and Julia Beil. Oncolytic viruses: challenges and considerations in an evolving clinical land- scape. Future Oncology, 18(24):2713–2732, 2022. [8] Malin Peter and Florian Kühnel. Oncolytic adenovirus in cancer immunotherapy. Cancers, 12(11):3354, 2020. [9] Sarah Di Somma, Carmelina Antonella Iannuzzi, Carmela Passaro, Iris Maria Forte, Raffaella Iannone, Vin- cenzo Gigantino, Paola Indovina, Gerardo Botti, Antonio Giordano, Pietro Formisano, et al. The oncolytic virus dl 922-947 triggers immunogenic cell death in mesothelioma and reduces xenograft growth. Frontiers in oncology, 9:564, 2019. [10] Thavasyappan Thambi, JinWoo Hong, A-Rum Yoon, and Chae-Ok Yun. Challenges and progress toward tumor- targeted therapy by systemic delivery of polymer-complexed oncolytic adenoviruses. Cancer Gene Therapy, 29(10):1321–1331, 2022. [11] Joung-Woo Choi, Jung-Sun Lee, Sung Wan Kim, and Chae-Ok Yun. Evolution of oncolytic adenovirus for cancer treatment. Advanced drug delivery reviews, 64(8):720–729, 2012. [12] Joung-Woo Choi, Young Sook Lee, Chae-Ok Yun, and Sung Wan Kim. Polymeric oncolytic adenovirus for cancer gene therapy. Journal of Controlled Release, 219:181–191, 2015. [13] Dominik Wodarz. Viruses as antitumor weapons: defining conditions for tumor remission. Cancer research, 61(8):3501–3507, 2001. [14] Dominik Wodarz. Gene therapy for killing p53-negative cancer cells: use of replicating versus nonreplicating agents. Human gene therapy, 14(2):153–159, 2003. [15] Natalia L Komarova and Dominik Wodarz. Ode models for oncolytic virus dynamics. Journal of theoretical biology, 263(4):530–543, 2010. [16] Yujie Wang, Jianjun Paul Tian, and Junjie Wei. Lytic cycle: a defining process in oncolytic virotherapy. Applied Mathematical Modelling, 37(8):5962–5978, 2013. [17] Akram Ashyani, HM Mohammadinejad, and Omid RabieiMotlagh. Hopf bifurcation analysis in a system for cancer virotherapy with effect of the immune system. Jordan J. Math. Stat., 9:93–115, 2016. [18] Tuan Anh Phan, Jianjun Paul Tian, et al. The role of the innate immune system in oncolytic virotherapy. Computational and mathematical methods in medicine, 2017, 2017. [19] Salma M Al-Tuwairqi, Najwa O Al-Johani, and Eman A Simbawa. Modeling dynamics of cancer virotherapy with immune response. Advances in Difference Equations, 2020:1–26, 2020. [20] Martial Kabong Nono, Elie Bertrand Megam Ngouonkadi, Samuel Bowong, and Hilaire Bertrand Fotsin. Hopf and backward bifurcations induced by immune effectors in a cancer oncolytic virotherapy dynamics. Interna- tional Journal of Dynamics and Control, 9(3):840–861, 2021. [21] Avner Friedman, Jianjun Paul Tian, Giulia Fulci, E Antonio Chiocca, and Jin Wang. Glioma virotherapy: effects of innate immune suppression and increased viral replication capacity. Cancer research, 66(4):2314–2319, 2006. A. T. Alshammari, N. Maan, M. A. M. Abdelaziz / Eur. J. Pure Appl. Math, 18 (4) (2025), 6804 24 of 24 [22] Najwa Al-Johani, Eman Simbawa, and Salma Al-Tuwairqi. Modeling the spatiotemporal dynamics of virother- apy and immune response as a treatment for cancer. Commun. Math. Biol. Neurosci., 2019:Article–ID, 2019. [23] Noufe H Aljahdaly and Nouf A Almushaity. A diffusive cancer model with virotherapy: Studying the immune response and its analytical simulation. AIMS Mathematics, 8(5):10905–10928, 2023. [24] Amal Theyab Alshammari, Normah Maan, and Mahmoud AM Abdelaziz. A new mathematical model approach on the oncolytic virotherapy potency. J. Math. Computer Sci, 36(1):99–120, 2025. [25] JF De Graaf, Lisanne de Vor, RAM Fouchier, and BG Van Den Hoogen. Armed oncolytic viruses: A kick-start for anti-tumor immunity. Cytokine & growth factor reviews, 41:28–39, 2018. [26] Nadishka Jayawardena, Laura N Burga, John T Poirier, and Mihnea Bostina. Virus–receptor interactions: Structural insights for oncolytic virus development. Oncolytic virotherapy, pages 39–56, 2019. [27] Fabrice Le Boeuf, Simon Gebremeskel, Nichole McMullen, Han He, Anna L Greenshields, David W Hoskin, John C Bell, Brent Johnston, Chungen Pan, and Roy Duncan. Reovirus fast protein enhances vesicular stomatitis virus oncolytic virotherapy in primary and metastatic tumor models. Molecular Therapy-Oncolytics, 6:80–89, 2017. [28] Lisette G de Pillis and Ami Radunskaya. A mathematical model of immune response to tumor invasion. In Computational fluid and solid mechanics 2003, pages 1661–1668. Elsevier, 2003. [29] Anyue Yin, Dirk Jan AR Moes, Johan GC van Hasselt, Jesse J Swen, and Henk-Jan Guchelaar. A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors. CPT: pharma- cometrics & systems pharmacology, 8(10):720–737, 2019. [30] Khaphetsi Joseph Mahasa, Amina Eladdadi, Lisette De Pillis, and Rachid Ouifki. Oncolytic potency and reduced virus tumor-specificity in oncolytic virotherapy. a mathematical modelling approach. PLoS One, 12(9):e0184347, 2017. [31] Chipo Mufudza, Walter Sorofa, Edward T Chiyaka, et al. Assessing the effects of estrogen on the dynamics of breast cancer. Computational and mathematical methods in medicine, 2012, 2012. [32] Abeer Hamdan Alblowy, Normah Maan, and Sana Abdulkream Alharbi. Role of glucose risk factors on human breast cancer: A nonlinear dynamical model evaluation. Mathematics, 10(19):3640, 2022. [33] Howard L Kaufman, Carl E Ruby, Tasha Hughes, and Craig L Slingluff. Current status of granulocyte– macrophage colony-stimulating factor in the immunotherapy of melanoma. Journal for immunotherapy of cancer, 2(1):1–13, 2014. [34] Lisette G de Pillis, Ami E Radunskaya, and Charles L Wiseman. A validated mathematical model of cell- mediated immune response to tumor growth. Cancer research, 65(17):7950–7958, 2005. [35] Mahmoud AM Abdelaziz, Ahmad Izani Ismail, Farah A Abdullah, and Mohd Hafiz Mohd. Codimension one and two bifurcations of a discrete-time fractional-order seir measles epidemic model with constant vaccination. Chaos, Solitons & Fractals, 140:110104, 2020. [36] Anuraj Singh, Abdelalim A Elsadany, and Amr Elsonbaty. Complex dynamics of a discrete fractional-order leslie-gower predator-prey model. Mathematical Methods in the Applied Sciences, 42(11):3992–4007, 2019. [37] Senol Kartal. Multiple bifurcations in an early brain tumor model with piecewise constant arguments. Inter- national Journal of Biomathematics, 11(04):1850055, 2018. [38] AMA El-Sayed and SM Salman. On a discretization process of fractional-order riccati differential equation. J. Fract. Calc. Appl, 4(2):251–259, 2013. [39] Purity M Ngina, Rachel Waema Mbogo, Livingstone S Luboobi, et al. Mathematical modelling of in-vivo dynamics of hiv subject to the influence of the cd8+ t-cells. Applied Mathematics, 8(08):1153, 2017. [40] Paul Van den Driessche and James Watmough. Further notes on the basic reproduction number. Mathematical epidemiology, pages 159–178, 2008. [41] Albert CJ Luo. Regularity and complexity in dynamical systems. Springer, 2012. [42] Xiaoliang Li, Chenqi Mou, Wei Niu, and Dongming Wang. Stability analysis for discrete biological models using algebraic methods. Mathematics in Computer Science, 5:247–262, 2011. [43] Guilin Wen. Criterion to identify hopf bifurcations in maps of arbitrary dimension. Physical Review E, 72(2):026201, 2005. [44] Guilin Wen, Shijian Chen, and Qiutan Jin. A new criterion of period-doubling bifurcation in maps and its application to an inertial impact shaker. Journal of sound and vibration, 311(1-2):212–223, 2008. [45] Wei Niu, Jian Shi, and Chenqi Mou. Analysis of codimension 2 bifurcations for high-dimensional discrete systems using symbolic computation methods. Applied Mathematics and Computation, 273:934–947, 2016. [46] Anuraj Singh and Vijay Shankar Sharma. Bifurcations and chaos control in a discrete-time prey–predator model with holling type-ii functional response and prey refuge. Journal of Computational and Applied Mathematics, 418:114666, 2023. [47] Qamar Din, Tzanko Donchev, and Dimitar Kolev. Stability, bifurcation analysis and chaos control in chlorine dioxide–iodine–malonic acid reaction. MATCH Commun. Math. Comput. Chem, 79(3):577–606, 2018.