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

: In this paper, we take the resting T cells into account and interpret the progression and regression of tumors by a predator-prey like tumor-immune system. First, we construct an appropriate Lyapunov function to prove the existence and uniqueness of the global positive solution to the system. Then, by utilizing the stochastic comparison theorem, we prove the moment boundedness of tumor cells and two types of T cells. Furthermore, we analyze the impact of stochastic perturbations on the extinction and persistence of tumor cells and obtain the stationary probability density of the tumor cells in the persistent state. The results indicate that when the noise intensity of tumor perturbation is low, tumor cells remain in a persistent state. As this intensity gradually increases, the population of tumors moves towards a lower level, and the stochastic bifurcation phenomena occurs. When it reaches a certain threshold, instead the number of tumor cells eventually enter into an extinct state, and further increasing of the noise intensity will accelerate this process.


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 antitumor 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.

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

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)) ∈ 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)) ∈ R 3 + , the system has an unique local positive solution (T (t), H(t), R(t)) for t ∈ [0, τ e ) a.s., where τ e represents the explosion time.In order to confirm the local solutions to be global ones, it is necessary to prove that τ e = ∞.
Let m 0 ≥ 1 be sufficiently large such that where we set inf ∅ = ∞.By the definition of τ m we observe that τ m increases with m approaching to the infinity and τ m < τ e , therefore lim Now, we prove that τ ∞ = ∞ 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 ≥ 0 and another constant ε ∈ (0, 1) such that Thus, there exists an integer m 1 ≥ m 0 such that for all m ≥ m 1 , we have Let us define a Lyapunov function V : Applying It ô's formula, we can differentiate the function V where then we can rewrite the above expression as: Integrating on both sides from 0 to τ m ∧ T 1 of the inequality (3.4) and then taking the expectation for every terms, we have: where . Utilizing the Gronwall inequality [26], we can establish the following result: Let , based on the Eq (3. 3), we have: Observe that for every ω in the set ω ∈ Ω m , it holds that at least one of Combining Eqs (3.5) and (3.6), we obtain: On the other hand, Let m → ∞, we have ∞ ≤ Ke 2bT < ∞, which implies that τ ∞ must be equal to infinity.Considering τ ∞ ≤ τ e , we conclude that τ e = ∞ is indeed equal to infinity.□

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 |x 0 | ≤ 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 ∈ (0, 1 + 2δ 1 /ξ 2 2 ), we can find a positive constant L 1 such that (3.9) Proof.we are going to prove the moment boundedness of T (t) first.Consider the following auxiliary procedure φ(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 ≤ T (t) ≤ φ(t) (t ≥ 0).By applying It ô's formula, the drift coefficient in the differential formula of φ p (t) can be expressed as (3.11) 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 ψ(t) where B 3 (t) is the standard Brownian motions defined in system (2.1).Then, we get 0 ≤ R(t) ≤ ψ(t) for t ≥ 0. Notice that ψ(t) follows the similar structure as φ(t), so we can obtain: where 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 ô's formula, we get The condition q ∈ (0, 1 > 0, then we can find a positive constant κ sufficiently small such that Define U(t) = e κt (y(t) + z(t)) q , we get where As Combining the continuity of q(y + z) (q−2) W(y, z), we can find a positive constant L satisfies:

Therefore,
L[e κt (y + z) q ] ≤ e κt L. (3.17) By integrating both sides of the inequality (3.17) from 0 to t and taking the expectation, we have lim By applying the comparison theorem [27], we have Extinction and persistence

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.Here, δ = min {δ 1 , δ 2 } , ξ 2 = min ξ 2 2 , ξ 2 3 .Proof.We denote Ψ(t) = H(t) + R(t).Taking the logarithm of Ψ(t) and applying It ô's formula yields By integrating both sides of the inequality (4.1) from 0 to t and dividing by t, we obtain log Combining the strong law of large numbers [26] Integrating the inequality (4.3) from 0 to t and dividing by t, results in As t approaches ∞, by combining the strong law of large numbers [26], we obtain that This implies that lim t→∞ T (t) = 0 a.s.(4.5) In other words, as t approaches ∞, 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, ∆t > 0 represents the time increment, ξ 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 × 10 9 (cells), H(0) = 4 × 10 6 (cells), R(0) = 3 × 10 8 (cells), respectively.First, verify the extinction of resting T cells R(t) with According to Theorem 4.1, resting T cells will eventually go extinct, which is consistent with the results shown in Figure 1.When ξ 1 = 0.8, ξ 2 = 2, ξ 3 = 2, the calculation yields α 1 − 1 2 ξ 2 1 = −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 ξ 2 = 2, ξ 3 = 2, and varying the noise intensity ξ 1 while ensuring the extinction condition for T (t).We take ξ 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 ξ 1 accelerates the extinction of tumor cells.This implies that stochastic environmental perturbations suppress the growth of tumor cells.

Persistence
Lemma 1 ( [31]).Assuming Z ∈ [Ω × [0, +∞) , R + ], if there exist positive constants η, t 0 and η 0 such that then we have 32]).Let T (t), H(t), R(t) be the solution of stochastic system (2.1) initial conditions 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] Using the It ô'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 ⟨x(t)⟩ = 1 t t 0 x(s)ds.Then, we can obtain According to the strong law of large numbers [26], we get For sufficiently small ε > 0, when α 1 − 1 2 ξ 2 1 > 0 , according to Lemma 1, it can be deduced that substituting Eq (4.15) into Eq (4.13) and combining with the Eq (4.14), it yields lim Similarly, when δ + 1 4 ξ 2 − α 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 lim Using the same method as in Theorem 4.2, the condition substituting Eq (4.19) into Eq (4.13) and combining with Eq (4.14), it yields lim According to formula Eqs (4.16), (4.18), (4.20) and using the comparison theorem [26], we get lim the tumor cell population T (t) is weakly persistently bounded, with lim sup t→∞ ⟨T (t)⟩ > 0.Moreover, it can be deduced that sup t→∞ ⟨T (t Proof.According to the It ô's differential rule, we have By integrating both sides from 0 to t and dividing by t, and taking the upper limit, we obtain lim sup Applying Lemma 2, we have By using the strong law of large numbers [26], we get lim t→∞ 1 t t 0 ξ 1 dB 1 (t) = 0 a.s., and when δ − ξ 2 4 − α 2 − k > 0, we have lim t→∞ H(t) = 0. Thus, lim sup t→∞ ⟨H(t)⟩ = 0 a.s.. Substituting it into Eq (4.23), we obtain Therefore, lim sup t→∞ ⟨T (t)⟩ > 0, which further indicates that the tumor cell population T is weakly persistent.In addition, log Applying Lemma 1, when 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.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 ξ 1 = 0.1, 0.3, 0.5 and using a time step of ∆t = 0.01, we sample N = 1×10 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 ξ 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.

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 ξ 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.

