THE REPLICABILITY OF ONCOLYTIC VIRUS : DEFINING CONDITIONS IN TUMOR VIROTHERAPY

The replicability of an oncolytic virus is measured by its burst size. The burst size is the number of new viruses coming out from a lysis of an infected tumor cell. Some clinical evidences show that the burst size of an oncolytic virus is a defining parameter for the success of virotherapy. This article analyzes a basic mathematical model that includes burst size for oncolytic virotherapy. The analysis of the model shows that 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. The study suggests that the tumor load can drop to a undetectable level either during the oscillation or when the burst size is large enough.

1. Introduction.Oncolytic viruses are genetically altered viruses that can infect and replicate in cancer cells but leave healthy normal cells unharmed.When oncolytic viruses are inoculated into a cancer patient or directly injected into a tumor, they spread throughout the tumor, and infect tumor cells.The viruses that are in the infected tumor cells replicate themselves.Upon a lysis of an infected tumor cell, a swarm of new viruses burst out of the infected cell and infect neighboring tumor cells.Over the last decade, great progress in understanding of the molecular mechanisms of viral cytotoxicity of oncolytic viruses has been providing a fascinating possible alternative therapeutic approach to cancer patients.This alternative approach could be especially beneficial in the case of malignant brain tumor, glioma, since the standard therapy of surgery-radiation-chemotherapy does not typically destroy all the tumor cells; survival rate in high grade glioma is measured in months.Recent experiments in animal brain tumors using genetically engineered viral strains, such as adenovirus, ONYX-15 and CV706, herpes simplex virus 1 and wild-type Newcastle disease virus show these viruses to be relatively non-toxic and tumor specific [4].Additionally, a wide array of viruses is being tested for potential as oncolytic viruses.
Although oncolytic viruses have shown promising outcomes in some clinical trials [1], viral oncolytic therapy has not yet lived up to its expectations.One reason is that once inside a cell the oncolytic virus, such as herpes simplex virus type 1 hrR3, replicate poorly; this difficulty may be overcome with advanced technology.Another reason for the failure of efficacy of virotherapy is because of limited understanding of the dynamics that underlie the spread of oncolytic viruses through tumors.Without such an understanding, much of the clinical research is based on trial and error.In this aspect, mathematical modeling can allow us to see the whole spectrum of possible outcomes, and provide the rationale to optimize treatments.Several attempts have been made to understand and characterize viral dynamics by mathematical models [6,18,12,14,15,16,2,11].These studies are largely of qualitative simulation in nature and examine how variation in viral and host parameters influences the outcome of treatment.The outcome of virotherapy depends in a complex way on interactions between viruses and tumor cells.A clear picture of the dynamics of virotherapy seems hard to obtain.
Experimental results [10] have shown that a novel glioma-selective HSV-1 mutant (rQNestin34.5)has a high replication ability and doubles the life span of glioma mice.Our work [6] has suggested that viral replication ability is a major factor for the success of virotherapy.Our model in [6] describes interactions among glioma cells, infected glioma cells, viruses, and immune cells.It also includes the immunosuppressive agent, cyclophosphamide.It is a free boundary problem with five nonlinear parabolic partial differential equations.Our study of a space-free virotherapy model [13] suggests that the oscillation is an intrinsic property of virotherapy dynamics and relates to virus replication ability.This paper aims to analyze the dynamics of oncolytic viruses related to replication ability in a relatively simpler setting.We will only consider tumor cells and viruses in a space-free fashion.
This paper is organized as follows.In section 2, we review models of the dynamics of oncolytic viruses in terms of ordinary differential equations, and introduce a common basic model.In section 3, we mathematically analyze our basic model.In section 4, using published data, we demonstrate the typical dynamic behaviors of our basic model by numerical simulation, and give a brief discussion and several open problems of virotherapy dynamics.
2. Mathematical models.Within the framework of modeling interacting populations by systems of ordinary differential equations, Wodarz and Komarova [16] proposed a general model (1) based on the law of mass action, where x stands for the uninfected tumor cell population and y the infected tumor cell population.The function F describes the growth properties of the uninfected tumor cells, and the function G describes the rate at which tumor cells become infected by the virus.These two functions could take several different forms, depending on how much detail of the biology is incorporated into the model.The coefficient β represents the infectivity of the virus.The infected tumor cells die with a rate ay.This is a basic framework to describe the dynamics of virotherapy.From our study [6,8], the free virus population itself is very important for virotherapy dynamics.
Actually, Wodarz [17] proposed a three population model (2), Here, v stands for the free virus population.
Here, ε is a constant.In order to consider the recombinant effect of measles virus, Dingli et al [2] proposed the following model ( 4) where ρxy describes cell-cell fusion effect during the virotherapy.This model includes an important fact that tumor cells and free virus both decrease due to infection.
In our PDE model [6] and ODE model [13], we model virus replication ability by its burst size.The burst size of a virus is the number of new viruses released from a lysis of an infected cell.Different types of viruses have different burst sizes.Viruses of the same type have almost the same burst size.The burst size is an important parameter of virus replicability.Including the burst size b, we propose a common basic model for virotherapy: Since tumor cells infected by most types of viruses do not combine each other [4], we do not include recombinant effect into this basic model.Logistic growth is a common growth pattern.We therefore take it as uninfected tumor growth instead of a generalized one.

