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.


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] 2 Computational and Mathematical Methods in Medicine 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 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 , 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 , 1 , 2 , , and 0 , 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 1 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].

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.
As (0) ≤ (0) + (0) ≤ 1, by the comparison theorem we get ( ) ≤ 1. On the other hand, since again by the comparison theorem we also have 0 ≤ ( ) + ( ) ≤ 1. It follows that 0 ≤ ( ) ≤ 1. 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.

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. In what follows we will analyze the existence of equilibrium solutions. First, let = ( , , V, ) and Then system (4) can be simply written as the autonomous system / = ( ). We assume that ( , , V, ) ∈ . The equilibrium points are solutions of the equation ( ) = 0; that is, If = 0, then, from the second and the third equation of (9), ( + 1) = 0 and = V( + ). Since + 1 > 0, then = 0. It leads to V( + ) = 0, which implies V = 0. The last equation of (9) gives − = 0, which implies = 0. Thus Notice that in order for the first three coordinates of 2 to be positive, it is enough that − − > 0; that is, > / + 1.
We summarize the details of the analysis above as follows.
(ii) When ̸ = 0, we have the following cases.
(a) If = 0, then (1) when V = 0, we have equilibrium solution , and V satisfies the following cubic equation: Computational and Mathematical Methods in Medicine 5 In this case, we can conclude the following.
= / = / 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 0 . If we do not consider the effect of the immune system, and the viruses are not powerful, that is, the burst size is smaller than a critical value, then the system has the equilibrium 1 with only tumor cells; if the viruses are powerful, that is, the burst size is greater than a critical value, then the system has the equilibrium 2 with the coexistence of tumor cells, infected tumor cells, and viruses. When we consider the immune effect, if the burst size is another critical value 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.

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 0 is locally unstable, 1 is locally asymptotically stable if < 1 + / and unstable if > 1 + / , and 2 is locally asymptotically stable if is in some range, while 3 , 4 , and 5 are locally unstable when , , and satisfy some conditions. We use Lyapunov functions to show 1 is globally asymptotically stable if < 1 + / and > 1. We apply the center manifold theorem to show 1 is locally asymptotically stable if = 1+ / . We summarize the main results in Lemma 3. For the combined parameter values, , = 1, 2, , , = 1, 2, 3, they will be defined in the following context. For methods we applied in this subsection, we refer to Allen [20] for basic knowledge of stability analysis and Carr [21] for center manifolds.
We look at the stability of trivial equilibrium solutions first. The variational matrix of system (4) is given by At the critical point 0 , the variational matrix is The corresponding eigenvalues are , −1, − , and − . We know that the local stability of 0 of system (4) is the same as that of the linearized system at 0 . Since > 0, 0 is locally unstable. For system (4), the local stable invariant manifold is tangent to the -V-subspace, and the unstable invariant manifold is tangent to the -axis. Biologically, the tumor will grow from an initial small value around 0 without viruses and infected tumor cells.
Proof. At the equilibrium point 1 , the variational matrix is The characteristic polynomial of this matrix is ( + )( + ) (  In fact, we can show that 1 is globally asymptotically stable in the positive invariant domain when < 1+ / and < by constructing two Lyapunov functions according to different ranges of the parameter . For convenience, we translate 1 into the origin by setting = 1 − , = , V = V, and = . After dropping all the bars over variables, system (4) becomes while the domain is translated to 1 and ( ) ≥ 0. Now we construct two Lyapunov functions corresponding to the values of the parameter to prove ( ) and V( ) approach 0, then we show ( ) and ( ) also tend to 0.
Biologically, Proposition 4 is easy to understand. When the viral burst size is smaller than a critical value which is = 1 + / , 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 1 when = 1 + / . This is a critical case, since the linearized system at 1 has three negative eigenvalues and one zero eigenvalue. we have to reduce the system to its local center manifold. We actually have the following proposition. Proof. Consider = 1+ / , which implies that = + . The linearized system at 1 has three negative eigenvalues and one zero eigenvalue. In order to determine the stability of 1 , we use the center manifold theorem to reduce system (22) into a center manifold, and then we study the reduced system. So we separate it into two parts, one with zero eigenvalue and the ) . Set = ( 1 , 2 , 3 , 4 ) and = ; then and = ( 1 , 2 , 3 , 4 ) is determined by Denote −1 = ( 1 , 2 , 3 , 4 ) ; then we can express , = 1, 2, 3, 4, in terms of : where , , , and are coefficients that can be easily determined. The transformed system can be expressed as It is easy to check that each , = 1, 2, 3, 4, is twice differentiable function, (0, 0, 0, 0) = 0 and (0, 0, 0, 0) = 0, where is the Jacobian matrix of the function . By the center manifold theorem, there exists a center manifold: with ℎ(0) = 0 and ℎ(0) = 0, and it satisfies Since ℎ(0) = 0 and ℎ(0) = 0, we can assume that ) ; here we use the variable instead of 3 for simplicity. Then By substituting 's into (34) and equating both sides of the equation, we can get 33 = − 2 ( +1)( −1)/( +1) < 0, since = 1 + / > 1. The asymptotically behavior of zero solution of system (31) is governed by that of the equation 3 Since 33 < 0, 3 = 0 is locally asymptotically stable. Therefore, the trivial solution of system (22) is locally asymptotically stable.
We next consider the stability of the equilibrium solution For lately defined and , we have a proposition as follows.
We show this proposition as follows. The variational matrix at 2 is given by ) .
(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 / .
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 3 , 4 , and 5 , when they exist, we derive conditions under which they are unstable.  ) .
Lastly, when (V) = 0 has two distinct positive real roots 0 < V 2 < V * 2 = < V 3 < * , the variational matrices at 4 and 5 are the same as the variational matrix at 3 except that is replaced by V 2 and V 3 , respectively. We obtain the corresponding characteristic polynomials of those matrices which are the same as the characteristic polynomial of except for replacing by V 2 and V 3 . By the same argument as above, when < (( + V 2,3 ) 2 / 2 )((2 2 / )V 2,3 − + ) − / + 1 š 2,3 , 4 and 5 are locally unstable.
It is interesting to notice that our model system has 3 positive equilibria when the viral burst size is not too large and the relative immune killing rate falls into some intervals. Proposition 8 gives conditions that ensure these equilibrium solutions are unstable. Biologically, when the relative immune killing rate and relative immune response decay rate fall into some ranges, we may genetically change the viral burst size to avoid coexistent equilibria.

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 1 and 2 . When < 1 , 2 is outside of the positive domain and 1 is locally asymptotically stable. As increases to 1 = 1 + / , 2 moves into and it coalesces with 1 which is still locally asymptotically stable. But when 1 < < 2 and ∈ , the stability of these equilibrium points interchanges, which means that 2 is locally asymptotically stable while 1 is unstable. Notice that when > 0 and ∈ , 2 is locally unstable.
Consider each coefficient of ( ) as a function of the parameter . Then The following theorem is our main result for occurring a Hopf bifurcation around 0 , which appears in [16] as Theorem 3.12. For completion, we restate this theorem and related lemmas and corollary and give slight different proofs. 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.
We now can prove our main Theorem 9.
Combining Proposition 6 and applying this theorem, we can obtain a statement about Hopf bifurcations of our model. Because we cannot find explicit algebraic expression for 0 , it is very hard to study the nature of periodical solutions that occur around 2 when is close to 0 such as the amplitudes, periods, and their stability. However, we can make some statements about the general properties of these periodical solutions as follows.
(i) If 2 is stable but not asymptotically stable, then any solution of system (3) in is periodical in a surface.
(ii) If 2 is asymptotically stable, then there is an asymptotically stable periodical solution ( ) in when is close to 0 . Any solution inside will spiral into (iii) If 2 is unstable, then there is an asymptotically stable periodical solutioñ( ) in when is close to 0 . Any solution starting at nearby 2 will spiral out and emerge intõwhen time is increasing, and any solution in nearby outsidẽ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 4 and 5 because their formulas are extremely cumbersome and therefore we will treat this by numerical simulations in the next section.
We close Section 2 with the following theorem that summarizes the main results about the our model and its biological implications.

