Analysis of a deterministic-stochastic oncolytic M1 model involving immune response via crossover behaviour: ergodic stationary distribution and extinction

1.   Introduction

Cancer is a category of infections characterized by malevolent tumor cell proliferation and transmission. It can result in mortality if not monitored. In 2018, the International Agency for Research on Cancer (IARC) predicted 18.1 million innovative scenarios of tumors and 9.6 deaths worldwide ^[1]. Numerous treatments, such as brachytherapy and chemotherapy, have been employed to handle leukaemia ^[2]. Even so, these treatments are ineffective and can infiltrate both tumors and healthy tissue. As a result, they may induce a diverse range of side effects, including exhaustion and baldness ^[3]. Oncolytic virotherapy is a new cancer treatment that has recently dominated diagnostic and mathematical research ^[4]. It employs preferential oncolytic pathogens that are designed to attack tumors while avoiding normal tissue. Oncolytic virions invade and replicate within tumors, leading to cell death. After a malignant microbe dies, a significant percentage of unique pathogens are distributed that can scatter and decapitate other tumor cells ^[5]. As a result, the capacity of oncolytic immunotherapy to entirely eliminate the tumor is dependent on its effectiveness.

Many oncogenes have been established and conducted clinical trials, with promising effects ^[6,7]. Nonetheless, there are considerable obstacles to this therapy method, which may lessen its effectiveness and necessitate additional research ^[8,9]. One of these difficulties is phagocytes against cancerous cells, which can constrain oncolytic pathogen growth and therefore diminish the quantity of pathogens. One solution is to create preferential pathogens that can swiftly expand and induce apoptosis ^[9]. Nevertheless, the relationship between oncolytic immunotherapy and malignant cell immunity is immensely complicated, and studies in this field are ongoing. Immune systems are often energised to endorse oncolytic drawbacks associated with and concise tumors ^[10].

In an attempt to formulate improved and more efficacious cancer therapies, mathematical models have been employed to assist in clarifying the intricate complexities of oncolytic immunotherapy. Several of these concepts are similar to those used to study HBV and HIV infectious diseases ^{[11,12,13,14]}. Wang et al. ^[6], for instance, lengthened the basic oncogenic virotherapy framework ^[15] and investigated the influence of pathogen interruption volume on highly contagious therapy, which reflects the number of newly acquired adenovirus infections constituted within a tumor cell. They discovered that the pathogen burst shape is an important factor to consider when creating innovative pathogens for viral therapeutic interventions. Okamoto et al. ^[8] developed a concept to assess the potential of judicious and non-selective cytotoxic agents to eradicate tumors and strengthen therapeutic approaches. Malinzi et al. ^[3] proposed a framework to investigate the impact of incorporating oncolytic virotherapy with pharmacotherapy. Their approach integrates tumors, immunological, oncolytic virotherapy, and gemcitabine. They demonstrated that virotherapy can improve treatments if the proper beneficial intake is utilized. Kim et al. ^[16] investigated the environments for the presence of a Hopf bifurcation in the connections between cancerous cells, the cytotoxic T lymphocyte (CTL) immune system, and virotherapy using a mathematical framework.

When searching through the publications, it is clear that numerous authors have offered so many proposals for finding kernels that could indeed be employed to achieve fractional derivative formulations ^{[17,18,19,20,21,22]}. The goal here is to employ similar mathematical techniques to encapsulate the mechanisms evidenced by serious challenges well with an index law kernel. It has surprisingly strengthened its extensive implementation in autocatalytic interactions, irregular exothermic reactions, deformable dissipative, Maxwell fluid, epidemiology, convection-diffusion methodology, thermoelectric disciplines, heat transfer, photonics, aerodynamics, and noncommunicable diseases; various fractional calculus compositions are incorporated in FDEs to appropriately construe and analyze memory ^{[23,24,25,26]}. Individuals including Coimbra, Davison, and Essex, Riesz, Riemann-Liouville, Hadamard, Weyl, Jumarie, Caputo and Fabrizio, Atangana and Baleanu, Grünwald and Letnikov, Liuville and Caputo have elaborated on various types of interpretations and conceptualizations of fractional formulations ^{[27,28,29,30]}. However, the Liouville-Caputo ^[31], Caputo-Fabrizio ^[32] and Atangana-Baleanu ^[33] operators are the ideal fractional classifiers ^[34,35,36]. While each derivative is assigned in a different way, the concept of the piecewise derivative proposed by Atangana and Araz ^[37], which can simulate the crossover mechanisms of these derivative instruments in a separate procedure, has recently gained popularity, see ^{[38,39,40,41,42]}.

Lin et al. ^[43] discovered an instinctually precipitating alphavirus M1 as a precise oncolytic pathogen that targets tumor cells lacking zinc-finger antibiotic nutrients. The pathogen was highly effective at assassinating cancerous cells while causing no detriment to immune tissues. Wang et al. ^[1] developed an ordinary differential modelling approach with competitive rivalry between healthy and cancerous molecules on a constricted feed ingredient to comprehend the function of this pathogen in oncogenic immunotherapy. They examined the impact of the oncolytic M1 pathogen on cancerous and regular tissue regeneration as well as the required virotherapy potency mandated to efficaciously negate the tumor. Elaiw et al. ^[44] described the dynamic behaviour of reaction-diffusion oncolytic M1 virotherapy involving immune reaction. Rashid et al. ^[45] expounded the numerical computations of the fractional oncolytic efficacy system with the M1 virus via a generalized fractional derivative with optimal criterion. Arik and Araz ^[41] presented the crossover impacts of a system of tumor growth and its response to radiotherapy. For further investigation of the oncolytic viruses and fractional calculus techniques on epidemiology we refer the readers, see ^{[46,47,48,49,50,51]} and the references cited therein.