3.
Analysis of the basic model.The main results can be summarized as follows.
Theorem 3.1.The system (5) has three equilibria, E 0 = (0, 0, 0), ).We will use Lyapunov functions, center manifold theorem to study equilibria.To study Hopf bifurcations, we will develop new properties of coefficient parameterized polynomials and properties of Routh-Hurwitz determinants since the bifurcation values are not algebraically expressible.Since the bust size is a measure of capability of a virus, the main results have the following possible biological implications.When viruses are not powerful, the burst size is less than b s1 , virotherapy fails.When viruses are powerful but not very powerful, the burst size is in the range b s1 < b < b s2 , virotherapy can be partially successful.When viruses are very powerful, virotherapy successes in the sense of tumor load dropping to a undetectable level.
In the rest of this section, we shall present preliminary results, equilibrium analysis, and Hopf bifurcation analysis.We will also highlight biological implications of mathematical statements when available.
3.1.Preliminary results.For simplicity, we non-dimensionalize the system (5) by setting τ = δt, x = K x, y = K ȳ, v = K v, and rename parameters r = λ δ , a = βK δ and c = γ δ .The system (5) becomes For convenience, dropping all bars over the variables and write τ as t, we have It is assumed that all parameters are nonnegative.
Proof.If the conclusion x(t) ≥ 0, y(t) ≥ 0, and v(t) ≥ 0 for t ≥ 0 is not true, there must be a time t 1 , such that there is at least one component that will be zero first.
We check each possible case.
If two components are zero simultaneously at t 1 , it is easy to check that the third component will be nonnegative when t > t 1 .
If the three components are zero simultaneously at t 1 , then from the uniqueness of the solution, x(t) = 0, y(t) = 0, and v(t) = 0 for t ≥ 0.
Therefore, x(t) ≥ 0, y(t) ≥ 0, and v(t) ≥ 0 for t ≥ 0. Since the initial value of each component is nonnegative, the solution is nonnegative.So, Under the given conditions, the solution components x(t) and y(t) are both nonnegative and smaller than 1. v (t) = by − axv − cv ≤ by − cv ≤ b − cv, by the comparison theorem, v(t) ≤ b c + v(0)e −ct .So, sup v(t) ≤ b c + v(0)e −ct , and lim t→+∞ sup v(t) ≤ b c .Define the domain D = {(x, y, v) : x ≥ 0, y ≥ 0, v ≥ 0, 0 ≤ x + y ≤ 1}.From Lemma 3.2 and Lemma 3.3, D is a positive invariant domain for the system (6).It is also a biological meaningful range for variables.We then refer the whole domain D as a "global"domain.
3.2.Equilibrium analysis.The system (6) always has two equilibrium solutions in the positive invariant domain D, E 0 = (0, 0, 0) and a , it has three equilibria, E 0 , E 1 , and a unique positive equilibrium solution The variational matrix of the system ( 6) is given by At the equilibrium point E 0 , the variational matrix is The eigenvalues are r, -1, and -c.Since r is positive, E 0 is unstable.For the linearization of the system ( 6) at E 0 , the stable invariant subspace is the y-v plane, and the unstable invariant subspace is the x axis.It happens that for the system (6), the local stable invariant manifold is in the y-v plane, and the unstable invariant manifold is in the x axis.A possible interpretation is that without viruses and infected tumor cells the tumor will grow from an initial small value around E 0 .If all tumor cells are infected, then the tumor will shrink until it disappears into E 0 .
At the equilibrium point E 1 , the variational matrix is The eigenvalues are λ 1 = −r and Proof.At the equilibrium E 1 , the eigenvalue λ 1 = −r and ) are both negative for all non-negative parameter values.The eigenvalue Hence, when b < 1 + c a , all three eigenvalues are negative.This implies E 1 is locally asymptotically stable.Similarly, when b > 1 + c a , we have Actually, we can prove that the equilibrium solution E 1 is globally asymptotically stable when b < 1 + c a .Here by the word "globally", we mean the whole positive invariant domain D. We will construct two Lyapunov functions according to different ranges of the parameter b.For convenience, we make a translation of variables, x = 1 − x, y = y, and v = v.After dropping the bars over variables, the system becomes while the domain D is translated to Proof.For any initial condition (x 0 , y 0 , v 0 ) in the Domain D 1 , we know from Lemma 3.2 and Lemma 3.3, the solution satisfies 0 ≤ x(t) ≤ 1, 0 ≤ y(t) ≤ 1 and v(t) ≥ 0. We will construct two Lyapunov functions according to the values of the parameter b to prove y(t) and v(t) approach 0, and then prove x(t) also approaches 0.
It is easy to see V 2 (x, y, v) > 0. The derivative along a solution is given by V2 Since 1 ≤ b < 1 + c a , that is, ab − a − c < 0 and 1 − b < 0, hence, V2 (x, y, v) < 0. Therefore, combining these two Lyapunov functions, we have y(t) → 0 and v(t) → 0 as t → +∞ when b < 1 + c a .Considering the component x(t), we see that where k is the initial value for x(t).The limit ke −rt → 0 as t → +∞, while the limit = 0.
Thus, x(t) → 0 as t → +∞.Such, (x(t), y(t), v(t)) approaches the origin wherever the initial value is in D 1 .Therefore, the original system (6) has a global attractor E 1 .Theorem 3.5 has an important biological implication.When the burst size of the oncolytic virus is smaller than 1 + c a , the virotherapy always fails.This failure is determined by the replication ability of the oncolytic virus, and it does matter what the initial tumor size, the initial infected portion of the tumor, and the initial amount of injected virus are as long as they are in the domain D. The burst size does not depends on initial conditions.The number 1 + c a is a threshold value for the burst size.We denote it by b s1 as the first threshold.Namely, b s1 = 1 + c a .Consider the critical case b = b s1 .It implies ab = a + c.In this case, the linearized system at E 1 has two negative eigenvalues and one zero eigenvalue.In order to determine the stability of the equilibrium solution E 1 , we will use the center manifold theorem to reduce the system (7) into a center manifold, and study the reduced system.We state the result first.
Theorem 3.6.E 1 is unstable when b = b s1 = 1 + c a .Proof.To reduce the system (7) to a center manifold, we first standardize the system, separate it into two parts, the part with zero eigenvalue and the part with negative eigenvalues.The system (7) does not have complex eigenvalues in this case.We start with the matrix corresponding to the linear part of the system (7), The eigenvalue −r has an associated eigenvector T , 0 has an associated eigenvector V 3 = (ar + a, ar, r) T .We set a transformation matrix to be T = (V 1 , V 2 , V 3 ).Denote X = (x, y, v) T .Then the system (7) can be written as , where A ij , B ij and C ij are coefficients that can be easily determined.
Then the transformed system can be written as where , A = (0), and Z = (y 1 , y 2 ).
It is easy to check that A has zero eigenvalue, B has negative eigenvalues, f 1 , f 2 and f 3 are C 2 differentiable functions, f k (0, 0, 0) = 0 and Df k (0, 0, 0) = 0, k = 1, 2, 3, where Df is the first derivative of the function f .Then from the Center Manifold Theorem [3], there is a center manifold given by Z = h(y 3 ) = , and it satisfies It is known that h(0) = 0 and h (0) = 0. So, we assume that h we use variable u instead of y 3 for simplicity.Then the equation for the center manifold is For simplicity, we only consider the order up to 5, and we will know it is enough late.Then, Similarly, we have expressions for f 2 (h(u), u) and f 3 (h(u), u).We substitute f k into the equation 9, and compare the coefficients on the both sides of the equation.We get Solving these equations for coefficients of h(u), we obtain ), and ).Now we reduce the system (7) to its center manifold, which is a single equation and is where It is easy to see that y 3 = 0 is a node regardless of the sign of A33C13 r + B33C23 1+ab .Therefore, for the system (7), the origin is unstable.
a , there is a third equilibrium solution E * = (x, y, v) while E 1 is unstable, where x = c a(b−1) , y = rc(ab−a−c) a(b−1)(ab−a+rc) , and v = r(ab−a−c) a(ab−a+rc) .It is the unique positive equilibrium solution.The variational matrix at this point is given by The characteristic polynomial is given by where , a 2 = rc(bc+b−1) a(b−1) 2 + rc(ab−a−c)(r−a) a(b−1)(ab−a+rc) , and a 3 = rc(ab−a−c) a(b−1) .By the Routh-Hurwitz Criterion [7], all roots of the polynomial (11) have negative real parts if and only if When b > b s1 , we have ab − a − c > 0. So, H 1 = a 1 > 0, and a 3 > 0. Since H 3 = H 2 a 3 , we only need to consider H 2 .
Define a function ϕ(b) = a(b−1)(ab−a+rc) ab−a+rc+abc − (bc+b−1)(ab−a+rc) (b−1)(ab−a−c) .Since b > 1 + c a , we can conclude from the derivation above that H 2 > 0 if and only if ϕ(b) < r − a, H 2 < 0 if and only if ϕ(b) > r − a, and H 2 = 0 if and only if ϕ(b) = r − a.We will also consider H 2 as a function of the parameter burst size b, and denote H(b) = H 2 .We thus get the following theorem.Theorem 3.7.When ϕ(b) < r − a, the equilibrium solution E * is locally asymptotically stable.
We will refine this theorem in the next subsection, and explain its biological implication.We define a function Φ(x) by It is easy to see that a such that Φ(x 2 ) = 0. Similarly, since Φ(x) approaches −∞ as x approaches −∞, there exists at least one element x 1 < c a such that Φ(x 1 ) = 0. From the relation between functions H(x) and Φ(x), H(b) = 0 has at least two roots, 1 + x 2 > 1 + c a and 1 + x 1 < 1 + c a .For the function Φ(x), there are three different cases: Φ(x) has four distinct real roots, or three distinct real roots (one repeated), or two real roots and two conjugate complex roots.We will check each case.
If Φ(x) has four distinct real roots, then there will be three roots or one root that are greater than c a .Figure 1 shows the case where three roots are greater than c a .We take the greatest root to be x 0 , then x 0 > c a and Φ (x 0 ) < 0 since the polynomial Φ(x) is decreasing around x 0 .Actually Φ(x) is decreasing to −∞ as x goes to +∞ from any point that is slightly bigger than x 0 .

