Longtime evolution and stationary response of a stochastic tumor-immune system with resting T cells

1.   Introduction

A tumor is a new organism formed by the proliferation of local tissue in the body under the action of various oncogene factors. In recent years, an increasing interest has been focused on the role of the immune system in fighting a tumor growth ^{[1,2,3,4,5,6]}. T cells in the immune system are one of the important components in resisting tumor. There are two population of T cells; one is the helper T cells expressed by CD4+T proteins, they cannot kill tumor cells directly, but once they are stimulated by the macrophage or their cognate antigen, they can produce a kind of cytokine such as interleukin -2 (IL-2) and send biochemical signals to a special type of effector cell called cytotoxic T lymphocyte CD8+ T cells ^[7]. Cytotoxic T lymphocyte cells are the other important population in the immune response, they can eliminate or kill the infected cells by mounting a cytotoxic reaction after they are activated by the helper T cells ^[8,9]. Therefore, the roles of the helper T cells and the cytotoxic T cells are distinct, they perform the complementary functions to eliminate the tumor.

It is a key task to identify the components of the cells to understand how the immune system works to overcome the tumor. However, due to the complexity of cell dynamics, most work concentrates on the discussion about the interaction between the tumor cells and the cytotoxic T lymphocytes, but neglects the role played by the helper T cells in tumor elimination. In recent years, some researchers have proposed that the helper T cells are a vital and essential component in tumor-immune response ^{[10,11,12,13]}, and built a system of ODEs to describe the population in a system consist of the helper T cells, cytotoxic T cells and the tumor cells. The representative work can be found in Refs. ^{[7,14,15,16,17]}.

Furthermore, according to the biochemical features of the tumor cells and the immune cells, some models use a prey-predator interaction model. Immune cells are the predators, who attack and destroy the tumor cells. That is, the helper T cells in a resting state interact with antigens to release cytokines and stimulate the cytotoxic T cells. As follows, the cytotoxic T cells in a hunting state attack or destroy the tumor cells. While tumor cells are the prey, which are eliminated or killed by the immune cells. Predator hunting cells and predator resting cells provide a useful model to learn the dynamical behaviors such as equilibration, stable solutions, longtime evolution and so on. In fact, a system with a component of the resting T cells provide insight to learn the tumor growth. For example, Sarkar ^[18] in 2005 considered a prey-predator like system with hunting and resting cells, and investigated the stability of the positive equilibrium from the analytical and numerical perspective, they obtained a threshold value for the rate of predation of the tumor cells by the hunting cells and discussed the control of malignant tumor growth in deterministic case. In 2014, Kaur ^[14] analyzed tumor growth progression and regression in a prey-predator system, their analysis indicated that the tumor growth will persist from a recurring state to a dormant state with the increase of the resting cells. In 2018, Dritschel ^[7] developed a mathematical model to examine the role of the helper T cells in an anti-tumor immune response, their work revealed that the immunoediting depends highly on the infiltration rates of the helper and cytotoxic T cells. In 2023, Dehingia ^[16] proposed a modified system consisting of tumor cells, resting helper T cells and hunting cytotoxic T cells, he discussed the stability and Hopf bifurcation of this system.

The works mentioned above regarding the prey-predator like systems considering two kinds of T cells and the tumor cells are mainly investigated in the deterministic condition. However, as is widely recognized, a variety of random factors in the cell environment should be taken into account. That is, the biochemical proliferation of the cells may undergo stochastic variations influenced by factors such as nutrient availability, temperature, radiation, oxygen supply, and gene expression ^{[19,20,21,22,23,24,25]}. Therefore, an idealized tumor-immune system should be modeled stochastically.

In this paper, the novelty of our efforts will focus on the special tumor-immune system including hunting and resting T cells, at the same time, stochastic fluctuations on the proliferation parameters will be explored, in order to understand the impact of the noise on the cells dynamics, as well as give insight on how to eliminate tumors and how to improve cancer treatment.

2.   Stochastic prey-predator like tumor-immune system

