The Role of the Innate Immune System in Oncolytic Virotherapy.

The complexity of the immune responses is a major challenge in current virotherapy. This study incorporates the innate immune response into our basic model for virotherapy and investigates how the innate immunity affects the outcome of virotherapy. The viral therapeutic dynamics is largely determined by the viral burst size, relative innate immune killing rate, and relative innate immunity decay rate. The innate immunity may complicate virotherapy in the way of creating more equilibria when the viral burst size is not too big, while the dynamics is similar to the system without innate immunity when the viral burst size is big.

1. Introduction

Oncolytic 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: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:where λ is tumor growth rate, C is the carrying capacity of tumor cell population, β is the infected rate of the virus, μ is immune killing rate of infected tumor cells, δ 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, γ is clearance rate of viruses, s is the stimulation rate of the innate immune system, and ρ 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 μ/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 ρ/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 b, b, b, , and b0, where is a function of the relative innate immune killing rate K and relative innate immune response decay rate N. When b is smaller than 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 . When b is greater than b or 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 , 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 τ = δt, , , , and rename parameters r = λ/δ, a = Cβ/δ, c = μC/δ, d = kC/δ, e = γ/δ, m = sC/δ, and n = ρ/δ. Then system (2) becomesFor convenience, dropping all the bars and writing τ as t, we obtainWe assume that all parameters are nonnegative.

Lemma 1 .

If x(0) ≥ 0, y(0) ≥ 0, v(0) ≥ 0, and z(0) ≥ 0, then the solution of system (4) x(t) ≥ 0, y(t) ≥ 0, v(t) ≥ 0, and z(t) ≥ 0 for t ≥ 0. Furthermore, if 0 ≤ x(0) + y(0) ≤ 1, then 0 ≤ x(t) + y(t) ≤ 1 for t ≥ 0, and lim supv(t) ≤ 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 ≤ x(0) + y(0) ≤ 1, then 0 ≤ x(t) + y(t) ≤ 1 and lim supv(t) ≤ b/e. In fact, since x(t), y(t), and v(t) are nonnegative for all t ≥ 0, As x(0) ≤ x(0) + y(0) ≤ 1, by the comparison theorem we get x(t) ≤ 1. On the other hand, since again by the comparison theorem we also have 0 ≤ x(t) + y(t) ≤ 1. It follows that 0 ≤ y(t) ≤ 1. Then v′(t) = by − axv − dvz − ev ≤ b − ev, so by the comparison theorem v(t) ≤ b/e + v(0)exp⁡(−et). Taking lim sup both sides yield lim supv(t) ≤ b/e.

We then conclude that the positive invariant domain of system (4) is 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), , and with b > 1 + e/a, system (4) has 3 equilibrium solutions: E0, E1, and E2. When either K ≤ (r/a)(1/N − 1) or K > 1/(a + e − aN) and with g(K) > 0, system (4) has a unique positive equilibrium solution: E3. When either K ≤ (r/a)(1/N − 1) or K > 1/(a + e − aN) and with g(K) > 0, system (4) has two distinct positive equilibrium solutions: E4 and E5.

In what follows we will analyze the existence of equilibrium solutions. First, let X = (x, y, v, z) and Then system (4) can be simply written as the autonomous system dX/dt = F(X). We assume that (x, y, v, z) ∈ D. The equilibrium points are solutions of the equation F(X) = 0; that is,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 ≠ 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 E1 = (1,0, 0,0) is an equilibrium point.

Now if v ≠ 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 Since axv = y, we have v = y/ax = r(ab − a − e)/a(ab − a + re). Thus we get the third critical point Notice that in order for the first three coordinates of E2 to be positive, it is enough that ab − a − e > 0; that is, b > e/a + 1.

Next, if z ≠ 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)≕v3 + a2v2 + a1v + a0 = 0, where Since f(0) = a0 > 0 and limf(v) = −∞, f(v) has at least one negative root, say v1 < 0. Assume that v1 is the least root. If a1 > 0 and a2 > 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 a1 > 0 and a2 > 0 are equivalent to (r/a)(1/N − 1) < c/d < 1/(a + e − aN) and . In this case, the system only has 3 equilibrium points: E0, E1, and E2.

