Mathematical Modeling and Analysis of Tumor Therapy with Oncolytic Virus

In this paper, we have proposed and analyzed a nonlinear mathematical model for the study of interaction between tumor cells and oncolytic viruses. The model is analyzed using stability theory of differential equations. Positive equilibrium points of the system are investigated and their stability analysis is carried out. Moreover, the numerical simulation of the proposed model is also performed by using fourth order Runge-Kutta method which supports the theoretical findings. It is found that both infected and uninfected tumor cells and hence tumor load can be eliminated with time, and complete recovery is possible because of virus therapy, if certain conditions are satisfied. It is further found that the system appears to exhibit periodic limit cycles and chaotic attractors for some ranges of the system parameters.


Introduction
A connection between cancer and viruses has long been theorized and case reports of cancer regression (cervical cancer, Burkitt lymphoma, Hodgkin lymphoma) after immunization or infection with an unrelated virus appeared at the beginning of the 20th century [1].Efforts to treat cancer through immunization or deliberate infection with a virus began in the mid 20th century [1,2].As the technology for creating a custom virus did not exist, all early efforts focused on finding natural oncolytic viruses.During the 1960s, promising research involved in using poliovirus [3], adenovirus [1], Coxsackie virus [4], and others [2].The early complications were occasional cases of uncontrolled infection, resulting in significant morbidity and mortality; the very frequent development of an immune response, while harmless to the patient [1], destroyed the virus and thus prevented it from destroying the cancer [3].In a number of cases, cancer cells exposed to viruses have experienced widespread necrosis, which cannot be entirely accounted for by viral replication alone.Cytotoxic T-cell responses directed against virusinfected cells have been identified as an important factor in tumor necrosis.However, since viruses are normal human pathogens, they induce an immune response, which reduces the effectiveness of viruses.For example, increased antibody could deactivate viruses before the tumor has been destroyed.This can be overcome by using parental viruses that are not normal human pathogens, thereby avoiding any pre-existing immunity.However, this does not avoid subsequent antibody generation.Alternatively, the viral vector can be coated with a polymer such as polyethylene glycol, shielding it from antibodies, but this also prevents viral coat proteins adhering to host cells.Deactivation of the immune system is not desirable, since it has a positive effect on tumor necrosis.Even when a response was seen, these responses were neither complete nor durable [1].
With the later development of advanced genetic engineering techniques, researchers gained the ability to deliberately modify existing viruses, or to create new ones.All modern research on oncolytic viruses involves viruses that have been modified to be less susceptible to immune suppression, to more specifically target particular classes of cancer cells, or to express desired cancersuppressing genes.The first Oncolytic virus to be approved by a regulatory agency was a genetically modified adenovirus named H101 by Shanghai Sunway Biotech.It gained regulatory approval in 2005 from China's State Food and Drug Administration (SFDA) for the treatment of head and neck cancer [5,6].Sunway's H101 and the very similar Onyx-15 have been engineered to remove a viral defense mechanism that interacts with a normal human gene p53, which is very frequently deregulated in cancer cells [6].Onyx-015 is a adenovirus that was developed in 1987 with the function of the E1B gene knocked out, meaning cells infected with Onyx-015 are incapable of blocking p53's function.If Onyx-015 infects a normal cell, with a functioning p53 gene, it will be prevented from multiplying by the action of the p53 transcription factor.However if Onyx-015 infects a p53 deficient cell it should be able to survive and replicate, resulting in selective destruction of cancer cells [7,8].
The interaction between the growing tumor and the replicating oncolytic virus are highly complex and nonlinear.Thus to precisely define the conditions that are required for successful therapy by this approach, mathematical models are needed.Several mathematical models that describe the evolution of tumors under viral injection were recently developed [9][10][11].Other mathematical models for tumor-virus dynamics are, mainly, spatially explicit models, described by systems of partial differential equations (which is an obvious and necessary extension of ordinary differential equations models in as much as most solid tumors have distinct spatial structure), the local dynamics is usually modeled by systems of ordinary differential equations that bear close resemblance to a basic model of virus dynamics [12].Wu et al. modeled and compared the evolution of a tumor under different initial conditions [13].Friedman and Tao (2003) presented a rigorous mathematical analysis of somewhat different model [14].Our model builds upon the model of Artem S. Novozhilov [9] with a modified functional response between the cells.Novozhilov presented a mathematical model that describes the interaction between two types of tumor cells (the cells that are not infected but are susceptible so far as they have the cancer phenotype) with ratio dependent functional response between them.We consider a more realistic type of functional response with saturation effect of virus-cell interaction as even when a virus is oncolytic and it punches a hole in a tumor, the immune response of the individual to the virus occurs so fast that the effects are quickly wiped out and the tumor continues to grow.We discuss the linear stability analysis of the biologically feasible equilibrium states of this model.The ranges of the system parameters for which the system has chaotic behavior is found.Some authors discussed the problem of chaos and stability analysis of some biological models such as cancer and tumor model, genital herpes epidemic, stochastic lattice gas prey-predator model and many other models, see, for example, [15][16][17][18][19].The problem of optimal control of the unstable equilibrium states of cancer self-remission and tumor system using a feedback control approach is studied by Sarkar and Banerjee and El-Gohary [15,19].The objective of present work is to study the interaction between growing tumor and the replicating Oncolytic virus with a functional response with saturation effect.The saturation effect accounts for the fact that the number of contacts an individual cell reaches some maxi-mal value as our immune system evolves to stop virus just as the virus evolves to enter cells and replicate.Our model exhibits that complete elimination of tumor is possible with the help of oncolytic virus therapy in the treatment of tumor.
This paper is organized as follows: In Section 2, we outline our model.In Section 3, Boundedness of solutions of the system is studied.Section 4, gives a review of equilibrium points of the system.Sections 5, deals with stability analysis of equilibrium points.In Section 6, we have determined conditions for global stability of interior equilibrium point.It is further followed by numerical simulation in Section 7. Lastly, a short discussion is represented in Section 8.

