Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations

Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations

Journal Pre-proof

Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations Ying Deng, Meng Liu PII: DOI: Reference:

S0307-904X(19)30598-0 https://doi.org/10.1016/j.apm.2019.10.010 APM 13071

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

18 April 2019 19 September 2019 2 October 2019

Please cite this article as: Ying Deng, Meng Liu, Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.10.010

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights • A stochastic hybrid tumor-immune model with impulse is developed. • Threshold between persistence and extinction is obtained. • White noises on growth rate and destruction rate of neoplasm are conducive to the elimination of neoplasm. • If neoplasm stays longer in states with negative “stochastic growth rate”, then it will be eliminated. • Positive impulse could make neoplasm become persistent, negative impulse could make neoplasm become extinct.

1

Analysis of a stochastic tumor-immune model with regime switching and impulsive perturbations Ying Deng, Meng Liu∗ School of Mathematical Science, Huaiyin Normal University, Huaian 223300, PR China.

Abstract In this article, an impulsive stochastic tumor-immune model with regime switching is formulated and explored. Firstly, it is proven that the model has a unique global positive solution. Then sufficient criteria for extinction, non-persistence in the mean, weak persistence and stochastic permanence are provided. The threshold value between extinction and weak persistence is gained. In addition, the lower- and the upper-growth rates of tumor cells are estimated. The results demonstrate that the dynamics of the model are intimately associated with the random perturbations and impulsive perturbations. Finally, biological implications of the results are addressed with the help of real data and numerical simulations. Keywords: tumor-immune model, random noises, impulse, persistence, extinction 1. Introduction Tumors are a serious threat to human lives and health. For example in 2018, more than 9.6 million people died from cancer, and more than 18.1 million people became new cancer cases [1]. According to a report released by the World Health Organization in 2016, for people under the age of 70, malignant tumor is the top 2 leading cause of death in 91 countries [2]. This motivates scholars to explore the inherent law of neoplasm evolution. These studies are of theoretical and clinical values for the prevention and treatment of neoplasm [3]. Mathematical models have been proved to be powerful in oncology not only in expounding some real phenomena [4, 5], but also in assisting doctors to make better and more scientific treatments [6, 7]. In recent decades, many deterministic mathematical models for tumor growth have been proposed, for example, the logistic model, the Gompertz model and the von Berta∗

Corresponding author. Tel:+86051783525318 Email address: [email protected] (Meng Liu)

Preprint submitted to Elsevier

October 14, 2019

lanffy model [8]. However, these deterministic models ignored the effects of immune system and environmental perturbations. Studies [9] demonstrated that the immune system can enhance T cell response against malignant tumors. In addition, the growth of malignant tumors relies heavily on the synthesis of proteins which is very sensitive to the environmental perturbations [10]. Therefore, researches on stochastic tumor-immune systems have received considerable concern in recent decades, for example, to study the properties of the stationary probability distribution by solving the corresponding Fokker-Planck equation ([11, 12, 13]), to explore the stochastic resonance by using the approach of potential functions ([14]), to investigate the persistence and extinction by applying Lyapunov’s second method ([3, 15, 16, 17]), to numerically show the dynamics of the model ([18, 19]). On the other hand, many patients with malignancy receive treatments, for example, surgery, radiation, chemotherapy, targeted therapy. Generally, these treatments could reduce the number of tumor cells in a short time [20]. In addition, metastasis is an important characteristic of malignancy, which is the main cause of mortality (about 90%) [21]. Three primary routes of metastasis are: through the blood system [22], through the lymphatic system [22], and transcoelomic metastasis [23]. Therefore, in order to describe the law of tumor evolution much better, one might need to consider stochastic tumor-immune systems with impulsive perturbations. However, the study of this topic still in its infancy. Motivated by these, this paper explores a stochastic tumor-immune model with impulsive perturbations. We formulate the model in Section 2 and give some preliminaries in Section 3. In Section 4, we establish sufficient criteria for extinction, non-persistence in the mean, weak persistence and permanence of the model. The critical value between extinction and weak persistence is gained under a hypothesis. In Section 5, we estimate the lower- and the uppergrowth rates of tumor cells. In Section 6, we discuss the effects of random perturbations and impulsive perturbations on the dynamics of the model, and illustrate these effects with the help of real data and numerical simulations. Finally, we give some concluding remarks in Section 7. 2. Model formulation 2.1. The deterministic model Let Ψ(t) stand for the number of tumor cells in a volume element of a tissue (in units: cells/ml) at time t. Without the immune reaction, it is assumed that the tumor evolution follows the 3

logistic growth, i.e.,

  dΨ(t) = Ψ(t) r − ρΨ(t) , dt

(1)

where r and ρ represent the growth rate (in units: day−1 ) and the intra-specific competition rate (in units: ml/(cells·day)), respectively. Experimental data of B cell lymphoma in mice obtained by deactivating the immune system of the mice support the validity of this assumption [12, 24]. Now we consider the immune response. In general, the immune response has a saturation effect and relies on the number of the tumor cells [12]. Several authors [3, 11, 16, 18] have used

λΨ(t) 1+Γ−1 Ψ(t)

to describe the immune response, where λ > 0 represents the destruction rate of

neoplastic cells (in units: day−1 )), Γ > 0 characterizes the upper limit of the tumoral population per volume (in units: cells/ml). Adopting this idea, model (1) is replaced by:   λ dΨ(t) = Ψ(t) r − ρΨ(t) − . dt 1 + Γ−1 Ψ(t)

(2)

