Stochastic tumor-immune interaction model with external treatments and time delays: An optimal control problem

1.   Introduction

Cancer is one of the most dangerous diseases in the world, and it affects both developing and developed nations. According to the World Health Organization, cancer is a prominent contributor to global mortality, with approximately 10 million deaths recorded in the year 2020 ^[1]. Through a multi-step or multifaceted process, cancer arises when normal cells turn into cancerous cells. It usually progresses from a pre-cancerous lesion to a fully cancerous cell. In this multistep mechanism, various genetic and molecular changes accumulate over time, ultimately leading to cancer. These transformations are the consequence of the interaction among a person's genetic circumstances and three classifications of external forces, such as physical, chemical and biological carcinogens ^[2]. An early diagnosis of cancer increases the likelihood of a favorable response to effective treatment, leading to a higher survival rate, reduced morbidity, and a more cost-effective treatment option. One of the most important and challenging questions is to understand cancer biology and how to treat it effectively. Therefore, many researchers have put time and efforts in creating and improving an efficient medication. Additionally, identifying the methods to support the patient improve their natural immunity that implement them to fight against tumor cells.

Sophisticated mathematical models are necessary to accurately depict the intricate interplay between tumor cells and the immune system. Different mathematical models have been investigated and analyzed to recognize tumor-immune dynamics ^{[3,4,5,6,7,8,9,10,11,12,13]}. In ^[14], Kirschner and Panetta introduces the tumor-immune interaction model of three types of cell populations in vitro. The interaction between effector cells and tumor cells was mathematically modeled by Kuznetsov et al. ^[15] in one of the earliest proposals of such a model. An optimal control problem (OC-P) is typically formulated as a mathematical optimization problem, in which the objective is to find the control trajectory that minimizes or maximizes the cost function while satisfying any constraints on the state variables or control inputs of the system ^[16,17]. According to Rihan et al. ^[18] a study was carried out using an OC-P model for a delay differential model to examine tumor-immune interactions under immuno-chemotherapy.

In actuality, deterministic models are not adequate to describe the dynamic process of the proliferation of cancer cells. Research has been conducted in this instance to amplify the deterministic models to stochastic counterpart ^{[19,20,21,22,23,24]}. Tumor-immune interactions are extremely complex, so it makes sense to include noise on the tumor-immune system in order to take into account a variety of relevant phenomena. An example of this may be variations in the intensity of neoantigens that stimulate the immune response, or changes in the expression of molecules that activate T cells; see ^[22,23,25]. As a result of the significant uncertainty inherent in the treatment process, authors in ^[26,27,28] modified the deterministic model to incorporate stochastic factors, such as variations in cellular reproduction and death rates, changes in the immune system's ability to fight tumor cells, and fluctuations in chemotherapy efficacy. In addition, the tumor-immune interaction model illustrates a delay between an immune response and a subsequent response by cancer cells ^[29]. A delay may affect the stability of the system, resulting in instability or bifurcation. In order to formulate effective therapies, it is crucial to understand how time-delays play a role in tumor-interaction models ^[29].

Up to the best of authors' knowledge, few studies have been conducted so far in tumor-immune interactions with stochastic noises and time delays involving immunological boosts and chemotherapeutic treatment therapies. There is very little research relating stochastic models to interacting cells in cancer dynamics ^[30,31]. Drawing motivation from the aforementioned research studies, this paper attempts to investigate the impact of two time-delays and stochastic white noises on both the dynamics and optimal control performance of a model representing the interaction between tumors and the immune system. We provide a stochastic optimal control problem to reduce tumor growth and reduce the load tumors and increase the size of effector cells.

This paper is structured as follows: Section 2 introduces a stochastic tumor-immune model with external treatments and time delay. Section 3 discusse the time-delayed stochastic model's global positive solution. The stationary distribution and extinction of tumor cells are discussed in Section 4. Section 5 examines the OC-P governed by the stochastic model using the stochastic maximal criterion. Section 6 presents numerical simulations incorporating various stochastic perturbations in order to validate the theoretical findings. Our conclusion is presented in Section 7.

Preliminaries

Here, we provide some preliminaries for the discussion that follows. Assume that $\left(\Omega, {\cal M}, \left\lbrace {\cal M}\right\rbrace _{t\geq0}, {\cal P}\right)$ is a complete probability space with a filtration $\{ {\cal M}_t\}_{t\geq0}$ that meets the usual criteria. It is continuous with the right-hand side and $\mathcal{M}_0$ contains all $\mathcal{P}$-null sets.

Assume $Z(t)$ is a regular time-homogeneous Markov process in $C([-\tau, 0];\mathbb{R}_+^m)$ and satisfies the following stochastic delay differential equation:

with the initial value $Z(s) = Z_0\in C([-\tau, 0];\mathbb{R}_+^m).$ $W(t)$ stands for the m-dimensional standard Brownian motion defined on the complete probability space $\left(\Omega, {\cal M}, \left\lbrace {\cal M}\right\rbrace _{t\geq0}, {\cal P}\right)$. Diffusion matrix $Z(t)$ is then defined as follows:

In order to have a nonnegative Lyapunov function $V(Z_t, t)$, it must be continuously twice differentiable in $C^{2, 1}(C([-\tau, 0];\mathbb{R}_+^m)\times[-\tau, \infty); \mathbb{R}_+)$ and once differentiable in $t$. Reference [21] defines the differentiable operator ${\cal L}$ of (1.1).

2.   Stochastic model for tumor-immune interactions with external treatments