Based on the deterministic prey-predator like tumor-immune model in ^[14], we study a modified model to take the impact of random perturbations on the cell evolution into account, and use standard Brownian motions to characterize it in the tumor-immune model. The cells'dynamics is possible to be changed in the case of stochastic perturbation, so it is necessary to have a relevant research from the point of random variables. On the other hand, there are plenty of valuable results about this model in the deterministic case, we expect to give some new insights by studying a stochastic one. The model is given by a set of differential equations in a form of the follow

where $T(t), H(t)$ and $R(t)$ are the population of tumor cells, predator hunting T cells and predator resting T cells in the given physiologic space, respectively. $\alpha_{1}$ and $\alpha_{2}$ represent the inherent growth rates of tumor cells and resting T cells. $K_{1}, K_{2}$ indicate the maximum carrying capacity of the environment. $\beta_{1}$ and $\beta_{2}$ represent the rate at which hunting T cells kill tumor cells and the rate at which tumor cells deplete hunting T cells. $\delta_{1}$ and $\delta_{2}$ indicate the apoptosis rates of hunting T cells and resting T cells. $\gamma$ represents the activation rate of resting T cells to hunting T cells, $k$ is the proliferation rate of resting T cells, $\eta$ is the half-saturation coefficient of proliferation term. $B_{1}(t), B_{2}(t), B_{3}(t)$ are mutually independent standard Brownian motions, and $\xi_{1}, \xi_{2}, \xi_{3}$ represent noise intensity of them.

3.   Dynamical evolution and properties

3.1.   Existence and uniqueness of the positive solution

Since $T(t), H(t), R(t)$ represent the quantities of tumor cells, hunting T cells, and resting T cells respectively, it is necessary to verify that the solutions of the system (2.1) are non-negative. Therefore, in this section, we begin with the following theorem on the global existence and non-negativity of solutions to the system (2.1).

Theorem 3.1. For any given initial values $(T(0), H(0), R(0))\in R_{+}^{3}$, there exists a global solution $(T(t), H(t), R(t))$ of system (2.1) such that

Proof. Due to the coefficients of the system (2.1) satisfying the local Lipschitz condition, for a given initial value $(T(0), H(0), R(0))\in R_{+}^{3}$, the system has an unique local positive solution $(T(t), H(t), R(t))$ for $t\in [0, \tau_e)$ a.s., where $\tau_e$ represents the explosion time. In order to confirm the local solutions to be global ones, it is necessary to prove that $\tau_e = \infty$.

Let $m_{0} \ge 1$ be sufficiently large such that $T(0)\in [\frac{1}{m_{0} }, m_{0}], H(0)\in [\frac{1}{m_{0} }, m_{0}], R(0)\in [\frac{1}{m_{0} }, m_{0}]$. For $m > m_{0}$, we define

where we set $\inf\emptyset = \infty$. By the definition of $\tau _{m}$ we observe that $\tau_m$ increases with $m$ approaching to the infinity and $\tau_m < \tau_e$, therefore

Now, we prove that $\tau_\infty = \infty$ is true almost surely first. Moreover, a proof by contradiction is applied, that is, assume the condition does not hold. If so, There are a constant $T_{1}\ge 0$ and another constant $\varepsilon \in(0, 1)$ such that

Thus, there exists an integer $m_{1}\ge m_{0}$ such that for all $m\ge m_{1}$, we have

Let us define a Lyapunov function $V: Int(R_{+}^3)\rightarrow R_{+}$:

Applying It$\hat{o}$'s formula, we can differentiate the function $V$

where $\mathcal{L}[\cdot]$ is the drift term of the $[\cdot]$ derivative. Denote

then we can rewrite the above expression as:

Integrating on both sides from $0$ to $\tau_{m}\wedge T_{1}$ of the inequality (3.4) and then taking the expectation for every terms, we have:

where $I_{A}(\cdot)$ represents the indicator function of set $A$, $\tau_{m}\wedge T_{1} = \min(\tau_{m}, T_1)$. Utilizing the Gronwall inequality ^[26], we can establish the following result:

Let $\Omega_{m} = \{\omega\in\Omega{:}\tau_{m} = \tau_{m}(\omega)\leq T_1\}{:}(m\geq m_{1})$, based on the Eq (3.3), we have:

Observe that for every $\omega$ in the set $\omega\in \Omega_m$, it holds that at least one of $T(\tau_m, \omega), H(\tau_m, \omega), R(\tau_m, \omega)$ is equal to $m$ or $\dfrac{1}{m}$. Therefore,

Combining Eqs (3.5) and (3.6), we obtain:

On the other hand,

Let $m\to \infty$, we have $\infty \le Ke^{2bT} < \infty$, which implies that $\tau_\infty$ must be equal to infinity. Considering $\tau_\infty\le \tau_e$, we conclude that $\tau_e = \infty$ is indeed equal to infinity. □

3.2.   Moment boundedness

Based on the existence of the positive solution for every cell population, this section focuses on the asymptotic estimation of the moments of $T(t), H(t)$, and $R(t)$.

Definition 1. $p$-order statistical moments of the solutions to the system (2.1) are bounded. there exists a function $r = r(x_{0}, t_{0})$ such that when $\left | x_{0} \right |\le m$, if for any $m > 0$, it holds that

If the function $r(x_{0}, t_{0})$ mentioned here is independent of $t_{0}$, then the solution to the system (2.1) is said to have uniformly bounded $p$-order moments.

Theorem 3.2. For any $p > 1$, we have

For any $p\in (0, 1+2\delta_{1}/\xi_{2}^{2})$, we can find a positive constant $L_{1}$ such that

Proof. we are going to prove the moment boundedness of $T(t)$ first. Consider the following auxiliary procedure $\varphi (t)$

where $B_1(t)$ is the standard Brownian motions defined in system (2.1). Using the comparison theorem ^[27], we can obtain that $0\le T(t)\le \varphi(t)$ $(t\ge 0)$. By applying It$\hat{o}$'s formula, the drift coefficient in the differential formula of $\varphi ^{p}(t)$ can be expressed as

Taking expectations on both sides of Eq (3.11) and applying the Holder's inequality ^[26], we have

Therefore,

Thus, for any $p > 1$, we have

Next, we will demonstrate the moment boundedness of $R(t)$. Similarly, consider the following auxiliary process $\psi (t)$

where $B_3(t)$ is the standard Brownian motions defined in system (2.1). Then, we get $0\le R(t) \le \psi(t)$ for $t\ge 0$. Notice that $\psi(t)$ follows the similar structure as $\varphi(t)$, so we can obtain:

where $A_{1} = \alpha_{2}+k-\delta_{2}$. Therefore, for any $p > 1$, we have

Finally, we will establish the moment boundedness of $H(t)$. Consider the following two auxiliary processes, $y(t)$ and $z(t)$

where $B_{2}(t), B_{3}(t)$ are mutually independent standard Brownian motions defined in system (2.1). Consider the function $V_{1}(t)$ defined as $V_{1}(t) = (y(t)+z(t))^{q}$. According to the It$\hat{o}$'s formula, we get

The condition $q\in(0, 1+2\delta_{1}/\xi_2^{2})$ implies $\delta_{1}-\dfrac{q-1}{2} \xi_{2}^{2} > 0$, then we can find a positive constant $\kappa$ sufficiently small such that

Define $U(t) = e^{\kappa t}(y(t)+z(t))^{q}$, we get

where

As $\delta_{1}-\dfrac{q-1}{2}\xi_{2}^{2}-\dfrac{\kappa}{q} > 0$ and $\dfrac{\alpha_{2}}{K_{2}} > 0$, it implies that

Combining the continuity of $q(y+z)^{(q-2)} W(y, z)$, we can find a positive constant $L$ satisfies:

Therefore,

By integrating both sides of the inequality (3.17) from $0$ to $t$ and taking the expectation, we have

By applying the comparison theorem ^[27], we have $H(t)\le y(t)$ and $R(t)\le z(t)$. Let $L_{1} = \dfrac{L}{\kappa}$ such that

□

4.   Extinction and persistence

4.1.   Extinction

In this section, we will analyze the extinction of the resting T cells $R(t)$ and tumor cells $T(t)$ in system (2.1), and provide numerical examples for verification.