Mathematical Model
The model considers two types of tumor cells x and y growing in logistic fashion.x is the size of the uninfected tumor cell population and y is the size of infected tumor cell population.It is assumed that oncolytic viruses enter tumor cells and replicate.These tumor cells, infected with oncolytic viruses, further cause infection in other tumor cells.Oncolytic virus preferentially infects and lysis cancer cells both by direct destruction of the tumor cells, and, if modified, as vectors enabling genes expressing anticancer proteins to be delivered specifically to the tumor site.Based on these assumption model takes the following form: Here 1 r and 2 r are the maximum per capita growth rates of uninfected and infected cells correspondingly; K is the carrying capacity , b is the transmission rate (this parameter also includes the replication rate of the viruses) ; The expression by x y a   , displays a saturation effect accounting for the fact that the number of contacts an individual cell reaches some maximal value as our immune system evolves to stop virus just as the viruses evolve to enter cells and replicate.a is the measure of the immune response of the individual to the viruses which prevents it from destroying the cancer and  is the rate of infected cell killing by the viruses (cy-totoxicity).All the parameters of the model are supposed to be nonnegative.

Boundedness of Solutions
Boundedness may be interpreted as a natural restriction to grow because of limited resources.To establish the biological validity of the model system, we have to show that the solutions of system (2.1) are bounded for this we find the region of attraction in the following lemma.Lemma 1.: All the solutions of (2.1) starting in the R  denote the non-negative cone of 2 R including its lower dimensional faces.
Proof: From system (1) we get: then by usual comparison theorem, we get the following expression as Thus, it suffices to consider solutions in the region  .Solutions of the initial value problem starting in  and defined by (2.1) exist and are unique on a maximal interval [20].Since solutions remain bounded in the positively invariant region  , the maximal interval is well posed both mathematically and epidemiologically.

Existence and Uniqueness of Equilibrium Points
An equilibrium point is a point at which variables of a system remain unchanged over time.System (2.1) possesses the following equilibria: are positive solutions of the system of algebraic equations given below: Now substituting the value of y  in Equation (4.1), we get To show existence of x  , it suffices to show that Equation (4.3) has a unique positive solution. Taking Thus, there exist a x  in the interval 0 Corresponding value of y  is given by  .Thus if all the three conditions (4.5)-(4.7)hold an unique interior equilibrium always exists.However, these conditions depend upon the parameter values so they do not always hold.We will show further that if Equation (4.5) and Equation (4.6) does not hold other equilibrium points of the system become locally asymptotically stable.

Local Stability Analysis of Equilibrium Points
An equilibrium point is locally asymptotically stable if all solutions of the system approaches it as t   .To discuss the local stability of equilibrium points we compute the variational matrix of system (2.1).The signs of the real parts of the eigenvalues of the variational matrix evaluated at a given equilibria determine its stability.The entries of general variational matrix are given by differentiating the right hand side of system (2.1) with respect to , x y .The matrix is given by We denote the variational matrix corresponding to i E by   i V E , 0,1, 2,3 i 

  0 0,0 E