Adopting the aforesaid proclivity, we intend to broaden the oncolytic M1 infection model by incorporating the influence of CTL immune reactions on oncolytic immunotherapy by convolving the white noise. The entire procedure from the start to the end of the virotherpy is taken into consideration within a piecewise scenario where the distinctive derivative formulations can be employed. Besides that, the specifications in the aforesaid model are clearly not constant and always change drastically around certain average values due to regular volatility in the environment. Additionally, the qualitative aspects of the model are presented in terms of the global positive recurrence, extinction and persistence, ergodic and stationary distribution by using appropriate Lyapunov candidates. It was also indicated that the framework be inferred using piecewise derivatives, which we will accomplish in this document with an oncolytic M1 infection model.

The key findings are presented in the following section. Section 2 presents the fundamental characteristics of the fractional calculus and solutions for the proposed model are investigated. Section 3 demonstrates a sufficient threshold for ensuring the existence and ergodicity of a stationary distribution. Meanwhile, the model's extinction conditions for all biological assumptions and viral infections are being debated. Then, we accomplish numerical computations for each scenario where procedures are taken into account in Section 4. In a nutshell, we summarize the main findings in the conclusion section.

2.   Model description and preliminaries

Elaiw et al. ^[44] discovered a naturally produced alphavirus M1 as a preferential oncolytic pathogen that targets cancerous cells lacking zinc-finger antiviral protein (ZAP). The pathogen was highly effective at assassinating cancerous cells while causing no detriment to normal tissue. Their model has the appearance of

where $\tilde{\mathbf{H}}(\xi), \; \tilde{\mathbf{N}}(\xi), \; \tilde{\mathbf{Y}}(\xi), \; \tilde{\mathbf{V}}(\xi)$ and $\tilde{\mathbf{Z}}(\xi)$ constitute nutrient, normal cell, tumor cell, free M1 virus and CTL cells, respectively. The component $\chi$ symbolizes the nutritional recruitment rate as well as the M1 virus's minimum efficacious potency. The nutrient is consumed by normal and tumor cells at rates of $\rho_{1}\tilde{\mathbf{N}}\tilde{\mathbf{H}}$ and $\rho_{2}\tilde{\mathbf{H}}\tilde{\mathbf{Y}}$, respectively. The growth rate of normal cells resulting of nutritional consumption is provided by $\lambda_{1}\rho_{1}\tilde{\mathbf{H}}\tilde{\mathbf{N}}$, while a rate of expansion of tumor cells is provided by $\lambda_{2}\rho_{2}\tilde{\mathbf{H}}\tilde{\mathbf{Y}}$. At a rate of $\rho_{3}\tilde{\mathbf{Y}}\tilde{\mathbf{V}}$, the pathogen infiltrates and destroys tumor cells, and it reproduces at a rate of $\lambda_{3}\rho_{3}\tilde{\mathbf{Y}}\tilde{\mathbf{V}}$. The configuration $\zeta$ represents the nutrition and microbes soaking constant rate. CTLs intrusion and kill tumor cells at a rate of $\rho_{4}\tilde{\mathbf{Y}}\tilde{\mathbf{Z}}$ and are energized at a rate of $\lambda_{4}\rho_{4}\tilde{\mathbf{Y}}\tilde{\mathbf{Z}}$. The characteristics $\delta_{1}$, $\delta_{2}$, $\delta_{3}$ and $\delta_{4}$ depict the natural mortality rate constants of normal cells, tumor cells, the M1 virus and the natural death rate constant of CTLs, respectively.

The deterministic model accurately describes the propagation of specific pathogens. However, it does not apply to pathogens used in network warfare, and no comprehensive studies on the effect of environmental noise on virotherpy or tumor cell killing have been conducted. These may be referred to as stochastic perturbation occurrences. In reality, stochastic DEs involving standard Brownian motion have received a lot of attention in recent decades, including COVID-19 models ^[52], functional stochastic HBV models ^[53], stochastic tuberculosis models ^[54], stochastic SVI models ^[55] and two-stage model of social insects with egg cannibalism ^[56]. Wei and Fangxiang ^[57] contemplated the stochastic permanence of an SIQS epidemic model with saturated incidence and independent random perturbations. Thus, the stochastic oncolytic M1 virotherapy with immune reaction model under stochastic perturbations takes the following form:

where $B_{\iota}(\xi)\; (\iota = 1, \cdots, 5)$ are two independent standard Brownian motions described on a complete filtered probability space $\{\Omega, \mathfrak{F}, \{\mathfrak{F}_{\xi}\}_{\xi\geq0}, \mathbf{P}\}$ having a $\wp$-filtration $\{\mathfrak{F}_{\xi}\}_{\xi\geq0}$ ^[58]. Also, $\wp_{\iota}, \; \iota = 1, \cdots, 5$ is the intensity of the white noise. we present $\mathbb{R}_{+}^{5} = \{x_{1}\in\mathbb{R}^{3}:{x_{1}}_{\iota} > 0\; for\; all \; 1\leq\iota\leq5\}$ and $\tilde{\mathbb{R}}_{+}^{5} = \{x_{1}\in\mathbb{R}^{3}:{x_{1}}_{\iota}\geq0\; for\; all \; 1\leq\iota\leq5\}$. For an integral function $f_{1}(\xi)$ defined on $[0, \infty)$, we define $f_{1}^{u} = \sup\{f_{1}(\xi)\setminus \xi\geq0\}$ and $f_{1}^{l} = \inf\{f_{1}(\xi)\setminus \xi\geq0\}$. The Itô's strategy is described in the following formula:

The stochastic DE in $\mathfrak{d}$-dimensions is described as follows:

where $\mathbf{f}:\mathbb{R}^{\mathfrak{d}}\times[\xi_{0}, \mathbf{T}]\mapsto\mathbb{R}^{\mathfrak{d}}$ and $\mathbf{g}:\mathbb{R}^{\mathfrak{d}}\times[\xi_{0}, \mathbf{T}]\mapsto\mathbb{R}^{\mathfrak{d}\times m_{1}}$ are Borel measurable with $\mathbf{W} = \{\mathbf{W}(\xi)\}_{\xi\geq \xi_{0}}$ is an $\mathbb{R}^{m_{1}}$-valued Wiener technique, and $\mathbf{u}_{0}$ is an $\mathbb{R}^{\mathfrak{d}}$-valued random variable stated on $\mho$.

Furthermore, $\mathbb{C}^{2, 1}(\mathbb{R}^{\mathfrak{d}}\times[\xi_{0}, \infty); \mathbb{R}_{+})$ is regarded as the collection of all positive mappings $\mathbf{S}(\mathbf{u}, \xi)$ on $\mathbb{R}^{\mathfrak{d}}\times[\xi_{0}, \infty)$ that are continuously twice differentiable in $\mathbf{u}\in\mathbb{R}^{\mathfrak{d}}$ and once in $\xi\in[\xi_{0}, \infty)$. The differential operator $\mathbb{L}$ for the Stochastic DE (2.3) is provided by