2.2. Introducing the stochastic perturbations We are in the position to take the stochastic perturbations into account. In the first place let us consider the small environmental perturbations which often exist in real tumor-immune reactions. For example, both the growth of malignant tumors and the immune response depends heavily upon the synthesis of proteins ([10]) which is quite sensitive to some environmental factors such as the supply of nutrients, temperature, pH value, drug concentration, and so on. R.M. May [25] have pointed out that the small environmental perturbations could be described by white noise. There are various ways to incorporate with the white noise. A widely adopted approach is to hypothesize that the white noise mainly influences one or more parameters in a system (see, e.g., ([3, 11, 17, 18],[26]-[37]). Following this approach, r → r + ξ1 B˙ 1 (t), − ρ → −ρ + ξ2 B˙ 2 (t), − λ → −λ + ξ3 B˙ 3 (t), where B˙ i (t) is a white noise and ξi2 is its intensity, i = 1, 2, 3. Then model (2) is replaced by   λ ξ3 Ψ(t) dΨ(t) = Ψ(t) r −ρΨ(t)− dt+ξ1 Ψ(t)dB1 (t)+ξ2 Ψ2 (t)dB2 (t)+ dB3 (t), −1 1 + Γ Ψ(t) 1 + Γ−1 Ψ(t) (3) where B(t) = (B1 (t), B2 (t), B3 (t))T is a three-dimensional Brownian motion defined on a given complete probability space (Ω, F, {Ft }t≥0 , P). Additional, the growth rate of tumors, the strength of the immune response and other parameters in a tumor-immune model often change from the current values to others because 4

of environmental variation. For instance, at 28 mg/kg, paclitaxel could effectively inhibit the growth of tumors, but at 10 mg/kg, the impacts could be neglected ([38]). The shifts of the values are memoryless, and the waiting time between two different shifts has an exponential distribution ([15, 27, 39, 40]). Thereby, one could utilize a continuous-time finite-state Markov chain to characterize these regime shifts [15, 27, 39, 40, 41]. As a consequence, from model (3), we obtain  λ(β(t)) dt + ξ1 (β(t))Ψ(t)dB1 (t) 1 + Γ−1 (β(t))Ψ(t) ξ3 (β(t))Ψ(t) + ξ2 (β(t))Ψ2 (t)dB2 (t) + dB3 (t), 1 + Γ−1 (β(t))Ψ(t)

 dΨ(t) = Ψ(t) r(β(t)) − ρ(β(t))Ψ(t) −

(4)

where β(t) is a continuous-time finite-state Markov chain. 2.3. Introducing impulsive perturbations Finally, taking the treatments and the metastasis of tumor cells into account, we introduce the impulse into model (4), and consider the following stochastic regime-switching model with impulsive perturbations:    λ(β(t))   dt + ξ1 (β(t))Ψ(t)dB1 (t) dΨ(t) = Ψ(t) r(β(t)) − ρ(β(t))Ψ(t) −    1 + Γ−1 (β(t))Ψ(t)  ξ3 (β(t))Ψ(t) dB3 (t), t 6= tk , k ∈ N. + ξ2 (β(t))Ψ2 (t)dB2 (t) +   1 + Γ−1 (β(t))Ψ(t)     Ψ(t+ ) = (1 + µ )Ψ(t ), k ∈ N, k k k

(5)

where N represent the set of all positive integers, 0 < t1 < t2 < ..., lim tk = +∞. When k→+∞

µk > 0, the impulsive perturbations imply that the tumor cells metastasize from other tissues to the tissue under consideration; When µk < 0, the impulsive perturbations indicate that some tumor cells are killed by treatments. Due to the fact that all treatments can not eliminate tumor cells completely, in this article we suppose that 1 + µk > 0 for ∀k ∈ N.

The evolution of system (5) is as follows. Hypothesize that at time t = 0, β(0) = i ∈ M,

where M = {1, ..., m} is the state space of β(t), then (5) obeys the following system    λ(i)   dΨ(t) = Ψ(t) r(i) − ρ(i)Ψ(t) − dt + ξ1 (i)Ψ(t)dB1 (t)    1 + Γ−1 (i)Ψ(t)  ξ3 (i)Ψ(t) + ξ2 (i)Ψ2 (t)dB2 (t) + dB3 (t), t 6= tk , k ∈ N.   1 + Γ−1 (i)Ψ(t)     Ψ(t+ ) = (1 + µ )Ψ(t ), k ∈ N, k k k 5

until β(t) jumps to another state, for instance, j ∈ M. Let t¯ stand for the jump-time. If t¯ does

not coincide with any impulse-time, then from the time t¯, (5) obeys the following system    λ(j)   dΨ(t) = Ψ(t) r(j) − ρ(j)Ψ(t) − dt + ξ1 (j)Ψ(t)dB1 (t)    1 + Γ−1 (j)Ψ(t)  ξ3 (j)Ψ(t) dB3 (t), t 6= tk , k ∈ N. + ξ2 (j)Ψ2 (t)dB2 (t) +   1 + Γ−1 (j)Ψ(t)     Ψ(t+ ) = (1 + µ )Ψ(t ), k ∈ N, k k k

(6)

until β(t) jumps again. If t¯ coincides with an impulse-time, for instance, tk∗ , then from the time t¯ = tk∗ , (5) obeys the following system    λ(j)   dΨ(t) = Ψ(t) r(j) − ρ(j)Ψ(t) − dt + ξ1 (j)Ψ(t)dB1 (t)    1 + Γ−1 (j)Ψ(t)  ξ3 (j)Ψ(t) dB3 (t), t > tk∗ , t 6= tk , k ∈ N. + ξ2 (j)Ψ2 (t)dB2 (t) +   1 + Γ−1 (j)Ψ(t)     Ψ(t+ ) = (1 + µ )Ψ(t ), k ≥ k ∗ , k ∈ N, k k k

until β(t) jumps again. In this paper, we always suppose that the jump-time of the Markovchain does not coincides with any impulse-time because the probability of coincidence is zero. However, our theoretical results given in the following sections still hold even though the jumptimes of the Markov-chain coincide with impulse-times. Model (5) is called a hybrid system. For each j ∈ M, (6) is called a subsystem of (5). Remark 1. We want to explain why we adopted Ψ(t+ k ) = (1 + µk )Ψ(tk ) not Ψ(tk ) = (1 + µk )Ψ(t− k ). Mathematically, the former means the solution is left continuous, the latter implies that the solution is right continuous. Biologically, these two assumptions are acceptable. The former is a traditional assumption and has been widely used in the literature, for example, the theory of impulsive differential equations [42], deterministic population models with impulse [43, 44, 45, 46], stochastic models with impulse [28, 29, 30, 31, 32, 33, 34, 37, 47]. In order to directly use some previous results (for example, the definition of the solution for stochastic differential equations with impulse [31]) and to be in accord with the traditional assumption, we adopted the former one. However, if one adopt the latter one, some similar results can be obtained by using the approaches given below.

6

3. Preliminaries Denote by H = (hij )m×m the generator of β(t), that is to say,     hij ∆t + o(∆t), P{β(t + ∆t) = j|β(t) = i} =    1 + h ∆t + o(∆t), ii

where hij ≥ 0 for i, j ∈ M with j 6= i, and

X j∈M

if j 6= i; if j = i,

hij = 0 for i ∈ M.

Throughout this article, it is hypothesized that a product equals 1 if the number of factors is 0. In addition, hypothesize that • β(t) is independent of B(t); min{ρ(i), λ(i)} > 0; i∈M

• β(t) is irreducible, which means that β(t) could jump from any state to any other state. As a result, β(t) has a stationary distribution, π = (π1 , π2 , ..., πm ), which is the solution of the following equation πH = 0 subject to

X i∈M

πi = 1 and πi > 0, i ∈ M.

(7)

For simplicity, define
au = max{a(i)}, al = min{a(i)}, b(i) = r(i) − λ(i) − ξ12 (i) + ξ32 (i) /2, i ∈ M, i∈M

i∈M

b=

X i∈M

πi b(i), bb = lim sup t−1 t→+∞

 X

0
 ln(1 + µk ) + b.

System (5) is a model for tumor growth, thus we need to give some criteria under which model (5) has a unique global positive solution. To this end, we recall a result. Consider a stochastic differential equation with Markovian switching of the form: dx(t) = f (x(t), β(t), t)dt + g(x(t), β(t), t)dB(t),

(8)

where f : Rn1 × M × [0, +∞) → Rn1 , g : Rn1 × M × [0, +∞) → Rn1 ×n2 . Lemma 1. (Theorem 3.15 and Theorem 3.17 in Mao and Yuan [48]) Assume that ( locally Lipschitz condition) for every T > 0 and integer k ≥ 1, there is a positive constant KT,k such that for all t ∈ [0, T ], i ∈ M and all x, y ∈ Rn1 with |x| ∨ |y| ≤ k,

|f (x, i, t) − f (y, i, t)|2 ∨ |g(x, i, t) − g(y, i, t)|2 ≤ KT,k |x − y|2 , 7

where x ∨ y = max{x, y}. Then (8) admits a unique maximum local solution. Theorem 1. For any initial data (Ψ(0), β(0)) ∈ (0, +∞) × M, model (5) has a unique solution

(Ψ(t), β(t)) ∈ (0, +∞) × M for all t almost surely (a.s.).

Proof. First of all, let us pay attention to the following equation  Y ξ 2 (β) λ(β) Q dφ(t) = r(β) − 1 − ρ(β) (1 + µk )eφ(t) − −1 2 1 + Γ (β) 0
1 + Γ−1 (β)

ξ (β) Q3

(9)

0
0
+ µk )eφ(t)

dB3 (t)

with initial data φ(0) = ln Ψ(0). Here we have dropped t from β(t) for simplicity. Since the coefficients in Eq.(9) obey the locally Lipschitz conditions, by virtue of Lemma 1, Eq.(9) admits a unique maximum local solution on [0, τe ), where τe represents the explosion time. Define ψ(t) = eφ(t) . Then Itˆo’s formula (see, e.g., [48]) implies that for t ∈ [0, τe ), (eφ(t) , β(t)) is the unique maximum local solution of the following equation:   Y λ(β) Q dψ(t) = ψ(t) r(β) − ρ(β) (1 + µk )ψ(t) − dt + ξ1 (β)ψ(t)dB1 (t) 1 + Γ−1 (β) 0
(10)

with ψ(0) = Ψ(0). Next we claim that τe = +∞. The proof is motivated by Theorem 2.1 in [39] and Theorem 2.1 in [49]. Choose a sufficiently large l0 such that ψ(0) < l0 . For every integer l, define τl = inf{t ∈ [0, τe ) : ψ(t) ≥ l}. Let τ∞ = lim τl . It is easy to see that τ∞ ≤ τe . If τe < +∞, then there are constants T > 0 l→+∞

and ε ∈ (0, 1) such that P{τ∞ ≤ T } > ε. Therefore, there exists an integer l1 ≥ l0 such that P{Ωl } ≥ ε, ∀ l ≥ l1 , 8

(11)

where Ωl = {ω : τl ≤ T }. Define [49] V (ψ) = ψ α , ψ > 0, 0 < α < 1. Then Itˆo’s formula means that dV (ψ) = αψ

α



r(β) − ρ(β)

Y

(1 + µk )ψ −

Γ−1 (β)

λ(β) Q

1+ 0
Hence

EV (ψ(τl ∧ T )) ≤ V (ψ(0)) + α|r|

u

≤ V (ψ(0)) + α|r|u It then follows from Gronwall’s inequality that

Z

τl ∧T

Z0 T 0

EV (ψ(s))ds

EV (ψ(τl ∧ s))ds.

