**Applied Mathematics** Vol.3 No.9(2012), Article ID:23020,8 pages DOI:10.4236/am.2012.39160

Non-Linear Mathematical Model of the Interaction between Tumor and Oncolytic Viruses

^{1}Department of Mathematics, The Madura College, Madurai, India^{}

^{2}Department of Mathematics, V V College of Engineering, Tisaiyanvilai, India

Email: ^{*}raj_sms@rediffmail.com

Received June 23, 2012; revised July 25, 2012; accepted August 2, 2012

**Keywords:** Mathematical Modeling; Non-Linear Differential Equations; Numerical Simulation; Homotopy Analysis Method; Tumor Cells; Oncolytic Viruses

ABSTRACT

A mathematical modeling of tumor therapy with oncolytic viruses is discussed. The model consists of two coupled, deterministic differential equations allowing for cell reproduction and death, and cell infection. The model is one of the conceptual mathematical models of tumor growth that treat a tumor as a dynamic society of interacting cells. In this paper, we obtain an approximate analytical expression of uninfected and infected cell population by solving the nonlinear equations using Homotopy analysis method (HAM). Furthermore, the results are compared with the numerical simulation of the problem using Matlab program. The obtained results are valid for the whole solution domain.

1. Introduction

Oncolytic viruses are viruses that infect and kill cancer cells but not normal cells [1-4]. Oncolytic virus therapy originated early in the last century upon the observation of occasional tumor regressions in cancer patients suffering from virus infections or those receiving vaccinations. Many types of oncolytic viruses have been studied as candidate therapeutic agents including adenoviruses, herpes viruses, reoviruses, paramyxoviruses, retroviruses, and others [2,4]. Probably, the best-characterized oncolytic virus, that has drawn a lot of attention, is ONYX- 015, an attenuated adenovirus that selectively infects tumor cells with a defect in the p53 gene [3]. This virus has been shown to possess significant antitumor activity and has proven relatively effective at reducing or eliminating tumors in clinical trials [5-7]. Furthermore, a small number of patients who were treated with the oncolytic virus showed regression of metastases [2]. Although safety and efficacy remain substantial concerns, several other oncolytic viruses acting on different principles, including tumor-specific transcription of the viral genome, have been developed, and some of these viruses have entered in trials [2,8-10].

The oncolytic effect has several possible mechanisms that yield complex results. The first such mechanism is the result of viral replication within the cell and rupture out of the cell [11,12]. The third mechanism involves virus infection of cancer cells that induces antitumoral immunity. Surviving mice acquired a resistance to rechallenge with tumor cells [13]. Host immune response maximizes antitumor immunity but also interferes with virus propagation. Wakimoto et al. [14] studied the limitation of virus propagation caused by host immune response in the central nervous system. Ikeda et al. [15] showed that the viral survival term was prolonged and that virus propagation was increased by the anti-immune drug, cyclophosphamide.

Several mathematical models that describe the evolution of tumors under viral injection were recently developed. Wodarz [13,16] presented a mathematical model that describes interaction between two types of tumor cells (the cells that are infected by the virus and the cells that are not infected by the virus) and the immune system. However, to the best of our knowledge, till date no general analytical expressions for the mathematical modeling of two populations of cells namely uninfected tumor cells and infected cells [17]. The purpose of this paper is to derive approximate analytical expression of two types of cells growing namely uninfected and infected tumor cells by solving the non-linear differential equations using Homotopy analysis method [18-20].

2. Mathematical Models

The model, which considers two types of cells growing in the logistic fashion, has the following form [17]:

, (1)

(2)

The equation can be solved subject to the following initial conditions:

(3)

(4)

where X is the size of the uninfected cell population; Y is the size of the infected cell population; and are the maximum per capita growth rates of uninfected and infected cells; K is the carrying capacity; b is the transmission coefficient; and a is the rate of infected cell killing by virus. We introduce the following set of dimensionless variables as follows [17]:

(5)

