# The Role of the Innate Immune System in Oncolytic Virotherapy.

1. IntroductionOncolytic virotherapy is a promising therapeutic strategy to destroy tumors [1]. Oncolytic viruses are viruses that selectively infect and replicate in tumor cells but spare normal cells, which have two types: wild-type oncolytic viruses which preferentially infect human cancer cells, and gene-modified viruses engineered to achieve selective oncolysis. In oncolytic viral therapy, oncolytic viruses infect tumor cells and replicate themselves in tumor cells; upon lysis of infected tumor cells, new virion particles burst out and proceed to infect additional tumor cells. This idea was initially tested in the middle of the last century and merged with renewed ones over the last 30 years due to the technological advances in virology and in the use of viruses as vectors for gene transfer (for the history of oncolytic viruses, see [2]). Oncolytic viruses have shown efficacy in clinical trials [3]. However, the immune response presents a challenge in maximizing efficacy. The major problem is the complexity of the innate and adaptive immune responses in the process of oncolytic viral therapy [4].

Mathematical models have been applied to the understanding of oncolytic virotherapy since fifteen years ago. Wu et al. [5] and Wein et al. [6] proposed and analyzed a system of partial differential equations that is essentially a radially symmetric epidemic model embedded in a Stefan problem to describe some aspect of cancer virotherapy. They were interested in three alternative virus-injection strategies: a fixed fraction of cells preinfected with the virus is introduced throughout the entire tumor volume, within the tumor core, or within the tumor rim. Wodarz [7] and his review paper [8] formulated a simple model with three ordinary differential equations. He studied three hypothetical situations: viral cytotoxicity alone kills tumor cells, a virus-specific lytic CTL response contributes to killing of infected tumor cells, and the virus elicits immunostimulatory signals within the tumor, which promote the development of tumor-specific CTL. Komarova and Wodarz [9] and Wodarz and Komarova [10] analyzed several possible mathematical formulations of oncolytic virus infection in terms of two ordinary differential equations, while Novozhilov et al. [11] analyzed ratio based oncolytic virus infection models. Bajzer et al. [12] used three ordinary differential equations to model specific cancer virotherapy with measles virus, and then they considered optimization of viral doses, number of doses, and timing with a simple mathematical model of three ordinary differential equations for cancer virotherapy [13].

Friedman et al. [14] proposed a free boundary problem with nonlinear partial differential equations to study brain tumor glioma with mutant herpes simplex virus therapy. The model incorporated immunosuppressive agent cyclophosphamide to reduce the effect of the innate immune response. This model reveals the oscillation of cell populations in the process of oncolytic viral therapy. Vasiliu and Tian [15] proposed a simplified model to study the cell population oscillation in oncolytic virotherapy, which may be caused by interaction between infected tumor cells and innate immune cells. To obtain a basic dynamic picture of oncolytic viral therapy, Tian [16] proposed a simple model with three ordinary differential equations to represent interaction among tumor cells, infected tumor cells, and oncolytic viruses and concluded that the viral therapeutic dynamics is largely determined by the viral burst size. To further understand how the viral burst size is affected, Wang et al. [17] and Tian et al. [18] incorporated virus lytic cycle as delay parameter into the basic model. These delay differential equation models gave another explanation of cell population oscillation and revealed a functional relation between the virus burst size and lytic cycle. In a recent paper [19], Choudhury and Nasipuri considered a simple model of three ordinary differential equations for the dynamics of oncolytic virotherapy in the presence of immune response. However, this model did not include the free virus population, and it may not give a complete picture of dynamics of viral therapy with innate immune response.

All proposed mathematical models have given some insights into oncolytic virotherapy. However, there is a considerable need to understand the dynamics of oncolytic virotherapy in the presence of immune responses [4], particularly, to understand the different effects of the innate immune system and adaptive immune system on virotherapy. In response to this call in [4], we plan to construct a comprehensive mathematical model for oncolytic virotherapy with both innate and adaptive immune responses. Toward this end, we will first build a mathematical model for oncolytic virotherapy with the innate immune system based on our basic model proposed in [16]. There are several types of cells that are involved in the innate immune response in virotherapy. So far, the experiments show that natural killer cells, macrophages, and neutrophils have significant effects in viral therapy [4]. For the sake of simplicity, we lump all these innate immune cell types to one variable, the innate immune cell population, in our mathematical model.

Our basic dynamical model for oncolytic virotherapy studied in [16] is as follows:

dx/dt = [lambda]x(1 - x + y/K) - [beta]xv, dy/dt = [beta]xv - [delta]y, dv/dt = b[delta]y - [beta]xv - [gamma]v, (1)

where x stands for uninfected tumor cells, y for infected tumor cells, and v for free viruses. For the details of explanations and results, the reader is referred to [16]. The innate immune response reduces infected tumor cells and viruses [4, 14]. We incorporate these effects into the basic model.

Denoting the innate immune cell population by z, we have the following system:

dx/dt = [lambda]x(1 - x + y/C) - [beta]xv, dy/dt = [beta]xv - [mu]yz - [delta]y, dv/dt = b[delta]y - [beta]xv - kvz - [gamma]v, dz/dt = zyz - [rho]z, (2)

where [lambda] is tumor growth rate, C is the carrying capacity of tumor cell population, [beta] is the infected rate of the virus, [mu] is immune killing rate of infected tumor cells, [delta] is death rate of infected tumor cells, b is the burst size of oncolytic viruses (i.e., the number of new viruses coming out from a lysis of an infected cell), k is immune killing rate of viruses, [gamma] is clearance rate of viruses, s is the stimulation rate of the innate immune system, and [rho] is immune clearance rate.

In this study, we analyze this four-dimensional system (2). Our analysis and numerical study show that the dynamics of the model is largely determined by the viral burst size b and parameters related to the innate immune response. We can denote the dynamical behaviors of the model by b, the ratio of the innate immune killing rate of infected tumor cells over the innate immune killing rate of free viruses by [mu]/k, the relative innate immune killing rate of viral therapy by K, the ratio of the innate immune clearance rate over the stimulation rate of the innate immune system by [rho]/Cs, and the relative innate immune response decay rate by N. These two combined parameters are related to the innate immune response. Comparing with the basic model in [16], there are several critical values of the oncolytic viral burst size [mathematical expression not reproducible], where [bar.b] is a function of the relative innate immune killing rate K and relative innate immune response decayrate N. When b is smaller than [bar.b] and the two relative rates are constrained in some intervals, the system has 5 equilibrium solutions and 2 of them are positive, while these two positive equilibrium points coalesce when b = [bar.b]. When b is greater than [mathematical expression not reproducible] or [bar.b] and the two relative rates are in the complement intervals, the system has at most 3 equilibrium solutions with 0 innate immune components. An interesting fact is that the equilibrium points where Hopf bifurcations arise for both models, (2) and the one in [16], are corresponding to each other. Therefore, we may conclude that the innate immune response complicates the oncolytic virotherapy in the way of creating more equilibrium solutions when the oncolytic viral burst size is not too big, say less than [bar.b], while the dynamics is similar to the system without the presence of the innate immune response when the oncolytic viral burst size is big.

The rest of this article is organized as follows. Section 2 presents analysis of model (2) in 4 subsections. Section 2.1 gives preliminary results about the model, Section 2.2 calculates equilibrium solutions, Section 2.3 studies stability of equilibrium solutions, and Section 2.4 provides bifurcation analysis of the model and the main Theorem 14 to summarize the dynamical behaviors of the model (2). Section 3 provides a numerical study and discussion, where we numerically compute some dynamical characteristics and simulate the model, and we also compare the dynamics of our model with the basic model of oncolytic virotherapy in [16].

2. Analysis of the Mathematical Model

We conduct a detailed analytical study of our proposed model in this section. The major properties of dynamical behaviors of our model are summarized in Theorem 14. For each analysis result, we also provide biological interpretations or implications as appropriate.

2.1. Positive Invariant Domain. In order to simplify system (2), we apply nondimensionalization by setting [tau] = [delta]t, x = C[bar.x], y = C[bar.y], v = C[bar.v], z = C[bar.z] and rename parameters r = [lambda]/[delta], a = C[beta]/[delta], c = [mu]C/[delta], d = kC/[delta], e = [gamma]/[delta], m = sC/[delta], and n = [rho]/[delta]. Then system (2) becomes

d[bar.x]/d[tau] = r[bar.x](1 - [bar.x] - [bar.y]) - a[bar.x][bar.y], d[bar.y]/d[tau] = a[bar.x][bar.v] - c[bar.y][bar.z] - [bar.y], d[bar.v]/d[tau] = b[bar.y] - a[bar.x][bar.v] - d[bar.v][bar.z] - e[bar.v], d[bar.z]/d[tau] = [bar.m][bar.z] - n[bar.z]. (3)

For convenience, dropping all the bars and writing [tau] as t, we obtain

dx/dt = rx(1 - x - y) - axv, dy/dt = axv - cyz - y, dv/dt = by - axv - dvz - ev, dz/dt = myz - nz. (4)

We assume that all parameters are nonnegative.

Lemma 1. If x(0) [greater than or equal to] 0, y(0) [greater than or equal to] 0, v(0) [greater than or equal to] 0, and z(0) [greater than or equal to] 0, then the solution of system (4) x(t) [greater than or equal to] 0, y(t) [greater than or equal to] 0, v(t) [greater than or equal to] 0, and z(t) [greater than or equal to] 0 for t [greater than or equal to] 0. Furthermore, if 0 [less than or equal to] x(0) + y(0) [less than or equal to] 1, then 0 [less than or equal to] x(t) + y(t) [less than or equal to] 1 for t [greater than or equal to] 0, and lim [sup.sub.t[right arrow]+[infinity]] v(t) [less than or equal to] b/e.