u

EV (ψ(τl ∧ T )) ≤ V (ψ(0))eα|r| T .

(12)

For ω ∈ Ωl , V (ψ(τl , ω)) ≥ lα . This, together with (11) and (12), means that V (ψ(0))eα|r|

uT

≥ E[1Ωl (ω)V (ψ(τl , ω)] ≥ εlα .

Letting l → +∞ gives the contradiction ∞ > V (ψ(0))eα|r|

uT

= ∞.

Consequently, τe = +∞. Therefore, Eq. (10) has a unique solution (ψ(t), β(t)) ∈ (0, +∞) × M for all t. Let Ψ(t) =

Y

(1 + µk )ψ(t).

0
Hence Ψ(t) > 0. Notice that for every k ∈ N and tk ∈ [0, +∞), Y Y Ψ(t+ ) = lim (1 + µ )ψ(t) = (1 + µj )ψ(t+ j k k) + t→tk

0
= (1 + µk )

Y

0
(1 + µj )ψ(tk ) = (1 + µk )Ψ(tk ).

0
9

(13)

Furthermore, Ψ(t− k ) = lim− t→tk

Y

(1 + µj )ψ(t) =

0
Y

(1 + µj )ψ(t− k) =

0
Y

(1 + µj )ψ(tk ) = Ψ(tk ).

0
In addition, for ∀t 6= tk ,  Y  Y dΨ(t) = d (1 + µk )ψ(t) = (1 + µk )dψ(t) 0
λ(β) Q

1 + Γ−1 (β) 0
0
which completes the proof. 4. Persistence and extinction First of all, we recall several concepts. • If lim Ψ(t) = 0, then Ψ(t) is called extinction; t→+∞ Rt • If lim t−1 0 Ψ(s)ds = 0, then Ψ(t) is called non-persistent in the mean. t→+∞

• If lim sup Ψ(t) > 0, Ψ(t) is called weak persistent. t→+∞

• If for ∀ε ∈ (0, 1), there are two positive constants %1 = %1 (ε), %2 = %2 (ε) such that   lim inf P Ψ(t) ≥ %1 ≥ 1 − ε, lim inf P Ψ(t) ≤ %2 ≥ 1 − ε, t→+∞

t→+∞

then Ψ(t) is called stochastic permanent.

Assumption 1. mini∈M {Γ(i)ρ(i) − λ(i) − ξ32 (i)} > 0. Theorem 2. Let Assumption 1 hold. If bb < 0, then Ψ(t) will become extinct a.s..

10

 dt

Proof. An application of Itˆo’s formula to Eq.(10) results in  Y ξ 2 (β) λ(β) Q d ln ψ(t) = r(β) − 1 − ρ(β) (1 + µk )ψ(t) − −1 2 1 + Γ (β) 0
ξ (β) Q3

1 + Γ−1 (β)

0
0
As a result,

+ µk )ψ(t)

dB3 (t).

 Z t ξ12 (β(s)) r(β(s)) − ln ψ(t) = ln ψ(0) + ds 2 0 Z t Y λ(β(s)) Q − (1 + µk )ψ(s) + ρ(β(s)) −1 1 + Γ (β(s)) 0
(14)

k

where L1 (t) =

Z

t

ξ1 (β(s))dB1 (s), L2 (t) =

0

L3 (t) =

Z

Y

t

ξ2 (β(s))

0

Z

t

1+

0

ξ3 (β(s)) Q

Γ−1 (β(s))

(1 + µk )ψ(s)dB2 (s),

0
0
+ µk )ψ(s)

dB3 (s).

Notice that the quadratic variations of L1 (t), L2 (t) and L3 (t) are hL1 (t), L1 (t)i =

Z

0

t

ξ12 (β(s))ds

and hL3 (t), L3 (t)i =

Z

≤

(ξ12 )u t,

hL2 (t), L2 (t)i =

t

0