Figure 1. Φ(x) has four roots
If Φ(x) has three distinct real roots, then there is one repeated root.Figure 2 shows cases that Φ(x) has three or two roots.We can not choose this repeated root since the graph of Φ(x) is tangent to the x−axis at the point of the repeated root.Hence Φ(x) does not change sign around the point of the repeated root, and is not monotonic around this point.If the repeated root is in the between of these roots, x 1 and x 2 , we choose the greater one as x 0 .If the repeated root is the smallest one, we choose the greatest root to be x 0 .If the repeated root is the greatest root, we choose the middle one to be x 0 .We claim that the middle root is bigger than c a .If it is smaller than c a , then c a must be in the between of the middle root and the repeated root.Since Φ( c a ) > 0, the graph of Φ(x) can not be tangent to x−axis at the next zero which is the point of the repeated root.This is a contradiction.Moreover, Φ(x) is decreasing around the point of the middle root in this case.If Φ(x) has only two distinct real roots, we take the greater one to be x 0 .It is also true that Φ (x 0 ) < 0.
As we choose x 0 for all cases, x 0 > c a , Φ(x 0 ) = 0 and Φ (x 0 ) < 0. Since we then have We rename the variable x as b, and set b 0 = 1 + x 0 .Then, H (b 0 ) < 0. Since H (b) is continuous, so there exists a positive number δ 1 that can be made smaller Since the tumor load at the equilibrium point E * is less than 1, this equilibrium point represents a partial success of virotherapy at a modest value of the burst size.Although we are not able to prove it is globally asymptotically stable (this is one open problem), it is locally asymptotically stable.This also implies that a permanent reduction of tumor load can be reached.Lemma 3.10.A cubic polynomial λ 3 + a 1 λ 2 + a 2 λ + a 3 = 0 with real coefficients has a pair of pure imaginary roots if and only if a 2 > 0 and a 3 = a 1 a 2 .When it has pure imaginary roots, the pure imaginary roots are given by ±i √ a 2 , the real root is given by −a 1 , and a 1 a 3 > 0.
Differentiating the given polynomial with respect to the parameter τ , we have that 3λ 2 λ + a 1 λ 2 + 2a 1 λλ + a 2 λ + a 2 λ + a 3 = 0. Evaluate the left-hand side of this expression at τ = τ 0 and notice that α(τ 0 ) = 0.If α (τ 0 ) = 0, we have Separate the real part and the imaginary part, we have Solve these two equations for β (τ 0 ), we get respectively.Therefore, . That is, From the proofs of Lemma 3.8, Lemma 3.10 and Lemma 3.11 and Routh-Hurwitz Criterion, we have the following corollary for the characteristic polynomial (11).We consider each coefficient of p(λ) to be a function of the parameter b.Actually, where Proof.
Since the exact value of b 0 is not algebraically expressible, the stability, amplitudes, and periods of periodical solutions that occur around b 0 and E * are difficult to analyze.We will simulate the system in the next section to demonstrate some typical dynamics.Here, we simply state a corollary about Hopf bifurcation [9].Corollary 2. If the equilibrium solution E * is stable, but not asymptotically stable at b = b 0 , then all solutions of the system (6) in a neighborhood of E * are periodical (in a surface).If the equilibrium solution E * is asymptotically stable or unstable at b = b 0 , then there is an asymptotically stable periodical solution in a neighborhood of E * when b is close to b 0 .
It is clear that for all cases in Figure 2, we only can have one Hopf bifurcation at b 0 = 1 + x 0 .However, for the case in Figure 1, we can have three Hopf bifurcations at b = 1 + x 0 , b = 1 + x 2 and b = 1 + x 3 .There will be three families of periodic solutions rising from these bifurcation.The similar bifurcation analysis can be done for the parameter values b = 1 + x 2 and b = 1 + x 3 .The only downside is that we can not study their nature sine we don't have explicit algebraic expressions for these values.
From the proof of Theorem 3.12, when the burst size b is big enough, there is at least one root of the characteristic polynomial (14) with positive real part.Actually, there exists b = b g > b 0 such that for every b > b g , the characteristic polynomial ( 14) has a root with positive real part.Therefore, E * is unstable for any large burst size.However, the tumor load at this equilibrium T (E * ) = x + y can be made as small as possible by increasing the burst size b.We state this in the following theorem.
For these b, H 1 > 0 and H 3 = H 2 a 3 > 0, hence, all real parts of roots of ( 14) are negative.Hence, E * is locally asymptotically stable.When b < b s1 and b is in a neighborhood of 1 + c a , the coefficient a 3 in ( 14) is negative.So, the polynomial ( 14) has at least one root with positive real part.Therefore, E * is unstable.Combining Theorem 3.4 and Theorem 3.6, there is a type of transcritical bifurcation at b = b s1 = 1 + c a .We then obtain the main Theorem 3.1.4. Numerical simulation and discussion.We will present some numerical results on the dynamics of oncolytic virus by using published data, and some open problems in this section.4.1.Numerical simulation.In order to demonstrate the analytical results about dynamics of virotherapy, we use some data from our previous research [6,10] to simulate our basic model.The tumor was considered as a solid sphere.When the tumor radius is 2 mm in mice brain, it is considered to be visible and a surgery can be performed.Once its radius reaches 6 mm, the mouse is regarded to be dead from the tumor.Since our model is space-free, we convert tumor size to cell numbers by using the constant of cell density 10 6 per cubical millimeter as all ODE models do.After non-dimensionlization, the parameter values are r = 0.36, a = 0.11, and c = 0.44.We compute b s1 = 5.Solving H(b) = 0, we have two conjugate complex zeroes b = 0.6321 ± 0.4062i and two real zeroes b = −0.347and b = 27.766.This happens to be the most simple case, where I 0 has only one element I 0 = {27.766},I p = (5, 27.766),I n = (27.766,+∞), and so b s2 = b 0 = 27.766.There is one family of periodic solutions when b is greater than b s0 .When b < 5, virotherapy always fails. Figure 3 shows the dynamics of virotherapy when b = 4.
When 5 < b < 27.766, E * is locally asymptotically stable while E 1 is unstable.Somehow this partial success of virotherapy can be always achieved.Figure 4 shows the treatment will eventually reach the equilibrium point E * after a damped oscillation when b = 27.For each value of the burst size b between 5 and 27, we get the similar behavior of the dynamics of viral treatments.The difference is that for different burst sizes the system reaches different equilibrium values.Since the equilibrium tumor load is a decreasing function of the burst size b, the treatment reaches a possible minimum of the equilibrium tumor load when b = 27.When b > 27.766, there exists a family of periodic solutions that arise from the Hopf bifurcation at the bifurcation value b = 27.766.For b = 28, there is a periodic solution.Figure 5 shows an eventually periodic solution when b = 28.Figure 5 also shows that the population of the tumor cells has a large amplitude and the population of the infected tumor cells has a small amplitude.The population of viruses also oscillates (not shown in this figure).This periodic solution is Lyapunov stable.It represents a predator-prey dynamics among tumor cells and viruses.The tumor cells and viruses coexist in a dynamical way.For even big burst size b, the solution still shares a similar behavior.
If the burst size b is taken to be very large, say b = 2000, the solution will behave differently.Figure 6 shows the profile of the uninfected tumor cells when b = 2000.It is a pulsating oscillation.The minimum of the tumor cell population can reach a very small value, say 10 −7 .That means there are about 100 tumor cells since the maximum tumor size is 6mm [6].If a tumor with 100 cells is considered as an    undetectable tumor, we then consider the tumor is eradicated.The minimum of the pulsating oscillation solution is decreasing as the burst size b increases.We can make it even smaller by increasing b. 4.2.Discussion.In this article we presented the analysis of a common basic model for virotherapy.It highlights key values of the burst size of a virus in oncolytic virus treatment.There are two threshold values for the burst size.When the burst size is smaller than the first threshold value, virotherapy always fails.When the burst size is in the between of the two threshold values, we have a partial success of virotherapy represented by the stable positive equilibrium solution.Since the tumor load is a decreasing function of the burst size, the minimum tumor load can be reached by genetically increasing the burst size of the virus up to the second threshold value.If the set in which the positive equilibrium solution is stable has more than one open intervals, we can increase the burst size up to the supreme value of this set, and still have stable partial therapeutic success with even lower tumor load.Once the burst size is greater than the second threshold value, there are one or three families of stable periodic solutions to the system of virotherapy dynamics.This indicates that there is a predator-prey type interaction among three populations, 0, 0), and the positive equilibrium E * .E 0 is always unstable for all positive values of b.E 1 is globally asymptotically stable when 0 < b < b s1 , and it is unstable when b ≥ b s1 .At b = b s1 , the positive equilibrium E * moves into the domain D, a type of transcritical bifurcation occurs with E 1 and E * .As b increases, while b s1 < b < b s2 and b ∈ I p , E * is locally asymptotically stable; when b > b 0 and b ∈ I n , E * is unstable.Hopf bifurcations occur for some b ≥ b s2 , and these bifurcations give rise to one or three families of periodic solutions.As b becomes large enough,