Proof. The proof is similar to that of Lemma 3.2 in [16]. Here we only show the second part of this lemma by using comparison theorem for ODEs; that is, if 0 [less than or equal to] x(0) + y(0) [less than or equal to] 1, then 0 [less than or equal to] x(t) + y(t) [less than or equal to] 1 and lim [sup.sub.t[right arrow]+[infinity]] v(t) [less than or equal to] b/e. In fact, since x(t), y(t), and v(t) are nonnegative for all t [greater than or equal to] 0,

x'(t) = rx(1 - x - y) - axv [less than or equal to] rx(1 - x - y) [less than or equal to] rx(1 - x). (5)

As x(0) [less than or equal to] x(0) + y(0) [less than or equal to] 1, by the comparison theorem we get x(t) [less than or equal to] 1. On the other hand, since

x'(t) + y'(t) = rx(1 - x - y) - cyz - y [less than or equal to] rx(1 - x - y) [less than or equal to] r(1 - x - y), (6)

again by the comparison theorem we also have 0 [less than or equal to] x(t) + y(t) [less than or equal to] 1. It follows that 0 [less than or equal to] y(t) [less than or equal to] 1. Then v'(t) = by - ax v - dvz - ev [less than or equal to] b - ev, so by the comparison theorem v(t) [less than or equal to] b/e + v(0) exp(-et). Taking lim sup both sides yield lim [sup.sub.t[right arrow]+[infinity]] v(t) [less than or equal to] b/e.

We then conclude that the positive invariant domain of system (4) is

D = {(x, y, v, z) : x [greater than or equal to] 0, y [greater than or equal to] 0, v [greater than or equal to] 0, z [greater than or equal to] 0, 0 [less than or equal to] x + y [less than or equal to] 1}. (7)

This is also a biological meaningful range for the variables. We will regard the whole domain D as a global domain.

2.2. Equilibrium Solutions. We know that the dynamics of oncolytic viral therapy without the presence of the immune response can be characterized by the burst size b [16]. The effects of the innate immune system on the virotherapy in our model are encoded in the parameters c, d, and m. To understand how the innate immune system affects the dynamics of oncolytic viral therapy, we use three combined parameters, the viral burst size b, the relative immune killing rate K = c/d, and the relative immune response decay rate N = n/m, to describe the solution behaviors of our model. We summarize the possible equilibrium solutions in the following lemma.

Lemma 2. When (r/a)(1/N - 1) < K < 1/(a + e- aN), N < 1 + (1/2)(e/a + 1/r - [square root of [(e/a + 1/r).sup.2] + 4/r))], and b > [bar.b] with b > 1 + e/a, system (4) has 3 equilibrium solutions: [E.sub.0], [E.sub.1], and [E.sub.2]. When either K [less than or equal to] (r/a)(1/N - 1) or K > 1/(a + e - aN) and b = [bar.b] with g(K) > 0, system (4) has a unique positive equilibrium solution: [E.sub.3]. When either K [less than or equal to] (r/a)(1/N - 1) or K > 1/(a + e - aN) and b < [bar.b] with g(K) > 0, system (4) has two distinct positive equilibrium solutions: [E.sub.4] and [E.sub.5].

In what follows we will analyze the existence of equilibrium solutions. First, let X = [(x, y, v, z).sup.T] and

F(X) = [(rx (1 - x - y) - axv, axv - cyz - y, by - axv - dvz - ev, myz - nz).sup.T].

Then system (4) can be simply written as the autonomous system dX/dt = F(X). We assume that (x, y, v, z) [member of] D. The equilibrium points are solutions of the equation F(X) = 0; that is,

x(r(1 - x - y) - av) = 0, axv = y (cz + 1), by = v (ax + dz + e), z (my - n) = 0. (9)

If x = 0, then, from the second and the third equation of (9), y(cz + 1) = 0 and by = v(dz + e). Since cz + 1 > 0, then y = 0. It leads to v(dz + e) = 0, which implies v = 0. The last equation of (9) gives -nz = 0, which implies z = 0. Thus E0 = (0, 0, 0, 0) is an equilibrium point.

If x [not equal to] 0, the first equation of (9) implies r(1 - x - y) = av. Consider the last one z(my - n) = 0. If z = 0, from the second and the third equation of (9), we get axv = y and by = v(ax + e). These follow abxv = v(ax + e) and hence v(abx - ax - e) = 0. If v = 0, then y = 0 and r(1 - x) = 0, which implies that x = 1. So [E.sub.1] = (1, 0, 0, 0) is an equilibrium point.

Now if v [not equal to] 0, then a(b - 1)x = e. Since we want to find positive equilibrium points, we assume b > 1. Then x = e/a(b - 1). From the equation r(1 - x - y) = av, we have rx(1 - x) = rxy + axv = rxy + y = y(1 + rx). It follows that

y = rx(1 - x)/1 + rx = re(ab - a - e)/a(b - 1)(ab - a + re). (10)

Since axv = y, we have v = y/ax = r(ab - a - e)/a(ab - a + re). Thus we get the third critical point

[E.sub.2] = (e/a(b - 1), re(ab - a - e)/a(b - 1)(ab - a + re), r(ab - a - e)/a(ab - a +re), 0). (11)

Notice that in order for the first three coordinates of [E.sub.2] to be positive, it is enough that ab - a - e > 0; that is, b > e/a + 1.

Next, if z [not equal to] 0, then y = n/m. Set N = n/m, then y = N. From the equation r(1 - x - y) = av, we can derive x = 1 - N - av/r. By the third equation of (9), z = ((b - 1)N - ev)/(cN + dv). Plugging these expressions into the second equation gives f(v) = [v.sup.3] + [a.sub.2][v.sup.2] + [a.sub.1]v + [a.sub.0] = 0, where

[a.sub.2] = c/d N + r/a (N - 1), [a.sub.1] = rN/[a.sup.2] x c/d (-a + aN - e + d/c), [a.sub.0] = b x r[N.sup.2]/[a.sup.2] x c/d. (12)

Since f (0) = [a.sub.0] > 0 and [lim.sub.v[right arrow]--[infinity]] f(v) = - [infinity], f(v) has at least one negative root, say [v.sub.1][sigma] < 0. Assume that [v.sub.1] is the least root. If [a.sub.1] > 0 and [a.sub.2] > 0, then f has no sign changes at all and hence the other two roots of f are either both negative or both complex. Notice that [a.sub.1] > 0 and [a.sub.2] > 0 are equivalent to (r/a)(1/N - 1) < c/d < 1/(a + e - aN) and N < 1 + (1/2)(e/a + 1/r - [square root of (([e/a + 1/r).sup.2] + 4/r))]. In this case, the system only has 3 equilibrium points: [E.sub.0], [E.sub.1], and [E.sub.2].

We now look at other situations of f(v). Taking derivative, f'(v) = 3[v.sup.2] + 2[a.sub.2]v + [a.sub.1]. Bylong division,

f(v) = (1/3 v + [a.sub.2]/9) f'(v) + 2/9 (3[a.sub.1] - [a.sup.2.sub.2]) v - 1/9 [a.sub.1][a.sub.2] + [a.sub.0]. (13)

Suppose that [a.sup.2.sub.2] - 3[a.sub.1] > 0, then f' has 2 distinct roots, [v.sup.*.sub.1] and [v.sup.*.sub.2], where [v.sup.*.sub.2] = (-[a.sub.2] + [square root of ([a.sup.2.sub.2] - 3[a.sub.1])/3 =: A)]. If [v.sup.*.sub.2] > 0 and f([v.sup.*.sub.2]) = 0, then f has one negative root, [v.sub.1], and one positive root, [v.sub.2] = [v.sub.3] = [v.sup.*.sub.2] = (9[a.sub.0] - [a.sub.1][a.sub.2])/2([a.sup.2.sub.2] - 3[a.sub.1]). In this case, in addition to the 3 equilibrium points [E.sub.0], [E.sub.1], and [E.sub.2], system (4) has one positive equilibrium point: [E.sub.3] = (1 - N - aA/r, N, A, ((b - 1)N - eA)/(cN + dA)). To guarantee all four coordinates of [E.sub.3] are positive, we need to impose 1 - N - aA/r > 0 and (b - 1)N - eA > 0, which imply that [v.sup.*.sub.2] = A < u* = min{(1 - N)(r/a), ((b - 1)/e)N}.