We now look at other situations of f(v). Taking derivative, f′(v) = 3v2 + 2a2v + a1. By long division, Suppose that a22 − 3a1 > 0, then f′ has 2 distinct roots, v1 and v2, where . If v2 > 0 and f(v2) = 0, then f has one negative root, v1, and one positive root, v2 = v3 = v2 = (9a0 − a1a2)/2(a22 − 3a1). In this case, in addition to the 3 equilibrium points E0, E1, and E2, system (4) has one positive equilibrium point: E3 = (1 − N − aA/r, N, A, ((b − 1)N − eA)/(cN + dA)). To guarantee all four coordinates of E3 are positive, we need to impose 1 − N − aA/r > 0 and (b − 1)N − eA > 0, which imply that v2 = A < u = min⁡{(1 − N)(r/a), ((b − 1)/e)N}.

On the other hand, if v2 > 0 and f(v2) < 0, then f has one negative root, v1, and 2 distinct positive roots, 0 < v2 < v2 < v3 < u. Hence system (4) has two positive equilibrium points: E4,5 = (1 − N − av2,3/r, N, v2,3, ((b − 1)N − ev2,3)/(cN + dv2,3)). Notice that v2 = A > 0 is equivalent to either a2 ≤ 0, a22 − 3a1 > 0 or a2 > 0 > a1. Furthermore, a2 ≤ 0, a22 − 3a1 > 0 are equivalent to c/d ≤ (r/a)(1/N − 1) and g(c/d) > 0, where g(x) = x2 + (r/aN − r/a + 3re/a2N)x + (r2/a2)(1/N − 1)2 − 3r/a2N; a2 > 0 > a1 is equivalent to c/d > max⁡{(r/a)(1/N − 1), 1/(a + e − aN)}. The condition f(v2) ≤ 0 is equivalent to v2 ≥ (9a0 − a1a2)/2(a22 − 3a1). That is, , and we have . The critical value is important for describing the dynamics of system (4).

We summarize the details of the analysis above as follows.

When x = 0, we have equilibrium solution E0 = (0,0, 0,0).

When x ≠ 0, we have the following cases.

If z = 0, then

when v = 0, we have equilibrium solution E1 = (1,0, 0,0).

when v ≠ 0 and b > 1 + e/a, we have

If z ≠ 0, then x = 1 − N − (a/r)v, y = N, z = ((b − 1)N − ev)/(cN + dv), and v satisfies the following cubic equation:

where a2 = (c/d)N + (r/a)(N − 1), a1 = rN/a2 · (c/d)(−a + aN − e + d/c), a0 = b · rN2/a2 · c/d.

In this case, we can conclude the following.

If K ∈ ((r/a)(1/N − 1), 1/(a + e − aN)), , and b > 1 + e/a, then system (4) has three equilibrium points: E0, E1, and E2.

If either K ≤ (r/a)(1/N − 1), g(K) > 0, or K > 1/(a + e − aN) and , then system (4) has a unique positive equilibrium point:

where v2 = A < u≔min⁡{(1 − N)(r/a), (b − 1)N/e} and

If either K ≤ (r/a)(1/N − 1), g(K) > 0, or K > 1/(a + e − aN) and , then system (4) has two distinct positive equilibrium points:

where v2 and v3 are two distinct positive real roots of the above cubic equation that satisfy 0 < v2 < v2 < v3 < u, in which .

In order to interpret our mathematical conditions biologically, we need to understand the combined parameters N and K first. N = ρ/sC can be considered as a relative immune response decay rate since ρ is innate immune cell death rate, s is innate immune stimulating rate by infection, and C is tumor carrying capacity. K = c/d = μ/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 E0. 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 E1 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 E2 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 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 E0 is locally unstable, E1 is locally asymptotically stable if b < 1 + e/a and unstable if b > 1 + e/a, and E2 is locally asymptotically stable if b is in some range, while E3, E4, and E5 are locally unstable when b, K, and N satisfy some conditions. We use Lyapunov functions to show E1 is globally asymptotically stable if b < 1 + e/a and N > 1. We apply the center manifold theorem to show E1 is locally asymptotically stable if b = 1 + e/a. We summarize the main results in Lemma 3. For the combined parameter values, b, i = 1,2, J, , 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 0 is unstable. E1 is globally asymptotically stable when b < 1 + e/a and N > 1 and unstable when b > 1 + e/a. E2 is locally asymptotically stable when b ∈ (b, b)∩J. E3 is unstable when . E4 and E5 are unstable when .

