DYNAMICS OF A MODEL OF TUMOR-IMMUNE INTERACTION WITH TIME DELAY AND NOISE

. We propose a model of tumor-immune interaction with time delay in immune reaction and noise in tumor cell reproduction. Immune response is modeled as a non-monotonic function of tumor burden, for which the tumor is immunogenic at nascent stage but starts inhibiting immune system as it grows large. Without time delay and noise, this system demonstrates bistability. The eﬀects of response time of the immune system and uncertainty in the tumor innate proliferation rate are studied by including delay and noise in the appropriate model terms. Stability, persistence and extinction of the tumor are analyzed. We ﬁnd that delay and noise can both induce the transition from low tumor burden equilibrium to high tumor equilibrium. Moreover, our result suggests that the elimination of cancer depends on the basal level of the immune system rather than on its response speed to tumor growth.


1.
Introduction. The immune system is a host defense system that distinguishes pathogens from one's own healthy cells and destroy them. Pathogens include bacteria and viruses as well as abnormal cells. There are evidence that immune systems can detect and eliminate cancer [27]. Cancer immunology has seen renewed interests due to recent achievements in immunotherapy [22]. Unlike common pathogens, cancer cells are not as distinguishable. Moreover, cancer can evade the control of immune systems by developing ways to achieve immune suppression [30]. Recent progress in immunotherapy has focused on removing immune suppression so that cancer cells will be under attack of immune systems. Further development of immunotherapy hinges on a better understanding of tumor-immune system interaction.
The immune system consists of immune cells of several types and complex signaling network. For example, immune cells include cytotoxic T cells, helper cells, natural killer cells, and signaling networks involve cytokines, e.g., interlukin-2. The immune response falls into two categories: innate and adaptive. Though the adaptive response enables fast reaction to certain antigens, immune responses oftentimes involve recruiting immune cells from bone marrows and further training or activation. So it is inherently a delayed process. On the other hand, cancer is well-known for its heterogeneity and characterized by fast mutation, which can render it the ability to suppress immune response [6,30]. In this process, stochasticity plays an important role.
Given the complexity of tumor-immune system interaction, mathematical modeling naturally offers some insights by capturing the key mechanisms. Literature on modeling non-spatial tumor growth under immune surveillance is abundant (see [5] for a review). In [17], the authors derived a system of five equations from a kinetic scheme and further reduced it to two equations tracking only effector cells and tumor cells. In their model, immune suppression is represented as annihilation by mass action of tumor cells and effector cells. Their paper inspired a number of later work. In [3], the author summarized models of this type into a family as follows where x is the size or density of tumor population, y is the size or density of immnue effector population and θ(t) represents treatment . With biological-relevant assumptions on f (x), φ(x) and q(x), some general conclusions of were drew in [3].
In [13], the authors took into account cytokines in their model in order to gain further insight into immunotherapy. In most cases, introducing more equations does not introduce new dynamics but is necessary to shed light on cancer treatments [5]. Signaling is arguably a fast process, so sometimes a quasi-steady state approximation can reduce more complicated models into the family of models studied by [3]. For example, in a recent study on immune check point inhibitor [25], a variation of (1) is derived from a more complicated model consisting of a system of 14 partial differential equations [18].
Building on aforementioned earlier work, the effect of time delay in tumorimmune interactions are also extensively studied [1,4,8,28]. On the other hand, stochasticity in tumor-immune interaction is less often studied, and mostly focused on single-equation models of tumor population growth ignoring explicit interaction with immune system [20,2,21]. Models that take into account both time delay and stochasticity are even rarer. To the best of our knowledge, the work of [11] is the only exception. In [11], a single equation was considered and noise was introduced in an ad-hoc way. To fill this gap, we will introduce time delay and noise one at a time to a variation of (1) with the novelty in its immune response term being non-monotonic.
The paper is organized as follows. In section 2, we introduce our model formulation and state general properties of the model in absence of noise and time delay. In section 3, we study the model with time delay where persistence, local stability and global stability results are presented. Section 4 is devoted to the stochastic version of the model where tumor extinction and persistence conditions are presented. In section 5, the paper comes to the culmination where both noise and time delay are present in the model, for which we study the effects of time delay on the stationary distribution. We conclude the paper with a discussion of future work in section 6.
2. ODE model formulation and properties. We consider the following equations where x denotes the tumor burden and y denotes the level of immune response. This is a special case of (1) for which we specify σq(x) =β +α x κ+x 2 to represent the immune response to tumor burden in a non-monotonic fashion. In spite of being phenomenological, it accounts for the fact that the cancer cells acquire mutations that can down-regulate immune response as the tumor grows bigger while the tumor is immuogenic at its nascent stage. Death term is assumed to be independent of tumor burden and thus we have Ψ(x)y = −μy. Note that without tumor, immune system activity is maintained at the basal level β/µ. We assume logistic growth of tumor and that killing of tumor cells by the immune system is represented as a mass action term so that we have f (x) = 1− x K and φ(x) =γx. Both are common choices in modeling literature [8,25]. Even though Gompertz growth is also commonly used in modeling tumor growth and backed up by a lot of experimental data, it predicts unbound growth rate as tumor burden approaches zero [5]. Hence it is not suitable to our purpose.
Note that we represent the two populations generically as tumor burden and level of immune response, and intentionally avoid using specific units of cell density or volume. In some sense, it justifies our choice of Ψ(x)y = −μy since we are not tracking explicitly effector cells. The rationale behind our generic representation is that our objective is to study interesting dynamics that can arise from a simple model with time delay and noise. It is our intention to bring attention to the possible key mechanisms underlying cancer-immune system dynamics without making the wrong impression that it is capable of clinical prediction.
By nondimensionlization, we can reduce the number of parameters. Let x = √ κu, y = α µ √ κ v and t = 1 µt . We get the dimensionless system where ρ =ρ µ , K =K √κ , γ =αγ µ 2 √κ and β =β √κ α are dimensionless (positive) parameters. It can be shown that there is always a tumor free equilibrium E 0 = {0, β} which is stable if and only if ρ < γβ by linear stability analysis. If ρ > γβ, depending on the parameters, there is additional either one or three interior equilibria. In this paper, we focus on the case where there are three interior equilibria, namely low tumor equilibrium E 1 , intermediate tumor equilibrium E 2 and high tumor equilibrium E 3 (see Figure 1). Since the intersections of nullclines are roots of a cubic polynomial, any results involving analytical expressions of parameters for stability conditions, if possible, would be unwieldy. Nevertheless, qualitative argument similar to that of in ( [24] Page 226-230) can be made regarding to the local stability . By doing so, we obtain the theorem below.
Theorem 2.1. Suppose that ρ > γβ and there are three interior equilibria E 1 , E 2 and E 3 of (4), then E 1 and E 3 are stable, and E 2 is unstable.
Proof. Consider the Jacobian matrix at the equilibrium where f u , f v , g u and g v are the partial derivatives of f (u, v) and g(u, v) evaluated at the equilibrium. The signs of f u , f v , g u and g v can be determined by geometric arguments. For example, to see the sign of g u at E 1 , we first note that g(u, v) < 0 in the region above g(u, v) = 0 and g(u, v) < 0 in the region below. Thus as it moves across g(u, v) = 0 at E 1 along the direction of u-axis from left to right, g(u, v) increases from being negative to being positive. Therefore, g u > 0 at E 1 . By the same argument, we find out the signs of the entries of Jacobian matrices J 1 , J 2 and J 3 at E 1 , E 2 and E 3 as follows To determine stability, we need to find signs of their determinants too. This can be done by noting the sign of the slopes of the nullclines. For example, at . Similarly, we have det(J 2 ) < 0 and det ( J 3 ) > 0. Therefore, we proved the claimed stability.  3. With delay. It is known that it takes time for the immune system to respond and take effect. So we introduce time delay into the term representing immune response. The model now reads as follows where τ ≥ 0 is the time delay of immune response. We are curious about any possible new dynamics due to the time delay in contrast to the ODE model, especially if the time delay in immune response can lead to oscillatory solutions that may represent cancer growth episodes.
3.1. Analysis. First we want to establish positivity and boundedness of the system (6) given appropriate initial values. For initial values, we assume v(0) ≥ 0 and The boundedness is straightforward.
The above proposition ensures that the solution of (6) is biological meaningful and will also be useful prerequisite for our later analysis. Next we present a theorem of the local stability of (6) as a counterpart to Theorem 2.1. In particular, we show that as τ increases there is a possible stability switching for low tumor equilibrium E 1 . Recalling the same definition for f u , f v , g u and g v as before, we have the following theorem.
Theorem 3.1. Suppose that ρ > γβ. Then the tumor free equilibrium E 0 is unstable for for all τ ≥ 0. In addition, suppose that there are three interior equilibria E 1 , E 2 and E 3 of (6). Then E 3 is stable, E 2 is unstable for all τ ≥ 0 and E 1 change its stability as τ increases if f u g v + f v g u < 0.
Proof. Linearizing (6) at the equilibrium gives Substituting u = Ae λt , v = Be λt into the above equations yields a linear system of A and B for which the existence of a solution entails the characteristic equation At E 0 , (7) becomes λ 2 − (ρ − γβ)λ − ρ + γβ = 0, which always has a positive root and thus E 0 remains unstable for all τ ≥ 0. Theorem 2.1 gives the stability results of E 1 , E 2 and E 3 at τ = 0. We want to study if their stability changes as time delay τ increases. If the equilibrium changes its stability as τ increases, there must exists ω > 0 such that λ = iω is a solution of (7) [14]. Substituting λ = iω into (7) gives Collecting real and imaginary parts gives Equating the real and imaginary parts to zeros respectively gives Squaring and adding the above equations gives From (5), we see that (9) is not satisfied for E 3 and hence it has no stability switching. Stability switching is possible for E 1 and E 2 . For E 1 , (9) entails an additional requirement which is stated in the theorem. The critical value of τ can be found by substituting ω = ω 0 to (8), which leads to To ensure stability switching, We need further confirm that the imaginary root crosses the imaginary axis as τ increases beyond τ 0 . To do so, we think λ as a function of τ and differentiating (7) gives It follows that So for E 1 which is stable at τ = 0, the imaginary root indeed crosses the imaginary axis and stability switching is ensured. But for E 2 which is unstable at τ = 0, the imaginary root does not cross the imaginary axis and there is no change in its stability.
The above theorem indicates that the time delay in immune response can lead to tumor escaping the control of the immune system as a consequence of loss of stability of the low tumor equilibrium. The condition (10) for this stability switching is not explicit in terms of parameters. Nevertheless, it has a geometric interpretation that the slope of the v-nullcine is larger than the magnitude of the slope of u-nullcline at their interception E 1 . In a biological sense, it means that there would be instability if the objective of immune response level is set to increase steeply with respect to tumor burden while immune response time is not fast enough to catch up.
Next we study the global stability of the tumor free equilibrium. Proof. We first shift v to make its equilibrium value to be zero by letting w = v − β.
If the stability condition of tumor-free equilibrium is not satisfied, we can prove that the tumor population is persistent. Proof. Suppose not. Then lim t→∞) u(t) = 0. From this, it can be shown that lim sup t→∞ v(t) ≤ β. There are two cases we need to consider. The first case is that u is eventually monotonically decreasing to zero. Then lim t→∞ v(t) = β. It follows that there is a t * > 0 such that u(t * ) is sufficiently small and v(t * ) is sufficiently close to β so that u (t * ) = u(ρ − ρu K − γv) > 0. So we have a contradiction. The second case is that u approaches zero in an oscillatory fashion. That is, there is a sequence of {t i } such that u(t i ) → 0 and u(t i ) is local minimum, i.e., u (t i ) = u(ρ − ρu K − γv) = 0. Then for any δ > 0 there is an N (δ) such that u(t i ) < δ for i > N (δ). Since ρ > γβ, there is a δ such that ρ − ρδ K − γ(β + δ) > 0. We see that for this δ, there is a T 1 such taht t > T 1 implies that v(t) < β + δ and an N such that for i > N , t i > T 1 . Hence for i > N , u (t i ) > 0 which is a contradiction to u (t i ) = 0.
Putting the conditions of the above two theorems in original units, we find an important dimensionless value R 0 =ρμ γβ . It is the ratio of the proliferation ability of the tumor verse the strength of the immune system, which can be viewed as the tumor reproduction number: if the tumor can be totally cleaned out by immune system. The biological meaning of (11) is more clear by recognizing that it is equivalent to saying that at the onset of tumor growth, the killing rate of tumor by the immune system (γβ/μ) is larger than the tumor growth rate (ρ). Theorem 3.3 tells us that otherwise (R 0 > 1) the tumor will always exist.
3.2. Numerical simulation. We studied (6) numerically using dde23 in Matlab. We confirmed the criteria for tumor persistence and extinction in our numerical experiment. The existence of a stability switching for low tumor equilibrium as τ increases is also demonstrated ( Figure 2). Moreover, we discovered that as τ is further increased, there is a transition to high tumor equilibrium (bottom right of Figure 2). Interestingly, the jerky transition to high tumor equilibrium is reminiscent of saltatory growth which is a commonly observed pattern in tumor growth [15]. This suggests that a weakened immune system response may be indicated by its slow action which takes longer time delay, which may result in increased tumor growth.
4. With noise. Tumor cell growth is known to lack of certain regulations compare to normal cells, we thus assume that the innate proliferation rate of tumor cells is subject to uncertainty. Adding white noise to it results in the following stochastic differential equations where B(t) is a scalar Brownian motion defined on the complete probability space (Ω, F, {F t } t≥0 , P ) with filtration {F t } t≥0 . We use a ∧ b to denote min(a, b), a ∨ b to denote max(a, b) and a.s. to denote almost surely.