On the other hand, if [v.sup.*.sub.2] > 0 and f([v.sup.*.sub.2]) < 0, then f has one negative root, [v.sub.1], and 2 distinct positive roots, 0< [v.sub.2] < [v.sup.*.sub.2] < [v.sub.3] < u*. Hence system (4) has two positive equilibrium points: [E.sub.4,5] = (1 -N - a[v.sub.2,3]/r, N, [v.sub.2,3], ((b - 1)N - e[v.sub.2,3])/(cN + d[v.sub.2,3])). Notice that [v.sup.*.sub.2] = A > 0 is equivalent to either [a.sub.2] [less than or equal to] 0, [a.sup.2.sub.2] - 3[a.sub.1] > 0 or [a.sub.2] > 0 > [a.sub.1]. Furthermore, [a.sub.2] [less than or equal to] 0, [a.sup.2.sub.2] - 3[a.sub.1] > 0 are equivalent to c/d [less than or equal to] (r/a)(1/N - 1) and g(c/d) > 0,where g(x) = [x.sup.2] + (r/aN - r/a + 3re/[a.sup.2]N)x + ([r.sup.2]/[a.sup.2])[(1/N - 1).sup.2] - 3r/[a.sup.2]N; [a.sub.2] > 0 > [a.sub.1] is equivalent to c/d > max{(r/a)(1/N - 1), 1/(a + e - aN)}. The condition f([v.sup.*.sub.2]) [less than or equal to] 0 is equivalent to [v.sup.*.sub.2] [greater than or equal to] (9[a.sub.0] - [a.sub.1][a.sub.2])/2([a.sup.2.sub.2] - 3[a.sub.1]). That is, (-[a.sub.2] + [square root of ([a.sup.2.sub.2] - 3[a.sub.1])/3 [greater than or equal to] (9[a.sub.0] - [a.sub.1][a.sub.2])/2([a.sup.2.sub.2] - 3[a.sub.1]), and we have [bar.b] := [a.sup.2]d/rc[N.sup.2] x (2[a.sup.3.sub.2] + 9[a.sub.1][a.sub.2] + 2[([a.sup.2.sub.2] - 3[a.sub.1]).sup.3/2])/27 [greater than or equal to] b. The critical value [bar.b] is important for describing the dynamics of system (4).

We summarize the details of the analysis above as follows.

(i) When x = 0, we have equilibrium solution [E.sub.0] = (0, 0, 0, 0).

(ii) When x [not equal to] 0, we have the following cases.

(a) If z = 0, then

(1) when v = 0, we have equilibrium solution [E.sub.1] = (1, 0, 0, 0).

(2) when v [not equal to] 0 and b > 1 + e/a, we have

[E.sub.2] = (e/a(b - 1), re(b - a - e)/a(b - 1)(ab - a + re), r(ab - a - e)/a(ab - a + re), 0). (14)

(b) If z [not equal to] 0, then x = 1 - N - (a/r)v, y = N, z = ((b - 1)N - ev)/(cN + dv), and v satisfies the following cubic equation:

[v.sup.3] + [a.sub.2][v.sup.2] + [a.sub.1]v + [a.sub.0] = 0, (15)

where [a.sub.2] = (c/d)N + (r/a)(N - 1), [a.sub.1] = rN/[a.sup.2] x (c/d)(-a + aN - e + d/c), [a.sub.0] = b x r[N.sup.2]/[a.sup.2] x c/d.

In this case, we can conclude the following.

(1) If K [member of] ((r/a)(1/N - 1), 1/(a + e - aN)), N < 1 + (1/2)(e/a + 1/r - [square root of ([(e/a + 1/r).sup.2] + 4/r))], and b > 1 + e/a, then system (4) has three equilibrium points: [E.sub.0], [E.sub.1], and [E.sub.2].

(2) If either K [less than or equal to] (r/a)(1/N - 1), g(K) > 0, or K > 1/(a + e - aN) and b = [bar.b], then system (4) has a unique positive equilibrium point:

[E.sub.3] = (1 - N - aA/r, N, A, (b - 1)N - eA/CN + dA), (16)

where [v.sup.*.sub.2] = A < u* := min{(1 - N)(r/a), (b - 1)N/e] and

g(x) = [x.sup.2] + (r/aN - r/a + 3re/[a.sup.2]N)x + [r.sup.2]/[a.sup.2][(1/N - 1).sup.2] - 3r/[a.sup.2]N. (17)

(3) If either K [less than or equal to] (r/a)(1/N - 1), g(K) > 0, or K > 1/(a + e - aN) and b < [bar.b], then system (4) has two distinct positive equilibrium points:

[E.sub.4] = (1 - N - [v.sub.2]/q, N, [v.sub.2], (b - 1) N - e[v.sub.2]/cN + d[v.sub.2]), [E.sub.5] = (1 - N - [v.sub.3]/q, N, [v.sub.3], (b - 1) N - e[v.sub.3]/cN + d[v.sub.3]), (18)

where [v.sub.2] and [v.sub.3] are two distinct positive real roots of the above cubic equation that satisfy 0< [v.sub.2] < [v.sup.*.sub.2] < [v.sub.3] < u*, in which [v.sup.*.sub.2] = A := (-[a.sub.2] [square root of ([a.sup.2.sub.2] - 3[a.sub.1])/3)].

In order to interpret our mathematical conditions biologically, we need to understand the combined parameters N and K first. N = [rho]/sC can be considered as a relative immune response decay rate since p is innate immune cell death rate, s is innate immune stimulating rate by infection, and C is tumor carrying capacity. K = c/d = [mu]/k is the ratio of the rate of immune killing infected tumor cells over the rate of immune killing viruses, which can be considered as a relative immune killing rate of viral therapy since it describes the ability of the innate immune system destroying infection versus destroying viruses. Biological interpretation of Lemma 2 is as follows. If there are no tumor cells, we have zero equilibrium [E.sub.0]. If we do not consider the effect of the immune system, and the viruses are not powerful, that is, the burst size is smaller than a critical value, then the system has the equilibrium [E.sub.1] with only tumor cells; if the viruses are powerful, that is, the burst size is greater than a critical value, then the system has the equilibrium [E.sub.2] with the coexistence of tumor cells, infected tumor cells, and viruses. When we consider the immune effect, if the burst size is another critical value [bar.b] and the relative immune killing rate satisfies some conditions, the system has a unique positive equilibrium; if the burst size is greater than that critical value and the relative immune killing rate satisfies certain similar conditions, the system has other two positive equilibria. Combining stability analysis, we can have more biological implications in the next two subsections.

2.3. Stability Analysis of Equilibrium Solutions. In this subsection, we apply various methods to analyze the asymptotically stable behaviors of equilibrium solutions. By finding the eigenvalues of the variational matrix of system (4) at the equilibrium points, we show [E.sub.0] is locally unstable, [E.sub.1] is locally asymptotically stable if b < 1 + e/a and unstable if b >1 + e/a, and [E.sub.2] is locally asymptotically stable if b is in some range, while [E.sub.3], [E.sub.4], and [E.sub.5] are locally unstable when b, K, and N satisfy some conditions. We use Lyapunov functions to show [E.sub.1] is globally asymptotically stable if b < 1 + e/a and N > 1. We apply the center manifold theorem to show [E.sub.1] is locally asymptotically stable if b = 1 + e/a. We summarize the main results in Lemma 3. For the combined parameter values, [mathematical expression not reproducible], i = 1, 2, J, [[bar.b].sub.j], j = 1, 2, 3, they will be defined in the following context. For methods we applied in this subsection, we refer to Allen [20] for basic knowledge of stability analysis and Carr [21] for center manifolds.

Lemma 3. [E.sub.0] is unstable. [E.sub.1] is globally asymptotically stable when b < 1 + e/a and N > 1 and unstable when b > 1 + e/a. [E.sub.2] is locally asymptotically stable when b [member of] ([mathematical expression not reproducible]) [intersection] J. [E.sub.3] is unstable when [bar.b] < [[bar.b].sub.1]. [E.sub.4] and [E.sub.5] are unstable when [bar.b] < [[bar.]sub.2,3].

We look at the stability of trivial equilibrium solutions first. The variational matrix of system (4) is given by

[mathematical expression not reproducible]. (19)

At the critical point [E.sub.0], the variational matrix is

[mathematical expression not reproducible]. (20)

The corresponding eigenvalues are r, -1, -e, and -n. We know that the local stability of [E.sub.0] of system (4) is the same as that of the linearized system at [E.sub.0]. Since r > 0, [E.sub.0] is locally unstable. For system (4), the local stable invariant manifold is tangent to the y-v-z subspace, and the unstable invariant manifold is tangent to the x-axis. Biologically, the tumor will grow from an initial small value around [E.sub.0] without viruses and infected tumor cells.

Proposition 4. The equilibrium solution [E.sub.1] is locally asymptotically stable when b < 1 + e/a, and it is locally unstable when b > 1 + e/a. [E.sub.1] is globally asymptotically stable when b < 1 + e/a and N > 1.

Proof. At the equilibrium point [E.sub.1], the variational matrix is

[mathematical expression not reproducible]. (21)

The characteristic polynomial of this matrix is ([lambda] + n)([lambda] + r)([[lambda].sup.2] + (1 + a + e)[lambda] + a + e - ab). The eigenvalues are [[lambda].sub.1] = -r, [[lambda].sub.2] = - n, and [[lambda].sub.3,4] = (1/2)(-1 - a - e [+ or -] [square root of ([(1 - a - e).sup.2] + 4ab))]. Since the eigenvalues [[lambda].sub.1] = -r, [[lambda].sub.2] = -n, and [[lambda].sub.3] = (1/2)(-1 - a - e - [square root of ([(1 - a - e).sup.2] + 4ab))] are negative for all positive parameters, [E.sub.1] is locally asymptotically stable if and only if [[lambda].sub.4] < 0. This is equivalent to [square root of ([(1 - a - e).sup.2] + 4ab)] < 1 + a + e, which is the same as b < 1 + e/a. Similarly, if b > 1 + e/a, then [[lambda].sub.4] = (1/2)(-1 - a - e + [square root of ([(1 - a - e).sup.2] + 4ab))] > 0, and hence [E.sub.1] is unstable.

In fact, we can show that [E.sub.1] is globally asymptotically stable in the positive invariant domain D when b < 1 + e/a and m < n by constructing two Lyapunov functions according to different ranges of the parameter b. For convenience, we translate [E.sub.1] into the origin by setting x = 1 - [bar.x], y = [bar.y], v = [bar.v], and z = [bar.z]. After dropping all the bars over variables, system (4) becomes

dx/dt = -rx + ry + av + r[x.sup.2] - rxy - axv, dy/dt = av - axv - cyz - y, dv/dt = by + axv - dvz - (a + e)v, dz/dt = myz - nz, (22)

while the domain D is translated to [D.sub.1] = {(x, y, v, z) : 0 [less than or equal to] x [less than or equal to] 1, y [greater than or equal to] 0, v [greater than or equal to] 0, z [greater than or equal to] 0, 0 [less than or equal to] x - y [less than or equal to] 1}. For any initial condition ([x.sub.0], [y.sub.0], [v.sub.0], [z.sub.0]) [member of] [D.sub.1], the solution of (22) satisfies 0 [less than or equal to] x(t) [less than or equal to] 1, 0 [less than or equal to] y(t) [less than or equal to] 1, v(t) [greater than or equal to] 0, and z(t) [greater than or equal to] 0. Now we construct two Lyapunov functions corresponding to the values of the parameter b to prove y(t) and v(t) approach 0, then we show x(t) and z(t) also tend to 0.

When 0 < b < 1, we define the Lyapunov function [V.sub.1] (x, y, v, z) = y + v. It is clear that [V.sub.1] (x, y, v, z) > 0, and the orbital derivative d[V.sub.1]/dt = dy/dt + dv/dt = (b - 1)y - cyz - dvz - ev <0. When 1 [less than or equal to] b < 1 + e/a, consider the Lyapunov function [V.sub.2](x, y, v, z) = (1/2)ab(a + e)[y.sup.2] + [a.sup.2]byv + (1/2)[a.sup.2][v.sup.2]. Obviously, [V.sub.2](x, y, v, z) > 0. The orbital derivative along a solution is given by

[mathematical expression not reproducible]. (23)

1 [less than or equal to] b < 1 + e/a; that is, ab - a - e < 0 and 1 - b [less than or equal to] 0; therefore d[V.sub.2]/dt < 0. Combining these two Lyapunov functions gives y(t) [right arrow] 0 and v(t) [right arrow] 0 as t [right arrow] [infinity] when b < 1 + e/a. Considering the component x(t), we have

dx/dt = -rx + ry + av + r[x.sup.2] - rxy - axv = (1 - x)(ry + av - rx) [less than or equal to] ry + av - rx. (24)

By the comparison theorem,

0 [less than or equal to] x(t) [less than or equal to] x(0)[e.sup.-rt] + [e.sup.-rt] [[integral].sup.t.sub.0] (ry(s) + av(s)) [e.sup.rs]ds. (25)

Taking limit of both sides as t [right arrow] [infinity] and using the L'Hospital's Rule and the fact that y(t) and v(t) approach 0 yield x(t) [right arrow] 0. Finally, since y(t) [less than or equal to] 1, we have dz/dt [less than or equal to] (m - n)z. By the comparison theorem, 0 [less than or equal to] z(t) [less than or equal to] z(0)[e.sup.(m-n)t] [right arrow] 0 as t [right arrow] [infinity], since m - n < 0. Therefore system (3) has a global attractor [E.sub.1].

Biologically, Proposition 4 is easy to understand. When the viral burst size is smaller than a critical value which is b = 1 + e/a, there will not be enough newly produced viruses to infect tumor cells. The therapy fails. The model system will be stable in the state of tumor cells and free of infected tumor cells, viruses, and immune cells. Proposition 5 ensures mathematically that this critical burst size is the smallest one that will make the virotherapy completely fail. One obvious medical implication is that we have to genetically increase the viral burst size in order to have effective virotherapy.

We next consider the stability of [E.sub.1] when b = 1 + e/a. This is a critical case, since the linearized system at [E.sub.1] has three negative eigenvalues and one zero eigenvalue. we have to reduce the system to its local center manifold. We actually have the following proposition.

Proposition 5. The equilibrium solution [E.sub.1] is locally asymptotically stable when b = 1 + e/a.

Proof. Consider b = 1 + e/a, which implies that ab = a + e. The linearized system at [E.sub.1] has three negative eigenvalues and one zero eigenvalue. In order to determine the stability of [E.sub.1], we use the center manifold theorem to reduce system (22) into a center manifold, and then we study the reduced system. So we separate it into two parts, one with zero eigenvalue and the other with negative eigenvalues. The matrix corresponding to the linear part of system (22) is

[mathematical expression not reproducible], (26)

which has eigenvalues -r, -(1 + ab), 0, and -n. The associated eigenvectors of these eigenvalues, respectively, are [V.sup.T.sub.1] = (1, 0, 0, 0), [V.sup.T.sub.2] = (ab - r, 1 + ab - r, -b(1 + ab - r), 0), [V.sup.T.sub.3] = (ar + a, ar, r, 0), and [V.sup.T.sub.4] = (0, 0, 0, 1). System (22) can be written as dX/dt = LX + F, where F = [(r[x.sup.2] - rxy - axv, -axv-cyz, axv - dvz, myz).sup.T]. Set T = ([V.sub.1], [V.sub.2], [V.sub.3], [V.sub.4]) and X = TY; then

dY/dt = [T.sup.-1]LTY + [T.sup.-1]F, (27)

where

[mathematical expression not reproducible], (28)

and Y = [([y.sub.1], [y.sub.2], [y.sub.3], [y.sub.4]).sup.T] is determined by

x = [y.sub.1] + (ab - r) [y.sub.2] + (ar + a) [y.sub.3], y = (1 + ab - r) [y.sub.2] + ar[y.sub.3], v = -b (1 + ab - r) [y.sub.2] + r[y.sub.3], z = [y.sub.4]. (29)

Denote [T.sup.-1]F = [([f.sub.1], [f.sub.2], [f.sub.3], [f.sub.4]).sup.T]; then we can express [f.sub.i], i = 1, 2, 3, 4, in terms of [y.sub.i]:

[mathematical expression not reproducible], (30)

where [A.sub.ij], [B.sub.ij], [C.sub.ij], and [D.sub.ij] are coefficients that can be easily determined. The transformed system can be expressed as

[mathematical expression not reproducible], (31)

where

[mathematical expression not reproducible], (32)

It is easy to check that each [f.sub.i], i = 1, 2, 3, 4, is twice differentiable function, [f.sub.i](0, 0, 0, 0) = 0 and D[f.sub.i](0, 0, 0, 0) = 0, where D[f.sub.i] is the Jacobian matrix of the function [f.sub.i]. By the center manifold theorem, there exists a center manifold:

[mathematical expression not reproducible], (33)

with h(0) = 0 and Dh(0) = 0, and it satisfies

[mathematical expression not reproducible]. (34)

Since h(0) = 0 and Dh(0) = 0, we can assume that

[mathematical expression not reproducible]; (35)

here we use the variable u instead of [y.sub.3] for simplicity. Then we can compute

[mathematical expression not reproducible]. (36)

By substituting [f.sub.i]'s into (34) and equating both sides of the equation, we can get [C.sub.33] = -[a.sub.2](r + 1)(b - 1)/(ab + 1) < 0, since b = 1 + e/a > 1. The asymptoticallybehavior of zero solution of system (31) is governed by that of the equation d[y.sub.3]/dt = [f.sub.3](h([y.sub.3]), [y.sub.3]) or d[y.sub.3]/dt = [C.sub.33][y.sup.2.sub.3] + O([y.sup.3.sub.3]). Since [C.sub.33] < 0, [y.sub.3] = 0 is locally asymptotically stable. Therefore, the trivial solution of system (22) is locally asymptotically stable.

We next consider the stability of the equilibrium solution [E.sub.2] = ([bar.x], [bar.y], [bar.v], [bar.z]), where [bar.x] = e/a(b - 1), [bar.y] = re(ab - a - e)/a(b - 1) (ab - a + re), [bar.v] = r(ab - a - e)/a(ab - a + re), and [bar.z] = 0. For lately defined [mathematical expression not reproducible] and J, we have a proposition as follows.

Proposition 6. When b [member of] ([mathematical expression not reproducible]) [intersection] J [not equal to] 0, [E.sub.2] is locally asymptotically stable.

We show this proposition as follows. The variational matrix at [E.sub.2] is given by

[mathematical expression not reproducible]. (37)

The characteristic polynomial of L is

p([lambda]) = [absolute value of ([lambda]I - L)] = [[lambda] - (m[bar.y] - n)] q([lambda]), (38)

where q([lambda]) = [[lambda].sup.3] + [a.sub.1][[lambda].sub.2] + [a.sub.2][lambda] + [a.sub.3], and

[a.sub.1] = re + a - a + abe/a(b - 1), [a.sub.2] = re(be + b - 1)/a[(b - 1).sup.2] + re(ab - a - e)(r - a)/a(b - 1)(ab - a +re), [a.sub.3] = re(ab - a - e)/a(b - 1). (39)

By Routh-Hurwitz Criterion, all roots of q([lambda]) have negative real parts if and only if

[mathematical expression not reproducible]. (40)

Since [mathematical expression not reproducible]. And those conditions in (40) are equivalent to [H.sub.2] = [a.sub.1][a.sub.2] - [a.sub.3] > 0. This inequality is the same as

[phi](b) = a(b - 1)(ab - a + re)/ab - a + re + abe - (be + b - 1)(ab - a + re)/(b - 1)(ab - a - e) < r - a. (41)

Therefore, we can conclude that if [bar.y] < N and [phi](b) < r - a; then [E.sub.2] is locally asymptotically stable. Now we can refine this result by considering H(b) = [H.sub.2], and

[PHI](x) = -[a.sup.3][x.sup.4] + [a.sup.2](3e + [e.sup.2] + r - a - ae + 1)[x.sup.3] + ae (3re + 3a + r[e.sup.2] + 3ae + r + [r.sup.2] - [a.sup.2])[x.sup.2] + [e.sup.2] (3ar + 2aer + [r.sup.2]e + 2[a.sup.2])x + r[e.sup.3] (r + a). (42)

It is easy to check that H(b) = re[PHI](b - 1)/[a.sup.2][(b - 1).sup.3](ab - a + re). Since b > 1 + e/a, H(b) and [PHI](b - 1) have the same roots. Since [PHI](e/a) = [e.sup.3](1 + r + e + a)(1 + e + a)((1 + r)/a) > 0 and [lim.sub.x[right arrow][+ or -][infinity]] [PHI](x) = -[infinity], there are at least one [x.sub.1] < e/a and one [x.sub.2] > e/a so that [PHI]([x.sub.1]) = [PHI]([x.sub.2]) = 0. Then H(b) has at least one root [mathematical expression not reproducible] and one root [mathematical expression not reproducible]. We consider three different cases as follows.

(i) If [PHI] has 4 distinct real roots, then either 3 roots are bigger than e/a or only 1 root is bigger than e/a.

(ii) If [PHI] has 3 distinct real roots, then one root must be repeated.

(iii) If [PHI] has 2 distinct real roots, then one root must be bigger than e/a.

In all cases, we always have at least one root bigger than e/a; denote it by [x.sub.0]. Then [x.sub.0] > e/a, [PHI]([x.sub.0]) = 0, and [PHI]'([x.sub.0]) < 0. Let [b.sub.0] = 1 + [x.sub.0]; since

H'(x) = re/[a.sup.2] x (x - 1)(ax - a + re) [PHI]'(x - 1) - (4ax - 4a + 3re) [PHI](x - 1)/[(x - 1).sup.4][(ax - a + re).sup.2], (43)

we get H'([b.sub.0]) = H'(1 + [x.sub.0]) = (re/[a.sup.2])([PHI]'([x.sub.0])/[x.sup.3.sub.0](a[x.sub.0] + re)) < 0. As H'(b) is continuous, there is [[delta].sub.1] >0 that can be made smaller than [mathematical expression not reproducible] so that H'(b) < 0 in ([b.sub.0] - [[delta].sub.1]], [b.sub.0] + [[delta].sub.1]). Thus H(b) is monotonically decreasing in this interval. Thus, we have proved the following property.

Property 7. The function H(b) has at least 2 real roots, one of which is bigger than [mathematical expression not reproducible], and the other is less than [mathematical expression not reproducible]. Among all roots bigger than [mathematical expression not reproducible], there is a root [b.sub.0] and a neighborhood of [b.sub.0], ([b.sub.0] - [[delta].sub.1], [b.sub.0] + [[delta].sub.1]) where [mathematical expression not reproducible] so that H'([b.sub.0]) < 0 and H(b) is decreasing in this interval.

Define [mathematical expression not reproducible]. All these sets are nonempty and [I.sub.0] has either at least 1 element or at most 3 elements. Let [mathematical expression not reproducible]. Note that [mathematical expression not reproducible]. It is easy to check that when [mathematical expression not reproducible]

Next, we refine the condition [bar.y] < N, where N = n/m. This inequality is equivalent to

N[(a(b - 1)).sup.2] + re(N - 1)(a(b - 1)) + r[e.sup.2] > 0. (44)

Considering it as a quadratic polynomial of a(b - 1),we have [DELTA] = r[e.sup.2][r[(N - 1).sup.2] - 4N]. If [DELTA] < 0, then this inequality is always true for all N, so [bar.y] < N is always true. If [DELTA] [greater than or equal to] 0, then

(i) when N [greater than or equal to] 1, the right-hand side of the inequality has 2 negative roots [N.sub.1] < [N.sub.2] < 0, but since N > 0, this inequality is obviously true and hence [bar.y] < N is always true;

(ii) when N < 1, the inequality is equivalent to b [member of] (1, 1 + [N.sub.1]/a)[union](1 + [N.sub.2]/a, +[infinity]), where [N.sub.1] and [N.sub.2] are 2 positive roots of the right-hand side of the inequality (they may be equal).

Let J = (1, 1 + [N.sub.1]/a) [union] (1 + [N.sub.2]/a, + [infinity]). Then we have the following result: if ([mathematical expression not reproducible] and b is in this intersection, then [E.sub.2] is locally asymptotically stable.

Biologically, when the viral burst size is becoming larger and between two critical values, Proposition 6 says that the virotherapy will reach a stable state which is free of innate immune cells. It is reasonable that these two critical burst sizes are related to the relative immune response decay rate. Actually, in order to have this equilibrium, it requires that the relative immune response decay rate is small. In other words, when the relative immune response decay rate is small and the viral burst size is becoming larger, the virotherapy can have a partial success where the innate immune system has no effects on the therapy, and tumor cells, infected tumor cells, and viruses coexist.

For positive equilibrium solutions [E.sub.3], [E.sub.4], and [E.sub.5], when they exist, we derive conditions under which they are unstable.

Proposition 8. [E.sub.3] is locally unstable when [bar.b] < [[bar.b].sub.1]. [E.sub.4] and [E.sub.5] are locally unstable when [bar.b] < [[bar.b].sub.2,3].

Proof. When f(v) = 0 has a unique positive root [v.sub.2] = [v.sub.3] = [v.sup.*.sub.2] = A < u*, the system has a unique positive equilibrium solution [E.sub.3] = (1 - N - aA/r, N, A, ((b - 1)N - eA)/(cN + dA)). The variational matrix at [E.sub.3] is

[mathematical expression not reproducible]. (45)

The characteristic polynomial of this matrix is computed as the quartic polynomial p([lambda]) = [[lambda].sup.4] + [b.sub.3][[lambda].sup.3] + [b.sub.2][[lambda].sup.2] + [b.sub.1][lambda] + [b.sub.0], where

[mathematical expression not reproducible]. (46)

Assume that p(0) = [b.sub.0] < 0. Since [lim.sub.[lambda][right arrow][infinity]] p([lambda]) = [infinity], p([lambda]) has at least one positive root. The fact [b.sub.0] < 0 is equivalent to b < ([(cN + dA).sup.2]/cd[N.sup.2])((2[a.sup.2]/r)A - a + aN)- ce/d + 1. On the other hand, we can compute b in terms of coefficients of f(v) and note that the coefficients [a.sub.1], [a.sub.2] do not depend on b. Since f([v.sup.*.sub.2]) = 0, we have [v.sup.*.sub.2] = (9[a.sub.0] - [a.sub.1][a.sub.2])/2([a.sup.2.sub.2] - 3[a.sub.1]). As [v.sup.*.sub.2] = A = (-[a.sub.2] + [square root of ([a.sup.2.sub.2] - 3[a.sub.1])/3)], so (9[a.sub.0] - [a.sub.1][a.sub.2])/2([a.sup.2.sub.2] - 3[a.sub.1]) = (-[a.sub.2] + [square root of ([a.sup.2.sub.2] - 3[a.sub.1])/3)], which implies that

b = [a.sup.2]d/cr[N.sup.2] x -2[a.sup.3.sub.2] + 9[a.sub.1][a.sub.2] + 2[([a.sup.2.sub.2] - 3[a.sub.1]).sup.3/2]/27 =: [bar.b]. (47)

Thus, when [bar.b] < ([(cN + dA).sup.2]/cd[N.sup.2])((2[a.sup.2]/r)A - a + aN) - ce/d + 1 = [[bar.b].sub.1], [E.sub.3] is locally unstable.

Lastly, when f(v) = 0 has two distinct positive real roots 0 < [v.sub.2] < [v.sup.*.sub.2] = A < [v.sub.3] < [u.sup.*], the variational matrices at [E.sub.4] and [E.sub.5] are the same as the variational matrix at [E.sub.3] except that A is replaced by [v.sub.2] and [v.sub.3], respectively. We obtain the corresponding characteristic polynomials of those matrices which are the same as the characteristic polynomial of L except for replacing A by [v.sub.2] and [v.sub.3]. By the same argument as above, when [bar.b] < ([(cN + d[v.sub.2,3]).sup.2]/cd[N.sup.2])((2[a.sup.2]/r)[v.sub.2,3] - a + aN) - ce/d + 1 = [[bar.b].sub.2,3,], [E.sub.4] and [E.sub.5] are locally unstable.

It is interesting to notice that our model system has 3 positive equilibria when the viral burst size is not too large and the relative immune killing rate falls into some intervals. Proposition 8 gives conditions that ensure these equilibrium solutions are unstable. Biologically, when the relative immune killing rate and relative immune response decay rate fall into some ranges, we may genetically change the viral burst size to avoid coexistent equilibria.

2.4. Bifurcation Analysis. The dramatic changes of solutions may occur at bifurcations of parameter values. It is important to study bifurcation phenomena for any mathematical models. For our model (4), there are two types of bifurcations, transcritical bifurcations and Hopf bifurcations. For basic knowledge of Hopf bifurcations, we refer Hassard et al. [22].

A transcritical bifurcation occurs with [E.sub.1] and [E.sub.2]. When [mathematical expression not reproducible], [E.sub.2] is outside of the positive domain D and [E.sub.1] is locally asymptotically stable. As b increases to [mathematical expression not reproducible], [E.sub.2] moves into D and it coalesces with [E.sub.1] which is still locally asymptotically stable. But when [mathematical expression not reproducible], the stability of these equilibrium points interchanges, which means that [E.sub.2] is locally asymptotically stable while [E.sub.1] is unstable. Notice that when b > [b.sub.0] and b [member of] [I.sub.n], [E.sub.2] is locally unstable.

In order to study the Hopf bifurcation at b = [b.sub.0], we look at the characteristic polynomial (38):

p([lambda]) = [absolute value of ([lambda]I - L)] = [[lambda] -(m[bar.y] - n)]q([lambda]), (48)

where q([lambda]) = [[lambda].sup.3] + [a.sub.1][[lambda].sup.2] + [a.sub.2][lambda] + [a.sub.3]. For convenience, we assume that ([mathematical expression not reproducible]. From the derivation of Proposition 6, we know m[bar.y] - n < 0. That is, p([lambda]) has a negative root m[bar.y] - n < 0. Thus, the assumption [mathematical expression not reproducible] reduces the study of the quartic polynomial p([lambda]) to the cubic polynomial q([lambda]).

Consider each coefficient of q([lambda]) as a function of the parameter b. Then

q([lambda]) = [[lambda].sup.3] + [a.sub.1] (b)[[lambda].sup.2] + [a.sub.2] (b)[lambda] + [a.sub.3] (b), (49)

where [a.sub.1](b) = (re + ab - a + abe)/a(b - 1), [a.sub.2](b) = re(be + b - 1) /a[(b - 1).sup.2] + re(ab - a - e)(r - a)/a(b - 1)(ab - a + re), and [a.sup.3](b) = re(ab - a- e)/a(b - 1).

The following theorem is our main result for occurring a Hopf bifurcation around [b.sub.0], which appears in [16] as Theorem 3.12. For completion, we restate this theorem and related lemmas and corollary and give slight different proofs.

Theorem 9. There exists a neighborhood of [b.sub.0], ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]), such that for each b in this neighborhood q([lambda]) has a pair of complex conjugate eigenvalues X[lambda](b) = [alpha](b) [+ or -] i[beta](b), where [alpha](b) changes sign when b passes through [b.sub.0] and [beta](b) > 0. Moreover, when b = [b.sub.0], [alpha]([b.sub.0]) = 0 and [alpha]'([b.sub.0]) [not equal to] 0. Notice that [[delta].sub.0] can be made small enough so that ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]) [subset or equal to] J.

To prove this theorem, we need two lemmas about the properties of roots of cubic equations which appear in [16] as Lemmas 3.10 and 3.11.

Lemma 10. The cubic polynomial [[lambda].sup.3] + [a.sub.1] [[lambda].sup.2] + [a.sub.2] [lambda] + [a.sub.3] = 0 with real coefficients has a pair of pure imaginary roots if and only if [a.sub.2] > 0 and [a.sub.3] = [a.sub.1][a.sub.2]. When it has pure imaginary roots, these imaginary roots are given by [lambda] = [+ or -]i [square root of [a.sub.2]], the real root is given by [lambda] = -[a.sub.1], and [a.sub.1][a.sub.3] > 0.

Proof. Suppose the cubic polynomial has a pair of complex roots [lambda] = u[+ or -]vi and one real root [lambda] = [[lambda].sub.0]. Then ([[lambda].sub.2] - 2u[lambda] + [u.sup.2] + [v.sup.2])([lambda] - [[lambda].sub.0]) = [[lambda].sup.3] + [a.sub.1][[lambda].sup.2] + [a.sub.2][lambda] + [a.sub.3]. Expanding the left-hand side and then equating both sides give [a.sub.1] = -([[lambda].sub.0] + [2.sub.u]), [a.sub.2] = [u.sup.2] + [v.sup.2] + 2u[[lambda].sub.0], and [a.sub.3] = - ([u.sup.2] + [v.sup.2])[[lambda].sub.0]. This implies that [[lambda].sub.0] = -([a.sub.1] + 2u), [u.sup.2] + [v.sup.2] = [a.sub.3]/([a.sub.1] + 2u), and [a.sub.3]/([a.sub.1] + 2u) - 2u([a.sub.1] + 2u) = [a.sub.2]. The last equation yields 2([a.sub.2] + [([a.sub.1] + 2u).sup.2])u = [a.sub.3] - [a.sub.1][a.sub.2]. Thus, u = 0 if and only if [a.sub.2] > 0 and [a.sub.3] = [a.sub.1][a.sub.2]. If u = 0, then [[lambda].sub.0] = -[a.sub.1] and [v.sup.2] = [a.sub.2], which follows that [v.sup.2][a.sub.1] = [a.sub.3].

Lemma 11. Consider polynomial [[lambda].sup.3] + [a.sub.1]([tau])[[lambda].sup.2] + [a.sub.2]([tau])[lambda] + [a.sub.3]([tau]) = 0, where [a.sub.k]([tau]) [member of] [C.sup.1], k = 1, 2, 3. Let [lambda]([tau]) = [alpha]([tau]) + i[beta]([tau]) be the roots of the polynomial. Suppose there is [[tau].sub.0] such that a([[tau].sub.0]) = 0 and [beta]([[tau].sub.0]) [not equal to] 0. Moreover, if [mathematical expression not reproducible], then [a'.sub.2]([[tau].sub.0])[a.sub.3]([[tau].sub.0]) = [a.sub.2]([[tau].sub.0])[[a'.sup.3]([[tau].sub.0]) - [a.sub.2]([[tau].sub.0])[a'.sub.1]([[tau].sub.0])].

Proof. Differentiating the polynomial with respect to [tau] and evaluating the derivative at [[tau].sub.0], we notice that [alpha]([[tau].sub.0]) = [alpha]'([[tau].sub.0]) = 0, then we get, after equating the real part and the imaginary part, [beta]'([[tau].sub.0]) = [a'.sub.2]([[tau].sub.0])[beta]([[tau].sub.0]/3[[beta].sup.2]([[tau].sub.0]) - [a.sub.2]([[tau].sub.0])) = ([a'.sub.3]([[tau].sub.0]) - [[beta].sup.2]([[tau].sub.0])[a'.sub.1][([[tau].sub.0]))/2[a.sub.1]([[tau].sub.0 ])[beta]([[tau].sub.0]). By Lemma 10, since [[beta].sup.2]([[tau].sub.0]) = [a.sub.2]([[tau].sub.0]) = [a.sub.3]([[tau].sub.0])/[a.sub.1]([[tau].sub.0]), we obtain the desired result.

From the proofs of Lemmas 10 and 11, and Routh-Hurwitz Criterion, we have the following corollary about q([lambda]).

Corollary 12. There is a neighborhood of [b.sub.0], ([b.sub.0] - [[delta].sub.2], [b.sub.0] + [[delta].sub.2]), where [mathematical expression not reproducible], such that [a.sub.2](b) > 0 for all b [member of] ([b.sub.0] - [S.sub.2], [b.sub.0] + [[delta].sub.2]), and the real part [alpha](b) of the root [lambda](b) = [alpha](b) + i[beta](b) of q([lambda]) changes sign when b passes through [b.sub.0].

Proof. Since b > 1 + e/a, [a.sub.1](b) > 0 and [a.sub.3](b) > 0. As H([b.sub.0]) = [a.sub.1]([b.sub.0])[a.sub.2]([b.sub.0]) - [a.sub.3]([b.sub.0]) = 0, [a.sub.2]([b.sub.0]) = [a.sub.3]([b.sub.0])/[a.sub.1]([b.sub.0]) > 0. Since [a.sub.2](b) is continuous with respect to b, there is a neighborhood of [b.sub.0] such that [a.sub.2](b) > 0 in that neighborhood. Its radius [[delta].sub.2] can be made small enough so that [mathematical expression not reproducible]. We know that H(b) is decreasing in this neighborhood ([b.sub.0] - [[delta].sub.2], [b.sub.0] + [[delta].sub.2]). If b [member of] ([b.sub.0] - [[delta].sub.2], [b.sub.0]), then [H.sub.2] = H(b) > H([b.sub.0]) = 0. Since [H.sub.1] = [a.sub.1](b) > 0 and [H.sub.3] = [a.sub.3](b)H(b) > 0, by Routh-Hurwitz Criterion [alpha](b) < 0. If b [member of] ([b.sub.0], [b.sub.0] + [[delta].sub.2]), then [H.sub.2] = H(b) < H([b.sub.0]) = 0. Since [a.sub.1](b) > 0, [a.sub.2](b) > 0, and [a.sub.3](b) > 0, from the proof of Lemma 10 we have [alpha](b) = -H(b)/2([a.sub.2](b) + [([a.sub.1](b) + 2[alpha](b)).sup.2]) > 0 and [alpha]([b.sub.0]) = 0. Thus the sign of [alpha](b) changes when b passes through [b.sub.0].

We now can prove our main Theorem 9.

Proof. By Property 7 in previous section, [b.sub.0] > 1 + e/a and H([b.sub.0]) = 0. Then, for b > 1 + e/a, [a.sub.1](b) > 0 and [a.sub.3](b) > 0; hence [a.sub.2]([b.sub.0]) = [a.sub.3]([b.sub.0))/[a.sub.1]([b.sub.0]) > 0. By Lemma 10, p([lambda]) has a pair of purely imaginary roots, [lambda]([b.sub.0]) = [+ or -]i[beta]([b.sub.0]) = [+ or -]i[square root of ([a.sub.2]([b.sub.0]))], and the real root is -[a.sub.1]([b.sub.0]) < 0. Since [beta]([b.sub.0]) > 0 and [beta](b) is continuous, we can find a neighborhood of [b.sub.0] so that [beta](b) > 0 in this neighborhood and its radius can be chosen so that [[delta].sub.0] < min{[[delta].sub.1], [[delta].sub.2]} and ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]) [subset or equal to] J. By the above corollary, in the interval ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]), [alpha](b) changes sign when b passes through [b.sub.0]. Finally, we claim that [alpha]'([b.sub.0]) [not equal to] 0. By way of contradiction, if [alpha]'([b.sub.0]) = 0, then by Lemma 11 we have [a'.sub.2]([b.sub.0])[a.sub.3]([b.sub.0]) = [a.sub.2]([b.sub.0])([a'.sub.3]([b.sub.0]) - [a'.sub.1]([b.sub.0])[a.sub.2]([b.sub.0])). Then this yields H'([b.sub.0]) = H([b.sub.0])[a'.sub.2]([b.sub.0])/[a.sub.2]([b.sub.0]) = 0, a contradiction. This completes the proof.

Combining Proposition 6 and applying this theorem, we can obtain a statement about Hopf bifurcations of our model.

Theorem 13. Assuming [mathematical expression not reproducible], for system (4) [??] = f(X, b) f([E.sub.2], b) = 0 for all [mathematical expression not reproducible]. The variational matrix of f at [E.sub.2], L = ([partial derivative]f/[partial derivative]X)([E.sub.2], b), has 2 strictly negative roots and 2 conjugate complex roots [lambda](b) = [alpha](b) [+ or -] i[beta](b) in the neighborhood ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]) of [b.sub.0], in which [beta](b) > 0, [alpha](b) changes sign when b passes through [b.sub.0], and [alpha]'([b.sub.0]) [not equal to] 0. Consequently, in a neighborhood U of [E.sub.2] and for any b [member of] ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]), the system [??] = f(X, b) has nontrivial periodical solutions in U.

Because we cannot find explicit algebraic expression for [b.sub.0], it is very hard to study the nature of periodical solutions that occur around [E.sub.2] when b is close to [b.sub.0] such as the amplitudes, periods, and their stability. However, we can make some statements about the general properties of these periodical solutions as follows.

(i) If [E.sub.2] is stable but not asymptotically stable, then any solution of system (3) in U is periodical in a surface.

(ii) If [E.sub.2] is asymptotically stable, then there is an asymptotically stable periodical solution [bar.X](t) in U when b is close to [b.sub.0]. Any solution inside [bar.X] will spiral into [E.sub.2] when time is increasing and any solution in U outside [bar.X] will spiral and emerge into [bar.X] when time is increasing.

(iii) If [E.sub.2] is unstable, then there is an asymptotically stable periodical solution [??](t) in U when b is close to [b.sub.0]. Any solution starting at nearby [E.sub.2] will spiral out and emerge into [??] when time is increasing, and any solution in U nearby outside [??] will move away from [??] when time is increasing.

We will use numerical simulations to confirm some of these situations. Lastly, we will not conduct the analysis for the bifurcations arising around the positive equilibrium points [E.sub.4] and [E.sub.5] because their formulas are extremely cumbersome and therefore we will treat this by numerical simulations in the next section.

We close Section 2 with the following theorem that summarizes the main results about the our model and its biological implications.

Theorem 14. The dynamical behaviors of system (4) can be described as follows.

(i) When (r/a)(1/N - 1) < K < 1/(a + e - aN), N < 1 + (1/2)(e/a + 1/r - [square root of ([(e/a + 1/r).sup.2] + 4/r))], and b > [bar.b] with [mathematical expression not reproducible], system (4) has at most 3 equilibrium solutions: [E.sub.0], [E.sub.1], and [E.sub.2]. [E.sub.0] is unstable. [E.sub.1] is globally asymptotically stable if [mathematical expression not reproducible] and unstable if [mathematical expression not reproducible]. [E.sub.2] is locally asymptotically stable if b [member of] [mathematical expression not reproducible].

(ii) When either K [less than or equal to] (r/a)(1/N - 1) or K > 1/(a + e - aN) and b = [bar.b] with g(K) > 0, system (4) has a unique positive equilibrium solution: [E.sub.3]. [E.sub.3] is unstable if [bar.b] < [[bar.b].sub.1]

(iii) When either K [less than or equal to] (r/a)(1/N - 1) or K> 1/(a + e - aN) and b < [bar.b] with g(K) > 0, system (4) has two distinct positive equilibrium solutions: [E.sub.4] and [E.sub.5]. [E.sub.4] and [E.sub.5] are unstable if [bar.b] < [[bar.b].sub.2,3].

(iv) When [mathematical expression not reproducible], there exist a neighborhood ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]) of [b.sub.0] and a neighborhood U of [E.sub.2], such that for any b [member of] ([b.sub.0] - [[delta].sub.0], [b.sub.0] + [[delta].sub.0]), system (4) has nontrivial periodical solutions in U.