3. 3 .
Bifurcation analysis.Two types of bifurcations occur in the system (6) as the parameter b varies.A transcritical bifurcation at b = b s1 introduces the equilibrium point E * into the positive invariant domain D. The Hopf bifurcation at some value b > b s1 gives rise to the periodic solutions.In order to analyze these bifurcations, we study the function H(b).

When b = 1 Lemma 3 . 8 .
and b = 1 − rc a , H(b) and Φ(b − 1) have the same zeroes.We always assume all other parameters, a, c, and r are positive.The equation H(b) = 0 has at least two real roots, one is less than b s1 = 1 + c a , and the another one is greater than b s1 .Furthermore, among all roots that are greater than b s1 , there exists a root denoted by b 0 and a small neighborhood of b 0 , (b 0 − δ 1 , b 0 + δ 1 ), where δ 1 < b 0 − b s1 , such that H (b 0 ) = 0 and H(b) is monotonic in this interval.Proof.Consider the function Φ(x).It is easy to compute that Φ( c a ) = c 3 (1 + r + c + a)(1 + c + a) 1+r a > 0. Since Φ(x) approaches −∞ as x approaches +∞, from the intermediate value theorem of continuous functions, there exists at least one element x 2 > c

Figure 2 .
Figure 2. Φ(x) has three or two roots

