EFFECTS OF THE KILLING RATE ON GLOBAL BIFURCATION IN AN ONCOLYTIC-VIRUS SYSTEM WITH TUMORS∗

Oncologists and virologist are quite concerned about many kinds of issues related to tumor-virus dynamics in different virus models. Since the virus invasive behavior emerges from combined effects of tumor cell proliferation, migration and cell-microenvironment interactions, it has been recognized as a complex process and usually simulated by nonlinear differential systems. In this paper, a nonlinear differential model for tumor-virus dynamics is investigated mathematically. We first give a priori estimates for positive steadystates and analyze the stability of the positive constant solution. And then, based on these, we mainly discuss effects of the rate of killing infected cells on the bifurcation solution emanating from the positive constant solution by taking the killing rate as the bifurcation parameter.


Introduction
Mathematical modeling research of virus-cell interaction has a long history. Equally, extensive efforts have been dedicated to the mathematical modeling of cancer development by oncologists over many years. In this paper, we analyze a complex process that involves both virus-cell interaction and tumor growth, which is the interaction of the so called oncolytic viruses with tumors. The interaction between oncolytic virus and tumor cells is amenable to mathematical modeling using adaptations of techniques employed previously for modeling other types of virus-cell interaction. The corresponding model considers two types of cells growing in logistic fashion, which has the following form: where u and v are the sizes of uninfected and infected cells populations. The coefficients r 1 , r 2 , k, b, a are positive constants. r 1 , r 2 are the maximum per capita growth rates of uninfected and infected cells, k is the carrying capacity, b is the transmission rate, and a is the rate of infected cells killed by the virus. Model (1.1) is a version of the classical predator-prey model of a biological community first developed by Lotka [14] and Volterra [20]. The term buv describes the simplest correspondence between prey consumption and predator production similar to the law of mass action. The main result of the analysis of a more general model with respect to (1.1) in [22] consists in defining conditions required for maximum reduction of the tumor load, where whether the tumor cells would become extinct has not been stated clearly. In [16], a complete parametric analysis of dynamic regimes of a conceptual model (which contains model (1.1)) of anti-tumor virus therapy is presented. The authors described deterministic elimination of the tumor cells. Moreover, for different parameter values, they analyzed the asymptotic state to the system. As most solid tumors have distinct spatial structure, tumor-virus dynamics is usually described by systems of partial differential equations (such as [4-7, 10, 18, 24]), which is an obvious and necessary extension of ordinary differential equation models. In general, these partial differential equation systems are basically divided into two types: (I) The different populations of cells are continuously presented everywhere in the tumor at all times (see [6]); (II) The different populations of cells are separated by interfaces (see [7]).
In larger-scale of tumor-virus dynamical systems, where the models described by partial differential equations, the governing equations are typically of reactiondiffusion type since the distribution of the cells are spatially inhomogeneous. In recent decades, reaction-diffusion equations have been successfully used to model the dynamics of tumors, see [1,3,9,23] for examples. For a similar consideration, instead of (1.1), we are led to investigating the following reaction-diffusion system in a bounded domain Ω ⊂ R N with smooth boundary ∂Ω. Where ν is the outward unit normal vector on ∂Ω. With an appropriate change of variables, the model (1.2) can be simplified to the following system . Then m will be regarded as the rate of infected cells killed by the virus.
Bifurcation solution of partial differential equation systems is one of the principal concerns of many researchers, and many valuable results were obtained. For example, a diffusive predator-prey system with Holling-type functional response and homogenous Neumann boundary condition was considered in [8], where the authors considered the bifurcation solution emanating from a positive constant solution by taking the growth rate as a bifurcation parameter, and performed a detailed local and global bifurcation analysis to the diffusive system. In [17], the author surveyed some basic bifurcation results for equation in a setting of infinite dimensional Banach spaces. For a reaction-diffusion system, he found the bifurcating points and obtained a general bifurcation result by taking the diffusion coefficient as a bifurcation parameter. In [19], a semi-linear elliptic problem with nonlinear boundary conditions was studied. The author used the bifurcation results to investigate the set of positive solutions bifurcating from the trivial line. A generic Turing type reaction-diffusion system derived from the Taylor expansion near a constant equilibrium was considered in [11], where the authors studied the existence of Hopf bifurcations and steady-state bifurcations, analyzed the bifurcation direction and the stability of the bifurcating periodic obits. Moreover, numerical simulations were included to show the spatiotemporal dynamics. Regarded the growth rate of prey as the bifurcation parameter, a delayed predator-prey system with diffusion and Dirichlet boundary conditions was investigated in [12]. The authors showed that Hopf bifurcation occurs when the parameter varies. Furthermore, the explicit algorithm for determining the direction of the Hopf bifurcations and stability of the bifurcating periodic solutions were derived.
For a biological depletion model, the authors, in [21], discussed the stability of the positive constant steady-states and made a detailed description for the global bifurcation structure. Where they took the diffusion coefficient as the bifurcation parameter to investigate the bifurcation solution emanating from positive constant solutions for the model. It is worth noting that the bifurcation parameter sequence is non-monotonic (which means that the extra condition is needed) and the positive constant solutions both have no relationship with the bifurcation parameter.
Motivated by the techniques in [21], in this paper, we mainly study the effects of m on positive steady-states of (1.3), specifically, we consider the effects of m on bifurcations of (1.3), and our aim is to make a better description for the structure of non-constant steady-states. Different from [21], here, we will take m as the bifurcation parameter to discuss the bifurcation solution which emanates from the positive constant solution. The result shows that the local bifurcation emanating form the positive constant solution may be extended to global bifurcation if m is restricted in a suitable value range. Compared with [21], we find that the positive constant solution depends on the bifurcation parameter m. What's more, influenced by the model itself, in order to ensure the existence of the positive constant in our model, we must claim that m is in a bounded interval. So, the bifurcation points are in a bounded domain. Moreover, we see that the bifurcation points come in pairs, say, m − i and m + i . Besides that, {m − i } and {m + i } can be arranged orderly. Just for these reasons, it needs no extra conditions when we discuss the bifurcation solutions (While in [21], it does, see Remark 2.3 in Subsection 2.3). Our result shows that if the bifurcation curve Γ − j eventually meets some bifurcation points, and k is the biggest footnote of all these points at which Γ − j meets, then Γ − j must meet both (m − k , 0) and (m + k , 0).