Theorem 4.1. When $\delta+\dfrac{\xi^{2}}{4}-\alpha_{2}-k > 0$, the resting T cells will eventually become extinct, i.e.,

Here, $\delta = \min\left \{\delta_{1}, \delta_{2}\right \}, \xi^{2} = \min\left \{ \xi_{2}^{2}, \xi_{3}^{2} \right \}$.

Proof. We denote $\Psi(t) = H(t)+R(t)$. Taking the logarithm of $\Psi(t)$ and applying It$\hat{o}$'s formula yields

By integrating both sides of the inequality (4.1) from $0$ to $t$ and dividing by $t$, we obtain

Combining the strong law of large numbers ^[26], we have

Therefore, if $\delta+\frac{\xi^{2}}{4}-\alpha_{2}-k > 0$, we can derive

As a result, we conclude that $\lim_{t \to \infty} \Psi (t) = 0 \ a.s.$. Combining with the positivity of $R(t)$ and $H(t)$, we have $\lim_{t \to \infty} R(t) = 0 \ a.s.$. □

Theorem 4.2. If $\alpha_{1}-\frac{1}{2}\xi_{1}^{2} < 0$, tumor cells will eventually go extinct, such that

Proof. By combining the It$\hat{o}$' formula, we have

Integrating the inequality (4.3) from $0$ to $t$ and dividing by $t$, results in

As $t$ approaches $\infty$, by combining the strong law of large numbers ^[26], we obtain that

Therefore, if $\alpha_{1}-\frac{1}{2}\xi_{1}^{2} < 0$, we can derive

This implies that

In other words, as t approaches $\infty$, the tumor cell population $T(t)$ tends to $0$.

□

Next, we will use the Milstein method ^[28] to illustrate our main theoretical results. The discretized equations for the system (2.1) are given as follows:

Here, $\Delta t > 0$ represents the time increment, $\xi_{i}(i = 1, 2, 3)$ denotes the noise intensity, and $W_{i, k}(i = 1, 2, 3.\ k = 1, 2, ...)$ follows a standard normal distribution. The parameter values are shown in Table 1.

The initial values for tumor cells, hunting T cells, and resting T cells are given as $T(0) = 5\times10^{9} (cells), H(0) = 4\times10^{6} (cells), R(0) = 3\times10^{8} (cells)$, respectively. First, verify the extinction of resting T cells $R(t)$ with $\xi_1 = 2$, $\xi_2 = 2$, $\xi_3 = 1$. Calculate $\delta+\frac{\xi^{2}}{4}-\alpha_{2}-k = 0.103 > 0$. According to Theorem 4.1, resting T cells will eventually go extinct, which is consistent with the results shown in Figure 1. When $\xi_1 = 0.8$, $\xi_2 = 2$, $\xi_3 = 2$, the calculation yields $\alpha_{1}-\frac{1}{2}\xi_{1}^{2} = -0.14 < 0$. According to Theorem 4.2, tumor cells $T(t)$ will go extinct. As shown in Figure 2, tumor cells eventually converge to zero.

Next, let's discuss the effect of stochastic perturbation on tumor cell proliferation by fixing $\xi_2 = 2, \xi_3 = 2$, and varying the noise intensity $\xi_1$ while ensuring the extinction condition for $T(t)$. We take $\xi_1 = 0.8, 1.5, 2$ as examples, and the result are shown in Figure 3. From Figure 3, it can be observed that increasing the noise intensity $\xi_1$ accelerates the extinction of tumor cells. This implies that stochastic environmental perturbations suppress the growth of tumor cells.

4.2.   Persistence

Lemma 1 (^[31]). Assuming $Z\in \left [\Omega \times \left [0, +\infty \right), R^+ \right]$, if there exist positive constants $\eta, t_{0}$ and $\eta_{0}$ such that

then we have

Lemma 2 (^[32]). Let $T(t), H(t), R(t)$ be the solution of stochastic system (2.1) initial conditions $(T(0), H(0), R(0))$

1) if $\alpha_{1}-\dfrac{1}{2}\xi_{1}^2 > 0$, then

2) if $\delta_2+\frac{1}{2} \xi_{3}^{2}-\alpha_{2}-k > 0$, then

