DYNAMICAL BEHAVIOR OF A NEW ONCOLYTIC VIROTHERAPY MODEL BASED ON GENE VARIATION

. Oncolytic virotherapy is an experimental treatment of cancer patients. This method is based on the administration of replication-competent viruses that selectively destroy tumor cells but remain healthy tissue unaffected. In order to obtain optimal dosage for complete tumor eradication, we derive and analyze a new oncolytic virotherapy model with a ﬁxed time period τ and non-local infection which is caused by the diﬀusion of the target cells in a continuous bounded domain, where τ is assumed to be the duration that oncolytic viruses spend to destroy the target cells and to release new viruses since they enter into the target cells. This model is a delayed reaction diffusion system with nonlocal reaction term. By analyzing the global stability of tumor cell eradication equilibrium, we give diﬀerent treatment strategies for cancer therapy according to the diﬀerent genes mutations (oncogene and antioncogene).


1.
Introduction. In recent years, oncolytic virotherapy has got increasing attention since Martuza, et al [12] published an article in Science in which they claim that Transgene HSV plays an important role in malignant glioma therapy. Replicating oncolytic viruses is able to infect and lyse cancer cells, and spread through the tumor, while leaving normal cells largely unharmed. This makes them potentially useful in cancer therapy. A variety of viruses have shown promising results in clinical trials [9]. In 2014, Lin, et.al [10] identified a naturally occurring alphavirus (M1) as a novel selective killer targeting zinc-finger antiviral protein (ZAP)-deficient cancer cells. While some encouraging results were reported, oncolytic virotherapy never became a standard treatment. Thus many mathematical models are formulated by medical scientists and mathematicians. These mathematical models of oncolytic virus have provided useful insight into the understanding of underlying phenomena and offered helpful guidance to related public health decisions.
In order to get complete understanding of how virus and host characteristics influence the outcome of therapy, Wodarz [25] formulated an ODE model to describe the development of tumor cells where a virus-specific lytic CTL response is considered. In 2003, Wodarz [24] further improved his work, by introducing a free virus population and removing the growth term for infected cells. In our recent work [23], the minimum effective dosage of medication is explicitly found by means of an ODE model in which the competition between tumor cells and normal cells is considered to be purely exploitative, in the sense that the organisms simply consume the nutrient, thereby making it unavailable for a competitor. There are many other related models established in these years, see for example [2,5,9,18,22]. In most cases, logistic growth is employed to describe the growth of cells.
Motivated by previous researches, we use Lotka-Volterra competition model to describe relationship between normal cells and tumors cells, and give the proper biological interpretations.
Let u 1 denote normal cells, u 2 the tumor cells in system (1). Where a 1 (a 2 ) denotes the replication rate of normal (tumor) cells, b 1 (b 2 ) indicates that the normal cells have negative feedback on the growth of normal cells (tumor cells). Similarly, c 1 (c 2 ) indicates that the tumor cells have negative feedback on the growth of normal cells (tumor cells).
As we know, the dynamic behavior is described as Picture I-Picture IV. It is consistent with our intuition to name equilibrium E 0 = ( a1 b1 , 0), E 1 = (0, a2 c2 ) and E 3 = ( a2c1−a1c2 c1b2−b1c2 , a1b2−a2b1 c1b2−b1c2 ) as tumor eradication equilibrium, tumor cell equilibrium, and coexistence equilibrium respectively. Now let us introduce the concepts of oncogene and antioncogene. Oncogene is a gene that has the potential to cause cancer [3]. Most normal cells will suffer a programmed form of rapid cell death (apoptosis) when critical functions are altered. Activated oncogenes can cause those cells designated for apoptosis to survive and proliferate instead. By contrast, antioncogene is a gene that protects a cell from one step on the path to cancer. In a word, activated oncogenes will increase replication rate of tumor cells, and the loss of antioncogene will decrease the negative feedback on tumor cells. Now we denote the B i = ai bi , C i = ai ci , i = 1, 2. In particular, tumor replication rate a 2 times its carrying capacity c −1 2 is equal to C 2 = a2 c2 . Dividing tumor replication rate a 2 by the negative feedback effect of normal cells b 2 leads to B 2 = a2 b2 . Together with the concepts of antioncogene (oncogene), antioncogene is activated if B 2 is small enough (oncogene gene mutates if C 2 is large enough).
Together with Picture I-Picture IV,we see that (Picture I) C 1 ≤ C 2 and B 1 < B 2 indicates that both oncogene and antioncogene mutate which leads to the tumor cells increasing. (Picture II) C 1 > C 2 and B 2 ≥ B 1 implies that neither oncogene nor antioncogene mutates.
(Picture III) C 1 > C 2 , B 1 < B 2 shows that only antioncogene mutates which lead to the increase of tumor cells. (Picture IV) C 1 < C 2 and B 1 > B 2 implies that oncogene mutates but the antioncogene doesn't mutate.
We point out that E 4 = (0, 0) is always an unstable equilibrium. Here, we will regard case (III) and case (IV) as the early stage of cancer, because in such cases, there is only one gene mutates. However, in case (I), both oncogene and antioncogene mutate. Therefore, we regard this case as the later period of cancer.
Picture IV Actually, as early as 1996, Gatenby and Gawlinski [8] studied the model of cancer invasion which is similar to (1) except for the diffusion term. The diffusion phenomena of tumor cells are studied by many other researches, see. e.g. [1,7,13,21].
The system under consideration is as follows: where d 1 ≥ 0, d 2 > 0. Note that some researchers pointed out that the normal cells are well regulated and organized normally in an organ and therefore, there will be no spatial diffusions [8]. However, others think that the normal cells can move in the organ [14,21]. In view of these two situations, we assume that d 1 ≥ 0. Especially, in the case of d 1 = 0, it is assumed that no boundary conditions are specified.
In order to analyze the effect of oncolytic therapy, we introduce the oncolytic virus biomass density u 3 (t, x) into system (2). Let τ be the duration that viruses spend to destroy the target cells and to release new virus since they enter into the target cells. It is necessary to take account of movements of oncolytic virus during the replication process. Similar to [17,19], let P (t, a, x) be the generalised virus biomass under consideration at time t ≥ 0, replication age a ≥ 0 and location x ∈ Ω, and assume it satisfies ( ∂ ∂t + ∂ ∂a )P (t, a, x) = d 3 P (t, a, x) − dP (t, a, x), where µu 2 (t, x)u 3 (t, x) is the interaction rate which can be interpreted as a massaction law, analogous to mass-action principles in chemistry. Applying the method of characteristic curve and separation variables to (2), we have where Γ is the Green's function of operator d 3 −dI associated with certain boundary conditions. Then the oncolytic virus biomass is given by After a change of variables, u 3 (t, x) becomes Differentiation of above formulae leads to More generally, we can also regard Γ(x, y, τ ) as a probability distribution function which expresses the probability that an individual in location y at time t − τ moves to location x at time t.
Since main goal of our research is to study the effect of oncolytic therapy, we will consider the optimal coagulant dosage B in our model. As is formulated as follows where d is the mortality rate of virus, d i is diffusion coefficient (i = 1, 2, 3), and is the Laplacian operator. The term µu 2 (t, x)u 3 (t, x) in the second equation of system (3) reflects viral cytotoxicity causing tumor cells being killed. Ω is a relatively isolated bounded region in R n with smooth boundary ∂Ω. The boundary condition means that there is no individual across ∂Ω and therefore ∂Ω serves as a natural barrier to dispersal.
The rest of this paper is organized as follows. In the next section, we show that for any given initial values, there exists unique classical solution to (3), and the solution semiflow generated by (3) is point dissipative. In section 3, we study the stability of the equilibria of model (3) corresponding to tumor eradication. In the last section, we combine the mathematical results with biological background to give different treatment strategies for different situations.
2. Global existence and point dissipativity of solution semiflow of (3). Let X = C(Ω, R 3 ) be a Banach space with supremum norm · X and X + = C(Ω, R 3 + ), where Ω denotes the closure of set Ω. Then X + is a closed cone of X, which induces a partial ordering on X. Thus (X, · X ) is a Banach lattice.
First, we consider the linear system for i = 1, 2, 3, for φ = (φ 1 , φ 2 , φ 3 ) and x ∈ Ω. By using standard arguments [15], the system (3) can be rewritten as the following semilinear integral equation of the form Proof. By Martin and Smith [11], it suffices to show that f defined by (5) satisfies the Lipschitzian and subtangential conditions.
(i) For each R > 0 and any Let L 1 = (a 1 +2b 1 R+2c 1 R). By the similar method, it is easy to see that there exist constants L 2 , L 3 > 0 such that
Theorem 2.2. There exist constants U 1 , U 2 and U 3 > 0, and a positively invariant set Λ for system (3), where In order to prove Theorem 2.2, we need the following three preliminary lemmas. 0) on Ω, and either u(x, t) > u(x, t) or on ∂Ω × (0, T ], then either u ≡ u or u > u on Ω × (0, T ]. Lemma 2.4. The scalar reaction-diffusion equation bU (t, x)), x ∈ Ω, t > 0, where a, b, d > 0. Then system (6) has a unique positive equilibrium a b which is globally stable in R + .
Proof. Let U + be the unique solution to the ordinary differential equation subject to the initial condition U (0) = max x∈Ω U 0 (x). We set U − be the solution of equation (7) with initial condition U (0) = min x∈Ω U 0 (x). It is easy to see that x ∈ ∂Ω, t > 0, By Lemma 2.3, it is easy to see that U (t, x) ≤ U + (t). By a similar argument, we have U − (t) ≤ U (t, x). Fortunately, system (7) is well-known logistic model which has a globally stable equilibrium a b . It follows that system (6) has a unique positive equilibrium a b which is globally stable in R + .
Lemma 2.5. For the reaction-diffusion predator-prey model with time-delayed of the form where a, b, d 1 , d 2 , c, B, τ > 0, its solution semiflow Φ(t) : Proof. Using the Lemma 2.3 and Lemma 2.4, one proves that 0 ≤ Z(t, x) ≤ Z * , where we can choose Z * > a b as t → ∞. Thus, by making use of the boundedness of Z(t, x) and the property of Γ in the second equation in system (8), we see that where k = max x∈Ω,y∈Ω Γ(x, y, τ )Z * , V (t) = Ω V (x, t)dx. Now, we show that V (t) is bounded. By the Neumann boundary condition and the divergence theorem, we integrate the first equation in system (8), we get Similarly, integrating the second equation of (8) with respect to x ∈ Ω and making use of (10), we get where k 1 = mes(Ω) max x∈Ω,y∈Ω cΓ(x, y, τ ). Thus, By Theorem 2.1 and Lemma 2.4, there exists a positive number K * such that which implies that V (t) ≤ V * for some V * > 0. Now combining this with system (9) and Lemma 2.4, we conclude that there exists a positive number V * such that V (t, x) ≤ V * for t → ∞.
We expect to find conditions on the optimal coagulant dosage B such that the tumor cells will be extinct. Hence, we set u 2 = 0 in (3), the system (3) becomes It is easy to see that system (11) is uncoupled. From Lemma 2.4, it follows that system (11) admits a unique positive equilibrium ( a1 b1 , B d ) which is globally asymptotically stable. Linearizing system (3) at the tumor eradication equilibrium E = ( a1 b1 , 0, B d ), we get the following system for the tumor component u 2 : Consider the linear operator Suppose ϕ(x) is an eigenfunction of J corresponding to an eigenvalue σ. Then Since {φ n } is the completely orthogonal basis in L 2 (Ω), we can express ϕ as ϕ = n≥1 c n φ n .
Substituting this expression into (12), Thus, the eigenvalues of J are simply given by σ (n) = −d 2 λ n + a 2 − a1b2 b1 − µB d , n ≥ 1. Therefore we have the following results: , then the tumor eradication equilibrium E is locally asymptotically stable; If B < d µ (a 2 − a1b2 b1 ), then the tumor eradication equilibrium E is unstable.
Unfortunately, the condition B > d u (a 2 − a1b2 b1 ) does not suffice to ensure the global stability of the equilibrium E, see Figure 1. Actually, we can also conclude that the condition B ≥ d u (a 2 − a1b2 b1 ) is sufficient for the local stability of equilibrium E. Next, we discuss the global stability of the equilibrium E.
By Theorem 2.1, we have u i (t, x) ≥ 0 for every x ∈ Ω, t > 0, i = 1, 2, 3. Using the property of Γ together with the system (3), we have By Lemma 2.3 and Lemma 2.4, we know that for arbitrary ε > 0, there exists constant t 1 > 0 such that Theorem 3.2. Let (u 1 (t, x), u 2 (t, x), u 3 (t, x)) be the unique solution to system (3). Let (v 1 (t, x), v 2 (t, x)) be the unique solution to the diffusion equation: Then To prove Theorem 3.2, we need the following lemma.
Remark 1. This result is also true for some specified d j = 0, j = 1, 2, · · · , m which is a special case of Proposition 1 in [11] without time delay. Now, we give the proof of Theorem 3.2 Proof. Let u 2 = k − u 2 ≥ 0 for some constant k > a2 b2 . Then we have We denote It is easy to see that ∂g i /∂u j ≥ 0, for i = j. We also consider system with ∂ n u 1 = ∂ n u 2 = 0, t > t 1 , x ∈ ∂Ω and u 1 ( (15) can be rewritten as system (13). Hence, In order to get the condition on the global stability of equilibrium E, we consider the following system of ordinary differential equations . Furthermore, we have main result of this section as follows.
Proof. As we know, if Consequently, positive solutions of system (16) will tend to the equilibrium ( a1 b1 , 0). Besides, since ε is an arbitrary number, we see that equilibrium ( a1 b1 , 0) of system On the other hand, u 2 (t, x) ≤ v 2 (t) implies that u 2 (t, x) → 0 as t → ∞. Now, we consider the following system: System (3) is asymptotically autonomous -with limit equation (17). In other words, the semiflow Φ generated by system (3) is asymptotically autonomous with limit system Θ generated by system (17) (see [20]). As system (17)  Further the ω-Θ-limit set of any pre-compact Θ-orbit contains equilibrium ( a1 b1 , B d ). So the condition (E) in [20] is satisfied. By globally stable equilibrium ( a1 b1 , B d ), we also know that there is no Θ-cyclical chain of Θ-equilibria. Thus, E = ( a1 b1 , 0, B d ) is also globally stable in system (3)   In this paper, we derive a nonlocal delayed and diffusive oncolytic virotherapy model with bounded spatial domain. From the point of view of genetic variation, we use the Lotka-Volterra competition model to describe relationship between normal cells and tumors cells. When oncogenes mutate and express at high level or both oncogenes and antioncogenes mutate, the evolution of normal cells and tumor cells are depicted as Picture I. Picture II describes the evolution of cells in healthy tissue. When antioncogenes mutate to cause a loss or reduction in its function but oncogenes do not mutate, the relationship between tumor cells and normal cells are described as Picture III. For the last case, Picture IV is corresponding to the case that oncogenes are activated but the antioncogenes do not mutate.
Using the method in [17,19], we analyze the relationship among normal cells, tumor cells and oncolytic virus, and formulate a new oncolytic virotherapy model as (3). In order to find the correct treatment strategy, we focus on the stability of the tumor eradication equilibrium E. Then we find different treatment strategies for the tumors corresponding to different gene mutation (oncogene and antioncogene). Specifically, Theorem 3.1 states that the tumor eradication equilibrium E is locally stable with the treatment strategy B > d u (a 2 − a1b2 b1 ). Unfortunately, the condition B > d u (a 2 − a1b2 b1 ) can not guarantee the global stability of the tumor eradication equilibrium E, see Figure 1. In Figure 1, we see that B > d µ (a 2 − a1b2 b1 ) = 0, but the solution of system (3) with the initial values u 1 (0, x) = 1, u 2 (0, x) = 3, u 3 (0, x) = 0, x ∈ Ω will tend to the point (0, 0.5528, 0.4472) as t → ∞. Thus, such a dosage B may not lead to complete eradication of tumor cells.
Hence, we have to find conditions under which the tumor eradication equilibrium E is globally stable. Fortunately, Theorem 3.4 gives the sufficient conditions for the global stability of equilibrium E. By this theorem, the different treatment strategies are given below corresponding to different situations.
Case III. In such case, the positive equilibrium E 3 is globally stable. Antioncogenes mutate to cause a loss or reduction in its functions, but the oncogenes do not mutate. Using the treatment of oncolytic virus, we know that the correct dosage for complete tumor eradication is which is proved by Theorem 3.4.
Case IV. Since a1 c1 < a2 c2 , a2 b2 < a1 b1 , the oncogenes mutate but the antioncogenes do not mutate. In this case, there exist two possible cases: (i) the initial values locate in the stable region of normal cell equilibrium E 0 ; (ii) the initial values locate in the stable region of tumor cell equilibrium E 1 . In the former case, the tumor cells will extinct without any treatment. In the latter case, the oncolytic virotherapy is needed, and the correct dosage for complete tumor eradication is Case I. All positive solutions will tend to the tumor cell equilibrium (0, a2 c2 ) due to the mutations of both oncogene and antioncogene. Using the results of Theorem 3.4, the evolution of cells will be similar to case III without treatment.
In addition, if a1 c1 = a2 c2 , a2 b2 = a1 b1 , it is easy to see that little drug dose will lead to complete tumor eradication according to Theorem 3.4.