A priori estimates of non-constant positive solutions
As a preparation for bifurcation analysis, in this subsection, we give a priori upper bounds for non-constant positive solutions of the elliptic system (2.1) Here, to give a simple and direct proof for the a priori upper bounds, we prefer an elliptic maximum principle due to Lou and Ni.
is a non-constant positive solution to the elliptic system (2.1), then Proof. It is clear that u < 1 is an immediate result of Lemma 2.1. We only need to show the upper bound for v. If v attains its maximum at some x 0 ∈ Ω, then by the second equation of (2.1) and Lemma 2.1 we get The positivity of v implies that

Stability of the constant positive solution
It is easy to verify that the system (2.1) has a unique positive constant solution (u * , v * ) with if and only if one of the following two conditions is satisfied: (H1) β < m < αβ β+1 , α > β + 1; (H2) αβ β+1 < m < β, α < β + 1. We next prove the local stability of (u * , v * ). Let 0 = λ 0 < λ 1 < λ 2 < · · · be the sequence of eigenvalues for the elliptic operator −∆ subject to the Neumann boundary condition on Ω, where each λ i has multiplicity s i ≥ 1. Let φ ij , 1 ≤ j ≤ s i , be the normalized eigenfunctions corresponding to λ i . Then the set {φ ij }, i ≥ 0, 1 ≤ j ≤ s i , forms a complete orthogonal basis in L 2 (Ω). The local stability of (u * , v * ) is as follows.
Proof. Consider the linearized operator of (2.1) at (u * , v * ) where a ij , b ij are constants. We find that It follows that µ is an eigenvalue of L if and only if for some i ≥ 0 the determinant of the corresponding matrix is zero, that is, If α < β + 1, then by (2.2), Q i > 0 for all i ≥ 0. This implies that the real parts Reµ of µ are negative for all eigenvalues µ, and so the steady-state (u * , v * ) is asymptotically stable.
When (u * , v * ) is unstable, we may expect the existence of non-constant steadystates near (u * , v * ). In the following, we will discuss the bifurcation under (H1).