The governing non-linear differential equations (Equations (1) and (2)) expressed in the following non-dimensionless format:

(6)

(7)

An appropriate set of boundary condition are given by:

, (8)

3. Solution of Boundary Value Problem Using Homotopy Analysis Method

The Homotopy analysis method (HAM) is a powerful and easy-to-use analytic tool for nonlinear problems [21- 23]. It contains the auxiliary parameter h, which provides us with a simple way to adjust and control the convergence region of solution series. Furthermore, the obtained result is of high accuracy. Solving the Equations (6) and (7) using HAM (see Appendix A) and simultaneous equation method (see Appendix B and C), the steady state and transient contributions to the model are given by:

(9)

where

(10)

(11)

Similarly we can obtain as follows:

(12)

where

(13)

(14)

Here and represent a time independent steady state term and and denote the time dependent transient component. The steady state term will be important at long times as In contrast the transient term will be of important at short times as

4. Numerical Simulation

The non-linear equations (Equations (6) and (7)) for the boundary conditions (Equation (8)) are solved by numerically. The function ode45 in Scilab software is used to solve two-point boundary value problems (BVPs) for ordinary differential equations. The Matlab program is also given in Appendix C. The numerical results are also compared with the obtained analytical expressions (Equations (9) and (12)) for all values of parameters, , , and.

5. Results and Discussion

Equations (9) and (12) represent the simple approximate analytical expressions of the uninfected and infected tumor cells for all values of parameters, , , and. The two types of tumor cell growing using Equations (9) and (12) are represented in Figures 1-4. In Figures 1 and 2, the dimensionless uninfected tumor cells reach the constant value when for some fixed value of and different values of and. The dimensionless infected tumor cells reaches the

Figure 1. Size of the dimensionless uninfected cell population X^{*}(t^{*}) are plotted using Equation (9) for the values A_{0} = 10, B_{0} = 2, δ = 2, γ = 1, h = −0.1 (a) β = 5 () (b) β = 10 () (c) β = 15 () and (d) β = 20 ().

Figure 2. Size of the dimensionless uninfected cell population X^{*}(t^{*}) are plotted using Equation (9) for the values A_{0} = 10, B_{0} = 2, γ = 1, β = 6, h = −0.21 (a) δ = 2 () (b) δ = 6 () (c) δ = 10 () and (d) δ = 15 ().

steady state value when.

In Figures 3 and 4, the dimensionless infected tumor cells reach the constant value when for some fixed value of and different values of and. The dimensionless infected tumor cells reaches the steady state value when. Figures 5 and 6 give us the confirmation for the above discussion in three-dimensional graphs also.

Figure 3. Size of the dimensionless infected cell population Y^{*}(t^{*}) are plotted using Equation (12) for the values A_{0} = 10, B_{0} = 2, β = 2, γ = 1, h = −0.05, (a) δ = 5 () (b) δ = 10 () (c) δ = 15 () and (d) δ = 20 ().

Figure 4. Size of the dimensionless infected cell population Y^{*}(t^{*}) are plotted using Equation (12) for the values A_{0} = 10, B_{0} = 2, δ = 50, γ = 1, h = −0.75 (a) β = 1.5 () (b) β = 5 () (c) β = 10 () and (d) β = 20 ().

Figure 5. The normalized three-dimensional of dimensionless uninfected tumor cells X^{*}(t^{*}) (Equation (9)) for t^{*} = 0 to 1.

Figure 6. The normalized three-dimensional of dimensionless infected tumor cells Y^{*}(t^{*}) (Equation (12)) for t^{*} = 0 to 0.1.

6. Conclusion

In this work, the coupled system of time dependent differential equations for the two types of cells growing has been solved analytically using the Homotopy Analysis Method. Approximate analytical expressions for uninfected and infected cell population are derived for all values of parameters. Furthermore, on the basis of the outcome of this work, it is possible to calculate the approximate rate of the tumor cells growth. The extension of the procedure to other systems that include interaction between tumor cells and anticancer agents seem possible.