Biologically, we have interpreted most parts of this theorem. We may emphasize some biological implications here. If the burst size is smaller than the critical value [mathematical expression not reproducible] and the relative immune decay rate is greater than 1, the virotherapy always fails. If the burst size is greater than [??] and smaller than another critical value [??], the immune free equilibrium is stable; that is, the tumor cells, infected tumor cells, and viruses coexist. When the relative immune killing rate is too small or too big compared to a relative immune survival rate (1/N), according to the burst size, the model can have one more or two more positive equilibria, and these positive equilibria are unstable. This gives a chance for the model to have periodical solutions. That is, the virus cannot eradicate the tumor and the virus, tumor cells, and immune cells fight each other forever.

For positive equilibria, [E.sub.3], [E.sub.4], and [E.sub.5], [E.sub.3] is difficult to obtain in practice because it requires a particular threshold of the burst size. In rat experiments, the virus burst size can be genetically changed as we want, but usually, we can ensure a range of the burst size in the process of genetic modification. [E.sub.4] and [E.sub.5] are unstable if the burst size is smaller than a threshold. Biologically, these two equilibria are not important because of their instability. The immune free equilibrium [E.sub.2] can be stable. If the burst size is made big enough, the tumor cell portion will be small in [E.sub.2]. On the other hand, the model can have periodic solutions. This may provide an opportunity for combining surgery with the phase where the tumor cell portion is in the lowest state.