Define the mapping $\mathbf{S}\in\mathbb{C}^{2, 1}(\mathbb{R}^{\mathfrak{d}}\times[\xi_{0}, \infty)$, then

where $\mathbf{S}_{\xi}: = \frac{\partial\mathbf{S}}{\partial \xi}; \; \mathbf{S}_{\mathbf{s}_{1}} = (\mathbf{S}_{\mathbf{u}_{\mathbf{p}}}, ..., \mathbf{S}_{\mathbf{u}_{\mathfrak{d}}}), \; \; \mathbf{S}_{\mathbf{u}\mathbf{u}} = (\mathbf{S}_{\mathbf{u}_{\mathbf{p}}}, \mathbf{S}_{\mathbf{u}_{\mathbf{p}}})_{\mathfrak{d}\times \mathfrak{d}}$.

For $\mathbf{u}(\xi)\in\mathbb{R}^{\mathfrak{d}}$, then Itô's approach is defined as

Here, we provide the accompanying description to help readers who are acquainted with fractional calculus (see ^[31,32,33]).

where $\mathbf{Q}(\gamma)$ is defined to be normalized function having $\mathbf{Q}(0) = \mathbf{Q}(1) = 1$.

The formulation of the Atangana-Baleanu derivative is represented below:

where $ABC(\gamma) = 1-\gamma+\frac{\gamma}{\Gamma(\gamma)}$ signifies the normalization function.

3.   Qualitative analysis

Elaiw et al. ^[44] investigated the global characteristics of the nonlinear oncolytic M1 virotherapy having an immune reaction model, explaining that oncolytic virotherapy is a form of chemotherapy that targets and kills cancerous cells using recombination pathogens. Oncolytic alphavirus M1 is a naturally occurring pathogen that has been shown to have high sensitivity and potency in cancer progression in (2.1). All specifications and their interpretations are encapsulated in upcoming sections.

3.1.   Existence-uniqueness of non-negative solutions

The auxiliary outcome is proposed to describe the stochastic framework's (2.2) existence-uniqueness.

Theorem 3.1. Suppose there is a unique solution of the stochastic model (2.2) for $\xi\geq0$ with ICs $\tilde{\mathbf{X}}(0)\in\mathbb{R}_{+}^{5}$. In addition, the solution of $\tilde{\mathbf{X}}(\xi)$ will stay in $\mathbb{R}_{+}^{5}$ with probability $1$, i.e., $\tilde{\mathbf{X}}(0)\in\mathbb{R}_{+}^{5}, \; \forall\; \xi\geq0$, almost surely (a. s.).

Proof. The system's coefficients supposed for the initial values settings $(\tilde{\mathbf{X}}(\xi))\in\mathbb{R}_{+}^{5}$ are continuous and locally lipschtz. Consequently, the system $(\tilde{\mathbf{X}}(\xi))$ have only one solution for $\xi\in[0, \theta_{\varepsilon})$. For the explosive period $\theta_{\varepsilon}$ is thoroughly examined in ^[58]. In order to show the solution's diverse nature, we must prove that $\theta_{\varepsilon} = \infty$ (a.s.). Assume that we do have a somewhat large positive number $\psi_{0}$ such that every state's ICs fall inside the given interval $\big[\psi_{0}, \frac{1}{\psi_{0}}\big]$. Choosing $\psi\geq\psi_{0}$ be the terminal duration specification for each non-negative integer.

Throughout this investigation, we will employ $\inf \phi = \infty$, whilst $\phi$ refers empty set. The idea of $\psi$ compels us to claim that it increases as $\psi$ approaches $\infty$. Fixing $\theta_{\infty} = \lim\limits_{\psi\mapsto\infty}\theta_{\varepsilon}\geq\theta_{\infty}$ (a.s.). After verifying that $\theta_{\infty} = \infty$ (a.s.), We shall argue that $\theta_{\varepsilon} = \infty$ and thus $\tilde{\mathbf{X}}(\xi)$ stayed in $\mathbb{R}_{+}^{5}$ a.s $\forall\; \xi\geq0$. So, verifying that $\theta_{\infty} = \infty$ (a.s.). For this, we suppose two non-negative fixed values $\varepsilon\in(0, 1)$ and $\mathbf{T}$ must exist such that

So, the integer $\psi_{1}\geq\psi_{0}$ exists in the subsequent version

So that, we shall investigate a mapping $\mathcal{J}:\mathbb{R}_{+}^{5}\mapsto\mathbb{R}_{+}$ in the following manner:

The $\mathcal{J}$ is a positive function, which should be noticed and may be confirmed by the argument that $0\leq \mathbf{u}_{1}-\ln \mathbf{u}_{1}-1, \; \forall \mathbf{u}_{1} > 0$. Suppose the arbitrary terms $\psi_{0}\leq\psi$ and $\mathbf{T} > 0$.

Employing Itô's technique to (3.2) yields

In view of (3.4), let us introduce the subsequent functional $\mathbf{L}\mathcal{J}:\mathbb{R}_{+}^{5}\mapsto\mathbb{R}_{+}$ described as

Since $\omega$ is a non-negative constant and it is free of $\tilde{\mathbf{X}}$ as well as of $\xi$. Thus, we have

Therefore, we have

Letting $\Psi_{\psi} = \big\{\theta_{\psi}\leq\mathbf{T}\big\}$ for $\psi\geq\psi_{1}$ and (3.1) yields that $P(\Psi_{\kappa})\geq\epsilon$. Clearly, for every $\omega$ from $\Psi_{\psi}$, there exists at least one $\tilde{\mathbf{H}}(\theta_{\psi}, \omega)$, $\tilde{\mathbf{N}}(\theta_{\psi}, \omega)$, $\tilde{\mathbf{Y}}(\theta_{\psi}, \omega)$, $\tilde{\mathbf{V}}(\theta_{\psi}, \omega)$ and $\tilde{\mathbf{Z}}(\theta_{\psi}, \omega)$ which are equal to $\frac{1}{\psi}$ or $\psi$. Hence, $\mathcal{J}\big(\tilde{\mathbf{H}}(\theta_{\psi}), \tilde{\mathbf{N}}(\theta_{\psi}), \tilde{\mathbf{Y}}(\theta_{\psi}), \tilde{\mathbf{V}}(\theta_{\psi}), \tilde{\mathbf{Z}}(\theta_{\psi})\big)$ is no less than $\ln\psi-1+\frac{1}{\psi}$ or $\psi-1-\ln\psi$.

Consequently,

Using the fact of (3.1) and (3.5), we have

As seen, the indicator mapping of $\Psi$ is $1_{\Psi{(\omega)}}$. As $\psi\mapsto\infty,$ yields the contradiction $\infty > \mathcal{J}\big(\tilde{\mathbf{X}}(0)\big)+\Omega\mathbf{T} = \infty,$ which shows that $\theta_{\infty} = \infty$ a.s.

3.2.   Extinction and persistence

One of the critical focuses in disease transmission is how to control illness's complexities so that the illness becomes extinct and persists over time. In this section, we attempt to determine the significance level for ailment extinction and persistence.

This portion will examine the extinction and ergodic stationary distributions (ESD) of the framework in discussion. Let us summarize

Next, we will present the well-known result of the strong law of large numbers, which is mainly due to ^[59].

Lemma 3.1. ^[59] Suppose that there is a continuous and real-valued local martingale $\mathbf{Q} = \{\mathbf{Q}\}_{\xi\geq0}$, which disappears as $\xi\mapsto 0$, then

Let us classify one more threshold parameter for our upcoming requirements:

Theorem 3.2. For $\tau > \frac{1}{2}(\wp_{1}^{2}\vee\wp_{2}^{2}\vee\wp_{3}^{2}\vee\wp_{4}^{2}\vee\wp_{5}^{2})$ and let $\tilde{\mathbf{X}}(\xi)$ has a non-negative solution of the model (2.2) having the initial settings $\tilde{\mathbf{X}}(0)\in\mathbb{R}_{+}^{8}$, we find

(i) If $\mathbb{R}_{0}^{p} < 1$, then

which indicate that the disease will be exterminated in a long run.

(ii) If $\mathbb{R}_{0}^{p} > 1$, then

where $\Lambda_{1}: = \bigg(\frac{\rho_{3}\tau}{\zeta+\delta_{3}}+\frac{(\zeta+\delta_{1})^{2}}{\lambda_{2}}+\frac{\tau(\tau+\zeta+\delta_{1})}{\lambda_{3}}\bigg)$, which indicate that the disease will be persisted in a long run.

Furthermore,

Proof. (ⅰ) Implementing integration on (2.2), we have

Utilizing Itô's strategy on $\ln(\tilde{\mathbf{Y}}(\xi))$, we get

After integrating the aforesaid equation from 0 to $\xi$ on both sides, we have

In accordance with the strong law of large numbers ^[60], we have $\lim\limits_{\xi\mapsto\infty}\frac{1}{\xi}\int\limits_{0}^{\xi}dB_{3}(\mathbf{s})d\mathbf{s} = 0 \; (a.s.)$.

Attempting to take the superior limit and applying the stochastic comparison theorem, we get

Thus, it implies that $\lim\limits_{\xi\mapsto\infty}{\tilde{\mathbf{Y}}(\xi)} = 0 \; (a.s.)$.

After simplification, (3.8) reduces to

Thus, it implies that $\lim\limits_{\xi\mapsto\infty}{\tilde{\mathbf{H}}(\xi)} = \frac{\chi}{\zeta} \; (a.s.)$.

Similarly, we can show that $\lim\limits_{\xi\mapsto\infty}\tilde{\mathbf{N}}(\xi) = 0, \; a.s., \; \lim\limits_{\xi\mapsto\infty}\tilde{\mathbf{V}}(\xi) = \frac{\tau}{\zeta+\delta_{1}}, \; a.s., \; and\; \lim\limits_{\xi\mapsto\infty}\tilde{\mathbf{Z}}(\xi) = 0 \; a.s.$

This means that the disease extinction based on the value of the parameter $\mathbb{R}_{0}^{s}$, that is, if $\mathbb{R}_{0}^{s} < 1$, consequently the disease will extinct out in the long run.

(ⅱ) Define a $\mathbb{C}^{2}$-function $\mathcal{U}_{1}$ as

Then,

Taking $\Lambda_{1} = \bigg(\frac{\rho_{3}\tau}{\zeta+\delta_{3}}+\frac{(\zeta+\delta_{1})^{2}}{\lambda_{2}}+\frac{\tau(\tau+\zeta+\delta_{1})}{\lambda_{3}}\bigg)$. Consequently,

Integrating both sides of (3.9), we have

where $\mathbf{Q}(\xi) = \int\limits_{0}^{\xi}\wp_{3}dB_{3}(\mathbf{s})$ is a martingale. Utilizing Lemma 3.1, we have $\lim\limits_{\xi\mapsto\infty}\frac{\mathbf{Q}(\xi)}{\xi} = 0\; (a.s.)$.

Considering Lemma 3.1, we find from (3.10)

3.3.   Ergodicity and stationary distribution

Despite the absence of an endemic equilibrium point in the stochastic model (2.2), we intend to investigate the presence of an ergodic stationary distribution (ESD) (stochastic non-negative steady state), which further demonstrates ailment persistence. Initially, we present a few outcomes from Has'minskii's concept. Readers can find more information at ^[59].

Assume that $\tilde{Q}(\xi)$ be a homogeneous Markov procedure in $\zeta$ (the $\mathfrak{d}$-dimensional Euclidean space) that effectively addresses the stochastic DE below:

The diffusion matrix $A_{1}(\mathbf{y}) = (a_{\iota\Bbbk}(\mathbf{y}))$ and $a_{\iota\Bbbk}(\mathbf{y}) = \sum\limits_{\kappa = 1}^{n_{1}}g_{\kappa}^{(\iota)}(\mathbf{y})g_{\kappa}^{(\Bbbk)}(\mathbf{y})$.

Lemma 3.2. ^[59] Suppose there is a bounded domain $\mathcal{U}\subset\Lambda_{\mathfrak{d}}$ having a regular boundary $\Gamma$ such that

($Z_{1}$) There is a non-negative number $\mathbf{Q}$ such that $\sum\limits_{\iota, \Bbbk = 1}^{\mathfrak{d}}a_{\iota\Bbbk}(\mathbf{y})\xi_{\iota}\xi_{\Bbbk}\geq \mathbf{Q}\vert \xi\vert^{2}, \; \mathbf{y}\in\mathcal{U}, \; \xi\in\mathbb{R}^{\mathfrak{d}}$.

($Z_{2}$) There exists a positive $\mathbb{C}^{2}$-function $\mathcal{H}$ such that $\mathcal{L}\mathcal{H}$ is negative for any $\mathbf{y}\in\Lambda_{\mathfrak{d}}\setminus\mathcal{U}$ (particularly $\mathcal{L}\mathcal{H}\leq-1$, for every $\mathbf{y}\in\Lambda_{\mathfrak{d}}\setminus\mathcal{U}$), then the Markov technique $\tilde{Q}(\xi)$ has a unique ESD $\pi(.)$, and

satisfies $\forall\; \mathbf{y}\in\Lambda_{\mathfrak{d}}$, where $f_{1}(.)$ is an integrable mapping respecting to the measure $\pi$.

Define a parameter

Furthermore, we will prove prerequisites that ensure the existence of an ESD depending on Has'minskii's hypothesis.

Theorem 3.3. For $\mathbb{R}_{0}^{s} > 1$, then for an initial setting $\bar{\mathbf{y}(0)}\in\mathbb{R}_{+}^{8}$, the model (2.2) has a unique stationary distribution $\pi(.)$ and has the ergodic property.

Proof. To demonstrate Theorem 3.3, simply check assumptions ($Z_{1}$) and ($Z_{2}$) in Lemma 3.2. To begin, we assemble an appropriate Lyapunov function $\mho$ and consider a closed set $\mathcal{U}_{\epsilon}\setminus\mathbb{R}_{+}^{8}$ such that $\sup_{\tilde{\mathbf{X}}\in\mathcal{U}_{\epsilon}\setminus\mathbb{R}_{+}^{8}}\mathcal{L}\mho$ is negative to determine the viability of ($Z_{2}$) in Lemma 3.2.

Now we intend to define a positive $\mathbb{C}^{2}$-function $\mathcal{H}:\mathbb{R}_{+}^{5}\mapsto\mathbb{R}_{+}$ as

here, the non-negative constants $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$ must be computed later. We achieve the accompanying outcomes by using Ito's technique and the proposed model (2.2) as

Then we have

This yields that

Taking

Denoting

Thus, we have

Furthermore, we have

Here, $\eta_{4} > 0$ is a constant which will be determined later. Therefore, it is important to present

here, $\mathcal{U}_{\kappa} = (\frac{1}{\kappa}, \kappa)\times(\frac{1}{\kappa}, \kappa)\times(\frac{1}{\kappa}, \kappa)$. Furthermore, we illustrate that $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$ has unique minimum value $\mathcal{H}_{2}(\tilde{\mathbf{X}}(0))$.

The partial derivative of $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$ in respect to $\tilde{\mathbf{H}}, \; \tilde{\mathbf{N}}, \; \tilde{\mathbf{Y}}, \; \tilde{\mathbf{V}}, \; \tilde{\mathbf{Z}}$ is as follow:

It is easily demonstrated that $\mathcal{H}_{2}$ has a distinct stagnation point, which is as determined by the following:

Furthermore, the Hessian matrix of $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$ at $\tilde{\mathbf{X}}(0)$ is

It is clear that the aforesaid matrix is positive definite. Consequently, $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$ has a least value of $\mathcal{H}_{2}(\tilde{\mathbf{X}}(0))$.

Utilizing (3.22) and applying the continuity of $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$, we can observe that $\mathcal{H}_{2}(\tilde{\mathbf{X}}(\xi))$ has one or more value $\mathcal{H}_{2}(\tilde{\mathbf{X}}(0))$ stayed in $\mathbb{R}_{+}^{5}$.

Accordingly, we present a positive $\mathbb{C}^{2}$-function $\mathcal{H}_{3}:\mathbb{R}_{+}^{5}\mapsto\mathbb{R}_{+}$ as follows:

In view of Ito's formula and suggested system (2.2), we have

as a result, the previous hypothesis can be formulated as

where $\eta_{5} = 4\chi\big[(\mathbb{R}_{0}^{s})^{1/4}-1\big] > 0$. The definition of a set is provided by

where $\epsilon_{\iota}, \; \iota = 1, 2$ are constants, which are extremely small and will must be discovered later. We will split the whole $\mathbb{R}_{+}^{5}\setminus\mathcal{W}^{\diamond}$ into ten domains as follows:

Furthermore, we will prove that $\mathcal{L}\mathcal{H}_{3}((\tilde{\mathbf{X}}(\xi)) < 0$, which is analogous as showing it on the ten domains as particularized previously.

Case Ⅰ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{1}$, so using (3.28) yields

Selecting $\epsilon_{1} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{1}$.

Case Ⅱ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{2}$, so using (3.28) yields

Selecting $\epsilon_{1} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{2}$.

Case Ⅲ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{3}$, so using (3.28) yields

Selecting $\epsilon_{1}, \epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{3}$.

Case Ⅳ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{4}$, so using (3.28) yields

Selecting $\epsilon_{1} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{4}$.

Case Ⅴ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{5}$, so using (3.28) yields

Selecting $\epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{5}$.

Case Ⅵ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{6}$, so using (3.28) yields

Selecting $\epsilon_{1}, \epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{6}$.

Case Ⅶ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{7}$, so using (3.28) yields

Selecting $\epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{7}$.

Case Ⅷ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{8}$, so using (3.28) yields

Selecting $\epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{8}$.

Case Ⅸ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{9}$, so using (3.28) yields

Selecting $\epsilon_{1} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{9}$.

Case Ⅹ. If $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{10}$, so using (3.28) yields

Selecting $\epsilon_{2} > 0$, produces $\mathcal{L}\mathcal{H}_{3}(\tilde{\mathbf{X}}(\xi)) < 0$ for every $(\tilde{\mathbf{X}}(\xi))\in\mathcal{W}^{\diamond}_{10}$. As a result, we prove that a constant $\mathcal{P} > 0$ is one that guarantees

So that

Suppose that $((\tilde{\mathbf{X}}(\xi))(0)) = (\mathbf{u}_{1}, \mathbf{u}_{2}, \mathbf{u}_{3}, \mathbf{u}_{4}, \mathbf{u}_{5}) = \mathbf{u}\in\mathbb{R}_{+}^{5}\setminus\mathcal{W}^{\diamond}$ and $\theta^{\mathbf{u}}$ is the time period for which a path begin at $\mathbf{u}$ leads to the set $\mathcal{W}^{\diamond}$,

The subsequent result can be obtained by integrating the two hand sides of the variant (3.31) from $0$ to $\theta^{(n)}(\xi)$, attempting to take expectation, and applying Dynkin's computation.

Since $\mathcal{H}_{3}(\mathbf{u})$ is a positive number, hence

Therefore, we have $\mathbf{P}\{\theta_{\epsilon} = \infty\}$ as an immediate consequence of Theorem 3.3.

Conversely, the framework defined in (2.2) can be characterized as regular. In light of this, if we choose $\xi\mapsto\infty$ and $n\mapsto\infty$, we will almost surely get $\theta^{(n)}(\xi)\mapsto\theta^{\mathbf{u}}$.

Consequently, utilizing Fatou's lemma, we achieve

Obviously, $\sup\limits_{\mathbf{u}\in\mathcal{C}}\mathbf{E}\theta^{\mathbf{u}} < \infty$, here $\mathcal{C}\in\mathbb{R}_{+}^{5}$ is a compact subset. It validates assumption ($Z_{2}$) of Lemma 3.1.

Moreover, the diffusion matrix of the system (2.2) is

Selecting $\mathcal{M} = \min_{(\tilde{\mathbf{X}}(\xi))\in\bar{\mathcal{W}^{\diamond}}\in\mathbb{R}_{+}^{5}}\big\{\wp_{1}^{2}\tilde{\mathbf{H}}^{2}, \wp_{2}^{2}\tilde{\mathbf{N}}^{2}, \wp_{3}^{2}\tilde{\mathbf{Y}}^{2}, \wp_{4}^{2}\tilde{\mathbf{V}}^{2}, \wp_{5}^{2}\tilde{\mathbf{Z}}^{2}\big\}$, we find

where $r = (r_{1}, r_{2}, r_{3}, r_{4}, r_{5})\in\mathbb{R}_{+}^{5}$. This means that assumption ($Z_{1}$) of Lemma 3.1 is also valid.

According to the investigation that preceded before, Lemma 3.1 shows that the framework (2.2) is ergodic and has a single stationary distribution.

4.   Determination of the model with piecewise derivative

In this section, we introduce the piecewise derivative when the entangled derivatives are the classical and fractional differential operators considering singular and non-singular kernels. The fractional order will be between 0 and 1, that is, $0 < \gamma\leq1$.

4.1.   Caputo fractional derivative operator

In this part, we will look at the dynamics of the oncolytic M1 virotherpy model (2.1) and (2.2) that incorporate immune response, encompassing classical, index-law and ultimately stochastic procedures. If we describe $\mathbf{T}$ as the final time of transmission, that is, the final time, then the mathematical framework will be developed in the first round using the classical derivative formulation, then the index-law kernel in the second step, and finally the stochastic environment in the later phases. The mathematical formalism that explains this phenomenon is then presented as:

Here, we apply the technique described in ^[37] for the scenario of Caputo's derivative to analyze quantitatively the piecewise structure (4.1)–(4.3). We commence the technique as follows:

It follows that

where

and

4.2.   Caputo-Fabrizio fractional derivative operator

In this subsection, we will examine the system of DEs of the oncolytic M1 model involving immune response in the community, comprising classical, exponential decay law and stochastic mechanisms. If we describe $\mathbf{T}$ as the concluding time of transmission, that is, the final time, then the mathematical structure will be formed in the first round using the classical derivative implementation, then the exponential decay kernel in the second step, and eventually the stochastic environment in the subsequent periods. Regarding that, the mathematical approach used to illustrate this occurrence is presented as follows:

Here, we apply the technique described in ^[37] for the scenario of Caputo-Fabrizio derivative to analyze quantitatively the piecewise structure (4.7)–(4.9). We commence the technique as follows:

It follows that

4.3.   Atangana-Baleanu fractional derivative operator

Here, we will concentrate on the dynamics of the oncolytic M1 model in this portion, which demonstrates immue response and CTL cells efficacy in population, including classical, generalized Mittag-Leffler law, and lastly, stochastic causes. If we define $\mathbf{T}$ as the final time, the mathematical configuration will be constituted in the first round employing the classical derivative application, followed by the Mittag-Leffler kernel in the second step, and finally the stochastic environment in subsequent periods. In this regard, the mathematical model utilized to describe this phenomenon is as follows:

Here, we apply the technique described in ^[37] for the scenario of Atanagan-Baleanu-Caputo derivative to analyze quantitatively the piecewise structure (4.12)–(4.14). We commence the technique as follows:

It follows that

where $\Im_{1}, \Im_{2}$ and $\Im_{3}$ are stated before in (4.4)–(4.6).

4.4.   Mathematical formulation of framework

In this section, we will display numerical simulation models to validate the theoretical predictions using the Atangana and Araz techniques mentioned earlier ^[37]. Several numerical results are taken into account to demonstrate the appropriateness and usefulness of the suggested oncolytic M1 framework for the deterministic-stochastic context. MATLAB 21 software was used to accomplish all representational and numerical simulations.

Focused on the configuration values and taking the initial settings $(\tilde{\mathbf{H}}(0), \tilde{\mathbf{N}}(0), \tilde{\mathbf{Y}}(0), \tilde{\mathbf{V}}(0), \tilde{\mathbf{Z}}(0)) = \big(0.3(1+0.2\cos^{2}(\pi x_{1})),$$0.2(1+0.2\cos^{2}(\pi x_{1})), 0.1(1+0.2\cos^{2}(\pi x_{1})), 0.1(1+0.2\cos^{2}(\pi x_{1})),$$0.01(1+0.2\cos^{2}(\pi x_{1}))\big), x_{1}\in[0, 2]$. An insightful analysis based on the representations of $\mathbb{R}_{0}^{s}$ and $\mathbb{R}_{0}^{p}$ clearly demonstrates that $\mathbb{R}_{0}\geq\mathbb{R}_{0}^{p}\geq\mathbb{R}_{0}^{s}$, when merged with paper ^[44]. From Theorems 3.2 and 3.3, we can conclude that $\mathbb{R}_{0}^{s} > 1$ can be considered as a consolidated criterion for population persistence of oncolytic M1 virus in deterministic (2.1) and stochastic system (2.2). Analogously, there is no immune reaction, the oncolytic M1 virotherapy completely eradicates the tumor, repairing normal tissue and enhancing the patient's condition, such as $\mathbb{R}_{0}^{p} < 1$, a cohesive criterion for framework population extinction (2.1) can be considered (2.2). As a consequence, we classify Elaiw et al. ^[44] nonlinear dynamic results.

Besides that, it is potentially proven that model (2.2) has a unique stationary distribution when $\mathbb{R}_{0}^{s} > 1$ (see Theorem 3.2), but oncolytic M1 virus will become extinct if $\mathbb{R}_{0}^{p} < 1$ (see Theorem 3.3). At the conclusion of this part, we will demonstrate our simulation values using piecewise fractional differential operators numerical method ^[37]. For the sake of simplicity, we make the assumption that all white noises have the distinct intensities $\wp_{1} = 0.08$, $\wp_{2} = 0.082$, $\wp_{3} = 0.009$, $\wp_{4} = 0.0092$, $\wp_{5} = 0.1$. Then we derived $\mathbb{R}_{0}^{s} = \frac{\zeta\lambda_{2}\rho_{2}}{\big(\frac{\wp_{1}^{2}}{2}+\zeta\big)\big(\zeta+\delta_{1}+\frac{\wp_{2}^{2}}{2}\big)\big(\zeta+\delta_{2}+\frac{\wp_{3}^{2}}{2}\big)}$ and $\mathbb{R}_{0}^{p} = \frac{{\lambda_{2}\rho_{2}\chi(\zeta+\delta_{2})}-{\rho_{3}\tau \zeta}}{\zeta(\zeta+\delta_{2}+\frac{\wp_{3}^{2}}{2})(\zeta+\delta_{3})}$. We set the parameters $(\chi, \rho_{1}, \rho_{2}, \rho_{3}, \rho_{4}, \lambda_{1}, \lambda_{2}, \lambda_{3}, \lambda_{4}, \zeta, \tau, \delta_{1}, \delta_{2}, \delta_{3}) =$$(0.02, 0.03, 0.03, 0.1, 0.03, 0.8, 0.8, 0.5, 0.8, 0.02, 0.01, 0.04, 0.01, 0.008, 0.01)$. Figures 1–4 display the variability patterns of $\mathbb{R}_{0}^{s}$ and $\mathbb{R}_{0}^{p}$ with distinct random perturbations.

Example 4.1. For findings provided in Section 4, we computed $\mathbb{R}_{0}^{s} = 1.2301 > 1$. According to Theorem 3.3, the solution $\tilde{\mathbf{X}}(\xi)$ of system $(2.2)$ has a unique stationary distribution, as illustrated in Figures 1–4.

Following that, these frameworks are aggregated utilizing the core idea of piecewise derivative, allowing us to visualize the entire underlying virus of an individual who discovers he has an oncolytic M1 virus, continues therapies, and then describes the immune response workflow after intervention. The classical derivative is employed to simulate the patient's illness educational procedure, which occurs in the first time-frame.

Because the therapy practise is the disease's second practise, a reduction in dynamics of the model with the optimistic consequence of immune response may be an increase in the patient's innate immunity, so the stochastic derivative is employed in the second time interval. The final practise is analyzed in three instances: Caputo (4.1)–(4.3), Caputo-Fabrizio (4.7)–(4.9) and Atangana-Baleanu (4.12)–(4.14) fractional derivative. Due to the possibility of a marked reduction and cancellation being successfully treated, these features have been preferred because there will be a slowing following a dramatic reduction. Because no medication is administered in the first interval, there is occasionally a decrease as well as an increase during intervention.

According to the above investigation, when the immune reaction is absent, the oncolytic M1 virotherapy completely eradicates the tumor, repairing healthy tissue and striving to improve the patient's care. As a result, the M1 virus regulates the tumor with it being influenced by the innate immunity against tumor, calculating the lowest inhibitory therapeutic dose forced to separate the tumor ^[1].

Example 4.2. A straightforward computation for model (2.2) reveals that $\mathbb{R}_{0}^{p} = 0.7834 < 1$ and $\mathbb{R}_{0} = 2.8571 > 1$. According to Theorem 3.3 and work ^[44] model (2.1) has a global asymptotic stability non-negative steady-states, implying the population persistence of oncolytic M1 virus, but framework (2.2) has a extinction probability one, as represented by Figure 5.

5.   Conclusions

In this document, we analyzed the complexities of a deterministic-stochastic oncolytic virotherapy framework with a CTL immune system. The goal was to see how the CTL immune system to tumor cells affected the efficiency of oncolytic M1 virion therapeutic interventions. Following that, these frameworks are aggregated utilizing the core idea of piecewise derivative, allowing us to visualize the entire disease progression of an individual who discovers they have an oncolytic M1 virus, continues to receive therapies, and then describes the disease progression. In view of the computed findings, we procure a critical value for the extinction of tumor cells, free M1 virus and CTL cells, which is presented by $\mathbf{R}_{0}^{s} < 1$. Significant linear perturbations clearly result in $\mathbb{R}_{0}^{s} < 1$. As a result, the intermittent nature of untreated specific receptors has a significant impact on virus infection elimination. We contend that, while the generalized Mittag-Leffler kernel, exponential decay and power law have been revealed to be capable of depicting several crossover behaviours, their strengths to accomplish this may be restricted due to the enormity of nature.

We can conclude that variability can eliminate viral transmissions, whereas small random variations induce infectious diseases to persist. Furthermore, we propose using therapeutic interventions and immune responses based on $\mathbb{R}_{0}^{s}$ and $\mathbb{R}_{0}^{p}$ representations to increase the mortality rate of insidiously tumor cells, free M1 virus and CTL cells. Moreover, in a broad sense, the immune reaction cannot be considered hazardous as it can target cancer cells while protecting immune tissues from extermination. In actuality, one of the most energetic analysis areas is the impact of immune function on oncolytic viruses. Noticeably, despite impressive outcomes in drug development, oncogenic virotherapy is still viewed as a specific cancer treatment ^[50,61]. Introducing chemotaxis, Lévy noises, and time delays to the framework (2.2) may provide a profound understanding of the model, which can be actually achieved in future research.

Conflict of interest

The authors declare that they have no competing interests.