Stochastic perturbations can influence the activation of immune cells, and it can affect the immune system's ability to detect and eliminate tumor cells. In addition, environmental fluctuations can impact the response to anti-cancer treatments, including immunotherapies. These perturbations can affect treatment efficacy, resistance development and long-term outcomes. First, we modify Kirschner and Panetta's model ^[14] which describes the dynamics of the activated effector cells $E(t)$, tumor cells $T(t)$ and Interleukin-2 (IL-2) cells $I(t)$, with time delays $\tau_1$ and $\tau_2$. IL-2 plays a crucial role in regulating immune responses and promoting T cell proliferation and survival, it is also essential for generating effector and memory T cells. We consider the following assumptions: (ⅰ) The Lotka-Volterra form represents cell interaction. (ⅱ) The cancer cells are eradicated by the effector cells, while the IL-2 level increases as competition between tumor cells and effector cells increases. (ⅲ) Two distinct time-delays are incorporated, $\tau_1$ and $\tau_2$. The modified model takes the form

with the initial conditions $E(\theta) = \sigma_1(\theta), T(\theta) = \sigma_2(\theta), I(\theta) = \sigma_3(\theta),$ where $\sigma_p(\theta) \geq 0$ for $\theta \in [-\tau, 0), \tau = \max\left\lbrace \tau_1, \tau_2 \right\rbrace$ and $\sigma_p(0) > 0, p = 1, 2, 3.$ Tumor cell populations are reduced as a result of their interactions with effector cells $E(t)$, which occur at distinct rates represented by $-\xi E(t)T(t)$. By interacting with tumor effector cells $-\xi E(t)T(t)$, immune effector cells $E(t)$ reduce the population of tumor cells $T(t)$ at a rate denoted by $\xi$. Effector cells are activated by IL-2 by stimulating their proliferation and differentiation, with rate $\mu_1E(t)I(t)$. Effector cells are also stimulated by tumor cells, denoted by $\alpha E(t-\tau_1)T(t-\tau_1)$, where $\alpha$ represents the antigenicity rate of the tumor and $\tau_1$ is the duration over which IL-2 is induced in effector cells. $\eta_1 E(t)$ represents the natural degradation of effector cells, while $e_1$ represents the exogenous supply of effector cells. Proliferation and death of tumor cells are represented by the parameter $r$. $\beta$ represents the biological environment's maximum carrying capacity for tumor cells $T(t)$. $\eta_2E(t-\tau_2)T(t-\tau_2)$ represents the effector cells as the source of IL-2, which is induced by their interaction with tumor cells. Effectors and cancer cells compete at $\eta_2$. The time delay between stimulating effector cells and tumor cells is represented by $\tau_2$. The parameter $\mu_2$ represents the degradation of IL-2 inherent to the system, whereas the variable $e_2$ represents an exogenous IL-2 inflow. The parameters in model (2.1) are summarized in Table 1. See the Appendix, for the analysis of the deterministic model (2.1).

Most biological phenomena are characterized by random fluctuations, particularly variations in the intensity of neoantigens that trigger immune responses or the expression of molecules that activate T cells. Due to the complexity of the interaction between tumor cells and immune effectors, noise related to the tumor-immune system is justified. Since tumor-immune cells may have small populations, the impact of random fluctuations is significant. To investigate the impact of stochastic fluctuations within a deterministic model, there are primarily two methodologies ^{[32,33,34,35]}. One or more relevant parameters can be substituted with their stochastic counterparts in a deterministic model. In this process, white or colored noises are introduced to the deterministic parameters to introduce stochastic perturbations. Alternatively, it is possible to introduce a stochastic driving force directly into the deterministic model instead of modifying specific parameters. In the current study, we adopt the second alternative approach, wherein we assume the occurrence of stochastic perturbations in the variables, throughout the boundary and the interior equilibrium point. The perturbations are assumed to follow a white noise distribution, whose magnitude is proportional to the distances between equilibrium values $E(t)$, $T(t)$, and $I(t)$. The model structure is as follows:

with initial conditions

$W_i (i = 1, 2, 3)$ stand for the independent Brownian motions that are stated on a complete probability space $\left(\Omega, {\cal M}, \left\lbrace {\cal M}\right\rbrace _{t\geq0}, {\cal P}\right)$ with a filtration $\{ {\cal M}_t\}_{t\geq0}$ satisfying the usual conditions and $\mathcal{M}_0$ contains all $\mathcal{P}$-null sets. Let $\mathbb{R}_+ = (0, +\infty), \mathbb{R}_+^m = \left\lbrace (y_1, y_2, \ldots, y_m) \in \mathbb{R}^m | y_p > 0, p = 1, 2, \ldots, m\right\rbrace,$ $H_i = H_i(E, I, T)$, $(i = 1, 2, 3)$ are locally Lipschitz-continuous functions, such that $H_i$ has the following forms:

● In the first form $R_1$: $H_1 = E$, $H_2 = T$, and $H_3 = I$, i.e., the environmental impact on the cell described by stochastic perturbation ^[36].

● In the second form $R_2$: $H_1 = E-E^0$, $H_2 = T-T^0$, and $H_3 = I-I^0$, where $\mathcal{S} = (E^0, T^0, I^0)$ is the equilibrium point of (2.1), which investigates the behavior of stochastic perturbations around the equilibrium point ^[37].

When examining the dynamics of the proposed system, the positivity of the solution is of primary importance. The following section provides a detailed analysis of the solution's positivity.

3.   Global positive solution

Here, we examine the global positivity criteria pertaining to the model (2.2) of the interaction between tumor and immune cells. In order to determine the positivity of the solutions of model (2.2), it is necessary to consider the Banach space $\mathcal{C} = \mathcal{C}\left(\left[-\tau, 0 \right], \mathbb{R}^3_+ \right)$, which has continuous functions that map the interval $\left[-\tau, 0 \right]$ into $\mathbb{R}_+^3$ with topological uniform convergence, where, $\mathbb{R}_+^3 = \left\lbrace \left(E_0, T_0, I_0\right)|E_0\geq0;T_0\geq0;I_0\geq0 \right\rbrace$. The following two theorems demonstrate that the stochastic model can have a positive global solution (2.2) with forms $R_1$ and $R_2$.