7. Acknowledgements

This work was supported by the Council of Scientific and Industrial Research (CSIR No.: 01(2442)/10/EMR-II), Government of India. The authors also thank Mr. M.S. Meenakshisundaram, Secretary, The Madura College Board, Dr. R. Murali, The Principal, and Prof. S. Thigarajan, HOD, Department of Mathematics, The Madura College, Madurai, Tamilnadu, India for their constant encouragement. The authors S. Usha is very thankful to the Manonmaniam Sundaranar University, Tirunelveli for allowing the research work.

REFERENCES

- D. H. Kirn and F. McCormick, “Replicating Viruses as Selective Cancer Therapeutics,” Molecular Medicine Today, Vol. 2, No. 12, 1996, pp. 519-527. doi:10.1016/S1357-4310(97)81456-6
- K. A. Parato, D. Senger, P. A. Forsyth and J. C. Bell, “Recent Progress in the Battle between Oncolytic Viruses and Tumours,” Nature Reviews Cancer, Vol. 5, No. 12, 2005, pp. 965-976. doi:10.1038/nrc1750
- F. McCormick, “Cancer-Specific Viruses and the Development of ONYX-015,” Cancer Biology and Therapy, Vol. 2, No. 4, 2003, pp. S157-S160.
- H. Kasuya, S. Takeda, S. Nomoto and A. Nakao, “The Potential of Oncolytic Virus Therapy for Pancreatic Cancer,” Cancer Gene Therapy, Vol. 12, No. 9, 2005, pp. 725-736. doi:10.1038/sj.cgt.7700830
- D. Kirn, T. Hermiston and F. McCormick, “ONYX-015: Clinical Data Are Encouraging,” Nature Medicine, Vol. 4, No. 12, 1998, pp. 1341-1342. doi:10.1038/3902
- F. R. Khuri, J. Nemunaitis, I. Ganly, J. Arseneau, I. F. Tannock, L. Romel, M. Gore, J. Ironside, R. H. MacDougall, C. Heise, B. Randlev, A. M. Gillenwater, P. Bruso, S. B. Kaye, W. K. Hong and D. H. Kirn, “A Controlled Trial of Intratumoral ONYX-015, a SelectivelyReplicating Adenovirus, in Combination with Cisplatin and 5-Fluorouracil in Patients with Recurrent Head and Neck Cancer,” Nature Medicine, Vol. 6, No. 8, 2000, pp. 879-885. doi:10.1038/78638
- J. Nemunaitis, F. Khuri, I. Ganly, J. Arseneau, M. Posner, E. Vokes, J. Kuhn, T. McCarty, S. Landers, A. Blackburn, L. Romel, B. Randlev, S. Kaye and D. Kirn, “Phase II Trial of Intratumoral Administration of ONYX-015, a Replication-Selective Adenovirus, in Patients with Refractory Head and Neck Cancer,” Journal of Clinical Oncology, Vol. 19, No. 2, 2001, pp. 289-298.
- A. C. Shah, D. Benos, G. Y. Gillespie and J. M. Markert, “Oncolytic Viruses: Clinical Applications as Vectors for the Treatment of Malignant Gliomas,” Journal of Neurooncology, Vol. 65, No. 3, 2003, pp. 203-226. doi:10.1023/B:NEON.0000003651.97832.6c
- H. L. Kaufman, G. Deraffele, J. Mitcham, D. Moroziewicz, S. M. Cohen, K. S. Hurst-Wicker, K. Cheung, D. S. Lee, J. Divito, M. Voulo, J. Donovan, K. Dolan, K. Manson, D. Panicali, E. Wang, H. Horig and F. M. Marincola, “Targeting the Local Tumor Microenvironment with Vaccinia Virus Expressing B7.1 for the Treatment of Melanoma,” Journal of Clinical Investigation, Vol. 115, No. 7, 2005, pp. 1903-1912. doi:10.1172/JCI24624
- T. Reid, R. Warren and D. Kirn, “Intravascular Adenoviral Agents in Cancer Patients: Lessons from Clinical Trials,” Cancer Gene Therapy, Vol. 9, No. 12, 2002, pp. 979-986. doi:10.1038/sj.cgt.7700539
- D. L. Lichtenstein, K. Toth, K. Doronin, et al., “Functions and Mechanisms of Action of the Adenovirus E3 Proteins,” International Review of Immunology, Vol. 23, 2004, pp. 75-111. doi:10.1080/08830180490265556
- A. Zou, I. Atencio, W. M. Huang, et al., “Over Expression of Adenovirus E3-11.6K Protein Induces Cell Killing by Both Caspase-Dependent and Caspase-Independent Mechanisms,” Virology, Vol. 326, 2004, pp. 240- 249. doi:10.1016/j.virol.2004.06.007
- H. Takakuwa, F. Goshima, N. Nozawa, et al., “Oncolytic Viral Therapy Using a Spontaneously Generated Herpes Simplex Virus Type 1 Variant for Disseminated Peritoneal Tumor in Immunocompetent Mice,” Archives of Virology, Vol. 148, 2003, pp. 813-825. doi:10.1007/s00705-002-0944-x
- H. Wakimoto, P. R. Johnson, D. M. Knipe, et al., “Effects of Innate Immunityon Herpes Simplex Virus and Its Ability to Kill Tumor Cells,” Gene Therapy, Vol. 10, 2003, pp. 983-990. doi:10.1038/sj.gt.3302038
- K. Ikeda, T. Ichikawa, H. Wakimoto, et al., “Oncolytic Virus Therapy of Multiple Tumors in the Brain Requires Suppression of Innate and Elicited Antiviral Responses,” Nature Medicine, Vol. 5, No. 8, 1999, pp. 881-887. doi:10.1038/11320
- D. Wodarz and N. Komarova, “Computational Biology of Cancer,” World Scientific Publishing Company, 2005. http://www.worldscibooks.com/lifesci/5642.html
- G. P. Karev, A. S. Novozhilov and E. V. Koonin, “Mathematical Modeling of Tumor Therapy with Oncolytic Viruses: Effects of Parametric Heterogeneity on Cell Dynamics,” Biology Direct, Vol. 1, No. 30, 2006. doi:10.1186/1745-6150-1-30
- S. J. Liao, “The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems,” Shanghai Jiao Tong University, Shanghai, 1992.
- S. J. Liao, “Beyond Perturbation: Introduction to the Homotopy Analysis Method,” Chapman & Hall/CRC Press, Boca Raton, 2003. doi:10.1201/9780203491164
- S. J. Liao, “A Kind of Approximate Solution Technique Which Does Not Depend upon Small Parameters (II): An Application in Fluid Mechanics,” International Journal of Non-Linear Mechanics, Vol. 32, 1997, pp. 815-822. doi:10.1016/S0020-7462(96)00101-1
- S. J. Liao, “On the Homotopy Analysis Method for NonLinear Problems,” Applied Mathematics and Computation, Vol. 147, 2004, pp. 499-513. doi:10.1016/S0096-3003(02)00790-7
- S. J. Liao and Y. Tan, “A General Approach to Obtain Series Solutions of Nonlinear Differential Equations,” Studies in Applied Mathematics, Vol. 119, 2007, pp. 297-355. doi:10.1111/j.1467-9590.2007.00387.x
- S. J. Liao, “Beyond Perturbation: A Review on the Basic Ideas of the Homotopy Analysis Method and Its Applications,” Advanced Mechanics, Vol. 38, No. 1, 2008, pp. 1-34.