4.1.
Analysis. We first show existence and uniqueness of a global solution remains in D = (0, K)×(β, β +u m ) where u m = max u∈(0,K) { u 1+u 2 }, whenever it starts in D. Because the drift term and noise term in (12) do not satisfy linear growth condition, the general existence and uniqueness theorem (see [23] chapter 2) does not apply. The techniques employed here are standard and follow the similar lines as [29,9]. Let D n = ( 1 n , K − 1 n ) × (β + 1 n , β + u m − 1 n ) where n ∈ N + . Define the stopping time τ n = inf{t ∈ (0, τ e ) : (u(t), v(t)) / ∈ D n }. Clearly τ n increases as n. Let τ ∞ = lim n→∞ t n and we have τ e ≥ τ ∞ . If we can show that τ ∞ = ∞ a.s., then τ e = ∞ a.s. and we have the unique solution (u(t), v(t)) ∈ D for all t. Suppose on contrary that there is a pair of constants T > 0 and ∈ (0, 1) such that P (τ ∞ ≤ T ) > . Then there is a n 1 such that P (τ n ≤ T ) ≥ for n > n 1 .
where we note it is nonnegative for u, v ∈ D.
In general, given an Ito process dX(t) = b(X(t)dt + σ(X t )dB(t) and if f ∈ C 2 0 (R n ), then its generator A acting on f gives (see Theorem 7.3.3 on Page 126 of [26]).
Consider the generator A of (12) and let A act on V , i.e., AV : D → R. Some calculation shows that where C = ρ k ∨γ(β +u m )∨(γ(β +u m )+1) a positive constant. By Dynkin's formula (see, e.g., [26] Page 127), we have By Gronwall's inequality, On the other hand, consider the set Ω n = {τ n ≤ T } for n > n 1 for which we know P (Ω n ) ≥ . For ω ∈ Ω n , we know that V ∼ O(n). Thus if we let n → ∞. Hence we have a contradiction.
Proof. By Ito's formula, we have where f (y) = − 1 2 σ 2 y 2 + ρy − γβ. Note that f (y) is a concave-down quadratic function with axis of symmetry y = ρ σ 2 and y(0) < 0. There are two ways to ensure max y∈(0,1) f (y) < 0, i.e., which corresponds to the stated conditions i) or ii). Thus, lim sup where the term with integral is zero a.s. by law of large numbers of martingale (see, e.g., [23]).
To complement with the condition for tumor extinction, we present a theorem for tumor persistence in the following. The proof involves construction of a contradiction to show that lim sup t→∞ u(t) is bounded from below, and again Ito's formula is the workhorse behind this type of proof.
Proof. Assume ρ > σ 2 /2 + γ(β + u m ). We note that g(K) < 0 and g(0) = − 1 2 σ 2 + ρ − γ(β + u m ) > 0. Thus ξ exists and g (u) < 0 for u in a small neighborhood of ξ. Suppose on contrary that there is a small ∈ (0, 1) such that P (Ω 1 ) > where Ω 1 = {ω : lim sup t→∞ u(t, ω) < ξ − 2 }. Then there is T (ω) > 0 such that We choose small enough so that g(u(t, ω) > g(ξ − ) for t > T (ω). Moreover, by law of large numbers of martingale there is Ω 2 ⊂ Ω with P (Ω 2 ) = 1 such that for any ω ∈ Ω 2 , lim t→∞ Thus we have lim inf which implies that u(t, ω) → ∞, which contradicts (14).  Table 1. Parameter values used in Figure 3 We compare the conditions in the previous two theorems to their counterparts in section 3. In the original units, we find two noise magnitude regimes withρ µ as the threshold (we call σ 2 <ρ µ small noise regime and σ 2 >ρ µ big noise regime). In the small noise regime, the tumor extinction condition is weaker than the one for the deterministic system. The biological interpretation is that the small noise can possibly help eliminate the tumor. The condition for tumor persistence is stronger in the stochastic system. However, it is only a sufficient condition and there is a gap between between the persistence condition and the extinction condition.

4.2.
Numerical simulation. (12) is simulated using Milstein's method [12]. As in Figure 3, the simulation confirmed our extinction criteria. Also shown is that the tumor persists in the monostability where the solution simply fluctuates about the equilibrium and bistablity where there is a noise-induced transition to high tumor equilibrium from low tumor equilibrium. The parameter values used in the simulations are summarized in Table 1. 5. With delay and noise. In this section, we study the system when both delay and noise are present as follows 5.1. Analysis. The analysis in section 4 can be extended to (15). Because of the time delay τ , here we need to supply suitable initial data on [−τ, 0]. The extension of Theorem 4.1 is stated as follows there is a unique solution defined for all t ≥ −τ and remains in D a.s.
Proof. The proof follows similar lines as the one for Theorem 4.1 with an application of method of steps. We keep the same definition of τ n , τ ∞ , τ e and V (u, v). Same as before we want to show that τ e = ∞ by proving that τ ∞ = ∞. First we want to show that τ ∞ > τ . For any n ∈ N + and t ∈ [0, τ ∞ ), it can be shown by Ito formula that  Table 1 where where K 1 is a positive constant. For t 1 ∈ [0, τ ], integrating (16) from 0 to t 1 ∧ τ n and taking expectation gives , v(0))e τ < ∞. In particular, Suppose on contrary there is > 0 such that P (τ ∞ < τ ) > . Then there is n 1 such that for n > n 1 , P (Ω 1 ) ≥ where Ω 1 = {ω : τ n (ω) < τ }. Thus ≥ n → ∞ as n → ∞, which is a contradiction to 17. Thus for τ ∞ ≥ τ . By the same argument on t ∈ [τ, 2τ ],we have τ ∞ ≥ 2τ . Repeating this procedure gives τ ∞ = ∞.
It is easy to see that v(t) ∈ (β, β + 1 2 ) for all t > 0. This enables Theorem 4.2 and Theorem 4.3 be extended to the system (15) with similar arguments.

Numerical simulation.
From previous analysis, we notice that the system (15) bears a lot of similarity to the system (12). In this subsection, we instead focus on the effect of time delay on stationary distributions by numerical simulations. As seen in Figure 4, the stationary distribution is bimodal which is not a surprise since the underlying deterministic system is bistable. As the delay increases (larger τ ), the more density shifts to the high tumor stable state and the bimodality becomes indistinguishable at τ = 2. We also observe that the mean first passage time from the low tumor stable state to the high tumor stable state is reduced by increasing τ .
6. Discussion and conclusion. In this paper, we presented a simple model of tumor-immune system interactions, in which the immune response to tumor is modeled as a non-monotonic function of tumor burden. We studied the effects of time delay in the immune response and the uncertainty in the innate proliferating rate of cancer cells on the tumor growth dynamics. The conditions of tumor extinction and persistence were proved in presence of noise or time delay. We also performed numerical experiments to confirm analytical results and to guide future analytical work.
We showed that the magnitude of the tumor reproduction number R 0 =μρ γβ relative to 1 dictates the stability of tumor-free equilibrium. This condition does not depend on time delay. It suggests that the elimination of cancer depends on the basal level of the immune system rather than on its response speed to tumor growth. However, it maybe possible to have delay-induced stability switching for the low-tumor steady state. We also established the global stability of tumorfree equilibrium, but that of interior equilibria is challenging and remains an open problem.
For stochastic version of the model, we showed that the noise can help eliminate cancer in either case of big or small noise. Similar results in an epidemic model were obtained by [9]. We note that the criteria for persistence in Theorem 4.3 is not a necessary condition and in fact far from being a sharp result. Indeed, parameter sets (c) (d) do not satisfy the persistence condition but appear to be persistent in the numerical simulation. The improvement of the current result will be deferred to future work. Also, as seen in the simulation, there is a stationary distribution of (12). However, it is challenging to show this analytically because the diffusion matrix is degenerate, which makes the standard techniques employed in [29,9] not applicable. There have been recent studies on the stationary distribution resulting from a degenerate diffusion matrix [19]. It will also be the focus of our future work.
When including time delay with noise, we showed that the results we obtained for noise-only system (12) can be easily extended to (15). Moreover, we found numerically that the noise favors transition to high tumor stable state, which was also observed for an ecological model studied in [31]. The biological implication is that the less responsive the immune system is, the easier for the tumor to escape to high burden. Same as for the noise-only version of the model, the existence of a stationary distribution of (15) will be the focus of future work. There are some perturbation techniques applied to delayed Focker-Planck equation to study the effects of time delay on mean first passage time [16,10,7]. Those techniques are limited to a scalar equation. Possible extension to study a system of equations will be carried out in the future.
Time delay and stochasticity are two hallmarks of cancer dynamics. Serious modeling work shall not shy away from them. The model we studied was kept minimal but nevertheless exhibited interesting dynamics. It is hoped that this paper will stimulate interests in stochastic delayed differential equations in modeling cancer dynamics.