We look at the stability of trivial equilibrium solutions first. The variational matrix of system (4) is given by At the critical point E0, the variational matrix is The corresponding eigenvalues are r, −1, −e, and −n. We know that the local stability of E0 of system (4) is the same as that of the linearized system at E0. Since r > 0, E0 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 E0 without viruses and infected tumor cells.

Proposition 4 .

The equilibrium solution E1 is locally asymptotically stable when b < 1 + e/a, and it is locally unstable when b > 1 + e/a. E1 is globally asymptotically stable when b < 1 + e/a and N > 1.

At the equilibrium point E1, the variational matrix is The characteristic polynomial of this matrix is (λ + n)(λ + r)(λ2 + (1 + a + e)λ + a + e − ab). The eigenvalues are λ1 = −r, λ2 = −n, and . Since the eigenvalues λ1 = −r, λ2 = −n, and are negative for all positive parameters, E1 is locally asymptotically stable if and only if λ4 < 0. This is equivalent to , which is the same as b < 1 + e/a. Similarly, if b > 1 + e/a, then , and hence E1 is unstable.

In fact, we can show that E1 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 E1 into the origin by setting , , , and . After dropping all the bars over variables, system (4) becomeswhile the domain D is translated to D1 = {(x, y, v, z) : 0 ≤ x ≤ 1, y ≥ 0, v ≥ 0, z ≥ 0, 0 ≤ x − y ≤ 1}. For any initial condition (x0, y0, v0, z0) ∈ D1, the solution of (22) satisfies 0 ≤ x(t) ≤ 1, 0 ≤ y(t) ≤ 1, v(t) ≥ 0, and z(t) ≥ 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 V1(x, y, v, z) = y + v. It is clear that V1(x, y, v, z) > 0, and the orbital derivative dV1/dt = dy/dt + dv/dt = (b − 1)y − cyz − dvz − ev < 0. When 1 ≤ b < 1 + e/a, consider the Lyapunov function V2(x, y, v, z) = (1/2)ab(a + e)y2 + a2byv + (1/2)a2v2. Obviously, V2(x, y, v, z) > 0. The orbital derivative along a solution is given by 1 ≤ b < 1 + e/a; that is, ab − a − e < 0 and 1 − b ≤ 0; therefore dV2/dt < 0. Combining these two Lyapunov functions gives y(t) → 0 and v(t) → 0 as t → ∞ when b < 1 + e/a. Considering the component x(t), we have By the comparison theorem, Taking limit of both sides as t → ∞ and using the L'Hospital's Rule and the fact that y(t) and v(t) approach 0 yield x(t) → 0. Finally, since y(t) ≤ 1, we have dz/dt ≤ (m − n)z. By the comparison theorem, 0 ≤ z(t) ≤ z(0)e( → 0 as t → ∞, since m − n < 0. Therefore system (3) has a global attractor E1.

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 E1 when b = 1 + e/a. This is a critical case, since the linearized system at E1 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 E1 is locally asymptotically stable when b = 1 + e/a.

Consider b = 1 + e/a, which implies that ab = a + e. The linearized system at E1 has three negative eigenvalues and one zero eigenvalue. In order to determine the stability of E1, 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 which has eigenvalues −r, −(1 + ab), 0, and −n. The associated eigenvectors of these eigenvalues, respectively, are V1 = (1,0, 0,0), V2 = (ab − r, 1 + ab − r, −b(1 + ab − r), 0), V3 = (ar + a, ar, r, 0), and V4 = (0,0, 0,1). System (22) can be written as dX/dt = LX + F, where F = (rx2 − rxy − axv, −axv − cyz, axv − dvz, myz). Set T = (V1, V2, V3, V4) and X = TY; thenwhere and Y = (y1, y2, y3, y4) is determined by Denote T−1F = (f1, f2, f3, f4); then we can express f, i = 1,2, 3,4, in terms of y: where A, B, C, and D are coefficients that can be easily determined. The transformed system can be expressed aswhere It is easy to check that each f, i = 1,2, 3,4, is twice differentiable function, f(0,0, 0,0) = 0 and Df(0,0, 0,0) = 0, where Df is the Jacobian matrix of the function f. By the center manifold theorem, there exists a center manifold:with h(0) = 0 and Dh(0) = 0, and it satisfiesSince h(0) = 0 and Dh(0) = 0, we can assume that here we use the variable u instead of y3 for simplicity. Then we can compute By substituting f's into (34) and equating both sides of the equation, we can get C33 = −a2(r + 1)(b − 1)/(ab + 1) < 0, since b = 1 + e/a > 1. The asymptotically behavior of zero solution of system (31) is governed by that of the equation dy3/dt = f3(h(y3), y3) or dy3/dt = C33y32 + O(y33). Since C33 < 0, y3 = 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 , where , , , and . For lately defined b and J, we have a proposition as follows.

Proposition 6 .

When b ∈ (b, b)∩J ≠ ∅, E2 is locally asymptotically stable.

We show this proposition as follows. The variational matrix at E2 is given by The characteristic polynomial of L iswhere q(λ) = λ3 + a1λ2 + a2λ + a3, and By Routh-Hurwitz Criterion, all roots of q(λ) have negative real parts if and only ifSince b > b≕1 + e/a, a1 > 0 and a3 > 0. And those conditions in (40) are equivalent to H2 = a1a2 − a3 > 0. This inequality is the same as Therefore, we can conclude that if and φ(b) < r − a; then E2 is locally asymptotically stable. Now we can refine this result by considering H(b) = H2, and It is easy to check that H(b) = reΦ(b − 1)/a2(b − 1)3(ab − a + re). Since b > 1 + e/a, H(b) and Φ(b − 1) have the same roots. Since Φ(e/a) = e3(1 + r + e + a)(1 + e + a)((1 + r)/a) > 0 and limΦ(x) = −∞, there are at least one x1 < e/a and one x2 > e/a so that Φ(x1) = Φ(x2) = 0. Then H(b) has at least one root 1 + x1 < 1 + e/a = b and one root 1 + x2 > 1 + e/a = b. We consider three different cases as follows.

If Φ has 4 distinct real roots, then either 3 roots are bigger than e/a or only 1 root is bigger than e/a.

If Φ has 3 distinct real roots, then one root must be repeated.

If Φ 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 x0. Then x0 > e/a, Φ(x0) = 0, and Φ′(x0) < 0. Let b0 = 1 + x0; since we get H′(b0) = H′(1 + x0) = (re/a2)(Φ′(x0)/x03(ax0 + re)) < 0. As H′(b) is continuous, there is δ1 > 0 that can be made smaller than b0 − b so that H′(b) < 0 in (b0 − δ1, b0 + δ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 b = 1 + e/a, and the other is less than b. Among all roots bigger than b, there is a root b0 and a neighborhood of b0, (b0 − δ1, b0 + δ1) where δ1 < b0 − b so that H′(b0) < 0 and H(b) is decreasing in this interval.

Define I = {b > b : H(b) > 0}, I = {b > b : H(b) < 0}, and I0 = {b > b : H(b) = 0}. All these sets are nonempty and I0 has either at least 1 element or at most 3 elements. Let b = min⁡I0. Note that b0 ≥ b. It is easy to check that when b ∈ (b, b), H(b) > 0.

Next, we refine the condition , where N = n/m. This inequality is equivalent to Considering it as a quadratic polynomial of a(b − 1), we have Δ = re2[r(N − 1)2 − 4N]. If Δ < 0, then this inequality is always true for all N, so is always true. If Δ ≥ 0, then

when N ≥ 1, the right-hand side of the inequality has 2 negative roots N1 < N2 < 0, but since N > 0, this inequality is obviously true and hence is always true;

when N < 1, the inequality is equivalent to b ∈ (1,1 + N1/a)∪(1 + N2/a, +∞), where N1 and N2 are 2 positive roots of the right-hand side of the inequality (they may be equal).

Let J = (1,1 + N1/a)∪(1 + N2/a, +∞). Then we have the following result: if (b, b)∩J ≠ ∅ and b is in this intersection, then E2 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 E3, E4, and E5, when they exist, we derive conditions under which they are unstable.

Proposition 8 .

E 3 is locally unstable when . E4 and E5 are locally unstable when .

When f(v) = 0 has a unique positive root v2 = v3 = v2 = A < u, the system has a unique positive equilibrium solution E3 = (1 − N − aA/r, N, A, ((b − 1)N − eA)/(cN + dA)). The variational matrix at E3 is The characteristic polynomial of this matrix is computed as the quartic polynomial p(λ) = λ4 + b3λ3 + b2λ2 + b1λ + b0, where Assume that p(0) = b0 < 0. Since limp(λ) = ∞, p(λ) has at least one positive root. The fact b0 < 0 is equivalent to b < ((cN + dA)2/cdN2)((2a2/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 a1, a2 do not depend on b. Since f(v2) = 0, we have v2 = (9a0 − a1a2)/2(a22 − 3a1). As , so , which implies that Thus, when , E3 is locally unstable.

Lastly, when f(v) = 0 has two distinct positive real roots 0 < v2 < v2 = A < v3 < u, the variational matrices at E4 and E5 are the same as the variational matrix at E3 except that A is replaced by v2 and v3, 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 v2 and v3. By the same argument as above, when , E4 and E5 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 E1 and E2. When b < b, E2 is outside of the positive domain D and E1 is locally asymptotically stable. As b increases to b = 1 + e/a, E2 moves into D and it coalesces with E1 which is still locally asymptotically stable. But when b < b < b and b ∈ I, the stability of these equilibrium points interchanges, which means that E2 is locally asymptotically stable while E1 is unstable. Notice that when b > b0 and b ∈ I, E2 is locally unstable.

In order to study the Hopf bifurcation at b = b0, we look at the characteristic polynomial (38): where q(λ) = λ3 + a1λ2 + a2λ + a3. For convenience, we assume that (b, b0) ⊂ J. From the derivation of Proposition 6, we know . That is, p(λ) has a negative root . Thus, the assumption (b, b0) ⊂ J reduces the study of the quartic polynomial p(λ) to the cubic polynomial q(λ).

Consider each coefficient of q(λ) as a function of the parameter b. Then where a1(b) = (re + ab − a + abe)/a(b − 1), a2(b) = re(be + b − 1)/a(b − 1)2 + re(ab − a − e)(r − a)/a(b − 1)(ab − a + re), and a3(b) = re(ab − a − e)/a(b − 1).

The following theorem is our main result for occurring a Hopf bifurcation around b0, 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 b0, (b0 − δ0, b0 + δ0), such that for each b in this neighborhood q(λ) has a pair of complex conjugate eigenvalues λ(b) = α(b) ± iβ(b), where α(b) changes sign when b passes through b0 and β(b) > 0. Moreover, when b = b0, α(b0) = 0 and α′(b0) ≠ 0. Notice that δ0 can be made small enough so that (b0 − δ0, b0 + δ0)⊆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 λ3 + a1λ2 + a2λ + a3 = 0 with real coefficients has a pair of pure imaginary roots if and only if a2 > 0 and a3 = a1a2. When it has pure imaginary roots, these imaginary roots are given by , the real root is given by λ = −a1, and a1a3 > 0.

Suppose the cubic polynomial has a pair of complex roots λ = u ± vi and one real root λ = λ0. Then (λ2 − 2uλ + u2 + v2)(λ − λ0) = λ3 + a1λ2 + a2λ + a3. Expanding the left-hand side and then equating both sides give a1 = −(λ0 + 2u), a2 = u2 + v2 + 2uλ0, and a3 = −(u2 + v2)λ0. This implies that λ0 = −(a1 + 2u), u2 + v2 = a3/(a1 + 2u), and a3/(a1 + 2u) − 2u(a1 + 2u) = a2. The last equation yields 2(a2 + (a1 + 2u)2)u = a3 − a1a2. Thus, u = 0 if and only if a2 > 0 and a3 = a1a2. If u = 0, then λ0 = −a1 and v2 = a2, which follows that v2a1 = a3.

Lemma 11 .

Consider polynomial λ3 + a1(τ)λ2 + a2(τ)λ + a3(τ) = 0, where a(τ) ∈ C1, k = 1,2, 3. Let λ(τ) = α(τ) + iβ(τ) be the roots of the polynomial. Suppose there is τ0 such that α(τ0) = 0 and β(τ0) ≠ 0. Moreover, if (dα/dτ)| = 0, then a2′(τ0)a3(τ0) = a2(τ0)[a3′(τ0) − a2(τ0)a1′(τ0)].

Differentiating the polynomial with respect to τ and evaluating the derivative at τ0, we notice that α(τ0) = α′(τ0) = 0, then we get, after equating the real part and the imaginary part, β′(τ0) = a2′(τ0)β(τ0)/(3β2(τ0) − a2(τ0)) = (a3′(τ0) − β2(τ0)a1′(τ0))/2a1(τ0)β(τ0). By Lemma 10, since β2(τ0) = a2(τ0) = a3(τ0)/a1(τ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(λ).

Corollary 12 .

There is a neighborhood of b0, (b0 − δ2, b0 + δ2), where δ2 < b0 − b, such that a2(b) > 0 for all b ∈ (b0 − δ2, b0 + δ2), and the real part α(b) of the root λ(b) = α(b) + iβ(b) of q(λ) changes sign when b passes through b0.

Since b > 1 + e/a, a1(b) > 0 and a3(b) > 0. As H(b0) = a1(b0)a2(b0) − a3(b0) = 0, a2(b0) = a3(b0)/a1(b0) > 0. Since a2(b) is continuous with respect to b, there is a neighborhood of b0 such that a2(b) > 0 in that neighborhood. Its radius δ2 can be made small enough so that δ2 < b0 − b and δ2 < δ1. We know that H(b) is decreasing in this neighborhood (b0 − δ2, b0 + δ2). If b ∈ (b0 − δ2, b0), then H2 = H(b) > H(b0) = 0. Since H1 = a1(b) > 0 and H3 = a3(b)H(b) > 0, by Routh-Hurwitz Criterion α(b) < 0. If b ∈ (b0, b0 + δ2), then H2 = H(b) < H(b0) = 0. Since a1(b) > 0, a2(b) > 0, and a3(b) > 0, from the proof of Lemma 10 we have α(b) = −H(b)/2(a2(b)+(a1(b) + 2α(b))2) > 0 and α(b0) = 0. Thus the sign of α(b) changes when b passes through b0.

We now can prove our main Theorem 9.

By Property 7 in previous section, b0 > 1 + e/a and H(b0) = 0. Then, for b > 1 + e/a, a1(b) > 0 and a3(b) > 0; hence a2(b0) = a3(b0)/a1(b0) > 0. By Lemma 10, p(λ) has a pair of purely imaginary roots, , and the real root is −a1(b0) < 0. Since β(b0) > 0 and β(b) is continuous, we can find a neighborhood of b0 so that β(b) > 0 in this neighborhood and its radius can be chosen so that δ0 < min⁡{δ1, δ2} and (b0 − δ0, b0 + δ0)⊆J. By the above corollary, in the interval (b0 − δ0, b0 + δ0), α(b) changes sign when b passes through b0. Finally, we claim that α′(b0) ≠ 0. By way of contradiction, if α′(b0) = 0, then by Lemma 11 we have a2′(b0)a3(b0) = a2(b0)(a3′(b0) − a1′(b0)a2(b0)). Then this yields H′(b0) = H(b0)a2′(b0)/a2(b0) = 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 (b, b0) ⊂ J, for system (4) for all b > b. The variational matrix of f at E2, L = (∂f/∂X)(E2, b), has 2 strictly negative roots and 2 conjugate complex roots λ(b) = α(b) ± iβ(b) in the neighborhood (b0 − δ0, b0 + δ0) of b0, in which β(b) > 0, α(b) changes sign when b passes through b0, and α′(b0) ≠ 0. Consequently, in a neighborhood U of E2 and for any b ∈ (b0 − δ0, b0 + δ0), the system has nontrivial periodical solutions in U.

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

If E2 is stable but not asymptotically stable, then any solution of system (3) in U is periodical in a surface.

If E2 is asymptotically stable, then there is an asymptotically stable periodical solution in U when b is close to b0. Any solution inside will spiral into E2 when time is increasing and any solution in U outside will spiral and emerge into when time is increasing.

If E2 is unstable, then there is an asymptotically stable periodical solution in U when b is close to b0. Any solution starting at nearby E2 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 E4 and E5 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.

When (r/a)(1/N − 1) < K < 1/(a + e − aN), , and with b > b, system (4) has at most 3 equilibrium solutions: E0, E1, and E2. E0 is unstable. E1 is globally asymptotically stable if b < b and N > 1 and unstable if b > b. E2 is locally asymptotically stable if b ∈ (b, b)∩J.

When either K ≤ (r/a)(1/N − 1) or K > 1/(a + e − aN) and with g(K) > 0, system (4) has a unique positive equilibrium solution: E3. E3 is unstable if .

When either K ≤ (r/a)(1/N − 1) or K > 1/(a + e − aN) and with g(K) > 0, system (4) has two distinct positive equilibrium solutions: E4 and E5. E4 and E5 are unstable if .

When (b, b0) ⊂ J, there exist a neighborhood (b0 − δ0, b0 + δ0) of b0 and a neighborhood U of E2, such that for any b ∈ (b0 − δ0, b0 + δ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 b and the relative immune decay rate is greater than 1, the virotherapy always fails. If the burst size is greater than b and smaller than another critical value b, 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, E3, E4, and E5, E3 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. E4 and E5 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 E2 can be stable. If the burst size is made big enough, the tumor cell portion will be small in E2. 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.

Table 1

Parameters and their values.

Parameters	Description	Values	Dimensions
λ	Tumor growth rate	2 × 10−2	1/h
θ	Density of tumor cells	106	cells/mm3
β	Infection rate of the virus	7/10 × 10−9	mm3/h virus
μ	Immune killing rate of infected tumor cells	2 × 10−8	mm3/h immune cell
δ	Death rate of infected tumor cells	1/18	1/h
b	Burst size of free virus	50	viruses/cell
k	Immune killing rate of virus	10−8	mm3/h immune cell
γ	Clearance rate of the virus	2.5 × 10−2	1/h
s	Stimulation rate of the virus by infected cells	5.6 × 10−7	mm3/h infected cell
ρ	Immune clearance rate	20 × 10−8	mm3/h cell

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 b = 1 + e/a = 2.82. Solving H(b) = 0 gives 2 conjugate complex roots b = 0.8353 ± 0.2312i, and 2 real roots b = 0.299 and b = 19.012. Therefore, I0 = {19.012}, I = (2.82,19.012), I = (19.012, +∞), and hence b = b0 = 19.012. On the other hand, all coefficients of the cubic equation are a2 = −2.8964, a1 = 0.1603, and . Considering the case with the burst size b = 9, we find all equilibrium solutions of system (4): By the analysis of previous section, E0 is always unstable. The equilibrium point E1 is locally asymptotically stable if b < 2.82 and unstable if b > 2.82. For the equilibrium point E2, since and b = 9 ∈ (b, b) = (2.82,19.012), it is locally asymptotically stable. In order to check the stability of the equilibrium point E4, we need to compute the Jacobian matrix of F at this point, which is L = (∂F/∂X)(E4). Using the Matlab, we can calculate 4 eigenvalues of L, which are −1.69849, −0.00803, and −0.07654 ± 0.20995i. This guarantees the locally asymptotical stability of E4. For the last equilibrium point E5, we know that v3 = 2.258366, which is the largest root of f(v) = 0. Then, we can compute the quantity ((cN + dv3)2/cdN2)((2a2/r)v3 − a + aN) − ce/d + 1 = 27.052. As , E5 is unstable.

When b = 20, the bifurcation analysis gives us the appearance of periodical solutions around the equilibrium point E2 = (0.09569,0.03011,2.86098,0). However, in this case, two positive equilibrium points E4 and E5 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: E4 and E5. The eigenvalues of the Jacobian matrix at E5 are −1.26885, 0.00075, and −0.00003 ± 0.23144i, whereas when m = 1.18, we also have two positive equilibrium points: E4 and E5. The eigenvalues at E5 are 1.26023, 0.00045, and 0.000455 ± 0.23025i. This partially shows the existence of the Hopf bifurcation point 1.17 < m0 < 1.18. Using the Newton method for the computation of Hopf bifurcation points, we find m0 = 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 τ = δt, or simply, relative time.

Figure 1 shows that E1 is locally asymptotically stable.

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.

Figure 2 shows E1 is locally unstable since b = 9 > 1 + e/a = 2.82. Figure 3 shows E2 is locally asymptotically stable because b is between b and b and y = 0.058 < N = 0.06. Figure 4 shows E4 is locally asymptotically stable since all eigenvalues of the variational matrix at E4 are negative. Figure 5 shows E5 is locally unstable sine .

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.

Figure 3

Dynamics of the system when b = 9 and initial values are x = 0.227, y = 0.058, v = 2.337, and z = 0.01.

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.

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.

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

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.

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.

Figure 8

Periodic solutions from Hopf bifurcation.

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

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.

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, b, b, , and b0, where 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 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 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 b or and K and N 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 , 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.