Adv Syst Sci Appl 2024; 01:69–81 Published online at https://ijassa.ipu.ru. Combined Therapeutic Strategies for Cancer: Integrating Oncolytic Viruses and Inhibitors in a Mathematical Model Majda El Younoussi1*, Khalid Hattaf1,2, Noura Yousfi1 1Laboratory of Analysis, Modeling and Simulation (LAMS), Faculty of Sciences Ben M’Sick, Hassan II University of Casablanca, P.O Box 7955 Sidi Othman, Casablanca, Morocco 2Equipe de Recherche en Modélisation et Enseignement des Mathématiques (ERMEM), Centre Régional des Métiers de l’Education et de la Formation (CRMEF), 20340 Derb Ghalef, Casablanca, Morocco Abstract: The greatest cause of death worldwide continues to be cancer, a complicated set of diseases characterized by uncontrolled cell development. Although early identification and therapeutic approaches have improved, the incidence of the disease is still on the rise, demanding continued study into its underlying causes and cutting-edge treatment paradigms. For creative interventions and focused medicines, the variety of cancer kinds, which are influenced by genetics, way of life, and environmental variables, poses both difficulties and opportunities. In this paper, we present a mathematical model to treat cancer with combined therapies, oncolytic viruses and Mitogen-activated protein kinase inhibitors. We demonstrate that our model is both biologically and mathematically well-posed through the existence, the non-negativity and the boundedness of solutions. Furthermore, we study the equilibrium points as well as the stability of these equilibria. Finally, we use numerical simulations to illustrate the effect of this combined therapy on tumor cells. Keywords: MAPK inhibitors, oncolytic viruses, mathematical modeling, stability, Hopf bifurcation. 1. INTRODUCTION A type of biological therapy called oncolytic viruses is made to target and eliminate cancer cells while sparing normal cells. These viruses are produced naturally or genetically modified in a lab to target only cancer cells. The virus replicates inside the cancer cells after it has infected them, leading to the cells’ bursting and dying. This process also releases new viral particles, which can infect nearby cancer cells and continue the cycle of destruction [2]. Certain oncolytic adenoviruses are notably dependent on the Coxsackie-adenovirus receptor (CAR), and variations in CAR expression levels within target cells could potentially impact the efficacy of viral infection and the resulting therapeutic advantages. CAR has been linked to numerous facets of cancer biology, including cell adhesion, signaling, and migration, in addition to its function in promoting viral entry, making it a viable therapeutic target [9]. Additionally, Mitogen-Activated Protein Kinase, also referred to as MAPK or MEK, is a family of serine/threonine protein kinases that are important for cellular functions like cell growth, differentiation, proliferation, and promoting CAR expression. MAPK has the potential to exploit the complex interplay between CAR and oncolytic viruses to enhance their cancer-killing capabilities and stimulate immune responses against tumor cells [1]. ∗Corresponding author: majda.elyounoussi-etu@etu.univh2c.ma 70 The development of our knowledge of cancer biology and treatment, on the other hand, has been greatly aided by mathematical modeling, particularly the use of ordinary differential equations (ODEs). The intricate dynamics and interactions between various cellular processes, signaling pathways, and therapeutic interventions in cancer can be better understood using ODE models. Researchers can predict how different treatment modalities will affect cancer cells and their microenvironment by creating and analyzing ODE-based models. They can also find potential targets for novel therapies. Further, ODE models can help integrate experimental and clinical data, allowing for the quantitative assessment of cancer progression and treatment outcomes. In 2007, Zurakowski and Wodarz [11] used an ODE model to discuss the interactions between the populations of the average level of CAR expression on the surface of the cells, free virus populations, susceptible, uninfected and infected tumor cells. This model was generalized by Youshan and Qian in their work [10], the researchers constructed a mathematical model to simulate the impacts of both MEK inhibitors and viruses on tumor cells. This model is a free boundary problem, which means that it takes into account the growth and shrinkage of the tumor as the therapies are applied. The researchers used the model to explore how the combined therapies could reduce the tumor size. Recently, Nono et al. [7] recently expanded upon a prior model, applying it to brain cancer and introducing optimal control techniques to optimize the combination of oncolytic virotherapy and MEK inhibitors. Incorporating delays in mathematical models for cancer treatment can play a crucial role in capturing the realistic behavior of biological systems and improving the accuracy of the model predictions. Motivated by all that, we propose in this paper an ODE model to study the dynamics of oncolytic viruses and their interaction with the coxsackie-adenovirus receptor and MAPK inhibitors in the context of cancer therapy, taking into account the duration required by tumor cells that have been infected to generate fresh viruses following the entry of the virus. Our paper is organized into several key sections to present our research cohesively. We introduce our mathematical model in Section 2, and rigorously examine its well-posedness. Moving to Section 3, we delve into equilibrium points and their stability, including an exploration of the Hopf bifurcation, shedding light on the system’s dynamic behavior. Section 4 presents the results of our numerical simulations, providing practical insights into the behavior of the model under various conditions. Finally, in Section 5, we offer a comprehensive conclusion, summarizing our findings, discussing their implications, and highlighting the broader significance of our research. 2. PRESENTATION AND WELL-POSEDNESS OF THE MODEL Within this section, we present the subsequent ordinary differential equation (ODE) model: dS dt = r(1− u)S(t)(1− S(t)+I(t) K )− βW (t)S(t)V (t) 1+αV (t) − dS(t), dI dt = βW (t−τ)S(t−τ)V (t−τ)e−mτ 1+αV (t−τ) − δ(1− u)I(t)− dI(t), dV dt = Nδ(1− u)I(t)− βW (t)S(t)V (t) 1+αV (t) − µV (t), dW dt = ηu(γ −W (t))− hW (t), (2.1) where S(t), I(t), V (t) and W (t) are the concentration of uninfected tumor cells, infected tumor cells, free oncolytic virus particles and the average level of CAR molecules on cell surfaces at the time t, respectively. The factor r is the rate of tumor growth per individual in a population, slowed down by the value (1− u), where u represents the intensity of MAPK inhibitor and varies between 0 and 1. If u = 1, the MAPK inhibitor has the maximum possible effect. If u = 0, every cell in the first phase continue to grow and the production of CAR Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 71 molecule is stopped. Moreover, for biological and mathematical reasons, we will suppose as in [11] that d < r(1− u), where d is the rate of natural death of cells. The parameter K signifies the maximal tumor size. The term βSVW 1+αV models the rate of tumor cells infection by the virus in presence of CAR receptor and the interaction between them on the uninfected cell, where α measures the saturation effect, and β is the rate of infection process. The factor δ represents the virus induced death rate while the rate µ represents the decay of the virus. The parameter denoted as N represents the quantity of newly released viruses following the lysis of an infected tumor cell. Cells generate CAR molecules at a rate denoted as η and experience a loss of these molecules from their surface at a rate of h. The term γ − w characterizes the saturation of CAR expression. Furthermore, τ represents the time required for the transition from tumor cell infection to new virus production, wherem denotes the death rate for infected cells prior to virus production, and e−mτ signifies the probability of survival during the time interval [t− τ, t]. To prove that the model is mathematically well-posed and biologically meaningful, it is important to demonstrate the existence, the boundedness and the non-negativity of solutions as time evolves. Let C be the set of continuous functions from the interval [−τ, 0] to R4, with the supremum norm ||ϕ|| given by sup−τ≤ζ≤0 |ϕ(ζ)|, where ϕ ∈ C. Applying the fundamental theory of functional differential equations [4], we conclude that a single solution exists (S(t), I(t), V (t),W (t)), where the initial condition (S0, I0, V0,W0) are in C and we suppose that: S0(ζ) ≥ 0, I0(ζ) ≥ 0, V0(ζ) ≥ 0,W0(ζ) ≥ 0, ζ ∈ [−τ, 0]. (2.2) Theorem 2.1: Let’s suppose that the initial conditions fulfill (2.2). Then each solution of model (2.1) stays non-negative for all t ≥ 0. Proof Using (2.2), we derive the following: S(t) = S(0)e ∫ t 0 (1−u)(1−S(x)+I(x) K )−βW (x)V (x) 1+αV (x) −d dx, then for all t > 0, we get S(t) ≥ 0. The second equation of model (2.1) gives I(t) = I(0)e−αt + e−mτ−αt ∫ t 0 βS(x−τ)V (x−τ)W (x−τ) 1+αV (x−τ) eαxdx, where α = δ(1− u) + d. Then I(t) ≥ 0 for every t ≥ 0. From the third equation of model (2.1), we obtain V (t) = (V (0)e− ∫ t 0 βS(x) 1+αV (x) dx +Nδ(1− u) ∫ t 0 I(y)eµy− ∫ t y βS(x) 1+αV (x) dxdy)e−µt. Thus, V (t) ≥ 0 for every t ≥ 0. Utilizing the final equation from model (2.1), we acquire: W (t) = W (0)e ∫ t 0 ηuγ W (x) −ηu−h dx ≥ 0, for all t ≥ 0. Hence, every solution of model (2.1) is non-negative for all t ≥ 0. Theorem 2.2: Each solution of the model (2.1), given non-negative initial conditions (2.2), remains bounded for all t ≥ 0. Proof By the first equation of our model, we have Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 72 dS dt ≤ r(1− u)S(t)(1− S(t)−I(t) K ). Using the comparison principal, we get lim sup t→+∞ S(t) ≤ K. Therefore, S(t) is bounded. Let Z(t) = S(t− τ)e−mτ + I(t), hence, we have dZ dt = r(1− u)S(1− S + I K )e−mτ − dSe−mτ − δ(1− u)I − dI ≤ r(1− u)Ke−mτ − (r(1− u) + d)Se−mτ − (δ(1− u) + d)I ≤ r(1− u)Ke−mτ − cZ(t), where c = c ′ (1− u) + d and c′ = min{r, δ}. Then, lim sup t→+∞ Z(t) ≤ rK(1−u)e−mτ c , and we get lim sup t→+∞ I(t) ≤ rK(1−u)e−mτ c . Hence, I(t) is bounded. From the third equation we deduce dV dt ≤ Nδ(1− u)I − µV, by lim sup t→+∞ I(t) ≤ rK(1−u)e−mτ c , we obtain lim sup t→+∞ V (t) ≤ Nδ(1−u)2e−mτ µc . Thus, V (t) is bounded. By the fourth equation of (2.1) and u ∈ [0, 1], we get dW dt ≤ ηγ − (η + h)w, then lim sup t→+∞ W (t) ≤ ηγ η+h . Thus, W (t) is bounded. 3. EQUILIBRIA AND STABILITY ANALYSIS Within this section, we explore the three equilibrium points of model (2.1) along with their stability characteristics. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 73 3.1. Equilibrium Points of the Model When there is no virus, the model (2.1) admits two infection-free equilibrium. The equilibrium point E0 = (0, 0, 0,W0), which reflects the non-existence of cells and virus, and the equilibrium point E1 = (S1, 0, 0,W1) where S1 = K(1− d r(1−u) ) and W0 = W1 = ηuγ ηu+h . By r(1− u) > d, the equilibrium E1 exists. In the existence of the virus, there exists another equilibrium point called the endemic equilibrium E∗ = (S∗, I∗, V ∗,W ∗). We suppose that I∗ > 0, V ∗ > 0 and we put R0 = βNδW1S1(1− u)e−mτ (δ(1− u) + d)(µ+ βW1S1) . R0 is the reproduction number and represents the potential for the oncolytic virus to spread within the tumor cell population. By simple calculus, we show that the equilibrium E∗ exists if R0 > 1, this means that the virus can sustain its presence and spread within the tumor cells, leading to a persistent infection. (S∗, I∗, V ∗, W ∗) are the solution of this system: r(1− u)S∗(1− S∗ + I∗ K )− βW ∗S∗V ∗ 1 + αV ∗ − dS∗ = 0, (3.3) βW ∗S∗V ∗e−mτ 1 + αV ∗ − δ(1− u)I∗ − dI∗ = 0, (3.4) Nδ(1− u)I∗ − βW ∗S∗V ∗ 1 + αV ∗ − µV ∗ = 0, (3.5) ηu(γ −W ∗)− hW ∗ = 0. (3.6) By equation (3.6), we get W ∗ = ηuγ ηu+ h . (3.7) By adding the equation (3.4) to the equation (3.5), we get V ∗ = δ(1− u)(N − emτ )− demτ µ I∗. (3.8) By adding the equation (3.3) to the equation (3.4), we obtain I∗ = S∗(r(1− u)(K − S∗)− dK) r(1− u)S∗ +Kemτ (δ(1− u) + d) . (3.9) Obviously, by determining S∗ we will determine V ∗ and I∗, for that we will use (3.7), (3.8) and (3.9) to find S∗. For simplification, we put: P = δ(1−u)(N−emτ )−demτ µ and O = Nδ(1− u). By equation (3.5), we get (1− αV ∗)OI∗ − βW ∗S∗V ∗ − µ(1 + αV ∗)V ∗ = 0, then, O − βW ∗S∗P − Pµ+ (αPO − µαP 2)I∗ = 0. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 74 Using (3.9), we obtain aS2 + bS + c = 0, where, a = −βW ∗Pr(1− u)− αOPr(1− u) + αP 2µr(1− u), b = Or(1− u)− Pµr(1− u)− βW ∗P ∗K(1− u)emτ (δ + d) +αOPK(r(1− u)− d)− αP 2µr(1− u)K + αP 2µdK, c = (Oδ(1− u)− Pµδ(1− u) +Od− Pµd)Kemτ . Then S∗ = −b− √ ∆ 2a , where ∆ = b2 − 4ac. Thus S∗, I∗, V ∗ and W ∗ are defined. 3.2. Stability Analysis The characteristic equation at any equilibrium E = (S, I, V,W ) is given by∣∣∣∣∣∣∣∣∣∣∣∣∣∣ −r(1−u)S K − λ − r(1−u)S K − βWS (1+αV )2 − βSV 1+αV βWV e−(m+λ)τ 1+αV −δ(1− u)− d− λ βWSe−(m+λ)τ (1+αV )2 βV Se−(m+λ)τ 1+αV − βWV 1+αV Nδ(1− u) − βWS (1+αV )2 − µ− λ − βSV 1+αV 0 0 0 −ηu− h− λ ∣∣∣∣∣∣∣∣∣∣∣∣∣∣ = 0. (3.10) Theorem 3.1: The equilibrium state E0 = (0, 0, 0,W0) exhibits instability. Proof At E0 = (0, 0, 0,W0) we get the following equation: (r(1− u)− d− λ)(δ(1− u) + d+ λ)(µ+ λ)(ηu+ h+ λ) = 0, (3.11) then the roots are λ1 = r(1− u)− d, λ2 = −δ(1− u)− d, λ3 = −µ and λ4 = −ηu− h. Since we have r(1− u) > d, we get λ1 > 0. Thus, E0 is unstable. Theorem 3.2: The equilibrium point E1 = (S1, 0, 0,W1) is locally asymptotically stable for every τ ≥ 0 if R0 < 1 and unstable if R0 > 1. Proof At E1, (3.10) becomes( 2r(1− u) S1 K + d+ λ )( ηu+ h+ λ )( λ2 + (δ(1− u) + d+ βW1S1 + µ)λ +(δ(1− u) + d)(βW1S1 + µ)(1−R0e −λτ ) ) = 0. (3.12) Clearly, λ1 = −2r(1− u)S1 K − d and λ2 = −ηu− h represent two of the roots of the aforementioned equation, with the remaining roots arising from solutions to the subsequent Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 75 equation: λ2 + (δ(1− u) + d+ βW1S1 + µ)λ+ (δ(1− u) + d)(βW1S1 + µ)(1−R0e −λτ ) = 0. (3.13) If R0 > 1, let f(λ) = λ2 + (δ(1− u) + d+ βW1S1 + µ)λ+ (δ(1− u) + d)(βW1S1 + µ)(1−R0e −λτ ). We have f(0) = (δ(1− u) + d)(βW1S1 + µ)(1−R0) < 0 and lim λ→+∞ f(λ) = +∞. In this case, the equation f(λ) = 0 possesses at least one positive root. Hence, if R0 > 1 E1 is unstable. If R0 < 1, we discuss two cases: When τ = 0, we get δ(1− u) + d+ βW1S1 + µ > 0 and (δ(1− u) + d)(βW1S1 + µ)(1− R0) > 0. Then all the roots of (3.12) have negative real parts for τ = 0 and R0 < 1. When τ > 0, let iω be a purely imaginary root of (3.13) where ω > 0. Then,{ −ω2 + (δ(1− u) + d)(βW1S1 + µ) = (δ(1− u) + d)(βW1S1 + µ)R0cos(wτ), (δ(1− u) + d+ βW1S1 + µ)ω = −(δ(1− u) + d)(βW1S1 + µ)R0sin(ωτ), thus we get ω4 + ( (δ(1− u) + d)2 + (βW1S1 + µ)2 ) ω2 + (δ(1− u) + d)2(βW1S1 + µ)2(1−R2 0) = 0. Let x = ω2, then we obtain x2 + ( (δ(1− u) + d)2 + (βW1S1 + µ)2 ) x+ (δ(1− u) + d)2(βW1S1 + µ)2(1−R2 0) = 0. Hence if R0 < 1, there is no positive solution. Thus, E1 is locally asymptotically stable for R0 < 1. Theorem 3.3: If R0 < 1, the equilibrium E1 is globally asymptotically stable. Proof Take into consideration the presented Lyapunov function: L(t) = δ(1− u)emτI(t) + δ(1− u) + d N emτV (t) + δ(1− u) ∫ t t−τ βW (ξ)S(ξ)V (ξ) 1 + αV (ξ) dξ, then, we obtain dL dt = δ(1− u) βWSV 1 + αV − (δ(1− u) + d)emτ N ( βWSV 1 + αV + µV ) . We have V 1+αV ≤ V , lim sup t→∞ S(t) ≤ S1 and lim sup t→∞ W (t) ≤ W1. Then, we get dL dt ≤ (δ(1− u) + d)(βW1S1 + µ)(R0 − 1)V N . Hence, if R0 < 1 we obtain dL dt ≤ 0. Clearly, dL dt = 0 if and only if S = S0, I = 0, V = 0 and W = W0. Then the largest invariant set contained in {(S, I, V,W )|dL dt = 0} is the singleton {E1}. By LaSalle’s invariance principale [5], we deduce that E1 is globally asymptotically stable when R0 < 1. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 76 The characteristic equation at the equilibrium E∗ can be expressed in the following manner: (ηu+ h+ λ) ( λ3 + p1λ 2 + p2λ+ p3 + (q1λ+ q2)e −λτ ) = 0, (3.14) where, p1 = (1− u) ( rS∗ K + δ ) + βW ∗S∗ (1 + αV ∗)2 + d+ µ, p2 = (1− u) ( βW ∗S∗ (1 + αV ∗)2 + µ )( rS∗ K + δ + d 1− u ) − β2S∗V ∗W ∗2 (1 + αV ∗)3 + r(1− u)S∗ K ( δ(1− u) + d ) , p3 = (δ(1− u) + d)S∗ ( r(1− u) K ( βS∗W ∗ (1 + αV ∗)2 + µ ) − β2V ∗W ∗2 (1 + αV ∗)3 ) , q1 = (1− u)βS∗W ∗ 1 + αV ∗ ( rV ∗ K − δN 1 + αV ∗ ) e−mτ , q2 = (1− u)βS∗W ∗ 1 + αV ∗ ( µrV ∗ K + Nδ 1 + αV ∗ ( βV ∗W ∗ 1 + αV ∗ − r(1− u)S∗ K )) e−mτ . Since the root λ1 = −ηu− h is negative, it remains to determine the roots of the following equation: λ3 + p1λ 2 + p2λ+ p3 + (q1λ+ q2)e −λτ = 0, (3.15) When τ = 0, the equation (3.15) becomes λ3 + p1λ 2 + (p2 + q1)λ+ p3 + q2 = 0. (3.16) We have p1 > 0, and by simple calculation we can find that p3 + q2 > 0. By Routh-Hurwitz criterion, we conclude the result bellow: Lemma 3.4: Suppose that R0 > 1 and p1(p2 + q1)− (p3 + q2) > 0. Then in the absence of delay (τ = 0), all the roots of (3.15) have negative real parts. Hence, the equilibriumE∗ = (S∗, I∗, V ∗,W ∗) is locally asymptotically stable. When τ > 0, let iω (ω > 0) be a purely imaginary root of the equation (3.15). Thus,{ p1ω 2 − p3 = q1ω sin(ωτ) + q2 cos(ωτ), −ω3 + p2ω = −q1ω cos(ωτ) + q2 sin(ωτ), (3.17) then, ω6 + ( p21 − 2p2 ) ω4 + ( p22 − q21 − 2p1p3 ) ω2 + p23 − q22 = 0. (3.18) For x = ω2, the equation (3.18) is reduced to x3 + (p21 − 2p2)x 2 + (p22 − q21 − 2p1p3)x+ p23 − q22 = 0. (3.19) We consider the following function: f(x) = x3 + c2x 2 + c1x+ c0, (3.20) where c0 = p23 − q22 , c1 = p22 − q21 − 2p1p3 and c2 = p21 − 2p2. Obviously f ′(x) = 3x2 + 2c2x+ c1, and ∆′ = 4(c22 − 3c1) its discriminant. Hence, we get the following result: Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 77 Lemma 3.5: (i) If c0 < 0, then the equation f(x) = 0 has at least one positive root. (ii) If c0 ≥ 0 and ∆′ ≤ 0, then the equation f(x) = 0 has no positive roots. (iii) If c0 ≥ 0 and ∆′ > 0, then the equation f(x) = 0 has a positive root if x1 > 0 and f (x1) ≤ 0, where x1 = √ c22−3c1−c2 3 is a root of f ′(x) = 0. By Lemma 3.4 and the previous Lemma we deduce the following theorem: Theorem 3.6: Assume that R0 > 1, c0 ≥ 0 and p1 (p2 + q1)− (p3 + q2) > 0. If any of the subsequent conditions are met, • ∆′ ≤ 0, • ∆′ > 0 and x1 ≤ 0, • ∆′ > 0 and f (x1) > 0, then the equilibrium pointE∗ is locally asymptotically stable for any non-negative time delay. On the other hand, we study the Hopf bifurcation of model (2.1) at the equilibrium point E∗. The stability of the equilibrium point E∗ changes when the equation (3.15) has purely imaginary roots. So, we assume that ω1, ω2 and ω3 are these positive roots. Moreover, we consider τ as a parameter of bifurcation. Substituting ω = ωϵ and τ = τ ϵ in (3.17) where ϵ = 1, 2, 3, we get q2(p1ω 2 ϵ − p3)− q1ωϵ(−ω3 ϵ + p2ωϵ) = q22cos(ωϵτ ϵ) + q21ωϵ 2cos(ωϵτ ϵ), thus, we obtain τ ϵn = 1 ωϵ ( arccos ( q2(p1ω 2 ϵ − p3) + q1ω 2 ϵ (ω 2 ϵ − p2) q22 + q21ω 2 ϵ ) + 2πn ) , (3.21) where n ∈ N. Obviously, ±iωϵ are purely imaginary roots of (3.15) with τ = τ ϵn. Let τ0 = min ϵ∈{1,2,3} {τ ϵ0} and λ(τ) = ψ(τ) + iω(τ) be the root of the equation (3.15) where ψ (τ ϵn) = 0 and ω (τ ϵn) = ωϵ. Differentiating equation (3.15) with respect to τ , we get( dλ dτ )−1 = 3λ2 + 2p1λ+ p2 + q1e −λτ λ (q1λ+ q2) e−λτ − τ λ . Therefore, it is straightforward to deduce Re ( dλ dτ )−1 ∣∣∣∣∣ τ=τϵn = 3ω4 ϵ + 2 (p21 − 2p2)ω 2 ϵ + p22 − q21 − 2p1p3 q21ω 2 ϵ + q22 = f ′ (ω2 ϵ ) q21ω 2 ϵ + q22 . By f ′ (ω2 1) > 0, f ′ (ω2 2) < 0 and f ′ (ω2 3) > 0, the transversality condition is verified, and we suppose these conditions: (a) c0 < 0, (b) c0 ≥ 0,∆′ > 0, x1 > 0 and f (x1) ≤ 0, and we get the following result: Theorem 3.7: Suppose R0 > 1 and p1 (p2 + q1)− (p3 + q2) > 0. If one of conditions (a)-(b) is satisfied, the equilibrium point E∗ is locally asymptotically stable for all time delays τ ∈ [0, τ0). Furthermore, E∗ becomes unstable when τ > τ0. In addition, when τ = τ ϵn model (2.1) undergoes a Hopf bifurcation at E∗ where ϵ = 1, 2, 3 and n ∈ N. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 78 4. NUMERICAL SIMULATIONS This section presents numerical simulations of our system to demonstrate its dynamics and behavior. The system is evaluated with various initial conditions satisfying S0, I0, V0,W0 > 0, and the time interval is set from t = 0 to t = 2000. The following set of parameters, selected based on previous works [3, 6, 8, 11], are used: r = 0.5, u = 0.5, d = 0.1, K = 2× 109, α = 1.95× 10−10, β = 1.2× 10−10, δ = 0.5,N = 1000, µ = 20, η = 0.17, h = 0.07, γ = 7, τ = 2, m = 1. Using these values, we obtain R0 = 3.8481 > 1, p1(p2 + q1)− (p3 + q2) = 11.2223 > 0, c0 = 0.1873 > 0, ∆′ = 6.7108× 105 > 0 and f(x1) = 0.1224 > 0. According to Theorem 3.6, the equilibrium point E∗ is locally asymptotically stable for any τ ≥ 0. This is illustrated and validated in Figure 4.1. Fig. 4.1. Dynamical behavior of system (2.1) around the equilibrium point E∗ when τ = 2, µ = 20 and γ = 7. We then consider the same parameter values, except that γ is changed to 2 and µ to 0.65. In this case, we obtain c0 = −9.7104× 10−4 < 0, which means that the condition of Theorem 3.6 is not satisfied. Therefore, the equilibrium E∗ is unstable, as shown in Figure 4.2. When m = 0 and N = 100, we obtain R0 = 19.0795 > 1, p1(p2 + q1)− (p3 + q2) = 0.0054 > 0, and c0 = −9.0776× 10−4 < 0. Therefore, the conditions for Theorem 3.7 are satisfied. Using equation (3.19) and (3.21), we obtain ω = 0.4186 and τ0 = 0.8567. In Figure 4.3, when τ = 0.7 < τ0, the equilibrium point E∗ is locally asymptotically stable. However, when τ = τ0 = 0.8567, model (2.1) undergoes a Hopf bifurcation at the equilibrium E∗. Furthermore, when the parameter τ exceeds the critical threshold τ0 associated with the Hopf bifurcation, the equilibrium at the point E∗ transitions from stable steady state to stable limit cycle oscillation, as illustrated in Figure 4.4. Moreover, all these results confirm our theoretical results stated in Theorem 3.7 Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 79 Fig. 4.2. Dynamical behavior of system (2.1) around the equilibrium point E∗ when τ = 2, µ = 0.65 and γ = 2. Fig. 4.3. The stability of the equilibrium point E∗ when τ = 0.7 < τ0. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 80 Fig. 4.4. The instability of the equilibrium point E∗ when τ = 5 > τ0. 5. CONCLUSION This paper has presented a novel contribution by introducing an ordinary differential equation model with delay designed to comprehensively capture and analyze the intricate dynamics of tumor cells following the administration of oncolytic viruses. Our model has carefully considered the role of the Coxsackie adenovirus receptor (CAR), as well as the influence of mitogen-activated protein kinase inhibitors, factors that have played pivotal roles in the interaction between the virus and the tumor microenvironment. By incorporating these crucial elements, our research seeks to provide a deeper understanding of the underlying mechanisms governing the response of tumor cells to oncolytic virus treatment. Through this innovative model, our aim has been to provide valuable information that can inform the development of more effective therapeutic strategies to combat cancer. For that, we have proven that our proposed model is mathematically and biologically meaningful through its existence, non-negativity, and boundedness of solution. We have explored the equilibrium points of the model and have identified three possible steady states. The first, E0, represents a state where there are no cells (uninfected or infected) and no virus particles. This equilibrium may not have biological relevance, as it reflects the non-existence of tumor cells and virus particles. The second, E1, describes where there are uninfected tumor cells at a constant concentration S1, but no infected tumor cells or virus particles, and the average level of CAR molecules on the surface of the cells is at W1. This equilibrium could represent a state where the virus fails to infect and spread in the tumor. The third, E∗, was an endemic equilibrium in which uninfected tumor cells, infected tumor cells, virus particles, and CAR molecules on the surface of the cells reached constant concentrations over time: S∗, I∗, V ∗, and W ∗, respectively. This equilibrium could represent a coexistence between the tumor cells, infected cells, and the virus. Furthermore, an examination of the local stability of these three equilibrium points has been carried out employing the characteristic equation. Moreover, the global stability of the equilibrium point E1 has been established by employing Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) 81 an appropriate Lyapunov function. In addition to the stability analysis of the equilibria, we have also investigated the effects of time delay on the dynamics of the system. Our analysis has revealed that the introduction of a time delay parameter can lead to a Hopf bifurcation at the equilibrium point E∗, inducing a shift in equilibrium from a stable steady state to a stable limit cycle oscillation. Our numerical experiments additionally have validated the theoretical findings, demonstrating the impact of time delay on the system’s behavior and the emergence of the Hopf bifurcation. Overall, our study has highlighted the importance of considering the time delay in modeling the dynamics of oncolytic virus therapy and has provided valuable information on the long-term behavior of the system. ACKNOWLEDGEMENTS We would like to extend our appreciation to everyone who contributed their knowledge and ideas to this research, regardless of how small a way. The editors and anonymous referees are thanked by the authors for their insightful criticism and recommendations, which significantly increased the caliber of this work. REFERENCES 1. Cuschieri, J. & Maier, R. V. (2005) Mitogen-activated protein kinase (MAPK), Critical care medicine, 33, S417–S419. 2. Everts, B. & van der Poel, H. G. (2005) Replication-selective oncolytic viruses in the treatment of cancer, Cancer gene therapy, 12, 141–161. 3. Friedman, A., Tian, J. P., Fulci, G., Chiocca, E. A. & Wang, J. (2006) Glioma virotherapy: The effects of innate immune suppression and increased viral replication capacity, Cancer Research, 66, 2314–2319. 4. Hale, J. & Verduyn Lunel, S.M. (1993) Introduction to Functional Differential Equations. New York, NY: Springer. 5. Huo, H., Zhao, H. & Zhu, L. (2015) The effect of vaccines on backward bifurcation in a fractional order HIV model, Nonlinear Anal., Real World Appl., 26, 289–305. 6. Linsenmann, T., Jawork, A., Westermaier, T., Homola, G., Monoranu, C. M. & Vince, G. H. (2019) Tumor growth under rhGM-CSF application in an orthotopic rodent glioma model, Oncolytic Letter, 17(6), 4843–4850. 7. Nono, M. K. & Ngouonkadi, E. M. (2020) Synergistic effects of oncolytic adenovirus and MEK inhibitors on glioma treatment dynamics: Analysis and optimal control, Applied Mathematical Sciences, 14(16), 781–800. 8. Okamoto, K. W., Priyanga, A. I. & Petty, T. D. (2014) Modeling oncolytic virotherapy: Is complete tumor-tropism too much of a good thing?, Journal of Theoretical Biology, 358, 166–178. 9. Wunder, T., Schmid, K., Wicklein, D., Groitl, P., Dobner, T., Lange, T., Anders, M. & Schumacher, U. (2013) Expression of the coxsackie adenovirus receptor in neuroendocrine lung cancers and its implications for oncolytic adenoviral infection, Cancer gene therapy, 20, 25–32. 10. YouShan, T. & Guo, Q. (2008) A mathematical model of combined therapies against cancer using viruses and inhibitors, Science in China Series A: Mathematics, 51(12), 2315–2329. 11. Zurakowski, R. & Wodarz, D. (2007) Model-driven approaches for in vitro combination therapy using ONYX-015 replicating oncolytic adenovirus, J Theoret Biol, 245, 1–8. Copyright © 2024 ASSA. Adv Syst Sci Appl (2024) Introduction Presentation and well-posedness of the model Equilibria and stability analysis Equilibrium Points of the Model Stability Analysis Numerical simulations Conclusion