Appendix A: Basic Concept of Liao’s Homotopy Analysis Method (HAM)

Consider the following nonlinear differential equation

(A1)

where N is a nonlinear operator, t denotes an independent variable, is an unknown function. For simplicity, we ignore all boundary or initial conditions, which can be treated in the similar way. By means of generalizing the conventional Homotopy method, Liao constructed the so-called zero-order deformation equation as:

(A2)

where is the embedding parameter, is a nonzero auxiliary parameter, is an auxiliary function, L is an auxiliary linear operator, is an initial guess of, is an unknown function. It is important, that one has great freedom to choose auxiliary unknowns in HAM. Obviously, when and, it holds:

and (A3)

respectively. Thus, as p increases from 0 to 1, the solution varies from the initial guess to the solution. Expanding in Taylor series with respect to p, we have:

(A4)

where

(A5)

If the auxiliary linear operator, the initial guess, the auxiliary parameter h, and the auxiliary function are so properly chosen, the series (A4) converges at p = 1 then we have:

. (A6)

Define the vector

(A7)

Differentiating Equation (A2) for m times with respect to the embedding parameter p, and then setting and finally dividing them by m!, we will have the socalled mth-order deformation equation as:

(A8)

where

(A9)

and

(A10)

Applying on both side of Equation (A8), we get