3. Numerical Study and Discussion

3.1. Numerical Study. In order to demonstrate our analytical results about dynamical behaviors of the model, we use some data of parameter values from our previous research [14] to conduct numerical computations for all dynamical characteristics and perform some numerical simulations by using Matlab. The data of parameter values we used from [14] is recorded in Table 1. We applied the algorithm of the Newton method for finding Hopf bifurcation points [23], and Matlab codes are available upon request.

After nondimensionalization, the parameter values are r = 0.36, a = 0.11, c = 0.48, d = 0.16, e = 0.2, m = 0.6, and n = 0.036. Then [mathematical expression not reproducible]. Solving H(b) = 0 gives 2 conjugate complex roots b = 0.8353 [+ or -] 0.2312i, and 2 real roots b = 0.299 and b = 19.012. Therefore, [I.sub.0] = {19.012}, [I.sub.p] = (2.82, 19.012), [I.sub.n] = (19.012, + [infinity]), and hence [mathematical expression not reproducible]. On the other hand, all coefficients of the cubic equation are [a.sub.2] = -2.8964, [a.sub.1] = 0.1603, and [bar.b] = 7.4455. Considering the case with the burst size b = 9, we find all equilibrium solutions of system (4):

[E.sub.0] = (0, 0, 0, 0),