To explore local stability of trivial equilibrium point, we compute variational matrix of 0 E .The variational matrix of equilibrium point 0 E is given by E is an unstable equilibrium point since both the eigenvalues of the matrix are positive.

Local Stability Analysis of
 

1
,0 E K Now, to study the stability behavior of 1 E , we compute the variational matrix   From Equation (5.2) we observe that eigenvalues of plies that equilibrium point 1 E is locally asymptotically stable if death rate of infected cells due to the virus attack is larger than the rate of transmission of infection from virus to the uninfected cells if measure of the immune response of the individual to the viruses is very less as compared to carrying capacity of the cells.Stability of this equilibrium point suggests that infected cells would die without having time to infect other cells and tumor grows unaffectedly.
Remark: 1 E is stable if interior equilibrium point does not exist.

  2 0, E y
The variational matrix of equilibrium point From Equation (5.3) we observe that eigenvalues of the matrix  

Biological interpretation:
Local asymptotic stability of 2 E implies that 2 1 br r   , Since a is a positive parameter and it cannot be less than a negative quantity, . Hence 2 E is locally asymptotically stable under any of the two situations: 1) If net growth rate of uninfected cells is less than the rate of transmission of virus infection than growth rate of infected cells must be greater than the death rate of infected cells caused by the viruses.
2) If net growth rate of uninfected cells is more than the rate of transmission of virus infection then the ratio of net growth rate of uninfected cells to the death rate of infected cells caused by the viruses is more than the ratio of net growth rate of uninfected cells to rate at which they become infective.
Remark: 2 E is stable if interior equilibrium point does not exist.

Local Stability Analysis
where The signs of the real parts of   and   are negative.This implies that 3 E is always locally asymptotically stable if it exists.

Global Stability of Interior Equilibrium Point
An equilibrium point is globally asymptotically stable if system always approaches it regardless of its initial position.We construct Lyapunov functions that enable us to find biologically realistic conditions sufficient to ensure existence and uniqueness of a globally asymptotically stable equilibrium state.Global stability of the interior equilibrium point of system (2.1) is determined in the theorem 6 given below: Theorem 6.: If the following inequality holds then 3 E is globally asymptotically stable with respect to the solutions initiating in the interior of the positive orthant  .
Proof: Consider the following positive definite function about 3 E ln ln x y W x x x y y y x y Computing the derivative of W with respect to t and after some algebraic manipulations, we get We note that dW dt can be made negative definite inside After maximizing left hand side and minimizing right hand side of above inequality we note that dW dt will be negative definite if following hold: This completes the proof of the theorem (6).

Numerical Simulation
To substantiate the above analytical findings, the model is studied numerically.The system of differential equation is integrated using fourth order Runge-Kutta method under the following set of compatible (hypothetical) parameters, which satisfy the stability conditions.It is found that all the conditions for local asymptotic stability of interior equilibrium point are satisfied for above parameter values.Further, to illustrate the global stability of the equilibrium point graphically, numerical simulation is performed for different initial starts and the result is displayed in Figure 1.It is found from the graph that all the trajectories starting from different initial starts, reach the endemic equilibrium 3 E .In Figure 2 densities of tumor cells for parameter values (7.1) is shown.It is observed from the figure that infected tumor cell population first rise and then attain a constant equilibrium value whereas uninfected tumor cell population first rise abruptly and then decrease due to virus infection and cytotoxicity to attain its equilibrium value.
In Figure 3, density of tumor cells is drawn for

 
and other parameters remaining same as in Equation (7.1).This figure implies that the system converges to equilibrium point Numerically it is found that tumor load decreases when per capita growth rate of infected cells ( 2 r ) is less than or equal to the death rate of infected cells due to virus (  ) for high replication rate of viruses in the cells.   and other parameters fixed as above.Figure 12 shows table limit cycle for the range 10 30    .Further, it is found that no limit cycles exist for 30   .