3) if $\alpha_{2}-\delta_{2}+k-\frac{1}{2} \xi_{3}^{2} > 0$, then

Proof. Construct the following auxiliary procedures:

where $B_{1}(t), B_{2}(t), B_{3}(t)$ are mutually independent standard Brownian motions defined in system (2.1). The comparison theorem ^[27] leads to $0\leq T(t)\leq T_1(t), 0\leq H(t)\leq H_1(t), 0 \leq R(t)\leq R_1(t)$. Using the It$\hat{o}$'s formula to $\log T_{1}(t), \log H_{1}(t)$, and $\log R_{1}(t)$, then integrating both sides from $0$ to $t$, we get

For convenience to write, we define the notation $\left \langle x(t) \right \rangle = \dfrac{1}{t}\int_{0}^{t}x(s)ds$. Then, we can obtain

According to the strong law of large numbers ^[26], we get

For sufficiently small $\varepsilon > 0$, when $\alpha_{1}-\frac{1}{2}\xi_{1}^2 > 0$, according to Lemma 1, it can be deduced that

which implies that

substituting Eq (4.15) into Eq (4.13) and combining with the Eq (4.14), it yields

Similarly, when $\delta+\frac{1}{4} \xi^{2}-\alpha_{2}-k > 0$, Lemma 1 allows us to deduce

substituting Eq (4.17) into Eq (4.13) and combining with Eq (4.14), it yields

Using the same method as in Theorem 4.2, the condition $\delta_2+\frac{1}{2} \xi_{3}^{2}-\alpha_{2}-k > 0$ leads to $\lim_{t \to \infty} R_{1}(t) = 0$, then we get

substituting Eq (4.19) into Eq (4.13) and combining with Eq (4.14), it yields

According to formula Eqs (4.16), (4.18), (4.20) and using the comparison theorem ^[26], we get

□

Theorem 4.3. When $\alpha_{1}-\dfrac{1}{2}\xi_{1}^{2} > 0$ and $\delta-\dfrac{\xi^{2}}{4}-\alpha_{2}-k > 0$, the tumor cell population $T(t)$ is weakly persistently bounded, with $\lim \sup_{t \to \infty}\left \langle T(t) \right \rangle > 0$. Moreover, it can be deduced that $\sup_{t\to\infty}\langle T(t)\rangle\leq\frac{K_1\left (\alpha_1-\frac{1}{2}\xi_1^2 \right) }{\alpha_1}$.

Proof. According to the It$\hat{o}$'s differential rule, we have

By integrating both sides from $0$ to $t$ and dividing by $t$, and taking the upper limit, we obtain

Applying Lemma 2, we have

By using the strong law of large numbers ^[26], we get $\lim_{t\to\infty}\dfrac{1}{t}\int_0^t\xi_1dB_1(t) = 0$ a.s., and when $\delta-\dfrac{\xi^{2}}{4}-\alpha_{2}-k > 0$, we have $\lim_{t \to \infty} H(t) = 0$. Thus, $\lim\sup_{t\to\infty}\langle H(t)\rangle = 0$ a.s.. Substituting it into Eq (4.23), we obtain

Therefore, $\lim\sup_{t\to\infty}\langle T(t)\rangle > 0$, which further indicates that the tumor cell population $T$ is weakly persistent. In addition,

Applying Lemma 1, when $\alpha_{1}-\frac{1}{2} \xi_{1}^{2} > 0$, we have

□

We verify the weak persistence of the tumor cell population below. With noise intensities set to $\xi_{1} = 0.3, \xi_{2} = 2, \xi_{3} = 2$, we calculate that

In accordance with Theorem 4.3, we conclude that $T(t)$ is weakly persistently average. As illustrated in Figure 4, the tumor cell population persists at a low level.

In order to further investigate the dynamical behavior of tumor cells in the persistent state, we control $\xi_{1}$ by selecting $\xi_{1} = 0.1, 0.3, 0.5$ respectively to indicate the persistent state of tumor cells under different environmental noise.