(A11)

In this way, it is easily to obtain for at order, we have

(A12)

When, we get an accurate approximation of the original Equation (A1). For the convergence of the above method we refer the reader to Liao. If Equation (A1) admits unique solution, then this method will produce the unique solution. If Equation (A1) does not possess unique solution, the HAM will give a solution among many other (possible) solutions.

Appendix B: Steady State Solution

For the case of steady-state condition, the Equations (6) and (7) becomes as follows:

(B1)

(B2)

Solving the above Equations (B1) and (B2), we get

(B3)

and

. (B4)

Thus we can obtain and as in the text (Equations (10) and (13)).

Appendix C: Non-Steady State Solution of the Equations Using the HAM

The given differential equations for the non-steady state condition are of the form as:

(C1)

(C2)

For the transient part, the initial conditions are redefined as

(C3)

(C4)

In order to solve the Equations (C1) and (C2) by means of the HAM, we first construct the Zeroth-order deformation equation by taking,

(C5)

(C6)

We have the solution series as

(C7)

and

(C8)

where

(C9)

Substituting Equations (C7) and (C8) into Equations (C5) and (C6) and comparing the coefficient of like powers of p, we get

(C10)

(C11)

(C12)

(C13)

and so on.

The initial conditions are redefined as

(C14)

(C15)

and

for (C16)

Solving the Equations (C10) and (C10) by using the boundary conditions given in Equations (C14) and (C15), we get

(C17)

and

(C18)

Substituting the Equations (C17) and (C18) in Equations (C12) and (C12) and by using the boundary conditions given in Equation (C16), we get

(C19)

and

(C20)

Adding the Equations (C17) and (C19) and the Equations (C18) and (C20), we get the Equations (9) and (12) as in the text.

Appendix D: MATLAB Program to Find the Numerical Solution of Non-Linear Equations (6) and (7)

function main1 options= odeset('RelTol',1e-6,'Stats','on');

Xo = [10; 2];

tspan = [0,10];

tic

[t,X] = ode45(@TestFunction,tspan,Xo,options);

toc figure hold on plot(t, X(:,1),'blue')

plot(t, X(:,2),'green')

return function [dx_dt]= TestFunction(t,x)

a=16;

b=3;

r=10;

dx_dt(1) =x(1)*(1-(x(1)+x(2)))-(b*x(1)*x(2))/(x(1)+x(2));

dx_dt(2)=r*x(2)*(1-(x(1)+x(2)))+(b*x(1)*x(2))/(x(1)+x(2))-a*x(2);

dx_dt = dx_dt';

Appendix E: Nomenclature

Symbol Meaning

Size of the uninfected cell population

Size of the infected cell population

Rate of infected cell killing by the virus

Transmission coefficient

Maximum per capita growth rates of uninfected cells

Maximum per capita growth rates of infected cells

Carrying capacity

Time

Size of the dimensionless uninfected cell population

Size of the dimensionless infected cell population

Dimensionless time