Theorem 3.1. For any given initial condition $(E_0, T_0, I_0)\in\mathbb{R}_+^3$, the model (2.2) with form $R_1$ for all $t\geq-\tau$, where $\tau = \max\{\tau_1, \tau_2\}$, has a solution which is almost surely unique and positive.

Proof. For any initial condition $(E_0, T_0, I_0)\in\mathbb{R}_+^3$, as the coefficients of system (2.2) satisfy the local Lipschitz condition, so system (2.2) has a unique local solution $(E(t), T(t), I(t)) on t\in[-\tau, \tau_e)$ almost surely, where $\tau_e$ stands for the explosion time. The target is to show that this solution is global i.e. $\tau_e = \infty$ with probability one. Assume $c_0\geq 1$ to be sufficiently large such that $E(\theta), T(\theta) ＆ I(\theta) (\theta\in[-\tau, 0])$ are lying in the interval $[\frac{1}{c_0}, c_0]$. For each $c\geq c_0, c\in\mathbb{N}$, define the stopping time

Assume $\inf\phi = \infty$. Therefore, $\tau_c$ is increasing as $c\rightarrow \infty$. Let $\tau_\infty = \lim_{c\rightarrow \infty}\tau_c$, then $\tau_\infty\leq\tau_e$. One needs to show that $\tau_\infty = \infty$ with probability one, then $\tau_e = \infty$ a.s. and $(E(t), T(t), I(t))\in\mathbb{R}_+^3$ for all $t\geq-\tau$. When it is incorrect, there is a pair $\iota\in(0, 1)$ and $\tilde{T} > 0$ such that $P\{\tau_\infty\leq\tilde{T}\} > \iota$. Hence, there is an integer $c_1\geq c_0$ such that

By employing the Lyapunov functional methodology, it can be observed that the system represented by Eq (2.2) exhibits the presence of a positive solution that is global in nature. Thus, we define a ${\cal C}^2$-function ${\cal G}(E, T, I):\mathbb{R} _+^3\rightarrow \mathbb{R}_+$ by

Applying It$\hat{o}$'s formula to ${\cal G}(E, T, I)$, we get

where

and ${\cal N} > 0$ is a constant which is independent of $E(t), T(t), I(t)$. Thus,

Integrating (3.4) from 0 to $\tau_c\wedge\tilde{T} = \min\{\tau_c, \tilde{T}\}$ and then taking the expectation $\mathbb{E}$ on both sides, we get

Let $\Omega_c = \{\tau_c\leq\tilde{T}\},$ for $c\geq c_1$ and in view of (3.1), we obtain $P(\Omega_c)\geq\iota$ such that, for every $\omega\in\Omega_c$, there is at least one of $E(\tau_c, \omega), T(\tau_c, \omega)$, or $I(\tau_c, \omega)$ equaling either $c$ or $\frac{1}{c}$, and then one obtains

In view of (3.5), we have

where $1_{\Omega_c}$ stands for the indicator function of $\Omega_c$. As $c\rightarrow \infty$, we have $\infty > \mathbb{E}{\cal G}(E(0), T(0), I(0))+{\cal N}\tilde{T} = \infty,$ which leads to a contradiction. Thus, it can be concluded that $\tau_\infty = \infty$ with probability one, which proves the theorem. □

Theorem 3.2. There exists a unique solution $(E(t), T(t), I(t))$ of tumor-immune model (2.2) with form $R_2$, for all $t\geq-\tau$, where $\tau = \max\{\tau_1, \tau_2\}$, and any initial condition $(E_0, T_0, I_0)\in\mathbb{R}_+^3$. The solution will remain in $\mathbb{R}_+^3$, that is $(E(t), T(t), I(t)) \in \mathbb{R}^3_+$ for all $t\geq-\tau$ almost surely.

Proof. Herein, we establish a ${\cal C}^2$-function ${\cal U}(E, T, I):\mathbb{R}_+^3\rightarrow \mathbb{R}_+$ as follows:

Applying It$\hat{o}$'s formula to ${\cal U}(E, T, I)$, one obtains

Therefore, ${\cal L}{\cal U}(E, T, I):\mathbb{R}_+^3\rightarrow \mathbb{R}_+$ is as follows:

The remainder of the proof closely resembles the proof of Theorem 3.1. □

4.   Stationary distribution and extinction

The objective of this section is to establish adequate criteria for the presence of a unique ergodic stationary distribution.

For each $t\geq0$, and probability measure $\mu$ on $(C([-\tau, 0];\mathbb{R}_+^m), M_{[-\tau, 0]})$, where $M_{[-\tau, 0]}$ is the associated Borel $\sigma$-algebra in $[-\tau, 0]$, consider the probability measure $\mu P_t$ on $(C([-\tau, 0];\mathbb{R}_+^m), M_{[-\tau, 0]})$ defined by

Definition 4.1. (Stationary Distribution ^[38]). A stationary distribution for (1.1) is a probability measure $\pi$ on $(C([-\tau, 0];\mathbb{R}_+^m), M_{[-\tau, 0]})$ such that $(\pi P_t)(\Delta) = \pi(\Delta)$ for all $t\geq0$ and $\Delta\in M_{[-\tau, 0]}$.

In order to analyze the existence of the stationary distribution of the SDDE system (2.2), it is enough to take $m = 3.$ For the stochastic model (2.2), the threshold parameter is defined as follows:

such that $\hat{\eta_1} = \eta_1+\frac{\nu_1^2}{2}$, $\hat{\mu_2} = \mu_2+\frac{\nu_3^2}{2}.$

Theorem 4.1. The tumor-immune model (2.2) admits a stationary distribution $\pi(.)$ if $\hat{{\cal T}}_0 > 1$ for any initial conditions (2.3).

Proof. According to Theorem 3.1, it has been established that there is a unique solution $(E(t), T(t), I(t)) \in \mathbb{R}_{+}^3$ on $t \geq -\tau$ for (2.2) for any given initial values (2.3). Thus,

Step 1: The diffusion matrix of stochastic delayed tumour immune model (2.2) is

Let ${\cal X}$ be any bounded domain in $\mathbb{R}_+^3$, then there exists a positive constant $d_0 = \min\{\nu_1^2E^2, \nu_2^2 T^2, \nu_3^2I^2, (E, T, I)\in \tilde{{\cal X}}\}$ such that $\sum_{i, j = 1}^3 \varsigma_{ij}(y)a_ia_j = \nu_1^2E^2a_1^2+ \nu_2^2 T^2a_2^2+\nu_3^2I^2a_3^2\geq d_0|a|^2$ for all $(E, T, I)\in\tilde{{\cal X}}, a\in\mathbb{R}_+^3$. This implies that the smallest eigenvalue of the diffusion matrix $\Pi(E, T, I)$ is bounded away from zero.

Step 2: We construct a nonnegative twice continuously differentiable function ${\cal F}: \mathbb{R}_+^3\rightarrow \mathbb{R}$ is introduced as follows:

where $c_1 = \frac{e_1e_2\alpha }{\hat{\eta_1}^2\hat{\mu_2}}$ and $c_2 = \frac{e_1e_2\alpha }{\hat{\eta_1}\hat{\mu_2}^2}$.

Select a constant ${\cal Q} > 0$ large enough such that

where $\Omega = (r+\frac{\alpha e_1e_2}{\hat{\eta_1}\hat{\mu_2}})-\frac{\nu_2^2}{2} > 0$ since ${\cal T}_0^s > 1$, and $\rho = \max\{\rho_1, \rho_2, \rho_3\}$,

By the application of It$\hat{o}$'s formula to ${\cal F}_1$, the resulting expression is obtained as:

In the same manner, we can get

From Eqs (4.6) and (4.7), one gets

Let us define a bounded closed set for any value of $\epsilon > 0$

To enhance intuitiveness, we divide $\mathbb{R}_+^3\setminus{\cal Z} = \cup_{i = 1}^{6}{\cal Z}_i$, into the following six regions:

To prove ${\cal L}\tilde{{\cal F}}\leq-1$ for any $(E, T, I)\in\mathbb{R}_+^3\setminus{\cal Z} = \cup_{i = 1}^{6}{\cal Z}_i$, we consider six cases as follows:

C.Ⅰ: For any $(E, T, I)\in{\cal Z}_1$

from Eq (4.4) and choosing $\epsilon\leq\frac{1}{{\cal Q}r\beta}.$

C.Ⅱ. For any $(E, T, I)\in{\cal Z}_2$

from Eq (4.4) and choosing $\epsilon\leq\frac{1}{{\cal Q}\xi}.$

C.Ⅲ. For any $(E, T, I)\in{\cal Z}_3$

by choosing $\epsilon\leq\frac{\eta_2}{\rho_3}.$

C.Ⅳ. For any $(E, T, I)\in{\cal Z}_4$

by choosing $\epsilon\leq\frac{\xi}{\rho_3}.$

C.Ⅴ. For any $(E, T, I)\in{\cal Z}_5$

by choosing $\epsilon\leq\frac{r\beta}{\rho_3}.$

C.Ⅵ. For any $(E, T, I)\in{\cal Z}_6$

by choosing $\epsilon\leq\sqrt[3]{\frac{(\mu_1+\mu_2)}{\rho_3}}.$

This implies that if the solution $Z(t) = (E(t), T(t), I(t))\in \mathbb{R}_+^3/{\cal Z}$ of the stochastic system (2.2), the mean times $\tau_1 ＆ \tau_2$ at which the path issuing from $Z(t)$ reaches the set ${\cal Z}$ is finite for every compact subset ${\cal X}\subset\mathbb{R}_+^3$. □

Definition 4.2. ^[21] If $\lim\limits_{t\rightarrow \infty}T(t) = 0$ a.s, the tumor $T(t)$ is said to extinct with probability one.

Theorem 4.2. Suppose that $(E(t), T(t), I(t))$ is the solution set of the tumor-immune interaction model system (2.2) along with the initial condition $(E_0, T_0, I_0) \in \mathbb{R}_+^3$. If $r < \frac{\nu_2^2}{2}$, then

In other terms, $T(t)$ approaches zero exponentially almost surely, such that the tumor will be completely eradicated from the community with unit probability.

Proof. In order to prove the theorem, we will utilize the method of direct integration on the stochastic tumor-immune model system (2.2). Firstly, we will utilize It$\hat{o}$'s formula to derive the second equation of model (2.2). This yields the following result:

On the integration of relation (4.18) over the interval $0$ to $t$ and subsequent division by $t$, the resulting expression is obtained as follows:

We can conclude that based on the strong law of large numbers of Brownian motion ^[21],

By computing the limit superior of each side of the equation and since $r < \frac{\nu_2^2}{2}$

It means that whenever $T(t)$ goes to zero exponentially with probability one, then

The proof is completed here. □

5.   Formulation of treatment as an optimal control-problem