Theorem 3 . 9 .
From the proof of Lemma 3.8, I p is not empty, and it is a finite open interval or a union of two finite open intervals.I n is also not empty, and is an infinite open interval or a union of a finite open interval and an infinite open interval.I 0 is a closed set.I 0 at least has one point, and at most has 3 points.We denote the smallest number in I 0 by b s2 .For the case of Figure 1, b s2 = 1 + x 2 , b 0 = 1 + x 0 and b s2 < b 0 ; for all cases of Figure 2, b s2 = b 0 = 1 + x 0 .It is also easy to see that H(b) > 0 when b s1 < b < b s2 .b s2 is the second threshold value for the burst size.We now can refine Theorem 3.7 as the following theorem 3.9.When b s1 < b < b s2 , the equilibrium solution E * is locally asymptotically stable.At this equilibrium point, the tumor load T (E * ) = x+y = c(1+r) ab−a+rc < 1. Proof.From the proof of Theorem 3.8, H(b s1 ) > 0. Since H(b) is a continuous function, H(b) is positive until it become zero.So, H(b) > 0 for b s1 < b < b s2 .Then, all roots of the characteristic polynomial (11) have negative real parts by the Routh-Hurwitz Criterion.The verification of the second part is straightforward.T (E * ) = x + y = c a(b−1) + rc(ab−a−c) a(b−1)(ab−a+rc) = c(1+r) ab−a+rc .With the condition b > 1 + c a , we have c < ab − a, and then c + rc < ab − a + rc.Hence, c(1+r) ab−a+rc < 1.