Local bifurcation and global bifurcation
In this subsection, under the condition (H1), we fix α and β, and take m as a bifurcation parameter to analyze the bifurcating solution of (2.1), which bifurcates from (u * , v * ). By using the local bifurcation theory, we give a precise description for the structure of positive solutions near the bifurcation points. Then we prolong the bifurcation curves by the global bifurcation theory.
With the condition (H1), we suppose holds. Moreover, we define i α = i α (α, β, Ω) the largest positive integer such that Clearly, if (2.3) or (2.4)-(2.5) is satisfied, then 1 ≤ i α < ∞. In this case, we let It is easy to see that these points satisfy It is easy to show that whether α − 4β 2 − 5β − 1 is larger than 0 or smaller than 0, λ * is always positive.
Remark 2.3. In [21], the bifurcation parameter sequence {d i } is non-monotonic, so there needs the extra condition d i = d j for i = j. Here {m − i } and {m + i } are monotonic, and we do not need such additional conditions. Now, let Y = L 2 (Ω) × L 2 (Ω) be the Hilbert space with the inner produc- Then E is a Banach space with usual C 2 norm. Define operator F : where U = (u, v). Thus U is the solution of (2.1) if and only if U satisfies F (m, U ) = 0. With U * = (u * , v * ), we have F (m, U * ) = 0 for all β < m < βα β+1 . Theorem 2.3. Suppose (H1) holds. If