Theorem 14.
The dynamical behaviors of system (4) can be described as follows.
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 1 and the relative immune decay rate is greater than 1, the virotherapy always fails. If the burst size is greater than 1 and smaller than another critical value 2 , 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/ ), 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, 3 , 4 , and 5 , 3 is difficult to obtain in practice because it requires a particular threshold of the burst size. In rat experiments, the virus burst size can be genetically changed as we want, but usually, we can ensure a range of the burst size in the process of genetic modification. 4 and 5 are unstable if the burst size is smaller than a threshold. Biologically, these two equilibria are not important because of their instability. The immune free equilibrium 2 can be stable. If the burst size is made big enough, the tumor cell portion will be small in 2 . On the other hand, the model can have periodic solutions. This may provide an opportunity for combining surgery with the phase where the tumor cell portion is in the lowest state.

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.  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 1 is locally asymptotically stable. Figure 2 shows 1 is locally unstable since = 9 > 1 + / = 2.82. Figure 3 shows 2 is locally asymptotically stable because is between 1 and 2 and = 0.058 < = 0.06. Figure 4 shows 4 is locally asymptotically stable since all eigenvalues of the variational matrix at 4 are negative. Figure 5 shows 5 is locally unstable sine < 3 . Figures 6-8 show periodic solutions rising from Hopf bifurcations. Figure 9 shows the tumor cell population when the burst size is 1000. 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, 1 , 2 , , and 0 , 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 1 or and and 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.

Conflicts of Interest
The authors declare that they have no conflicts of interest.