Theorem 3 . 13 .
At the equilibrium solution E * , the tumor load T (E * ) is a decreasing function of the burst size b, and T (E * ) b→+∞ = 0. Proof.From Theorem 3.9, we have T (E * ) = c(1+r) ab−a+rc .Because of the conditions b > 1 + c a and c(1+r) ab−a+rc < 1, it is obvious that T (E * ) is a decreasing function of b, and therefore, we have T (E * ) b→+∞ = 0. We now study the relation between equilibria E 1 and E * .When b < b s1 , the equilibrium E * is not in the positive invariant domain D. As b increases to b s1 = 1 + c a , the equilibrium E * moves into D, and it coalesces with the equilibrium point E 1 .The stability of these equilibrium solutions exchange at b s1 .For the characteristic polynomial (14), since H(b s1 ) > 0, so for b > b s1 and b

Figure 3 .
Figure 3. Dynamics of virotherapy when b = 4. Relative time unit used after non-dimensionalization.Each variable is the ratio between its population size and the maximum size of the tumor.

Figure 6 .
Figure 6.The tumor cell population when b = 2000.The initial condition is x = 0.5, y = 0.5 and v = 1.5.
[5] tumor growth is modeled by logistic growth, and C is maximal tumor size.The term αy models the release of virions by infected tumor cells, and γv is the clearance rate of free virus particles by various causes including non-specific binding and generation of defective interfering particles.The death rate of tumor cells dx seems redundant, since it is included in the logistic model.Dingli et al[5]proposed a similar model by replacing the first equation in (2) by a generalized logistic model,