Conclusions
It is well known that the cancer is one of the greatest killers in the world and the control of tumor growth requires great attention.Various efforts have been made for its treatment.Equally extensive efforts have been dedicated over many years to mathematical modeling of cancer development.Here we address a complex process that involves both virus-cell interaction and tumor growth.The positive equilibrium points of the system are investigated.The stability and instability of the equilib-  rium points of the system are studied using the linear stability approach.To substantiate the analytical findings, the model is studied numerically and for which the system of differential equation is integrated using fourth order Runge-Kutta method, which supports the theoreticcal findings.
It is found that virus therapy fails if death rate of infected cells due to virus attack exceeds the transmission rate of infection from virus to the uninfected tumor cells for small immune response of the individual to the virus action.It is observed from our analysis that tumor load can be reduced to a lower value under any one of the two conditions: If net growth rate of uninfected cells is less than the rate of transmission of virus infection than  growth rate of infected cells must be greater than the death rate of infected cells caused by the viruses and if net growth rate of uninfected cells is more than the rate of transmission of virus infection then the ratio of net growth rate of uninfected cells to the death rate of infected cells caused by the viruses must exceed the ratio of net growth rate of uninfected cells to rate at which they become infective.Further it is found from our analysis that interior equilibrium point of the system exists under certain condition which is always locally asymptotically stable if we assume that rate of growth of uninfected tumor cells is more than the growth rate of infected tumor cells.
From numerical simulation it is found that tumor load decreases when per capita growth rate of infected cells is less than or equal to the death rate of infected cells due to virus for high replication rate of viruses in the cells.However, tumor grows unaffectedly when death rate of infected cells due to the virus attack exceeds the rate of transmission of infection from virus to the uninfected cells.The system appears to exhibit different attractors and stable limit cycles for some ranges of the system parameters.

1 )
The equilibrium points for this set of parameters are: 5295, 89.4258 E

Figure 1 .
Figure 1.Variation of infected tumor cells with uninfected tumor cells for different initial starts.

Figure 2 .
Figure 2. The densities of tumor cells for all the parameter values given in Equation (7.1).
remaining same as in Equation (7.1).It demonstrates the convergence of the system to equilibrium point   2 0,99.85E for higher infectivity of the oncolytic virus.

Figure 5
Figure 5 shows variation of tumor load with time for different transmission rate of virus infection ( b ).It is

Figure 3 .
Figure 3. Convergence of the system to the equilibrium point   1 100,0 E , for 0.3   and other parameters remaining same as in Equation (7.1).

Figure 4 . 2 
Figure 4. Convergence of the system to the equilibrium point   2 0,99.85E for 0.06 b  and other parameters remaining same as in Equation (7.1).

Figure 7
Figure 7 shows stable limit cycle for 30  b and keeping other parameters fixed as above.However, when 40 b  , tumor load decreases and strange chaotic attractors are obtained.Figure 8 display chaotic attractor for this range.

Figure 8
display chaotic attractor for this range.

Figure 5 .Figure 6 .
Figure 5. Variation of tumor load with time for different transmission rate of infection from virus and parameters 1 40 r  , 100 K  , 2 r 2  , 0.05 a  , 2  

Figure 9
Figure 9 shows that for a particular value of rate of transmission of infection, 50 b  tumor load vanishes when per capita growth rate of infected cells ( 2 r ) is less than or equal to the death rate of infected cells due to virus (  ) but tumor load remains at its maximal size if 2 r   .Now keeping all the parameters fixed at 1 40 r  , 100, K  2 r 2,  0.05, a  50 b  and varying ,  we observe that tumor load increases with increase in alpha and acquires maximum possible load for   50 b    which implies that tumor grows unaffectedly when death rate of infected cells due to the virus attack exceeds the rate of transmission of infection from virus to the uninfected cells.However, tumor load decreases giving periodic oscillations for b   .Figure 10 shows variation of tumor load with time for different values of  .Fig-

Figure 10
shows variation of tumor load with time for different values of  .Fig- ure 11 display three dimensional attractors obtained for the range 3 10

Figure 9 .
Figure 9. Variation of tumor load with time for different growth rate of infected tumor population relative to their mortality rate due to virus with parameters 1 40, r  100 K  , 0.05 a  , 50 b  .

Figure 10 .
Figure 10.Variation of tumor load with time for different values of  and parameters 1 40, r  100, K  2 r 1, 

Figure 12 .
Figure 12.Stable limit cycle for the set of parameters: 1 40 r  , 100 K  , 2 r 2  , 0.05 a  , 50 b  , 20   and This equilibrium point implies the complete elimination of tumor.Biologically, it means that both infected and uninfected tumor cells can be eliminated with time, and complete recovery is possible because of the virus therapy.