[E.sub.1] = (1, 0, 0, 0),

[E.sub.2] = (0.2273, 0.0584, 2.3377, 0),

[E.sub.4] = (0.48328, 0.06, 1.49471, 0.67571),

[E.sub.5] = (0.24994, 0.06, 2.25837, 0.07261). (50)

By the analysis of previous section, [E.sub.0] is always unstable. The equilibrium point [E.sub.1] is locally asymptotically stable if b < 2.82 and unstable if b > 2.82. For the equilibrium point [E.sub.2], since [bar.y] = 0.0584 < N = n/m = 0.06 and [mathematical expression not reproducible], it is locally asymptotically stable. In order to check the stability of the equilibrium point [E.sub.4], we need to compute the Jacobian matrix of F at this point, which is L = ([partial derivative]F/[partial derivative]X)([E.sub.4]). Using the Matlab, we can calculate 4 eigenvalues of L, which are -1.69849, -0.00803,and -0.07654 [+ or -] 0.20995i. This guarantees the locally asymptotical stability of [E.sub.4]. For the last equilibrium point [E.sub.5], we know that [v.sub.3] = 2.258366, which is the largest root of f(v) = 0. Then, we can compute the quantity ([(cN + d[v.sub.3]).sup.2] /cd[N.sup.2])((2[a.sup.2]/r)[v.sub.3] - a + aN) - ce/d + 1 = 27.052. As [bar.b] < 27.052, [E.sub.5] is unstable.