(1 + Γ−1 (β(s))

Z

0

t

 Y

ξ22 (β(s))

ξ32 (β(s)) Q

2 0
0
2 (1 + µk )ψ(s) ds,

ds ≤ (ξ32 )u t,

respectively. According to the strong law of large numbers for martingales (see e.g. [48], P. 16), one gets lim t−1 Li (t) = 0,

t→+∞

a.s., i = 1, 3.

(15)

At the same time, according to the exponential martingale inequality (see e.g. [48], ones derives     1 P sup L2 (t) − hL2 (t), L2 (t)i > 2 ln n ≤ 1/n2 . 2 0≤t≤n 11

Then the Borel-Cantelli lemma (see e.g. [48]) implies that that for almost all ω ∈ Ω, there exists an integer n0 = n0 (ω) such that for n ≥ n0 ,   1 sup L2 (t) − hL2 (t), L2 (t)i ≤ 2 ln n. 2 0≤t≤n Thus, for ∀0 ≤ t ≤ n, n ≥ n0 , 1 L2 (t) ≤ 2 ln n + hL2 (t), L2 (t)i = 2 ln n + 0.5 2

Z

0

 Y

t

ξ22 (β(s))

0
2 (1 + µk )ψ(s) ds.

In addition, a direction calculation indicates that if Assumption 1 holds, then   λ(i) ξ32 (i) ξ32 (i) min ρ(i)x + + = λ(i) + . x≥0 1 + Γ−1 (i)x 2(1 + Γ−1 (i)x)2 2

(16)

(17)

Applying (16) and (17) to (14), one can observe that for ∀0 ≤ t ≤ n, n ≥ n0 ,  Z t ξ12 (β(s)) + ξ32 (β(s)) ln ψ(t) − ln ψ(0) ≤ r(β(s)) − λ(i) − ds 2 0Z t Y (1 + µk )ψ(s)ds + 2 ln n + L1 (t) + L3 (t) − ρ(β(s)) =

Z

0

0
t

0

b(β(s))ds −

ρ(β(s))Ψ(s)ds + 2 ln n + L1 (t) + L3 (t).

0

Then (13) means that ln Ψ(t) − ln Ψ(0) ≤

X

ln(1 + µk ) +

Z

t

b(β(s))ds

0

Z0
(18)

0

As a consequence, for 0 < n − 1 ≤ t ≤ n, n ≥ n0 ,   Z t X −1 −1 −1 t ln Ψ(t)−ln Ψ(0) ≤ t ln(1+µk )+t b(β(s))ds+2(n−1)−1 ln n+t−1 L1 (t)+t−1 L3 (t). 0

0
In view of the ergodicity of β(·) and (15), one can see that  X  X −1 −1 lim sup t ln Ψ(t) ≤ lim sup t ln(1 + µk ) + πi b(i) = bb < 0, t→+∞

t→+∞

0
i∈M

which is the required assertion.

Theorem 3. Let Assumption 1 hold. If bb ≥ 0, then we have Z t bb −1 Ψ(s)ds ≤ l , a.s.. lim sup t ρ t→+∞ 0

Especially , if b = 0, Ψ(t) will be non-persistent in the mean a.s.. 12

(19)

Proof. For ∀ε > 0, there exists T > 0 such that for 0 < T < n − 1 ≤ t ≤ n, n ≥ n0 ,   Z t X −1 b(β(s))ds + 2 ln n + L1 (t) + L3 (t) ≤ bb + ε. ln Ψ(0) + ln(1 + µk ) + t 0

0
Let $ = bb + ε. It then follows from (18) that for all 0 < T < n − 1 ≤ t ≤ n, n ≥ n0 , Z t Z t X ρ(β(s))Ψ(s)ds + 2 ln n + L1 (t) + L3 (t) b(β(s))ds − ln Ψ(t) ≤ ln Ψ(0) + ln(1 + µk ) + ≤ $t − ρl

Define Λ(t) =

Rt 0

Z

0

0

0
Ψ(s)ds.

0 l

Ψ(s)ds. Consequently, eρ Λ(t) (dΛ/dt) ≤ e$t , t ≥ T. It follows that l

l

eρ Λ(t) ≤ eρ Λ(T ) + ρl $−1 e$t − ρl $−1 e$T .

Taking logarithm yields l −1

Λ(t) ≤ (ρ )



−1 $t

l

ln ρ $ e

ρl Λ(T )

+e

 −ρ$ e . l

−1 εT

Therefore, −1

lim sup t t→+∞

Z

0

t

l −1

Ψ(s)ds ≤ (ρ )

   −1 l −1 $t ρl Λ(T ) l −1 $T lim sup t ln ρ $ e + e −ρ$ e . t→+∞

Then L’Hospital’s rule implies that lim sup t t→+∞

−1

Z

0

t

Ψ(s)ds ≤

Using the arbitrariness of ε results in (23).

$ bb + ε = . ρl ρl

Theorem 4. If bb > 0, then Ψ(t) is weakly persistent a.s.

 Proof. Define F = ω : lim Ψ(t, ω) = 0 . Suppose that P{F } > 0. By (14), t→+∞

ln Ψ(t) − ln Ψ(0)  Z t X ξ12 (β(s)) = ln(1 + µk ) + r(β(s)) − ds 2 0 0
(20)

It is easy to see that for ∀ω ∈ F , lim Ψ(t, ω) = 0. As a result, t→+∞

−1

lim sup t [ln Ψ(t, ω) − ln Ψ(0)] < 0, t→+∞

−1

lim t

t→+∞

Z t ρ(β(s))Ψ(s)+ 0

lim t

t→+∞

−1

Z

t

ξ2 (β(s))Ψ(s)dB2 (s) = 0,

0

  X  λ(β(s)) ξ32 (β(s)) ξ32 (i) + ds = . πi λ(i)+ 1 + Γ−1 (β(s))Ψ(s) 2(1 + Γ−1 (β(s))Ψ(s))2 2 i∈M

Applying the above inequalities into (20), and then utilizing (15), one can observe that 0 > lim sup t−1 ln Ψ(t, ω) = bb > 0. t→+∞

This is a contradiction.

Remark 2. Under Assumption 1, Theorems 2 and 4 demonstrate that bb is the threshold between extinction and weak persistence of Ψ(t).

Permanence is one of the most important issues in biomathematics [50]. In order to investigate the permanence of model (5), we recall a lemma. Lemma 2. ([51], Lemma 2.3) Let ς ∈ Rm be a vector. If H is irreducible and πς = 0, then equation Hx = ς has a solution. Assumption 2. There exist two positive constants c1 and c2 such that for ∀t > 0, c1 ≤

Y

(1 + µk ) ≤ c2 .

0
Theorem 5. Under Assumption 2, if b > 0, then Ψ(t) is stochastically permanent. Proof. Define V1 (ψ) = 1/ψ 2 , ψ > 0. Applying generalized Itˆo’s formula to (10), one gets dV1 (ψ)

  Y 2 Y 2 = 2V1 (ψ) 1.5ξ2 (β) (1 + µk ) ψ 2 + ρ(β) (1 + µk )ψ − r(β) + 1.5ξ12 (β) 0
According to Lemma 2, equation Hx = −2b + ¯b(2, ..., 2)T has a solution (x1 , ..., xm )T , where

b = (b(1), ..., b(m))T . Therefore

1X hij xj + b(i) = ¯b. 2 j∈M 14

(21)

Let θ ∈ (0, 1) be sufficiently small such that for every i ∈ M,

By (21), we have b(i) −


1 − xi θ > 0, ¯b − θ ξ12 (i) + ξ32 (i) +

X xi θ hij xj > 0. 2(1 − xi θ) j∈M

(22)

X X X 1 1 xi θ hij (1 − xj θ) = b(i) + hij xj = ¯b + hij xj . (23) 2(1 − xi θ)θ j∈M 2(1 − xi θ) j∈M 2(1 − xi θ) j∈M

Define V2 (ψ, i) = (1 − xi θ)(1 + V1 (ψ))θ . Then the generalized Itˆo’s formula implies that Z t LV2 (ψ(s), β(s))ds, EV2 (ψ(t), β(t)) = V2 (ψ(0), β(0)) + E 0

where LV2 (ψ, i)



ξ12 (i) λ(i) Q + −1 (i) 2 1 + Γ 0
= 2(1 − xi θ)θ(1 + V1 (ψ))

− r(i) +

It is straightforward to see that there is a sufficiently small η > 0 such that for ∀i ∈ M   X xi θ η 2 2 b − θ ξ1 (i) + ξ3 (i) + hij xj − > 0. 2(1 − xi θ) j∈M 2θ 15

(24)

Define V3 (ψ, i) = eηt V2 (ψ, i). In view of Itˆo’s formula, EV3 (ψ(t), β(t)) = V2 (ψ(0), β(0)) + E

Z

0

t

LV3 (ψ(s), β(s))ds,

where LV3 (ψ, i)

= ηeηt (1 − xi θ)(1 + V1 (ψ))θ + eηt LV2 (ψ, i)     θ−2 2 2 ηt − b − θ ξ1 (i) + ξ3 (i) + ≤ e 2(1 − xi θ)θ(1 + V1 (ψ))

 X xi θ η hij xj − V 2 (ψ) 2(1 − xi θ) j∈M 2θ 1 P   η j∈M hij (1 − xj θ) 1.5 2 2 2 2 + ρ(i)c2 V1 (ψ) + − r(i) + 1.5(ξ1 (i) + ξ3 (i) + ξ2 (i)c2 ) + λ(i) + + V1 (ψ) (1 − xi θ)θ θ  X η 1 hij (1 − xj θ) + + ρ(i)c2 V10.5 (ψ) + 1.5ξ22 (i)c22 + =: eηt f (ψ, i). 2(1 − xi θ)θ j∈M 2θ

Thanks to (24), f1 := sup ψ>0

It follows that,



 max f (ψ, i) < +∞. i∈M

(1 − xi θ)E[eηt (1 + V1 (ψ(t)))θ ] ≤ (1 − xi θ)(1 + ψ −2 (0))θ + η −1 f1 (eηt − 1). As a result, lim sup E[ψ −2θ (t)] = lim sup E[V1θ (ψ(t))] ≤ lim sup E[(1+V1 (ψ(t)))θ ] ≤ t→+∞

t→+∞

t→+∞

f1 =: f2 . η(1 − θ maxi∈M {xi }) (25)

By virtue of (13) and Assumption 2, we get  Y −2θ  −2θ −2θ lim sup E[Ψ (t)] = lim sup E (1 + µk ) ψ (t) ≤ c−2θ 1 f2 := f3 . t→+∞

t→+∞

For ∀ε > 0, choose %1 =

ε f3


2θ1

0
. According to Chebyshev’s inequality, one obtains

  E[Ψ−2θ (t)] −2θ P Ψ(t) < %1 = P Ψ−2θ (t) > %−2θ ≤ = %2θ (t)]. 1 1 E[Ψ −2θ %1  As a consequence, lim sup P Ψ(t) < %1 ≤ %2θ 1 f3 = ε. Therefore t→+∞

 lim inf P Ψ(t) ≥ %1 ≥ 1 − ε. t→+∞

16

We are now in the position to demonstrate that exists a positive constant %2 such that  lim inf P Ψ(t) ≤ %2 ≥ 1 − ε. According to Itˆo’s formula, t→+∞

d(et V (ψ))   Y t α α = e ψ + αψ r(β) − ρ(β) (1 + µk )ψ −

λ(β) Q

1 + Γ−1 (β) 0
where c3 > 0 is a constant. As a result,

lim sup E[ψ α (t)] ≤ c3 , t→+∞

which implies that, α   Y α (1 + µk ) ψ (t) ≤ cα2 c3 . lim sup E[Ψ (t)] = lim sup E α

t→+∞

t→+∞

0
An application of Chebyshev’s inequality yields the desired assertion. 5. Lower- and the upper-growth rates Now let us estimate the lower- and the upper-growth rates of tumor cells. Theorem 6. If Assumption 2 holds, then lim sup t→+∞

ln Ψ(t) ≤ 1, ln t

a.s.

(26)

Proof. Applying Itˆo’s formula to Eq.(10) results in d[et ln ψ]   Y 2 Y ξ12 (β) ξ22 (β) t = e ln ψ + r(β) − − ρ(β) (1 + µk )ψ − (1 + µk )ψ 2 2 0
17

Therefore et ln ψ(t) − ln ψ(0) Z t  Y ξ 2 (β(s)) es ln ψ(s) + r(β(s)) − 1 − ρ(β(s)) (1 + µk )ψ(s) = 2 0 0
where

K1 (t) =

Z

t s

e ξ1 (β(s))dB1 (s),

K2 (t) =

0

Z

t

es ξ2 (β(s))

0

K3 (t) =

Z

t

0

Define

Y

(27)

(1 + µk )ψ(s)dB2 (s),

0
es ξ3 (β(s)) Q dB3 (s). 1 + Γ−1 (β(s)) 0
3 X

Ki (t).

i=1

Then

  Y 2 2 2 e ξ1 (β(s)) + ξ2 (β(s)) (1 + µk )ψ(s) hN (t), N (t)i = 0 0
t

2s

Thanks to the exponential martingale inequality, one can see that for arbitrary σ > 1 and ν > 0,     e−νj νj P sup N (t) − hN (t), N (t)i > σe ln j ≤ j −σ . 2 0≤t≤νj Then the Borel-Cantelli lemma implies that for almost all ω ∈ Ω, there exists a j0 such that for ∀j ≥ j0 , N (t) ≤

e−νj hN (t), N (t)i + σeνj ln j, 2

0 ≤ t ≤ νj.

As a consequence, for j ≥ j0 , 0 ≤ t ≤ νj,  Y 2  Z e−νj t 2s 2 2 e ξ1 (β(s)) + ξ2 (β(s)) (1 + µk )ψ(s) N (t) ≤ 2 0 0
18

Applying this inequality to (27) yields that for j ≥ j0 , 0 ≤ t ≤ νj, et ln ψ(t) − ln ψ(0) Z t  Y ξ 2 (β(s)) − ρ(β(s)) (1 + µk )ψ(s) ≤ es ln ψ(s) + r(β(s)) − 1 2 0 0
≤ c4 (et − 1) + σeνj ln j,

where c4 is a positive constant. As a result, if 0 < ν(j − 1) ≤ t ≤ νj and j ≥ j0 , then e−t ln ψ(0) c4 (1 − e−t ) σe−ν(j−1) eνj ln j ln ψ(t) ≤ + + . ln t ln t ln t ln(ν(j − 1))

That is to say lim sup t→+∞

ln ψ(t) ≤ σeν . ln t

Letting σ ↓ 1 and ν ↓ 0 results in lim sup t→+∞

ln ψ(t) ≤ 1. ln t

It then follows from Assumption 2 that P ln Ψ(t) ln ψ(t) 0 0, then with probability one, the decay 1

rate of tumor cells will not be faster than t− 2θ −ε , where θ > 0 is defined in (22). Theorem 7. Let Assumption 2 hold. If b > 0, then lim inf t→+∞

ln Ψ(t) 1 ≥− . ln t 2θ 19

(28)

Proof. It follows from (25) that there is a positive constant c5 such that  θ  E 1 + V1 (ψ(t)) ≤ c5 , t ≥ 0.

(29)

In addition, according to Itˆo’s formula, we obtain d[(1 + V1 (ψ))θ ]



ξ12 (β) λ(β) Q + −1 2 1 + Γ (β) 0
= 2θ(1 + V1 (ψ))

− r(β) +

where



    ξ12 (i) + ξ32 (i) 2 2 Υ1 = max − r(i) + + λ(i) + θ ξ1 (i) + ξ3 (i) , Υ2 = c2 max ρ(i) , i∈M i∈M 2     2 2 2 2 2 2 Υ3 = max − r(β) + 1.5(ξ1 (β) + ξ3 (β)) + λ(β) + (θ + 0.5)ξ2 (β)c2 , Υ4 = 1.5c2 max ξ2 (i) . i∈M

i∈M

Let c6 be sufficiently large such that

2θ(Υ1 V12 (ψ) + Υ2 V11.5 (ψ) + Υ3 V1 (ψ) + Υ2 V10.5 (ψ) + Υ4 ) ≤ c6 (1 + V1 (ψ))2 . Then we have θ

θ

θ−1

 ξ1 (β)V1 (ψ)dB1 (t)

d[(1 + V1 (ψ)) ] ≤ (1 + V1 (ψ)) c6 dt − 2θ(1 + V1 (ψ))  Y ξ3 (β)V1 (ψ) −1 Q + ξ2 (β) (1 + µk )ψ dB2 (t) + dB3 (t) . 1 + Γ−1 (β) 0
20

(30)

Let δ > 0 be sufficiently small such that q  q q 1 0.5 2 u 2 u 2 u c6 δ + 12θδ (ξ1 ) + c2 (ξ2 ) + (ξ3 ) < . 2

(31)

Let l = 1, 2.... In light of (30),   θ E sup (1 + V1 (ψ(t)) (l−1)δ≤t≤lδ Z t   θ   θ  1 + V1 (ψ(s)) ds ≤ E 1 + V1 (ψ((l − 1)δ)) + c6 E sup (l−1)δ≤t≤lδ Z t  θ−1 (l−1)δ  1 + V1 (ψ(s)) ξ1 (β(s))V1 (ψ(s))dB1 (s) + 2θE sup (l−1)δ≤t≤lδ

  ξ3 (β(s))V1 (ψ(s)) Q (1 + µk )ψ (s)dB2 (s) + dB3 (s) . −1 1 + Γ (β(s)) 0
Y

+ ξ2 (β(s))

0
(l−1)δ

−1

(32)

According to the Burkholder-Davis-Gundy inequality, one derives Z t   θ−1  E sup 1 + V1 (ψ(s)) ξ1 (β(s))V1 (ψ(s))dB1 (s) (l−1)δ≤t≤lδ (l−1)δ   Y ξ3 (β(s))V1 (ψ(s)) −1 . Q dB (s) + ξ2 (β(s)) (1 + µk )ψ (s)dB2 (s) + 3 −1 (β(s)) (1 + µ )ψ(s) 1 + Γ k 0
0
Z

lδ

(l−1)δ≤t≤lδ

In addition,  E sup

(l−1)δ≤t≤lδ

Z

  θ   Z lδ  θ  ds 1 + V1 (ψ(s)) ds ≤E 1 + V (ψ(s)) 1 (l−1)δ (l−1)δ   θ  ≤ δE sup 1 + V1 (ψ(t)) .

(33)

t

(l−1)δ≤t≤lδ

21

(34)

When (33) and (34) are used to (32) leads to    θ  θ E sup (1 + V1 (ψ(t)) ≤ E 1 + V1 ψ((l − 1)δ) (l−1)δ≤t≤lδ q    q q 0.5 2 u 2 u 2 u E sup + c6 δ + 12θδ (ξ1 ) + c2 (ξ2 ) + (ξ3 )

(l−1)δ≤t≤lδ

In light of (29) and (31), we get  E

θ

sup

(1 + V1 (ψ(t))

(l−1)δ≤t≤lδ



 θ  1 + V1 (ψ(t)) .

≤ 2c5 .

For ∀ε > 0, Chebyshev’s inequality indicates that   θ  1+ε P sup 1 + V1 (ψ(t) > (lδ) ≤ (l−1)δ≤t≤lδ

2c5 , l = 1, 2, ... (lδ)1+ε

According to the Borel-Cantelli lemma, for almost all ω ∈ Ω, there exists an integer l0 such that for ∀l ≥ l0 and (l − 1)δ ≤ t ≤ lδ, (1 + ε) ln(lδ) ln(1 + V1 (ψ(t))θ ≤ . ln t ln((l − 1)δ) As a result, lim sup t→+∞

ln(1 + V1 (ψ(t))θ ≤ 1 + ε. ln t

Letting ε → 0 gives lim sup t→+∞

ln(1 + V1 (ψ(t))θ ≤ 1. ln t

Hence lim sup t→+∞

ln(ψ −2θ (t)) ≤ 1, ln t

or lim inf t→+∞

ln(ψ(t)) 1 ≥− . ln t 2θ

In view of (13) and Assumption 2, (28) follows. To end of this section, we apply Theorems 2, 3, 4, 5, 6 and 7 to impulse free model (4), and hence obtain the following corollary. Corollary 1. For Eq. (4), (i) if Assumption 1 holds and ¯b < 0, then Ψ(t) in Eq. (4) will become extinct a.s.; (ii) If Assumption 1 holds and ¯b = 0, then Ψ(t) will be non-persistent in the mean a.s.; 22

(iii) If ¯b > 0, then Ψ(t) in Eq. (4) will be weakly persistent a.s.; (iv) If ¯b > 0, Ψ(t) in Eq. (4) will be stochastically permanent; (v) Ψ(t) in Eq. (4) obeys lim sup ln Ψ(t)/ ln t ≤ 1, a.s.; In addition, if ¯b > 0, then t→+∞

lim inf ln Ψ(t)/ ln t ≥ − t→+∞

1 , a.s.. 2θ

6. Discussions and simulations 6.1. Impacts of regime switching on the dynamics of model (5) Our results indicate that the persistence-extinction properties of model (5) are intimately associated with the random perturbations and impulsive perturbations. Let us first look at the impacts of regime switching. In order to understand these impacts more clearly, we focus on model (4) with M = {1, 2} and hypothesize that Assumption 1 holds. As a result, model (4) jumps between the following two subsystems:  dΨ(t) = Ψ(t) r(1) − ρ(1)Ψ(t) −

 λ(1) dt + ξ1 (1)Ψ(t)dB1 (t) 1 + Γ−1 (1)Ψ(t) ξ3 (1)Ψ(t) + ξ2 (1)Ψ2 (t)dB2 (t) + dB3 (t), 1 + Γ−1 (1)Ψ(t)

and

 dΨ(t) = Ψ(t) r(2) − ρ(2)Ψ(t) −

 λ(2) dt + ξ1 (1)Ψ(t)dB1 (t) 1 + Γ−1 (2)Ψ(t) ξ3 (2)Ψ(t) dB3 (t). + ξ2 (2)Ψ2 (t)dB2 (t) + 1 + Γ−1 (2)Ψ(t)

(35)

(36)

For models (35) and (36), there are two possibilities. • Models (35) and (36) have the same persistence-extinction behaviours. In this case, according to Corollary 1, hybrid model (4) has the same persistence-extinction behaviours. For example, both models (35) and (36) become extinct, then hybrid model (4) still becomes extinct. In other words, regime switching will not change the persistence-extinction behaviours in this case. • Models (35) and (36) have the different persistence-extinction behaviours. For example, model (35) becomes extinct and (36) is stochastic permanent. In this case, as a consequence of regime shift, Corollary 1 shows that hybrid model (4) has three possible persistenceextinction behaviours. If ¯b < 0, then the tumor cells will be eliminated; If ¯b > 0, then the tumor cells will be permanent in the tissue; If ¯b = 0, then the tumor cells will be 23

non-persistent in the mean. By the definition of ¯b one can find out that in this case, the persistence-extinction properties of model (5) are intimately associated with the Markov chain β(t): if β(t) stay longer in states with negative “stochastic growth rate” b(·), the tumor cells will be eliminated; on the contrary, the tumor cells will be permanent. On order to illustrate the above influences numerically, we need to give the values of the parameters. Different values of the growth rate r and the intra-specific competition rate ρ have been used in literature, for example, r = 0.18 and ρ = 3.6 × 10−10 [24], r = 0.51 and ρ = 5.14 × 10−10 [52, 53], r = 0.57 and ρ = 4.2 × 10−6 [12]. The first two values come from the

tumor growth in mice, the latter one comes from tumor cell cultivation in vitro. In this paper we choose r(1) = 0.18, ρ(1) = 3.6 × 10−10 and r(2) = 0.51, ρ(2) = 5.14 × 10−10 . The values of

λ and Γ are variable, relaying on the nature of immune response. In general, λ ∈ (10−2 , 10) and

Γ ∈ (0.1×109 , 5×109 ) [11, 54]. In this paper, we choose λ(1) = 0.05, λ(2) = 0.1, Γ(1) = 0.5×109 , Γ(2) = 109 . For the sake of convenience, the values of the parameters are listed in Table 1. Table 1. Parameter values Parameters

Units

Value

Source

r

day−1

r(1) = 0.18, r(2) = 0.51

[52, 53]

ρ

ml/(cells·day) ρ(1) = 3.6 × 10−10 , ρ(2) = 5.14 × 10−10

[52, 53]

λ

day−1

λ(1) = 0.05, λ(2) = 0.1

[11, 54]

Γ

cells/ml

Γ(1) = 0.5 × 109 , Γ(2) = 109

[11, 54]

Ψ(0)

cells/ml

106

[11, 54]

ξ1

day−1

ξ1 (1) = 0.7, ξ1 (2) = 0.5

Hypothesized

ξ2

ml/(cells·day)

ξ2 (1) = ξ2 (2) = 10−10

Hypothesized

ξ3

day−1

ξ3 (1) = ξ3 (2) = 0.2

Hypothesized

For parameters given in Table 1, it is easy to see that Assumption 1 holds and b(1) = −0.135, b(2) = 0.265. According to Corollary 1, model (35) becomes extinct (see Fig.1(a)) and model (36) is stochastic permanent (see Fig.1(b)). Now we let the stationary distribution of the Markov chain, π, vary. (i) Let π = (0.7, 0.3). Hence b = π1 b(1) + π2 b(2) = −0.015. It then follows from Corollary 1 that Ψ(t) in hybrid model (5) becomes extinct (see Fig.1(c)). 24

Figure 1: Hybrid model (4) with M = {1, 2} and parameters given in Table 1. (a) is a sample path of the solution of subsystem (35), which shows that subsystem (35) becomes extinct; (b) is a sample path of the solution of subsystem (36), which illustrates that subsystem (36) is stochastically permanent; (c) is a sample path of the solution of hybrid model (4) with π = (0.7, 0.3), which indicates that hybrid model (4) becomes extinct; (d) is a sample path of the solution of hybrid model (4) with π = (0.3, 0.7), which indicates that hybrid model (4) is stochastically permanent.

Figure 2: A sample path of the solution of Eq. (36) with ξ1 (2) = 0.9, the values of other parameters are the same with those in Fig.1(b).

(ii) Let π = (0.3, 0.7). Hence b = π1 b(1) + π2 b(2) = 0.145. By virtue of Corollary 1, Ψ(t) in hybrid model (5) is stochastic permanent (see Fig.1(d)).

6.2. Impacts of the white noise on the dynamics of model (5) Now we see the impacts of the white noise on the persistence-extinction properties of model (5). As said above, under Assumption 1, bb is the threshold between extinction and weak persistence

of Ψ(t). Notice that

bb = lim sup t−1 t→+∞

 X

0


ln(1 + µk ) +

X i∈M



πi r(i) − λ(i) −



ξ12 (i)

+

  /2 .

ξ32 (i)

One can see that bb decreases with the increase of ξ12 (·) or ξ32 (·), in addition, bb has nothing to do

with ξ22 (·). As a consequence, both the white noise perturbation on the growth rate (i.e., ξ12 ) and the white noise perturbation on the destruction rate (i.e., ξ32 ) are conducive to the elimination of neoplasm, however, the persistence-extinction properties of model (5) have nothing to do with the white noise perturbation on the intra-specific competition rate (i.e., ξ22 ). Let us numerically illustrate the impacts of ξ12 (the impacts of ξ22 and ξ32 can be illustrated in the same way and hence are omitted). In Fig.2, we let the values of the parameters be the same with those in Fig.1(b) except ξ1 (2) = 0.9. Comparing Fig.1(b) with Fig.2, one can observe that as ξ12 increase, tumor cells tend to become extinct. 6.3. Impacts of impulsive perturbations on the dynamics of model (5) Finally, we explore the impacts of impulsive perturbations on the dynamical properties of model (5). 25

Figure 3: A sample path of the solution of Eq. (37) with parameters given in Table 1 and tk = 10k, µk = e2 − 1,

k ∈ N.

(I) On the one hand, one can observe that if Assumption 2 holds (i.e., the impulse is bounded), then bb = ¯b. Hence Theorems 2, 3, 4, 5, 6 and 7 indicate that if the impulse is bounded,

then the impulse does not affect the persistence, extinction, stochastic permanence, the

lower- and the upper-growth rates of the solutions. Hence, limited frequency of treatments and metastasis could not change the persistence-extinction behaviours of neoplasm. (II) One the other hand, if the impulse is unbounded, then Theorems 2, 3 and 4 imply that the persistence-extinction properties of model (5) rely on the values of tk and µk , k ∈ N. Therefore, the impulse can change the persistence-extinction properties of model (5). More precisely, positive impulse (for example, the metastasis of the tumor cells from other tissues to the tissue under consideration) could make the neoplasm become persistent (see Example 1 below), negative impulse (for example, the treatments) could make the neoplasm become extinct (see Example 2 below). We are now in the position to numerically illustrate the above influences. Here we ignore the regime switching and only present the figures of case (II) to see the influences of the impulse. Example 1. For model (35) with parameters given in Table 1, Fig.1a shows that the neoplasm represented by model (35) becomes extinct. Now we take the impulsive perturbations into account, then model (35) becomes    λ(1)   dΨ(t) = Ψ(t) r(1) − ρ(1)Ψ(t) − dt + ξ1 (1)Ψ(t)dB1 (t)   1 + Γ−1 (1)Ψ(t)    ξ3 (1)Ψ(t) + ξ2 (1)Ψ2 (t)dB2 (t) + dB3 (t), t 6= tk , k ∈ N. 1 + Γ−1 (1)Ψ(t)        Ψ(t+ ) = (1 + µ )Ψ(t ), k ∈ N, k

k

(37)

k

For parameters given in Table 1 and tk = 10k, µk = e2 − 1, k ∈ N, it is easy to obtain that

bb = 0.065 > 0. According to Theorem 4, the neoplasm represented by the impulsive model (37) is weakly persistent (see Fig.3).

26

Figure 4: A sample path of the solution of Eq. (38) with parameters given in Table 1 and tk = 10k, µk = e−3 − 1,

k ∈ N.

Figure 5: Solution of Eq. (36) with r(2) = 0.75 and ξ32 (2) = 0.84, the values of other parameters are the same with those in Fig.1(b). Clearly, (40) holds. (a) is with Ψ(0) = 102 , which indicates that Ψ(t) becomes extinct; (a) is with Ψ(0) = 106 , which shows that Ψ(t) is weakly persistent.

Example 2. For model (36) with parameters given in Table 1, Fig.1b indicates that the neoplasm represented by model (36) is weakly persistent. When the impulsive perturbations are taken into account, model (36) becomes    λ(2)   dΨ(t) = Ψ(t) r(2) − ρ(2)Ψ(t) − dt + ξ1 (2)Ψ(t)dB1 (t)   1 + Γ−1 (2)Ψ(t)    ξ3 (2)Ψ(t) + ξ2 (2)Ψ2 (t)dB2 (t) + dB3 (t), t 6= tk , k ∈ N. 1 + Γ−1 (2)Ψ(t)        Ψ(t+ ) = (1 + µ )Ψ(t ), k ∈ N, k

k

(38)

k

For parameters given in Table 1 and tk = 10k, µk = e−3 − 1, k ∈ N, we have bb = −0.035 < 0. In view of Theorem 2, the neoplasm represented by the impulsive model (38) become extinct (see Fig.4).

7. Conclusions Understanding the inherent law of neoplasm evolution is useful for the prevention and treatment of neoplasm [3, 4]. In this paper, by taking the white noise, regime switching and impulsive perturbations into account, we proposed a stochastic hybrid model to describe the evolution of neoplasm. Especially, our model considered both the treatment of neoplasm and the metastasis of tumor cells. The main theoretical results are six theorems. Under Assumption 1, Theorems 2, 3 and 4 establish the threshold between weak persistence and extinction of tumor cells; Under the assumption that the impulse is bounded (i.e., Assumption 2), Theorems 5 give sufficient conditions for stochastic permanence of tumor cells while Theorems 6 and 7 estimate the lowerand the upper-growth rates of tumor cells respectively. These theoretical results indicate that the persistence-extinction properties of the model are intimately associated with the random perturbations and impulsive perturbations: 27

I if the Markovian chain stays longer in states with negative “stochastic growth rate”, the tumor cells will be eliminated; on the contrary, the tumor cells will be permanent (see Fig.1); I The white noise perturbations on the growth rate and the destruction rate are conducive to the elimination of neoplasm (compare Fig.1(b) with Fig.2), however, the white noise perturbation on the the intra-specific competition rate does not affect the persistenceextinction properties of tumor cells; I Positive impulse could make the neoplasm become persistent (compare Fig.1(a) with Fig.3), negative impulse could make the neoplasm become extinct (compare Fig.1(b) with Fig.4). Therefore, in order to control tumors, some strategies are as follows. • To enhance the white noise perturbations on the growth rate and the destruction rate; • To make the Markovian chain stay longer in states with negative “stochastic growth rate”. A feasible way is to keep the immune system strong and healthy, for example, eat appropriate fruits and vegetables, exercise regularly, get adequate sleep, try to minimize stress [55, 56]; • To increase the negative impulse and to reduce the positive impulse. Some feasible strategies are: to strengthen doses or frequencies of treatments, to inhibit the metastasis of tumor cells (a recent study [21] showed that KBU2046, a chemical, can significantly decrease the metastasis of human prostate, breast, colon, and lung cancer cells by up to 92%). However, the white noise perturbation on the the intra-specific competition rate and limited frequency of treatments and metastasis do not affect the persistence-extinction behaviours of neoplasm. Some questions deserve further consideration. Firstly, compared with model (5), the following model is more reasonable:  dΨ(t) = Ψ(t) r(β(t)) − ρ(β(t))Ψ(t) −

 λ(β(t)) dt + ξ1 (β(t))Ψ(t)dB1 (t) 1 + Γ−1 (β(t))Ψ(t) X ξ3 (β(t))Ψ(t) + ξ2 (β(t))Ψ2 (t)dB2 (t) + dB (t) + Ψ(t) µk δtk (dt), 3 1 + Γ−1 (β(t))Ψ(t) k 28

(39)

where δtk is the Dirac measure concentrated at tk . In this paper, we studied the properties of model (5) with the help of impulse-free model (10) and (13). However, how to analyze model (39) is still unknown (the approaches in this paper are invalid, one needs to develop some new approaches). Secondly, from Theorems 2 and 4, an interesting issue is that what happens if Assumption 1 does not hold? In fact, similar to the proof of Theorem 2, one can prove that if  X  X ξ 2 (i)
−1 ln(1 + µk ) + lim sup t πi r(i) − 1 < 0, 2 t→+∞ 0
then Ψ(t) will become extinct a.s.. Fig.5 indicates that if  X  X  ξ 2 (i)
2 −1 ln(1 + µk ) + πi r(i) − 1 min Γ(i)ρ(i) − λ(i) − ξ3 (i) ≤ 0, lim sup t > 0 > bb, i∈M 2 t→+∞ 0
(40)

the persistence and extinction of Ψ(t) depend on the initial value. However, the theoretical analysis is still known. In addition, it is also interesting to consider the effects of reaction diffusion (see, e.g., [57]) and color noises (see, e.g., [58]). We leave these problems for further research. Declarations of interest None. Acknowledgements The authors thank the editors and the referees for their careful reading and valuable comments. This work was supported by the National Natural Science Foundation of P.R. China (Nos. 11771174), Natural Science Foundation of Jiangsu Province (No. BK20170067), “333 High-level Personnel Training Project” and Qinglan Project of Jiangsu Province, Young Talent Lifting Project of Jiangsu Province Science and Technology Association. REFERENCES References [1] F. Bray, J. Ferlay, I. Soerjomataram, et al., Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA: Cancer J. Clin. 68 (2018) 394-424. 29

[2] World Health Organization, World Health Statistics 2016: Monitoring Health for the SDGs Sustainable Development Goals, World Health Organization, 2016. [3] D. Li, F. Cheng, Threshold for extinction and survival in stochastic tumor immune system, Commun. Nonlinear Sci. Numer. Simulat. 51 (2017) 1-12. [4] T.E. Wheldon, Mathematical Models in Cancer Research, Hilger Publishing, BostonPhiladelphia, 1988. [5] A. Bertuzzi, A. d’Onofrio, A. Fasano, A. Gandolfi, Regression and regrowth of tumour cords following single-dose anticancer treatment, Bull. Math. Biol. 65 (2003) 903-931. [6] A. d’Onofrio, A. Gandolfi, Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by Hahnfeldt et al. (1999), Math. Biosci. 191 (2004) 159-184. [7] A. d’Onofrio, A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences, Physica D 208 (2005) 220-235. [8] M. Maruˇsi´c, S. Vuk-Pavlovic, J.P. Freyer, Tumor growth in vivo and as multicellular spheroids compared by mathematical models, Bull. Math. Biol. 56 (1994) 617-631. [9] M.H. Gee , A. Han , S.M. Lofgren, et al., Antigen identification for orphan T cell receptors expressed on tumor-infiltrating lymphocytes, Cell 172 (2018) 549-563. [10] G. Sulli, A. Rommel, X. Wang, et al. Pharmacological activation of REV-ERBs is lethal in cancer and oncogene-induced senescence, Nature 553 (2018) 351-355. [11] R. Lefever, W. Horsthemke, Bistability in fluctuating environments. implications in tumor immunology, Bull. Math. Biol. 41 (1979) 469-90. [12] T. Bose, S. Trimper, Stochastic model for tumor growth with immunization, Phys. Rev. E 79 (2009) 051903. [13] L. Jiang, X. Luo, D. Wu, S. Zhu, Stochastic properties of tumor growth with coupling between non-Gaussian and Gaussian noise terms, Chin. Phys. B 21 (2012) 090503.

30

[14] W. Zhong, Y. Shao, Z. He, Pure multiplicative stochastic resonance of a theoretical antitumor model with seasonal modulability, Phys. Rev. E 74 (2006) 011916. [15] M. Liu, M. Deng, Permanence and extinction of a stochastic hybrid model for tumor growth, Appl. Math. Lett. 94 (2019) 66-72. [16] D. Li, F. Cheng, The extinction and persistence of tumor evolution influenced by external fluctuations and periodic treatment, Qual. Theory Dyn. Syst. (2019), doi: 10.1007/s12346019-00317-9. [17] J. Yang, Y. Tan, R.A. Cheke, Thresholds for extinction and proliferation in a stochastic tumour-immune model with pulsed comprehensive therapy, Commun. Nonlinear Sci. Numer. Simulat. 73 (2019) 363-378. [18] A. Fiasconaro, A. Ochab-Marcinek, B. Spagnolo, E. Gudowska-Nowak, Monitoring noiseresonant effects in cancer growth influenced by external fluctuations and periodic treatment, Eur. Phys. J. B 65 (2008) 435-442. [19] A. d’Onofrio, Bounded-noise-induced transitions in a tumor-immune system interplay, Phys. Rev. E 81 (2010) 021923-7. [20] M. Vanneman, G. Dranoff, Combining immunotherapy and targeted therapies in cancer treatment, Nat. Rev. Cancer 22 (2012) 237-251. [21] L. Xu, R. Gordon, R. Farmer, et al. Precision therapeutic targeting of human cancer cell motility. Nature Communications 9 (2018) 2454. [22] R. Paduch, The role of lymphangiogenesis and angiogenesis in tumor metastasis, Cellular Oncology 39 (2016) 397-410. [23] D.S.P. Tan, S.B. Kaye, Transcoelomic Metastasis. In: Schwab M. (eds) Encyclopedia of Cancer. Springer, Berlin, Heidelberg, 2014. [24] V.A. Kuznetsov, I.A. Makalkin, M.A. Taylor, A.S. Perelson, Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis, Bull. Math. Biol. 56 (1994) 295-321.

31

[25] R.M. May, Stability and Complexity in Model Ecosystems, Princeton Univ. Press, 1973. [26] J.R. Beddington, R.M. May, Harvesting natural populations in a randomly fluctuating environment, Science 197 (1977) 463-465. [27] X. Li, A. Gray, D. Jiang, X. Mao, Sufficient and necessary conditions of stochastic permanence and extinction for stochastic logistic populations under regime switching, J. Math. Anal. Appl. 376 (2011) 11-28. [28] M. Liu, Dynamics and simulations of a logistic model with impulsive perturbations in a random environment, Math. Comput. Simulation 92 (2013) 53-75. [29] R. Wu, K. Wang, Asymptotic properties of a stochastic Lotka-Volterra cooperative system with impulsive perturbations, Nonlinear Dynam. 77 (2014) 807-817. [30] R. Wu, X. Zou, K. Wang, Asymptotic behavior of a stochastic non-autonomous predatorprey model with impulsive perturbations, Commun. Nonlinear Sci. Numer. Simulat. 20 (2015) 965-974. [31] S. Zhang, D. Tan, Dynamics of a stochastic predator-prey system in a polluted environment with pulse toxicant input and impulsive perturbations, Appl. Math. Model. 39 (2015) 63196331. [32] R. Tan, Z. Liu, S. Guo, H. Xiang, On a nonautonomous competitive system subject to stochastic and impulsive perturbations, Appl. Math. Comput. 256 (2015) 702-714. [33] C. Lu, X. Ding, Persistence and extinction of an impulsive stochastic logistic model with infinite delay, Osaka J. Math. 53 (2016) 1-31. [34] S. Zhang, X. Meng, T. Feng, T. Zhang, Dynamics analysis and numerical simulations of a stochastic non-autonomous predator-prey system with impulsive effects, Nonlinear Anal. Hybrid Syst. 26 (2017) 19-37. [35] A. Das, G.P. Samanta, Stochastic prey-predator model with additional food for predator, Physica A 512 (2018) 121-141.

32

[36] A. Das, G.P. Samanta, Modeling the fear effect on a stochastic prey-predator system with additional food for the predator, J. Phys. A 51 (2018), 465601. [37] C. Lu, X. Ding, Periodic solutions and stationary distribution for a stochastic predator-prey system with impulsive perturbations, Appl. Math. Comput. 350 (2019) 313-322. [38] D. Belotti, V. Vergani, T. Drudis, et al., The microtubule-affecting drug paclitaxel has antiangiogenic activity, Clin. Cancer Res. 2 (1996) 1843-1849. [39] X. Li, D. Jiang, X. Mao, Population dynamical behavior of Lotka-Volterra system under regime switching, J. Comput. Appl. Math. 232 (2009) 427-448. [40] C. Zhu, G. Yin, On hybrid competitive Lotka-Volterra ecosystems, Nonlinear Anal. 71 (2009) e1370-e1379. [41] M. Liu, M. Deng, Analysis of a stochastic hybrid population model with Allee effect, Appl. Math. Comput. 364 (2020) 124582. [42] V. Lakshmikantham, D.D. Bainov, P.S. Simeonov, Theory of Impulsive Differential Equations, World Scientific, London, 1989. [43] G. Ballinger, X. Liu, Permanence of population growth models with impulsive effects, Math. Comput. Modelling 26 (1997) 59-72. [44] S. Ahmad, G.T. Stamov, Almost periodic solutions of N-dimensional impulsive competitive systems, Nonlinear Anal. Real World Appl. 10 (2009) 1846-1853. [45] I. Stamova, G.T. Stamov, Applied Impulsive Mathematical Models, Springer, New York, 2016. [46] M. He, F. Chen, Extinction and stability of an impulsive system with pure delays, Appl. Math. Lett. 91 (2019) 128-136. [47] G.P. Samanta, S.P. Bera, Analysis of a Chlamydia epidemic model with pulse vaccination strategy in a random environment, Nonlinear Anal. Model. Control 23 (2018) 457-474. [48] X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, 2006. 33

[49] J. Bao, C. Yuan, Stochastic population dynamics driven by L´evy noise, J. Math. Anal. Appl. 391 (2012) 363-375. [50] H.L. Smith, H.R. Thieme, Dynamical Systems and Population Persistence, American Mathematical Society, NY, 2011. [51] R.Z. Khasminskii, C. Zhu, G. Yin, Stability of regime-switching diffusions, Stochastic Process. Appl. 117 (2007) 1037-1051. [52] A. Diefenbach, E. Jensen, A. Jamieson, D. Raulet, Rae1 and H60 ligands of the NKG2D receptor stimulate tumor immunity, Nature 413 (2001) 165-171. [53] L.G. DePillis, A.E. Radunskaya, C.L. Wiseman, A validated mathematical model of cellmediated immune response to tumor growth, Cancer Res. 65 (2005) 7950-7958. [54] R.P. Garay, R. Lefever, A kinetic approach to the immunology of cancer: stationary states properties of effector-target cell reactions, J. Theor. Biol. 73 (1978) 417-438. [55] G.J.V. Nossal, Life, death and the immune system. Scientific American 269 (1993) 52-62. [56] E.M.V. Reiche, S.O.V. Nunes, H.K. Morimoto, Stress, depression, the immune system, and cancer. The Lancet Oncology 5 (2004) 617-625. [57] S.C. Ferreira Jr, M.L. Martins, M.J. Vilela, Reaction-diffusion model for the growth of avascular tumor. Phys. Rev. E 65 (2002) 021907. [58] G.P. Samanta, A stochastic two species competition model: nonequilibrium fluctuation and stability, Int. J. Stoch. Anal. 2011 (2011) 489386.

34

ig1a

5

14

x 10

Ψ(t) 12

Tumor cell number

10

8

6

4

2

0

0

5

10

15 Day

20

25

30

ig1b

8

10

x 10

Ψ(t) 9 8

Tumor cell number

7 6 5 4 3 2 1 0

0

100

200

300

400

500 Day

600

700

800

900

ig1c

8

7

x 10

Ψ(t) 6

Tumor cell number

5

4

3

2

1

0

0

500

1000

1500

2000 Day

2500

3000

3500

4000

ig1d

8

6

x 10

Ψ(t)

Tumor cell number

5

4

3

2

1

0

0

500

1000

1500 Day

2000

2500

3000

ig2

6

3.5

x 10

Ψ(t) 3

Tumor cell number

2.5

2

1.5

1

0.5

0

0

50

100 Day

150

200

ig3

6

12

x 10

Ψ(t)

Tumor cell number

10

8

6

4

2

0

0

100

200

300 Day

400

500

ig4

6

11

x 10

Ψ(t)

10 9

Tumor cell number

8 7 6 5 4 3 2 1 0

0

20

40

60 Day

80

100

ig5a

250 Ψ(t)

Tumor cell number

200

150

100

50

0

0

50

100

150 Day

200

250

300

ig5b

8

4

x 10

Ψ(t) 3.5

Tumor cell number

3

2.5

2

1.5

1

0.5

0

0

200

400

600

800

1000 Day

1200

1400

1600

1800

2000