From Figure 5, it can be observed that as the intensity of noise increases, the fluctuation in the number of tumor cells becomes more pronounced, and the center of oscillation gradually shifts downwards. When $\xi_{1} = 0.1$, the number of tumor cells fluctuates around $4.8\times10^{8}$, mainly within the range of $[3\times10^{8}, 7\times10^{8}]$. For $\xi_{1} = 0.3$, the number of tumor cells is concentrated within the range of $[1\times10^{8}, 1\times10^{9}]$ with the increase of fluctuation range and the decrease in the probability peak value. When $\xi_{1}$ is further increases to 0.5, the number of tumor cells start to fluctuate around zero, which indicating weaker persistence.

Figure 6(a)–(c) depicts the distribution of persistent states of tumor cells under different noise intensity. We combine the populations of hunting T cells and resting T cells to obtain the two-dimensional stationary probability distribution of tumor cells and two types of T cells. It can be observed that the number of two types of T cells is primarily concentrated at position 0. In order to obtain the probability distribution function of tumor cells, we simplify system (2.1) by assuming $H(t) = R(t) = 0$. Thus, system (2.1) degenerates into system (4.27).

By solving the Fokker-Planck equation of Eq (4.27), we obtain an analytical expression for the stationary probability density.

Here, $N$ is the normalization coefficient. Then we validate the accuracy of the analytical expression for the probability density by utilizing the Milstein simulation method ^[28] to sample. By taking $\xi_{1} = 0.1, 0.3, 0.5$ and using a time step of $\Delta t = 0.01$, we sample $N = 1\times10^{7}$ sample points, respectively, and obtain the quantities of tumor cells under different noise intensity.

The histogram in Figure 7(a)–(c) represents the distribution of the number of sample points under different noise intensity, and the curve fitted on the histogram edges represents the probability distribution function based on Eq (4.28). It can be observed from Figure 7 that the function fits the distribution of sample points perfectly. By approximating the frequencies of tumor cells at different quantities in Figure 7 as probabilities, an accurate expression for $p(x)$ is obtained by substituting it into Eq (4.28). Figure 8 illustrates the steady-state probability distributions of tumor cells under different noise intensities. It can be observed that when the noise intensity is low, the probability density distribution exhibits a unimodal shape, indicating that the tumor maintains stable growth within a certain quantity range. As the noise intensity increases, the tumor cell population moves towards lower levels, and the peak value of the probability decreases, indicating that the increased noise intensity effectively suppresses the expansion of the tumor population. When $\xi_{1}$ continues to increase beyond a certain threshold, the probability density function transitions from a unimodal shape to a decreasing shape, and the tumor population concentrates around zero, indicating that tumor growth is significantly restricted, with lower invasiveness and metastatic potential.

5.   Conclusions

We investigate the dynamic behavior and stationary response of tumor cells in a predator-prey like system with resting T cells. First, we prove the existence and uniqueness of global positive solutions for the stochastic system which guarantees the biological significance of the system model. Second, we establish the boundedness of moments for $T(t), H(t)$ and $R(t)$ by constructing appropriate auxiliary equations. Finally, we provide sufficient conditions for the extinction of resting T cells, as well as the threshold for persistence and extinction of tumor cells and validate the results through numerical simulations. When $\xi_{1}$ is small, the tumor remains in a persistent state. From Figure 6, it can be seen that in the persistent state, the quantities of two types of T cells tend to cluster around zero. By simplifying the system (2.1) through the assumption $H(t)+R(t) = 0$, we obtain the system (4.27) and derive an analytical expression for the stationary probability density of the tumor. The stochastic perturbations in the environment play a crucial role in eliminating tumor cells, increasing the noise intensity in the persistent state leads to stochastic bifurcation, and the stationary probability density of tumors transitions from unimodal to decreasing. When the noise intensity reaches a certain threshold, tumor cells eventually transition from a persistent state to an extinct state, and further increasing the noise intensity accelerates tumor extinction.

In fact, due to the influence of factors such as radiation and viruses, parameters in the tumor immune system may undergo mutations. Therefore, it is crucial to investigate the parameter changes and corresponding responses of the tumor immune system in different environments. The further efforts can be made to study the dynamic behavior of tumor immune system under stochastic switching.