In this section, we explore the dynamics outlined in Eq (2.2) in the context of two distinct agents: an immune boost agent labeled $\vartheta_1$ and a chemotherapeutic agent labeled $\vartheta_2$. Like Paclitaxel, the first agent is commonly regarded as cytotoxic or apoptotic. In contrast, the second agent is viewed as a form of rudimentary immunotherapy, which involves the application of an interleukin-derived drug to stimulate the immune system. We aim to gain a deeper understanding of how different agents can be used to control and manage tumor growth within the immune system, to identify the most effective combination therapy strategies, and to gain insights into the efficacy of biotherapy and immune boosts. The tumor-immune model (2.2) can be described as follows:

We assume that the initial conditions for system (5.1) satisfy

To facilitate clarity in our explanation, we shall establish a vector

such that

with the initial conditions

where, $f(\cdot, \cdot, \cdot)$ and $g(\cdot, \cdot, \cdot)$ are two vectors, each of which has components such that

and $g_1(x(t)) = \nu_1E(t)$, $g_2(x(t)) = \nu_2T(t)$, $g_3(x(t)) = \nu_3I(t)$.

Our goal is to minimize the number of tumor cells. As a result, the quadratic cost function is proposed and defined as follows:

where $K_i$, $S_i$ and $b_i$, for $i = 1, 2, 3$ are positive constants. The aim of this section is to find an optimal control $\vartheta^*(t) = [\vartheta_1^*(t), \vartheta_2^*(t)]^{'}$ that possesses the following property:

where, the control set is denoted by the set $U$ and is defined as follows:

where, $\vartheta_i^{\max} \in \mathbb{R}_+$ are constants. Before applying the stochastic maximal criterion, it is necessary to determine the Hamiltonian $H_m(p, q, \vartheta, x)$ such that

where $\langle., .\rangle$ represents the inner product in the Euclidean space; the vector fields $\langle f(x, \vartheta), p\rangle, \langle g(x), q\rangle$ and $l(x, \vartheta)$ are continuous functions on the inner product and continuously differentiable with repect to the state variable and the time variable; $p = [p_1, p_2, p_3]'$ and $q = [q_1, q_2, q_3]'$ refers to the two different sets of adjoint variables that are independent of one another. By employing a similar maximum criterion approach, we find:

where, the state $x^*(t)$ denotes the optimal trajectory followed by $x(t)$. The initial condition at time $t = 0$ for the state variables and the final condition $(t = t_f)$ for the adjoint variables $p(t)$ of Eqs (5.10) and (5.11) are as follows:

respectively. The optimal value of the control variable $\vartheta^*$ can be expressed as a function of the adjoint variables $p, q,$ and state $x^*$, as demonstrated by Eq (5.12):

where the value of $\phi$ is determined by Eq (5.12). Thus, the Hamiltonian is represented by:

Thus, we have

where

The Pontryagin maximum principle yields the adjoint system:

where, $p_1(t_f) = -b_1E(t), p_2(t_f) = -b_2T(t), p_3(t_f) = -b_3I(t).$ In addition to the condition of $p_i(t_f) = 0$ at the final time $t_f$, it holds true for values of $i = 1, 2, 3$. Similarly, the supplementary initial and final conditions are: $E^*(0) = \hat{E}, T^*(0) = \hat{T}, I^*(0) = \hat{I}$, $p(t_f) = -\frac{\partial h(x^*(t_f))}{\partial x}$ such that

Thus, the Hamiltonian function calculates the partial derivatives of $\vartheta_1$ and $\vartheta_2$, respectively; we get the following characterization of the OC

Remark 5.1. In ODEs/DDEs, time is continuous and the control inputs are typically continuous functions of time. The presence of randomness in stochastic differential equations (SDEs) makes optimal control more challenging, as the control input is required to account for the inherent uncertainty. In SDEs, optimal control policies may be stochastic, i.e., probability distributions over control inputs, rather than deterministic functions of time.

The following section pertains to the numerical simulation of outcomes derived from a stochastic optimal control problem.

6.   Numerical simulations

In this section, numerical simulations are conducted to validate the main results derived previously. The results are validated through numerical simulations utilizing Milstein's higher order method ^[39], which exhibits a strong order of convergence one. This method is employed to numerically solve the stochastic model (2.2). Subsequently, the discretization system that corresponds to the aforementioned method is obtained.

where the mutually independent $N(0, 1)$ random variables are denoted by $\chi_{1, n}$, $\chi_{2, n}$ and $\chi_{3, n}$, the integers $m_1$ and $m_2$ make the time-delays $\tau_1$ and $\tau_2$ to be represented with the step size $h$, with $\tau_1 = m_1h$ and $\tau_2 = m_2h$, respectively.

Consider the initial conditions of the system (2.2) as $E(0) = 3.8, T(0) = 1.0, I(0) = 0.5$, and the parameters are presented in Table 1. For the following numerical results, the time unit is days.

Based on the parameter values taken from Table 1, Figure 1 displays the phase and time figures of the system (2.1). Theorem 7.1 states that system (2.1) has a unique equilibrium point that is locally asymptotically stable. Figure 1(a) illustrates that all trajectories converge towards the equilibrium point $\mathcal{S}^* = \left(E^{*}, T^{*}, I^{*}\right)$. Figure 1(b) illustrates that the sizes of three distinct cell types remain constant following a period of fluctuation.

Figure 2 shows the simulations results of $(E(t), T(t), I(t))$ for stochastic model (2.2) with its corresponding deterministic model (2.1). According to Theorem 4.1 and Figure 3, the tumor size $T(t)$ in the system (2.2) is persistent. It demonstrates that the model has a unique stationary distribution and that the disease is chronic, resulting in relatively low white noise intensities where $\hat{{\cal T}}_0 > 1$, where the white noise disturbance level of system states is $\nu_i = 0.001, i = 1, 2, 3$. The solution of system (2.2) approaches equilibrium $\mathcal{S}^*$ and is stochastic asymptotically stable, as can be seen in Figure 3.

The efficacy of cancer treatments may vary from patient to patient or even within a single patient over time. When this variability is taken into account, white noise can be introduced into the model. A stochastic model refers to white noise as random fluctuations or variability. Different factors can cause it, and increasing the white noise levels can have interesting effects on the system's dynamics, such as accelerating tumor cell death. Increasing white noise in a model can be caused by external factors that introduce randomness or fluctuations into the environment.

Tumor growth can be influenced by fluctuations in nutrient availability, temperature, or oxygen levels, for example in Figure 4, we simulate the stochastic model (2.2) with different white noise levels. The results in Theorem 4.2 are validated when the white noise level is increased (see Figure 4(b)). As shown in Figure 4(a), the effector cells increase with the level of white noise. As shown in Figure 5(a), the histogram represents the population of tumor cells under the influence of white noise disturbance intensities $\nu_2 = 0.02$. On the other hand, Figures 5(b), (c) illustrate that tumor cell population decreases as noise disturbance intensities increase, as $\nu_2 = 0.1$ and $\nu_3 = 0.15$, respectively.

The graphical findings in Figure 6 demonstrate the biological importance of delay parameters. It is possible to compare tumor immune interactions with different delays based on these results. There is a significant correlation between delays and stochastic tumor immune system as shown in Figure 6. Compared to the model with time delay $\tau_1 = 0.4, \tau = 0.38$ with a small value of delay ($\tau_1 = 0.05$, $\tau_2 = 0.1$), the plot of model (2.2) with small value of delay ($\tau_1 = 0.05, \tau_2 = 0.1$) displays a significant deviation. Thus, delays significantly impact the analysis of stochastic tumor immune interactions' dynamic characteristics.

OC is an essential mechanism in mathematical modeling that targets to figure out the best appropriate control procedure for a given system, based on particular constraints. Therefore, OC can be used to outline medication protocols that minimize the growth of tumors and maximize the stimulating of the immune system, while considering the random fluctuations that arise in biological processes. To numerically solve the stochastic OC-P:

● Discretize the stochastic model with time-delays using a numerical scheme such as Milstein's higher order method.

● Carry out the forward simulations which includes the discretized stochastic model, with a particular set of control variables.

● Include computing the objective function at every step of the forward approximation.

● Execute the backward approximation by including the adjoint system in order to find the gradient of the objective function with regard to the input variables.

● Applying an optimization method, such as the stochastic approximation scheme, to minimize the objective function based on the system solutions and control constraints, and then modify the control inputs accordingly.

The numerical simulations of the optimal control problem (5.1), along with the adjoint equation (5.17) and the characterization of the optimal control (5.19), (5.20), are presented graphically. According to Table 1, the graphic outcomes of an optimal control problem were obtained by comparing them with the results of a scenario without a control system. Figure 7 shows that the controls employed are highly effective in eliminating the disease. Figure 7(a)–(f) illustrate the time evolution of effector, tumor, and IL-2 cells without and with controls $\vartheta_1$ and $\vartheta_2$. Figure 7(a)–(f) demonstrate the efficacy of the implemented controls in reducing the number of tumor cells and maximizing the number of effector cells. In the absence of chemotherapy and an immune boost, the tumor cell population exhibits an increasing trend over time, whereas the presence of treatment helps the immune system's ability to regulate the proliferation of tumor cells. Evidently, as the immune boost and chemotherapeutic agent are enhanced, there is a corresponding rise in effector cells. In contrast, tumor cells lead to their decline and eventual extinction. Meanwhile, IL-2 cells tend to return to baseline levels.

7.   Concluding remarks

Tumor modeling requires understanding the effects of environmental changes such as changes in nutrient availability, temperature, and oxygen levels on tumor-immune interactions. By incorporating discrete time-delay parameters and multiplicative white noise terms, the current work investigates the dynamical behavior of tumor-immune interactions with external treatments. The influence of environmental noises on the persistence and possible extinction of tumor cells has been studied under circumstances where the intensity of randomly varying driving forces is extremely diverse. A combination of white noises and time delays greatly affected the dynamics of the tumor-immune interaction model. It has been seen that, the presence of stochastic perturbations with a relatively small scale of white noise, tumor cells oscillate within a wide range of values, whereas large noises can lead to the eradication of the tumor cells.

To minimize tumor cells and maximize effector cells and IL-2 concentration, control variables are included in the stochastic model. An optimal control problem has been investigated to manage tumor growth and identify the most effective combination therapy strategies, and gain insights into the efficacy of biotherapy and immune boosts. Numerical results demonstrate optimality in the control variables and the effectiveness of introducing additional controls to reduce the growth of tumor cells. As a result of the presence of white noise in the optimality problem, tumor growth can be reduced.

A color-coded noise, such as telegraph noise, can be described as a random switch between two or more regimes of environmental influence on tumor immune response ^[40]. Future studies will examine how regime switching affects tumor cells, effector cells, and IL-2 cells. A sensitivity analysis can also be used to analyze tumor-free steady-state stability under small parameter variations.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work has been generously supported by UAE-ZU project #12S107.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix.   Equilibrium states and dynamics of deterministic model (2.1)

The equilibrium states of the system (2.1) are:

$(a)$ The tumor-free steady state $\mathcal{S}^{**}$, refers to a condition in which the population of tumor cells is absent, while the normal cells stay healthy:

$(b)$ The tumor-infection steady state $\mathcal{S}^*$:

where the $T^*$ represents the positive solutions of the below equation:

In order to analyze the stability of the equilibrium states, it is common practice to linearize the system around the equilibrium states and determine the Jacobian matrices. Currently, our focus is only on examining the stability of the equilibrium point in the presence of a tumor.

The Jacobian matrix $\Delta_{\mathcal{S}^*}(E^*, T^*, I^*)$ is obtained in order to determine the eigenvalues of the system (2.1).

Remark 7.1. The stability analysis of the equilibrium states is examined by the eigenvalues of $\Delta_{\mathcal{S}^*}$. If the real parts of all the eigenvalues of the matrix $\Delta_{\mathcal{S}^*}$ are negative, then the equilibrium state $\mathcal{S}^*$ exhibits local asymptotic stability. Conversely, if at least one of the eigenvalues exhibit a positive real part, the equilibrium state is considered to be unstable.

Following some calculations, the characteristic equation of (7.1) is obtained:

where,

Case (ⅰ): $\tau_1 = 0, \tau_2 = 0,$ we have

Based on the Routh-Hurwitz criterion, all the eigenvalues of the characteristic equation Det$(\Delta_{\mathcal{S}^*}) = 0$ have negative real parts if the following conditions hold:

Therefore, it can be concluded that the asymptotic stability of the system (2.1) is achieved in the equilibrium state $\mathcal{S}^*$, when $\tau_1 = \tau_2 = 0$.

Case (ⅱ): $\tau_1 = \tau_2 > 0,$ we have

Assume that the pure imaginary root $\lambda = i\theta, \theta > 0$. Then applying the $\lambda = i\theta$ into (7.5) and separating its real and imaginary parts, we have

Through mathematical calculations, the following equations are derived:

Using the trigonometric formula $\cos^2(\theta\tau)+\sin^2(\theta\tau) = 1$, we have

where,

Therefor, the equilibrium states of system (2.1) are asymptotically stable if and only if (7.9) has no positive real roots. Assume that $\omega = \theta^2,$ we have

If $h_p < 0, p = 1, \ldots, 5$, then (7.10) has at least one positive root and if $h_p > 0, p = 1, \ldots, 5$, then (7.10) has all the roots must be negative.

Case (ⅲ): $\tau_1 = 0, \tau_2 > 0$. Then, (7.2) can be written as

Let suppose $\lambda = {\bf i}\omega$ is a pure imaginary root of (7.11), then

It can be deduced from Eq (7.12) that,

where,

Denote $\beta = \omega^2$, (7.13) becomes

Then let

If the condition $(\mathcal{C}_1) z_3 < 0$ is satisfied, given that $\frac{d\mathcal{F}(\omega)}{d\omega} > 0$ $\forall \omega > 0$, it can be concluded that Eq (7.13) possesses at least one positive real root. Therefore, Eq (7.11) has at least one pair of purely imaginary roots. According to the findings of Sun et al. ^[42], the subsequent outcome is as follows:

Lemma 7.1. The following holds true for (7.11):

(a) If condition (7.4) is satisfied and the values of $z_1, z_2, z_3 > 0$, then equation (7.11) does not has any roots with zero real parts for values of $\tau_2\geq0$.

(b) If $z_3 < 0$ and $z_1, z_2 > 0$ holds, then (7.11) has a pair of purely imaginary roots $\pm {\bf i} \omega_0$ when $\tau_2 = \tau_{2j}$ and for any $\omega_0$ (unique positive zero of the function $\mathcal{G}(\omega)$)

The proof of Lemma 7.1 has resemblance to the proof of Lemma 3.1 as presented in the work of Li et al. ^[43]. Here, it is left out.

$(\mathcal{C}_2)$ $\tilde{\mathcal{A}}_1\tilde{\mathcal{A}}_2+\tilde{\mathcal{B}}_1\tilde{\mathcal{B}}_2 > 0$

where,

Case (ⅳ): When $\tau_1 > 0, \tau_2 = 0$, the Eq (7.2) becomes

For pure imaginary root $\lambda = {\bf i}\omega$, Eq (7.17) can be written as

It can be deduced from Eq (7.18) that,

where,

Take $\alpha = \omega^2$, then (7.19) becomes

and let

If the condition $(\mathcal{C}_3) r_5 < 0$ is satisfied, given that $\frac{d\mathcal{F}(\omega)}{d\omega} > 0$ $\forall \omega > 0$, it can be concluded that Eq (7.20) possesses at least one positive real root. Therefore, Eq (7.17) has at least one pair of purely imaginary roots. According to the findings of Sun et al. ^[42], the subsequent outcome is as follows:

Lemma 7.2. The following holds true for (7.17):

(a) If condition (7.4) is satisfied and the values of $r_i > 0, i = 1, 2, 3, 4, 5$, then Eq (7.17) does not have any roots with zero real parts for values of $\tau_1\geq0$.

(b) If $r_5 < 0$ and $r_i > 0, i = 1, 2, 3, 4$ holds, then (7.17) has a pair of purely imaginary roots $\pm {\bf i} \omega_0$ when $\tau_1 = \tau_{1j}$ and for any $\omega_0$ (unique positive zero of the function $\mathcal{F}(\omega)$)

The proof of Lemma 7.2 has resemblance to the proof of Lemma 3.1 as presented in the work of Li et al. ^[43]. Here, it is left out.

$(\mathcal{C}_4)$ $\hat{\mathcal{A}}_1\hat{\mathcal{A}}_2+\hat{\mathcal{B}}_1\hat{\mathcal{B}}_2 > 0$

where,

Case (ⅴ): $\tau_1 > 0, \tau_2 > 0.$ Let the delay $\tau_2 \in [0, \tau_{20})$ in (7.2) and choose $\tau_1$ as a bifurcating parameter. Thus let pure imaginary root (7.2) as $\lambda = {\bf i}\omega = \omega\left(\cos\frac{\pi}{2}+\sin\frac{\pi}{2}\right).$ Then (7.2) becomes

By (7.23), we have

In view of $\cos^2\omega\tau_2+\sin^2\omega\tau_2 = 1,$ we have

where

Then let the polynomial function

When the condition $(\mathcal{C}_5)\ q_6 < 0$ holds, given that $\frac{d\mathcal{H}(\omega)}{d\omega} > 0$ $\forall \omega > 0$, it can be concluded that Eq (7.26) possesses at least one positive real root. Therefore, Eq (7.2) has at least one pair of purely imaginary roots. According to the findings of Sun et al. ^[42], the subsequent outcome is as follows:

Lemma 7.3. The following holds for (7.2):

(ⅰ) If $q_i > 0, i = 1, 2, \ldots, 6$, then Eq (7.2) has no roots with zero real parts.

(ⅱ) If the value of $q_6$ is less than zero, it can be concluded that Eq (7.2) possesses at least one pair of purely imaginary roots.

(ⅲ) When given that $q_6 > 0$ and there exists a constant $\delta > 0$ such that $\mathcal{H}^{'}(\delta) < 0$, it can be concluded that Eq (7.2) possesses at least two pairs of roots that are totally imaginary.

The proof of Lemma 7.3 has resemblance to the proof of Lemma 3.1 as presented in the work of Li et al. ^[44]. Here, it is left out. Now, it is assumed that Eq (7.26) possesses a positive real root $\omega$. By (7.24), we have

where, $j = 0, 1, 2, 3, \ldots.$ Then $\pm{\bf i}\omega$ is a pair of roots of (7.2) when $\tau_2 = \tau_{2j}^{\star}.$

$(\mathcal{C}_6)$ $\bar{\mathcal{A}}_1\bar{\mathcal{A}}_2+\bar{\mathcal{B}}_1\bar{\mathcal{B}}_2 > 0$

where,

Now, we reach the following results:

Theorem 7.1. If system (2.1) has an tumor-infection equilibrium state $\mathcal{S}^*$, then the following statements are hold:

(ⅰ) For $\tau_1 = \tau_2 = 0,$ the equilibrium state $\mathcal{S}^*$ of system (2.1) is asymptotically stable if the conditions in (7.4) are true.

(ⅱ) For $\tau_1 = \tau_2\geq0$, the tumor-infection equilibrium state $\mathcal{S}^*$ is asymptotically stable if all the roots of Det$(\Delta_{\mathcal{S}^*}) = 0$ are negative real parts.

(ⅲ) If $\tau_1 = 0$ and $\tau_2 > 0$ and (7.4), $(\mathcal{C}_1), (\mathcal{C}_2)$ hold, then the tumor-infection equilibrium state $\mathcal{S}^*$ is asymptotically stable when $\tau_2 \in [0, \tau_{20})$.

(ⅳ) If $\tau_2 = 0$ and $\tau_1 > 0$ and (7.4), $(\mathcal{C}_3), (\mathcal{C}_4)$ hold, then the tumor-infection equilibrium state $\mathcal{S}^*$ is asymptotically stable when $\tau_1 \in [0, \tau_{10})$.

(ⅴ) If $\tau_1 \in [0, \tau_{10}),$ and (7.4), $(\mathcal{C}_5)$ and $(\mathcal{C}_6)$ are satisfied, then the tumor-infection equilibrium state $\mathcal{S}^*$ is asymptotically stable when $\tau_2 \in [0, \tau_{20}^{\star})$.

Remark 7.2. The stability of the model (2.1) is completely determined by the zeros of its characteristic polynomial (7.2). To find the number of zeros of (7.2) as time delays $\tau_1$ and $\tau_2$ vary, we write Det$(\Delta_{\mathcal{S}^*})$ in (7.2) as $\Delta_{\mathcal{S}^*}(\lambda, \tau_1, \tau_2)$. In this case the solution of the system is stable, the coefficient polynomials $\Lambda_1(\lambda) = \lambda^3+w_1\lambda^2+w_2\lambda+w_3, \Lambda_2(\lambda) = \lambda^2w_4+\lambda+w_5+w_6, \Lambda_3(\lambda) = w_7$ satisfies the following conditions^[41]: (1) $\operatorname{deg}\left(\Lambda_1(\lambda)\right) \geqslant \max \left\{\operatorname{deg}\left(\Lambda_2(\lambda)\right), \operatorname{deg}\left(\Lambda_3(\lambda)\right)\right\}$, (2) $\Lambda_1(0)+\Lambda_2(0)+\Lambda_3(0) \neq 0$, (3) The polynomials $\Lambda_1(\lambda), \Lambda_2(\lambda)$ and $\Lambda_3(\lambda)$ do not have any common zeros, (4) $\lim _{\lambda \rightarrow \infty}\left(\left|\Lambda_2(\lambda) / \Lambda_1(\lambda)\right|+\left|\Lambda_3(\lambda) / \Lambda_1(s)\right|\right) < 1$. In the absence of condition (1), the polynomial Det($\Delta_{\mathcal{S}^*}$) is not stable for any positive delays. If condition (2) is not satisfied, then the polynomial Det($\Delta_{\mathcal{S}^*}$) has 0 as a root for any $\tau_1$ and $\tau_2$, thus, it can never be stable. Condition (3) is natural. (4) is automatically satisfied since its left-hand side is zero. In condition (5) the number of zeros of $Det (\Delta_{\mathcal{S}^*})$ can change only if a zero appears on the imaginary axis.