8)
holds and λ j is simple for some j. Then (m − j , U * ) and (m + j , U * ) are both bifurcation points of F (m, U ) = 0 with respect to the curve (m, U * ).
Proof. Note that the Fréchet derivative of F at U * is Let Φ = (φ, ψ) T ∈ kerL with kerL being the kernel space of L. And write as in the proof of Theorem 2.2. Then 0≤i≤∞,1≤j≤si So, there only holds detB j = 0. Moreover, the eigenfunction corresponding to λ j is φ j1 in view of the multiplicity of λ j is 1. Then < 0. Here we note that g 0 < 0, under the assumption (H1). Denote by R(L) the range space of L. Then codimR(L) = dim ker(L * ) = 1 in view of R(L) = (kerL * ) ⊥ .
Finally, since . Then the conditions of the local bifurcation theorem [2] are satisfied, and it follows that (m − j , U * ) is a bifurcation point. This completes the proof.
According The zero set of F consists of two curves (m, U * ) and Γ − j (s) in a neighborhood of the bifurcation point (m − j , U * ). What's more, in the same way as above we can verify that (m + j , U * ) is also a bifurcation point and we have the similar conclusion for the bifurcation point (m + j , U * ). Remark 2.4. Theorem 2.3 assumes that λ j is simple for some j ≥ 1. In one dimensional case, it is well known that all eigenvalues of −∆ are simple. However, in multi-dimensional case, in general, it is hard to know the multiplicity of λ j , j ≥ 1. Sometimes, λ j , j ≥ 1 may be simple for special Ω. For example, if Ω is a rectangle region, say, Ω = (0, a)×(0, b), then the set of eigenvalues of −∆ is {λ m,n = ( mπ a ) 2 + ( nπ b ) 2 } ∞ m,n=0 , and ϕ m,n (x, y) = cos mπ a x cos nπ b y, (x, y) ∈ (0, a) × (0, b), is the eigenfunction corresponding to λ m,n . For m = n, it is easy to verify that the corresponding eigenvalues have multiplicity one indeed.
The above local bifurcation result gives a relatively precise description for the structure of positive solutions near bifurcation points, but it provides no information on the bifurcation curve far from the equilibrium. Therefore, a further study is necessary in order to understand the bifurcation curve when it is far away from the bifurcation point. In the following, we investigate the coexistence of the steadystate solution by considering the global bifurcation. For simplicity, we suppose that Ω is one dimensional, say Ω = (0, l). By applying the global bifurcation theory and the Leray-Schauder degree theory to the compact linear operator, we can obtain the following global result. (ii) Γ − j meets some bifurcation points, which are in where f and g are higher-order terms of u and v. The constant steady-state (u * , v * ) of (2.1) is converted to (0, 0) of this new system.
In order to apply the global bifurcation theory in [13], firstly, we need to verify that 1 is an eigenvalue of K(m − j ) with algebraic multiplicity one. By the proof of We now compute ker( Then By the definition of G and G we obtain In the one dimensional case, let {φ i } be a complete orthogonal basis in L 2 ((0, l)).
By calculation we know that detA i = detB i , then detA i = 0 only for i = j, and ker(K * (m − j ) − I) is spanned by (1, f1 λj −f0 )φ j (In the case of one dimensional, we denote φ j1 by φ j ). Since which shows that the eigenvalue 1 has algebraic multiplicity one indeed. Secondly, we need to verify that the index(I − K(m − j ) − H, (m − j , 0)) changes when m crosses m − j . That is, for > 0 sufficient small, we should show that the index inequality where B is a sufficiently small ball centering at 0, and p is the sum of the algebraic multiplicities of eigenvalues of K(m) that are lager than 1.
Suppose that µ is an eigenvalue of K(m) and (φ, ψ) T is the corresponding eigenfunction. Then Thus µ satisfies Let µ(m), µ(m) be the roots of the equation (2.10). Then .
Take m = m − j . Then µ = 1 is a root of (2.10) if and only if i = j. And in this case These points satisfy m − j ∈ (β, m * j ), m + j ∈ (m * j , βα β+1 ). For > 0 sufficiently small, there results µ(m − j + ) > 1, and µ(m − j − ) < 1. Hence K(m − j + ) has exactly one more eigenvalues that are larger than 1 than K(m − j − ) dose. In the same way as above we can verify that this eigenvalue has algebraic multiplicity one, so the desired index inequality holds.
By the global bifurcation theorem, we know that either Γ − j is noncompact in (β, βα β+1 ) × E or Γ − j meets m + k and m − k where m − k = m − j for k = j. Then, we verify if the second alternative occurs, it must be the situation (ii) in Theorem 2.4. If Γ − j meets some other bifurcation points, say, Γ − j meets (m − k , 0), but dose not meet (m − i , 0) and (m + i , 0) for any i > k. Consider the system (2.1) on the interval (0, l/k) with x ∈ (0, l/k) in equations. If U is a solution of (2.11), then we construct a solution U to (2.1) by a reflective and periodic extension. Let x n = nl/k, n = 0, 1, · · · , k, and where λ i = (kπi/l) 2 . Then λ i = λ ik (λ ik is the ik-th eigenvalue for the elliptic operator −∆ under the one dimensional case). We obtain thatm − 1 = m − k is a bifurcation point of (2.1). Let Γ − 1 be the bifurcation branch of (2.11) that meets (m − 1 , 0). Then by the global bifurcation theorem, either it is noncompact in (β, βα β+1 ) × E or it meets case (1): (m − s , 0) or (m + s , 0) with s > 1; case (2): (m + 1 , 0). If case (1) occurs, then by the extension we know that Γ − j meets (m − sk , 0) or (m + sk , 0), where sk > k. This violates the definition of k. So, case (1) will not occur. Suppose that case (2) also will not occur. Since we consider the bifurcation from positive constant solution (u * , v * ), the situation αβ β+1 = m = β will not be considered. Otherwise, (u * , v * ) ≡ (0, 0). For m < β and m > αβ β+1 , we know that (u * , v * ) is stable in light of Theorem 2.2, so there has no nonconstant positive solution bifurcating from (u * , v * ) indeed. And then, Γ − 1 is noncompact in (β, βα β+1 )× E, and by the extension, Γ − j is noncompact in (β, βα β+1 ) × E (More precisely, Γ − 1 and Γ − j are noncompact in a maximal connected subset of (β, βα β+1 ) × E in view of the boundedness of positive solutions), which means that the projection of the bifurcation curve Γ − j on the m-axis contains (β, m − j ) or (m − j , βα β+1 ), which is the situation (i) in Theorem 2.4; If case (2) occurs, by extension we know that Γ − j meets (m + k , 0), which is the situation (ii) in Theorem 2.4. If Γ − j meets some other bifurcation points, say, Γ − j meets (m + k , 0), but does not meet (m − i , 0) and (m + i , 0) for any i > k. By the same argument, we can get the same conclusion.
Using the same way as above, we can obtain that the bifurcation curve Γ + j has similar properties with Γ − j . The proof is accomplished.

Discussion and conclusion
For tumor-virus models, from valuable results that the scholars have obtained, we know that the size of the tumor may be reduced in a short time if the injected viral density is suitable. This means that both infected and uninfected tumors can be eliminated with time, and complete recovery is possible. We investigate a mathematical model for cancer treatment by using oncolytic viruses. The viruses specifically infect and kill cancer cells. We establish the existence of the non-constant positive solution and we know the rate of killing the infected cell has effects on the global bifurcation.
In order to make good use of oncolytic virus, there is still much to be learned and discovered about oncolytic viruses and virus therapy. Because of possible adverse side effects, it is important to achieve reduction of the tumor cells with small dosage, which remains a problem. Searching for optimal schedules of injections is also a promising area for future research. The role of the immune reaction in controlling the extent of infection is still unclear. One may try to explore these problems and results established in this paper, although they are based on a number of simplifying assumptions, they still can be used as a rigorous tool. As more information becomes available, our model can be refined and expanded to incorporated more of the known biology. To this end, the development of biologically realistic mathematical models is an important tool to explore these problems.