Figure 6 .
Figure 6.Stationary probability distributions of tumor cells and two types of T cells (a) Simulation of stationary probability distributions of tumor cells and two types of T cells with values of ξ 1 = 0.1, ξ 2 = 2, ξ 3 = 2. (b)Simulation of stationary probability distributions of tumor cells and two types of T cells with values of ξ 1 = 0.3, ξ 2 = 2, ξ 3 = 2. (c)Simulation of stationary probability distributions of tumor cells and two types of T cells with values of ξ 1 = 0.5, ξ 2 = 2, ξ 3 = 2.
are the population of tumor cells, predator hunting T cells and predator resting T cells in the given physiologic space, respectively.α 1 and α 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.β 1 and β 2 represent the rate at which hunting T cells kill tumor cells and the rate at which tumor cells deplete hunting T cells.δ 1 and δ 2 indicate the apoptosis rates of hunting T cells and resting T cells.γ represents the activation rate of resting T cells to hunting T cells, k is the proliferation rate of resting T cells, η is the half-saturation coefficient of proliferation term.B 1 (t), B 2 (t), B 3 (t) are mutually independent standard Brownian motions, and ξ 1 , ξ 2 , ξ 3 represent noise intensity of them.
T (t), H(t) and R(t)

Table 1 .
Biological significance and values of parameters.