When b = 20, the bifurcation analysis gives us the appearance of periodical solutions around the equilibrium point [E.sub.2] = (0.09569, 0.03011, 2.86098, 0). However, in this case, two positive equilibrium points [E.sub.4] and [E.sub.5] do not exist. Now we fix the burst size b = 20 and other parameters, considering m as a variable parameter. Observe that when m = 1.17, we have two positive equilibrium points: [E.sub.4] and [E.sub.5]. The eigenvalues of the Jacobian matrix at [E.sub.5] are -1.26885, 0.00075, and -0.00003 [+ or -] 0.23144i, whereas when m = 1.18, we also have two positive equilibrium points: [E.sub.4] and [E.sub.5]. The eigenvalues at [E.sub.5] are 1.26023, 0.00045, and 0.000455 [+ or -] 0.23025i. This partially shows the existence of the Hopf bifurcation point 1.17 < [m.sub.0] < 1.18. Using the Newton method for the computation of Hopf bifurcation points, we find [m.sub.0] = 1.1706885699.

For the sake of demonstration and simplicity, we conduct numerical simulations based on nondimensionalized model. Therefore, the unit of the tumor cells, infected tumor cells, viruses, and innate immune cells are not absolute numbers and are only relative numbers. For example, the quantity of tumor cells in all figures is the portion of tumor cells over the tumor carrying capacity. Similarly, other quantities have the same meaning. We just indicate them as relative tumor cells and so on in the figures. For the time, it can be considered as runs of infected tumor bursting since [tau] = [delta]t, or simply, relative time.

Figure 1 shows that [E.sub.1] is locally asymptotically stable.

Figure 2 shows [E.sub.1] is locally unstable since b = 9 > 1 + e/a = 2.82. Figure 3 shows [E.sub.2] is locally asymptotically stable because b is between [mathematical expression not reproducible] and y = 0.058 < N = 0.06. Figure 4 shows [E.sub.4] is locally asymptotically stable since all eigenvalues of the variational matrix at [E.sub.4] are negative. Figure 5 shows [E.sub.5] is locally unstable sine [bar.b] < [[bar.b].sub.3].

Figures 6-8 show periodic solutions rising from Hopf bifurcations.

Figure 9 shows the tumor cell population when the burst size b is 1000.

