Abstract
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 stands for uninfected tumor cells, for infected tumor cells, and 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 , we have the following system:where is tumor growth rate, 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, is the burst size of oncolytic viruses (i.e., the number of new viruses coming out from a lysis of an infected cell), is immune killing rate of viruses, is clearance rate of viruses, 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 and parameters related to the innate immune response. We can denote the dynamical behaviors of the model by , the ratio of the innate immune killing rate of infected tumor cells over the innate immune killing rate of free viruses by , the relative innate immune killing rate of viral therapy by , the ratio of the innate immune clearance rate over the stimulation rate of the innate immune system by , and the relative innate immune response decay rate by . 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 , , , , and , where is a function of the relative innate immune killing rate and relative innate immune response decay rate . When 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 is greater than or and the two relative rates are in the complement intervals, the system has at most 3 equilibrium solutions with 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 , , , , and rename parameters , , , , , , and . Then system (2) becomesFor convenience, dropping all the bars and writing as , we obtainWe assume that all parameters are nonnegative.
Lemma 1. If , , , and , then the solution of system (4) , , , and for . Furthermore, if , then for , and .
Proof. The proof is similar to that of Lemma in [16]. Here we only show the second part of this lemma by using comparison theorem for ODEs; that is, if , then and . In fact, since , , and are nonnegative for all , As , by the comparison theorem we get . On the other hand, since again by the comparison theorem we also have . It follows that . Then , so by the comparison theorem . Taking both sides yield .
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 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 [16]. The effects of the innate immune system on the virotherapy in our model are encoded in the parameters , , and . To understand how the innate immune system affects the dynamics of oncolytic viral therapy, we use three combined parameters, the viral burst size , the relative immune killing rate , and the relative immune response decay rate , to describe the solution behaviors of our model. We summarize the possible equilibrium solutions in the following lemma.
Lemma 2. When , , and with , system (4) has 3 equilibrium solutions: , , and . When either or and with , system (4) has a unique positive equilibrium solution: . When either or and with , system (4) has two distinct positive equilibrium solutions: and .
In what follows we will analyze the existence of equilibrium solutions. First, let and Then system (4) can be simply written as the autonomous system . We assume that . The equilibrium points are solutions of the equation ; that is,If , then, from the second and the third equation of (9), and . Since , then . It leads to , which implies . The last equation of (9) gives , which implies . Thus is an equilibrium point.
If , the first equation of (9) implies . Consider the last one . If , from the second and the third equation of (9), we get and . These follow and hence . If , then and , which implies that . So is an equilibrium point.
Now if , then . Since we want to find positive equilibrium points, we assume . Then . From the equation , we have . It follows that Since , we have . Thus we get the third critical point Notice that in order for the first three coordinates of to be positive, it is enough that ; that is, .
Next, if , then . Set , then . From the equation , we can derive . By the third equation of (9), . Plugging these expressions into the second equation gives , where Since and , has at least one negative root, say . Assume that is the least root. If and , then has no sign changes at all and hence the other two roots of are either both negative or both complex. Notice that and are equivalent to and . In this case, the system only has equilibrium points: , , and .
We now look at other situations of . Taking derivative, . By long division, Suppose that , then has 2 distinct roots, and , where . If and , then has one negative root, , and one positive root, . In this case, in addition to the 3 equilibrium points , , and , system (4) has one positive equilibrium point: . To guarantee all four coordinates of are positive, we need to impose and , which imply that .
On the other hand, if and , then has one negative root, , and 2 distinct positive roots, . Hence system (4) has two positive equilibrium points: . Notice that is equivalent to either , or . Furthermore, , are equivalent to and , where ; is equivalent to . The condition is equivalent to . 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.(i)When , we have equilibrium solution .(ii)When , we have the following cases.(a)If , thenwhen , we have equilibrium solution .2when and , we have (b)If , then , , , and satisfies the following cubic equation: where , , In this case, we can conclude the following.If , , and , then system (4) has three equilibrium points: , , and .If either , , or and , then system (4) has a unique positive equilibrium point: where and If either , , or and , then system (4) has two distinct positive equilibrium points: where and are two distinct positive real roots of the above cubic equation that satisfy , in which .
In order to interpret our mathematical conditions biologically, we need to understand the combined parameters and first. can be considered as a relative immune response decay rate since is innate immune cell death rate, is innate immune stimulating rate by infection, and is tumor carrying capacity. 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 . 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 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 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 is locally unstable, is locally asymptotically stable if and unstable if , and is locally asymptotically stable if is in some range, while , , and are locally unstable when , , and satisfy some conditions. We use Lyapunov functions to show is globally asymptotically stable if and . We apply the center manifold theorem to show is locally asymptotically stable if . We summarize the main results in Lemma 3. For the combined parameter values, , , , 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. is unstable. is globally asymptotically stable when and and unstable when . is locally asymptotically stable when . is unstable when . and 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 , the variational matrix is The corresponding eigenvalues are , , , and . We know that the local stability of of system (4) is the same as that of the linearized system at . Since , is locally unstable. For system (4), the local stable invariant manifold is tangent to the -- subspace, and the unstable invariant manifold is tangent to the -axis. Biologically, the tumor will grow from an initial small value around without viruses and infected tumor cells.
Proposition 4. The equilibrium solution is locally asymptotically stable when , and it is locally unstable when . is globally asymptotically stable when and .
Proof. At the equilibrium point , the variational matrix is The characteristic polynomial of this matrix is . The eigenvalues are , , and . Since the eigenvalues , , and are negative for all positive parameters, is locally asymptotically stable if and only if . This is equivalent to , which is the same as . Similarly, if , then , and hence is unstable.
In fact, we can show that is globally asymptotically stable in the positive invariant domain when and by constructing two Lyapunov functions according to different ranges of the parameter . For convenience, we translate into the origin by setting , , , and . After dropping all the bars over variables, system (4) becomeswhile the domain is translated to . For any initial condition , the solution of (22) satisfies , , , and . Now we construct two Lyapunov functions corresponding to the values of the parameter to prove and approach , then we show and also tend to .
When , we define the Lyapunov function . It is clear that , and the orbital derivative . When , consider the Lyapunov function . Obviously, . The orbital derivative along a solution is given by ; that is, and ; therefore . Combining these two Lyapunov functions gives and as when . Considering the component , we have By the comparison theorem, Taking limit of both sides as and using the L’Hospital’s Rule and the fact that and approach yield . Finally, since , we have . By the comparison theorem, as , since . Therefore system (3) has a global attractor .
Biologically, Proposition 4 is easy to understand. When the viral burst size is smaller than a critical value which is , 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 when . This is a critical case, since the linearized system at 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 is locally asymptotically stable when .
Proof. Consider , which implies that . The linearized system at has three negative eigenvalues and one zero eigenvalue. In order to determine the stability of , 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 , , , and . The associated eigenvectors of these eigenvalues, respectively, are , and . System (22) can be written as , where . Set and ; thenwhere and is determined by Denote ; then we can express , in terms of : where , , , and are coefficients that can be easily determined. The transformed system can be expressed aswhere It is easy to check that each , is twice differentiable function, and , where is the Jacobian matrix of the function . By the center manifold theorem, there exists a center manifold:with and , and it satisfiesSince and , we can assume that here we use the variable instead of for simplicity. Then we can compute By substituting ’s into (34) and equating both sides of the equation, we can get , since . The asymptotically behavior of zero solution of system (31) is governed by that of the equation or . Since , 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 and , we have a proposition as follows.
Proposition 6. When , is locally asymptotically stable.
We show this proposition as follows. The variational matrix at is given by The characteristic polynomial of iswhere , and By Routh-Hurwitz Criterion, all roots of have negative real parts if and only ifSince , and . And those conditions in (40) are equivalent to . This inequality is the same as Therefore, we can conclude that if and ; then is locally asymptotically stable. Now we can refine this result by considering , and It is easy to check that . Since , and have the same roots. Since and , there are at least one and one so that . Then has at least one root and one root . We consider three different cases as follows. (i)If has 4 distinct real roots, then either 3 roots are bigger than or only 1 root is bigger than .(ii)If has 3 distinct real roots, then one root must be repeated.(iii)If has 2 distinct real roots, then one root must be bigger than .
In all cases, we always have at least one root bigger than ; denote it by . Then , , and . Let ; since we get . As is continuous, there is that can be made smaller than so that in . Thus is monotonically decreasing in this interval. Thus, we have proved the following property.
Property 7. The function has at least 2 real roots, one of which is bigger than , and the other is less than . Among all roots bigger than , there is a root and a neighborhood of , where so that and is decreasing in this interval.
Define , , and . All these sets are nonempty and has either at least 1 element or at most 3 elements. Let . Note that . It is easy to check that when , .
Next, we refine the condition , where . This inequality is equivalent to Considering it as a quadratic polynomial of , we have . If , then this inequality is always true for all , so is always true. If , then (i)when , the right-hand side of the inequality has 2 negative roots , but since , this inequality is obviously true and hence is always true;(ii)when , the inequality is equivalent to , where and are 2 positive roots of the right-hand side of the inequality (they may be equal).
Let . Then we have the following result: if and is in this intersection, then 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 , , and , when they exist, we derive conditions under which they are unstable.
Proposition 8. is locally unstable when . and are locally unstable when .
Proof. When has a unique positive root , the system has a unique positive equilibrium solution . The variational matrix at is The characteristic polynomial of this matrix is computed as the quartic polynomial , where Assume that . Since , has at least one positive root. The fact is equivalent to . On the other hand, we can compute in terms of coefficients of and note that the coefficients , do not depend on . Since , we have . As , so , which implies that Thus, when , is locally unstable.
Lastly, when has two distinct positive real roots , the variational matrices at and are the same as the variational matrix at except that is replaced by and , respectively. We obtain the corresponding characteristic polynomials of those matrices which are the same as the characteristic polynomial of except for replacing by and . By the same argument as above, when , and 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 and . When , is outside of the positive domain and is locally asymptotically stable. As increases to , moves into and it coalesces with which is still locally asymptotically stable. But when and , the stability of these equilibrium points interchanges, which means that is locally asymptotically stable while is unstable. Notice that when and , is locally unstable.
In order to study the Hopf bifurcation at , we look at the characteristic polynomial (38): where . For convenience, we assume that . From the derivation of Proposition 6, we know . That is, has a negative root . Thus, the assumption reduces the study of the quartic polynomial to the cubic polynomial .
Consider each coefficient of as a function of the parameter . Then where , , and .
The following theorem is our main result for occurring a Hopf bifurcation around , which appears in [16] as Theorem . For completion, we restate this theorem and related lemmas and corollary and give slight different proofs.
Theorem 9. There exists a neighborhood of , , such that for each in this neighborhood has a pair of complex conjugate eigenvalues , where changes sign when passes through and . Moreover, when , and . Notice that can be made small enough so that .
To prove this theorem, we need two lemmas about the properties of roots of cubic equations which appear in [16] as Lemmas and .
Lemma 10. The cubic polynomial with real coefficients has a pair of pure imaginary roots if and only if and . When it has pure imaginary roots, these imaginary roots are given by , the real root is given by , and .
Proof. Suppose the cubic polynomial has a pair of complex roots and one real root . Then . Expanding the left-hand side and then equating both sides give , , and . This implies that , , and . The last equation yields . Thus, if and only if and . If , then and , which follows that .
Lemma 11. Consider polynomial , where , . Let be the roots of the polynomial. Suppose there is such that and . Moreover, if , then .
Proof. Differentiating the polynomial with respect to and evaluating the derivative at , we notice that , then we get, after equating the real part and the imaginary part, . By Lemma 10, since , we obtain the desired result.
From the proofs of Lemmas 10 and 11, and Routh-Hurwitz Criterion, we have the following corollary about .
Corollary 12. There is a neighborhood of , , where , such that for all , and the real part of the root of changes sign when passes through .
Proof. Since , and . As , . Since is continuous with respect to , there is a neighborhood of such that in that neighborhood. Its radius can be made small enough so that and . We know that is decreasing in this neighborhood . If , then . Since and , by Routh-Hurwitz Criterion . If , then . Since , , and , from the proof of Lemma 10 we have and . Thus the sign of changes when passes through .
We now can prove our main Theorem 9.
Proof. By Property 7 in previous section, and . Then, for , and ; hence . By Lemma 10, has a pair of purely imaginary roots, , and the real root is . Since and is continuous, we can find a neighborhood of so that in this neighborhood and its radius can be chosen so that and . By the above corollary, in the interval , changes sign when passes through . Finally, we claim that . By way of contradiction, if , then by Lemma 11 we have . Then this yields , 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 , for system (4) for all . The variational matrix of at , , has 2 strictly negative roots and 2 conjugate complex roots in the neighborhood of , in which , changes sign when passes through , and . Consequently, in a neighborhood of and for any , the system has nontrivial periodical solutions in .
Because we cannot find explicit algebraic expression for , it is very hard to study the nature of periodical solutions that occur around when is close to 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 is stable but not asymptotically stable, then any solution of system (3) in is periodical in a surface.(ii)If is asymptotically stable, then there is an asymptotically stable periodical solution in when is close to . Any solution inside will spiral into when time is increasing and any solution in outside will spiral and emerge into when time is increasing.(iii)If is unstable, then there is an asymptotically stable periodical solution in when is close to . Any solution starting at nearby will spiral out and emerge into when time is increasing, and any solution in 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 and 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 , , and with , system (4) has at most 3 equilibrium solutions: , , and . is unstable. is globally asymptotically stable if and and unstable if . is locally asymptotically stable if .(ii)When either or and with , system (4) has a unique positive equilibrium solution: . is unstable if .(iii)When either or and with , system (4) has two distinct positive equilibrium solutions: and . and are unstable if .(iv)When , there exist a neighborhood of and a neighborhood of , such that for any , system (4) has nontrivial periodical solutions in .
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 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 , 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, , , and , 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. and 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 can be stable. If the burst size is made big enough, the tumor cell portion will be small in . 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 , , , , , , and . Then . Solving gives 2 conjugate complex roots , and 2 real roots and . Therefore, , , , and hence . On the other hand, all coefficients of the cubic equation are , , and . Considering the case with the burst size , we find all equilibrium solutions of system (4): By the analysis of previous section, is always unstable. The equilibrium point is locally asymptotically stable if and unstable if . For the equilibrium point , since and , it is locally asymptotically stable. In order to check the stability of the equilibrium point , we need to compute the Jacobian matrix of at this point, which is . Using the Matlab, we can calculate 4 eigenvalues of , which are , , and . This guarantees the locally asymptotical stability of . For the last equilibrium point , we know that , which is the largest root of . Then, we can compute the quantity . As , is unstable.
When , the bifurcation analysis gives us the appearance of periodical solutions around the equilibrium point . However, in this case, two positive equilibrium points and do not exist. Now we fix the burst size and other parameters, considering as a variable parameter. Observe that when , we have two positive equilibrium points: and . The eigenvalues of the Jacobian matrix at are , , and , whereas when , we also have two positive equilibrium points: and . The eigenvalues at are , , and . This partially shows the existence of the Hopf bifurcation point . Using the Newton method for the computation of Hopf bifurcation points, we find .
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 , or simply, relative time.
Figure 1 shows that is locally asymptotically stable.
Figure 2 shows is locally unstable since . Figure 3 shows is locally asymptotically stable because is between and and . Figure 4 shows is locally asymptotically stable since all eigenvalues of the variational matrix at are negative. Figure 5 shows is locally unstable sine .
Figures 6–8 show periodic solutions rising from Hopf bifurcations.
Figure 9 shows the tumor cell population when the burst size 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 , for example, , , , and , where is a function of the relative innate immune response killing rate and relative innate immune decay rate , which we combine with innate immune parameters and to describe the dynamics of our model (4). When is smaller than and and 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 and relative innate immune decay rate . When is greater than or and and fulfill the complement conditions, system (4) has at most 3 equilibrium solutions with 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.
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.