3.2. Discussion. The dynamics of oncolytic virotherapy without the presence of the innate immune response is largely determined by the oncolytic viral burst size as studied by Tian [16]. Specifically, there are two threshold values of the burst size: below the first threshold, the tumor always grows to its maximum (carrying capacity) size; while passing this threshold, there is a locally stable positive equilibrium solution appearing through transcritical bifurcation; while at or above the second threshold, there exits one or three families of periodic solutions arising from Hopf bifurcations. And it also suggests that the tumor load can drop to an undetectable level either during the oscillation or when the burst size is large enough. When the model for oncolytic virotherapy is with the presence of the innate immune response, the dynamics becomes more complicated. There are several critical values for the oncolytic viral burst size b, for example, [mathematical expression not reproducible], where [bar.b] is a function of the relative innate immune response killing rate K and relative innate immune decay rate N, which we combine with innate immune parameters K and N to describe the dynamics of our model (4). When b is smaller than [bar.b] and K and N satisfy some constraints, system (4) has 5 equilibrium solutions and 2 of them are positive, while these two positive equilibrium points coalesce when b = [bar.b] and there are some constraints for the relative innate immune killing rate K and relative innate immune decay rate N. When b is greater than [mathematical expression not reproducible] fulfill the complement conditions, system (4) has at most 3 equilibrium solutions with 0 innate immune components. An interesting fact is that the equilibrium points where Hopf bifurcations arise for both models (4) and in [16] are corresponding to each other. Therefore, we may conclude that the innate immune response complicates the oncolytic virotherapy in the way of creating more equilibrium solutions when the oncolytic viral burst size is not too big, say less than [bar.b], while the dynamics is similar to the system without the presence of the innate immune response when the oncolytic viral burst size is big.

As we mention in the Introduction, the major challenge in current medical practice of oncolytic viral therapy is the complexity of the immune responses [4]. The innate immune response tends to reduce the efficacy of oncolytic viral treatments by reducing new virus multiplication and blocking the spreading of infection. However, the stimulated adaptive immune response tends to reduce tumor cells. These two opposite functions increase the complexity of oncolytic viral therapy. A balance between two functions needs to be determined in order to improve the efficacy of oncolytic virotherapy. This is a very subtle question. There are not too much experimental or clinical data about this balance in the literature. Therefore, a mathematical study of this question is highly demanded. Our current mathematical model is only dealing with the innate immune system. The extension of our model to incorporate the adaptive immune system is expected. We plan to carry out such study in other articles.

https://doi.org/10.1155/2017/6587258

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

Jianjun Paul Tian would like to thank E. Antonio Chiocca for his suggestions and discussion about the model and results in this study. Both authors would like to acknowledge the support by NSF (USA) Grant DMS-1446139, NIH (USA) Grant U54CA132383, and NNSF (China) Grant 11371048.

References

[1] E. A. Chiocca, "Oncolytic viruses," Nature Reviews Cancer, vol. 2, no. 12, pp. 938-950, 2002.

[2] E. Kelly and S. J. Russell, "History of oncolytic viruses: genesis to genetic engineering," Molecular Therapy, vol. 15, no. 4, pp. 651-659, 2007.

[3] R. H. I. Andtbacka, H. L. Kaufman, F. Collichio et al., "Talimogene laherparepvec improves durable response rate in patients with advanced melanoma," Journal of Clinical Oncology, vol. 33, no. 25, pp. 2780-2788, 2015.

[4] E. A. Chiocca and S. D. Rabkin, "Oncolytic viruses and their application to cancer immunotherapy," Cancer Immunology Research, vol. 2, no. 4, pp. 295-300, 2014.

[5] J. T. Wu, H. M. Byrne, D. H. Kirn, and L. M. Wein, "Modeling and analysis of a virus that replicates selectively in tumor cells," Bulletin of Mathematical Biology, vol. 63, no. 4, pp. 731-768, 2001.

[6] L. M. Wein, J. T. Wu, and D. H. Kirn, "Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: Implications for virus design and delivery," Cancer Research, vol. 63, no. 6, pp. 1317-1324, 2003.

[7] D. Wodarz, "Viruses as antitumor weapons: defining conditions for tumor remission," Cancer Research, vol. 61, no. 8, pp. 3501-3507, 2001.

[8] D. Wodarz, "Computational approaches to study oncolytic virus therapy: insights and challenges," Gene TherMol Biol, vol. 8, pp. 137-146, 2004.

[9] N. L. Komarova and D. Wodarz, "ODE models for oncolytic virus dynamics," Journal of Theoretical Biology, vol. 263, no. 4, pp. 530-543, 2010.

[10] D. Wodarz and N. Komarova, "Towards predictive computational models of oncolytic virus therapy: Basis for experimental validation and model selection," PLoS ONE, vol. 4, no. 1, Article ID e4271, 2009.

[11] A. S. Novozhilov, F. S. Berezovskaya, E. V. Koonin, and G. P. Karev, "Mathematical modeling of tumor therapy with oncolytic viruses: regimes with complete tumor elimination within the framework of deterministic models," Biology Direct, vol. 1, article 6, pp. 1-18, 2006.

[12] Z. Bajzer, T. Carr, K. Josic, S. J. Russell, and D. Dingli, "Modeling of cancer virotherapy with recombinant measles viruses," Journal of Theoretical Biology, vol. 252, no. 1, pp. 109-122, 2008.

[13] M. Biesecker, J.-H. Kimn, H. Lu, D. Dingli, and Z. Bajzer, "Optimization of virotherapy for cancer," Bulletin of Mathematical Biology, vol. 72, no. 2, pp. 469-489, 2010.

[14] A. Friedman, J. P. Tian, G. Fulci, E. A. Chiocca, and J. Wang, "Glioma virotherapy: Effects of innate immune suppression and increased viral replication capacity," Cancer Research, vol. 66, no. 4, pp. 2314-2319, 2006.

[15] D. Vasiliu and J. P. Tian, "Periodic solutions of a model for tumor virotherapy," Discrete and Continuous Dynamical Systems - Series S, vol. 4, no. 6, pp. 1587-1597, 2011.

[16] J. P. Tian, "The replicability of oncolytic virus: defining conditions in tumor virotherapy," Mathematical Biosciences and Engineering, vol. 8, no. 3, pp. 841-860, 2011.

[17] Y. Wang, J. P. Tian, and J. Wei, "Lytic cycle: A defining process in oncolytic virotherapy," Applied Mathematical Modelling, vol. 37, no. 8, pp. 5962-5978, 2013.

[18] J. P. Tian, Y. Kuang, and H. Yang, "Intracellular viral life-cycle induced rich dynamics in tumor virotherapy," Canadian Appl. Math. Quart.

[19] B. S. Choudhury and B. Nasipuri, "Efficient virotherapy of cancer in the presence of immune response," International Journal of Dynamics and Control, vol. 2, no. 3, pp. 314-325, 2014.

[20] L. J. S. Allen, An Introduction to Mathematical Biology, Pearson, 2006.

[21] J. Carr, Applications of Centre Manifold Theory, Applied Mathematics Sciences, Springer, New York, NY, USA, 1981.

[22] B. D. Hassard, N. D. Kazarinoff, and Y.-H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, 1981.

[23] D. Roose and V. Hlavacek, "A direct method for the computation of Hopf bifurcation points," SIAM Journal on Applied Mathematics, vol. 45, no. 6, pp. 879-894, 1985.

Tuan Anh Phan and Jianjun Paul Tian

Department of Mathematical Sciences, New Mexico State University, Las Cruces, NM

88001, USA

Correspondence should be addressed to Jianjun Paul Tian; jtian@nmsu.edu

Received 4 July 2017; Revised 16 October 2017; Accepted 6 November 2017; Published 12 December 2017

Academic Editor: Jan Rychtar

Caption: FIGURE 1: Dynamics of the model when b = 2 and initial values are x = 0.9, y = 0.01, v = 0.01, and z = 0.01.

Caption: FIGURE 2: Dynamics of the system when b = 9 and initial values are x = 0.99, y = 0.01, v = 0.01, and z = 0.01.

Caption: FIGURE 3: Dynamics of the system when b = 9 and initial values are x = 0.227, y = 0.058, v = 2337, and z = 0.01.

Caption: FIGURE 4: Dynamics of the system when b = 9 and initial values are x = 0.483, y = 0.059, v = 1.494, and z = 0.675.

Caption: FIGURE 5: Dynamics of the system when b = 9 and initial values are x = 0.249, y = 0.059, v = 2.258, and z = 0.072.

Caption: FIGURE 6: Periodic solutions from Hopf bifurcation when m= 1.17. The initial values are x = 0.0998, y = 0.0307, v = 2.8451, and z = 0.0331.

Caption: FIGURE 7: Periodic solutions from Hopf bifurcation when m = 1.18. The initial values are x = 0.0981, y = 0.0305, v = 2.8515, and z = 0.0198.

Caption: FIGURE 8: Periodic solutions from Hopf bifurcation.

Caption: FIGURE 9: The tumor cell population when b = 1000. The initial values are x = 0.5, y = 0.5, v = 1.5, and z =1.

TABLE 1: Parameters and their values. Parameters Description [lambda] Tumor growth rate [theta] Density of tumor cells [beta] Infection rate of the virus [mu] Immune killing rate of infected tumor cells [delta] Death rate of infected tumor cells b Burst size of free virus k Immune killing rate of virus [gamma] Clearance rate of the virus s Stimulation rate of the virus by infected cells [rho] Immune clearance rate Parameters Values Dimensions [lambda] 2 x [10.sup.-2] 1/h [theta] [10.sup.6] cells/[mm.sup.3] [beta] 7/10 x [10.sup.-9] [mm.sup.3]/h virus [mu] 2 x [10.sup.-8] [mm.sup.3]/h immune cell [delta] 1/18 1/h b 50 viruses/cell k [10.sup.-8] [mm.sup.3]/h immune cell [gamma] 2.5 x [10.sup.-2] 1/h s 5.6 x [10.sup.-7] [mm.sup.3]/h infected cell [rho] 20 x [10.sup.-8] [mm.sup.3]/h cell

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Phan, Tuan Anh; Tian, Jianjun Paul |

Publication: | Computational and Mathematical Methods in Medicine |

Article Type: | Report |

Date: | Jan 1, 2017 |

Words: | 13554 |

Previous Article: | Fast Parabola Detection Using Estimation of Distribution Algorithms. |

Next Article: | A Fusion-Based Approach for Breast Ultrasound Image Classification Using Multiple-ROI Texture and Morphological Analyses. |

Topics: |