Modeling of tumor growth undergoing virotherapy

Computers in Biology and Medicine 41 (2011) 922–935

Contents lists available at SciVerse ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Modeling of tumor growth undergoing virotherapy Ricardo Dunia , Thomas F. Edgar Department of Chemical Engineering, The University of Texas at Austin, Austin, TX 78712, United States

a r t i c l e i n f o

a b s t r a c t

Article history: Received 17 May 2010 Accepted 10 August 2011

Tumor growth models subject to virotherapy treatment are analyzed and compared in this paper. Tumor growth conditions are obtained for each model type based on the virus infection rate and immune suppressive drug delivery. Equilibrium conditions resulted into quadratic functions for which the tumor radius remained constant during virotherapy. An irrigation tumor model for virotherapy treatment was also proposed. This model consists of irrigation layers distributed radially along the tumor and attached to a common blood circulation compartment. The irrigation model has similar dynamic and steady state characteristics to the diffusion model, which has been supported by experimental results. The irrigation model considers the immune system cell generation and consumption outside the tumor boundary but inside the blood circulation compartment. These characteristics provide a great potential for advanced cancer treatment applications because therapy dose delivery and immune system measurements can be made at the blood compartment level of the irrigation model. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Virotherapy Cancer treatment Tumor model Irrigation model Tumor growth

1. Introduction Viruses that selectively replicate in tumor cells have recently demonstrated their potential use in cancer treatment [1–3]. These viruses, known as oncolytic viruses, are genetically altered to infect and reproduce inside the confined tumor mass without affecting the surrounding organs. An oncolytic virus consists of tumor-tropic RNA, which will only permit the virus to grow inside cells with a defective antiviral response system, in particular cells that lack a functional interferon response. Among the oncolytic viruses with potential use for virotherapy are the adenovirus Onyx-015 [4], the herpes simplex virus HSV-1 [5] and the Newcastle disease virus, NDV [6]. Although the analysis of lab results indicates that virotherapy has a very promising future for cancer treatment, there are several factors that could diminish its effectiveness. Among these factors is the development of undesirable immune responses that attack the oncolytic virus inside the tumor [7]. The innate immune system can destroy not only free virus particles but also infected tumor cells, which enable the tumor growing process. Therefore, immune suppression drugs like cyclophosphamide (CPA) have been used in conjunction with Oncolytic viruses in order to maintain the effectiveness of virotherapy during cancer treatment [8]. The delivery and blood concentration of such drugs are also limited to avoid the eradication of the immune system and the development of diverse infections in treated patients [9].

 Corresponding author.

E-mail address: [email protected] (R. Dunia). 0010-4825/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2011.08.003

Virotherapy dynamic models have contributed to the development of therapy protocols, tumor growth predictions and cancer treatment [10,11]. The most simple virotherapy models consider a uniform distribution of the different type of tumor cells and virus particles along the solid tumor [12,13]. Such models simplify the dynamic simulation for tumor treatment because the tumor mass is confined and described by a reduced number of dynamic variables or states. However, because the distribution and migration of oncolytic viruses plays an important role in the success of virotherapy, the need to consider a spatial coordinate that determines the uneven distribution of such particles becomes imperative for reliable cancer growth simulation and treatment [14]. Therefore, virotherapy cancer models with spatial coordinates were developed within the last decade by considering the tumor as radial symmetric spheroids [15,16]. Two tumor spheroid models have been recently developed to study the virotherapy effectiveness in conjunction with immune suppression treatment. The first model developed by Friedman et al. [7] and analyzed by Wang and Tian [17] considers a convective radial migration of immune system cells along the solid tumor. This model is referred as convective immunization and it is based on the assumption that immune cells are large enough to be considered as a part of the tumor mass. A second model developed by Tao and Guo [18] considers the immune cells to be small enough to penetrate and migrate by diffusion inside the tumor mass. Migration by diffusion is governed by concentration gradients. The convection and diffusion models are analyzed in this paper and conditions for tumor growth are determined based on the amount of immune suppressive drug delivered and the infection rate. The analysis presented in this paper clarifies the conditions at which virotherapy can avoid cancer growth and

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Nomenclature

Zðr,tÞ Z

A,B,C Dv Dz F k k0 N Nðr,tÞ Pðr,tÞ r R(t) Re R_

b

s t U(r) Vðr,tÞ w Xðr,tÞ Yðr,tÞ

linear system matrices virus diffusivity coefficient immune system cell diffusivity coefficient tumor cell flow variable immune killing rate take-up rate of viruses number of radial elements volumetric fraction of tumor dead cells concentration of cyclophosphamide dimensionless radial coordinate tumor size radius real part radius rate of change in time stimulation rate by infected cells time radial tumor cell velocity density of free virus particles clearance rate of immune cells volumetric fraction of tumor uninfected cells volumetric fraction of tumor infected cells

diminish the risk for metastasis. The larger the immune suppression drug dose, the more effective is the virotherapy, which will diminish the cancer growth. Moreover, the immune suppressive drug delivery is treated as a manipulated variable while the replication infection rate is a predetermined virus parameter [10]. Therefore, controllability conditions based on linearized models are verified. In the development of convection and diffusion models, the immune system cells are generated and consumed inside the tumor mass. In a similar manner, the immune suppression drug is applied at the tumor level and not conveyed as part of the blood stream system. This paper proposes an immunization model based on tumor irrigation, where the immune cell dynamics is considered outside the tumor mass, and the immune suppression drug is applied at the blood stream level. The concept of immunization by irrigation is introduced in this paper to make realistic the application of the immunization model at the blood stream level. This paper is organized as follows: Section 2 provides the mathematical formulation and solution approach used for the different virotherapy models. Convection and diffusion models are analyzed in Sections 3 and 4, respectively. In Section 5 the irrigation model is presented and analyzed. Such a model provides treatment in a more convenient manner than the models analyzed in previous sections. Section 6 introduces a model comparison in view of cancer treatment effectiveness. Conclusions and future work are provided in Section 7 of this paper.

2. Mathematical formulation The mathematical formulation of the virotherapy treatment for spherical tumors results in a nonlinear PDE system, with partial derivatives along the radial coordinate r. The azimuthal and polar angle coordinates are not considered due to radial symmetry assumption, which simplifies the problem without affecting significantly the accuracy of the solution. The size of the tumor is given by R(t), where t represents time. Because the tumor changes its size during treatment, any evaluation of an external boundary condition also changes with time.

g d y l

m r te tt f

923

volumetric fraction or density of immune system cells immune cell system density infection rate clearance rate of viruses infected cell lysis rate density of tumor cells proliferation rate of tumor cells removal rate of necrotic cells radial coordinate time constant at the blood circulation level time constant at the tumor level virus initial distribution parameter

Subscripts c d e i n R

convection diffusion external irrigation nominal at the radius

Nevertheless, this inconvenience is eliminated by defining the radial coordinate as r  r=RðtÞ, where r ¼1 represents the tumor surface. The variables used to describe the tumor progress and its interaction with the oncolytic virus and the immune system are:

 The volumetric fraction of uninfected tumor cells, X: this     

fraction represents the major cause of tumor enlargement and cancer metastasis. The volumetric fraction of infected tumor cells, Y: this fraction represents the tumor cells attacked by the virus. Infected tumor cells eventually become dead cells. The volumetric fraction of dead tumor cells, N: the increase of this fraction tends to slow down the tumor enlargement. The immune system concentration, Z: this variable is considered a volumetric fraction or a density in the convection or diffusion model, respectively. The density of the free virus particles, V: the virus particles are injected in the tumor at the start of the virotherapy treatment. The radial cell velocity, U: the propagation velocity of the tumor cells determines the rate at which the tumor enlarges or shrinks.

No analytical solution is available for the partial differential equations obtained for the different type of models. The numerical solution approach considers radial discretization using shellshaped finite elements. Fig. 1 illustrates how these shells are defined and distributed. A node is placed at every shell boundary to assure radial profile smoothness and continuity, and the node distribution is made such that each shell element holds the same volume. Some of the radial boundary conditions are a function of the model type in study, and no virus particle is expected to leave the tumor surface. The discretization of the radial coordinate leads to a system of differential equations in time. Each inner node provides differential equations in X, Y, Z, N and V. Even though there exists a linear dependency among the tumor fraction time derivatives, the condition that all fractions add to one is verified at each integration step. Variable step Runge-Kutta 23 (RK23) is used to integrate such a system of differential equations, with integration steps no larger than a minute.

924

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Table 1 Model parameters.

l y k s w Dv Dz

d k0

g m Zn

Fig. 1. Spheroidal tumor model considered for virotherapy analysis. Twenty-one radial nodes define a total of 20 shells. Nodes are distributed such that shells contain identical volumes.

The tumor growth is based on the convective velocity evaluated at the tumor surface, U(R). A positive velocity at the boundary indicates tumor enlargement. In general, the tumor radius and the convective velocity are related by R_ ¼ UðRÞ. The next sections provide details about the type of assumptions made for the different published models and the irrigation model developed in this paper.

3. Convection model The convection model considers immune system cells as part of the tumor tissue, X þY þNþZ ¼ 1

ð1Þ

where all tumor cells migrate by convection at the same velocity rate, U. The virus migration is based on a diffusion mechanism due to the size and nature of the virus particles. The system of partial differential equations that govern the convection model were developed in [7] and are given below ! @X Ur R_ @X þ ¼ lXbVXFc X ð2Þ @t R @r @Y Ur R_ þ @t R

!

@N Ur R_ þ @t R @Z Ur R_ þ @t R

!

!

@Y ¼ bVXkYZdYFc Y @r

ð3Þ

Value

Units

Description

2.0 106 2.0 56 20 3.6 0.036 5.6 1.0 25 2.1 0.06

10  2 l/h cells/mm3 10  8 mm3/h immune cell 10  8 mm3/h infected cell mm3/h cell 10  2 mm2/h 10  2 mm2/h 10  2 l/h 10  8 mm3/h immune cell 10  8 l/h 10  2 l/h Immune cell/mm3

Proliferation rate of tumor cells Density of tumor cells Immune killing rate Stimulation rate by infected cells Clearance rate of immune cells Virus diffusivity Immune cells diffusivity Infected cells lysis rate Take-up rate of viruses Clearance rate of viruses Removal rate of necrotic cells Nominal immune cell density

and P(t) represents the immune suppression drug dose during virotherapy treatment. Table 1 provides the parameters required to solve this PDE system. The initial conditions for the tumor cell fractions are R(0)¼ 2 mm, Xðr,0Þ ¼ 0:84, Yðr,0Þ ¼ 0:10, Zðr,0Þ ¼ 0:06. The initial virus profile is given by Vð0,rÞ ¼ fe4r

2

=f

2

ð9Þ

where f is obtained from the following expression: Z R 2 2 f r 2 e4r =f dr ¼ 0:45 0

These initial conditions, which includes a Gaussian distribution for the initial virus profile, were chosen to match experimental results [17]. Note that the velocity profile follows a quasi-steady state condition as there is no implicit time dependency in Eq. (7), but U(r) changes in time due to variations in Fc. On the other hand, radial symmetry at r¼0 indicates that @V=@rð0,tÞ ¼ 0 and Uð0,tÞ ¼ 0. Finally, virus particles do not cross the tumor surface, which indicates that @V=@rð1,tÞ ¼ 0. Equilibrium points are calculated by setting the right-hand side of Eqs. (2) to (6) to zero, which results in the following system of nonlinear algebraic equations: V¼

lF c b

ð10Þ

Z¼

ðF c þ PÞdsV g k0 sV wd

ð11Þ

Y¼

X¼

V ðk0 Z þ gÞ

ð12Þ

d Y ðkZ þ d þ F c Þ

ð13Þ

bV

N ¼ 1X Y Z @N ¼ kYZ þ dYmNFc N @r

ð4Þ

@Z ¼ ½sYwZFc PðtÞZ @r

ð5Þ

@V r R_ 2Dv  þ 2 @t R R r

!

@V Dv @2 V  ¼ dYk0 VZgV @r R2 @r 2

ð6Þ

ð14Þ

where X , Y , Z , N and V denote equilibrium values for the respective variables. The amount of drug delivered P as well as the infection rate parameter b can be adjusted to provide desirable equilibrium fractions. It is important to notice that because all cell fractions are limited to the range [0 1], not any ðP, bÞ pair provides feasible equilibrium fractions. The quantity F c represents Fc evaluated at the equilibrium conditions: 2

1 @ 2 ðr UÞ ¼ Fc Rr2 @r

F c ¼ lX mN þ sY Z wZ PZ ð7Þ

where Fc ¼ lXmN þ sYZwZ 2 PðtÞZ

ð8Þ

ð15Þ

Note that F c is not necessarily zero even though the cell fractions have approached equilibrium conditions. The calculation of the velocity profile at equilibrium conditions determines if virotherapy has successfully stopped the growth of the tumor. The solution of

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

925

The general expressions for the coefficients ac, bc and cc are provided in Appendix A.1. Substitution of the parameter values given in Table 1 results in polynomial coefficients that only depend on the infection rate parameter b. Therefore, the condition F c ¼ 0 is converted into a relation of the form PðbÞ. This relation defines the Static Line (SL) of tumor growth, for which R_ ¼ 0. Fig. 2 illustrates the SL based on P and b. The SL indicates treatment points for which F c ¼ 0, which makes R_ ¼ 0. In general, a larger b increases virus replication, which requires a lesser amount of immune suppression drug for R_ ¼ 0. Any point to the left or below of SL, like point A, represents a condition for which the tumor grows because the amount of drug being delivered is not enough to attenuate the immune cells. The opposite occurs

for points to the right or above SL, as it is the case of point C. The tumor shrinks if at a given b the amount of drug to suppress the immune system surpasses the one specified by SL. Point B is at SL, which indicates that under such a treatment the tumor radius will eventually reach steady state, R. Notice also that in the absence of CPA, b Z 9:77 for F c r0. Because typical values for b are around seven [10], CPA is usually delivered during virotherapy. For this reason virotherapy is also applied in conjunction with other types of cancer treatments [19]. Fig. 3 illustrates the dynamic response of the tumor size R(t) at points fA,B,Cg, located at contour lines F c ¼ f þ 0:5,0,0:5g in Fig. 2. All three points share the same infection rate parameter b ¼ 7:45, with CPA doses of P ¼ f0:23,4:21,8:7g, respectively. The dynamic responses for points A and C in Fig. 3 show the enlargement and shrinking of the tumor, while the curve associated to point B illustrates the response toward a steady state radius of 1.25 mm. These time responses were expected due to the relative location of points fA,B,Cg with respect to SL. A plot that corroborates the effect of F c into R_ is the phase diagram for the tumor radius, R. Fig. 4 represents the phase diagram R_ vs. R for points fA,B,Cg, located at contour lines _ F c ¼ f þ0:5,0,0:5g, respectively. Note that for F c ¼ 0, R-0 after 500 h of treatment. For point C, where F c ¼ 0:5, the rate of change for the tumor radius is negative and steady after 500 h of treatment. Such an outcome is expected because R_ is proportional to F c after cell fractions approach equilibrium conditions. Finally for point A, the tumor radius increases because F c 4 0. The steady-state conditions for the tumor radius were developed in terms of F c . Nevertheless, to reach F c , the tumor cell fractions require that equilibrium conditions be achieved. Fig. 5 demonstrates the average cell fraction variations with respect to the tumor radius at treatment point C. This figure also shows the cell fraction progress for 100, 500, 1000 and 1500 h of treatment, and its horizontal axis R has been inverted to illustrate the treatment progress. Note that the tumor radius continues to decrease even though cell fractions are inside the equilibrium region represented by the shade portion of the plot. This is because point C is at the right of SL. Even though stability is verified for point C, local stability should be satisfied in the

Fig. 2. Static Line (SL) for the convection model. Point A is below SL, which indicates tumor enlargement during treatment. The contrary occurs at the region above SL (point C), where the tumor shrinks when undergoing virotherapy. Treatment points at the SL, like point B, keep the tumor size constant after days of treatment.

Fig. 3. Dynamic responses for the tumor radius at the conditions given by points fA,B,Cg in Fig. 2.

the velocity profile U(r) for a constant F c is denoted by U ðrÞ and is obtained from the integration of Eq. (7) with respect to r: U ðrÞ ¼ 13F c RðtÞr

ð16Þ

The evaluation of this last expression at the tumor surface (r¼1) provides an approximation of the rate of change of the tumor radius after the cell fractions and virus density approach equilibrium: _ ¼ U ð1Þ ¼ 1F c RðtÞ RðtÞ 3

ð17Þ

Because F c is not necessarily zero, R(t) will not necessarily reach steady state even though the tumor cell fractions are at equilibrium. In such a situation, the tumor will continuously grow or shrink but its volumetric fractions will remain steady. From Eq. (17), the tumor size reaches steady state (R_ ¼ 0) either because F c ¼ 0 or for F c o0, which makes Rðt-1Þ ¼ 0. Nevertheless, the first condition (F c ¼ 0) provides more insight to short-term virotherapy treatment than making the tumor vanish after long treatment periods. Therefore, the condition F c ¼ 0 imposes an additional relation to the equilibrium values, where 2

lX mN þ sY Z wZ PZ ¼ 0

ð18Þ

The substitution of the equilibrium expressions in Eq. (18) provides a second-order polynomial in P: ac P 2 þ bc P þ cc ¼ 0

ð19Þ

926

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Fig. 4. Phase diagram for points fA,B,Cg in Fig. 2. The diagram represents 1500 h of treatment. The responses start at R¼ 2 mm, and the symbol  indicates conditions at 500 h of treatment.

Fig. 5. Cell fraction average X^ , Y^ and N^ with respect to R at treatment point C in Fig. 2. Tumor radius continues to decrease even though cell fractions are around their equilibrium values.

equilibrium region for F c r0 in order to demonstrate that all equilibrium conditions correspond to attainable stable points. To demonstrate local stability at the equilibrium region, the right-hand sides (rhs) of Eqs. (2) to (6) are linearized around equilibrium points. Appendix A.2 provides the corresponding systems matrices Ac and Bc for the linearized convection model. Substitution of the parameters provided in Table 1 in the expressions presented in Appendix A.2 give expressions for matrices Ac and Bc as a function of the equilibrium fractions. In order to evaluate Ac and Bc, the feasible range of the equilibrium fractions needs to be considered based on F c . Fig. 6 shows the equilibrium values for the tumor cell fractions at the Steady Line (SL) as a function of b. Note that b should be above 5.77 to prevent infeasible equilibrium values for Z . An important feature in Fig. 6 is that X is always greater than N. Such a situation

Fig. 6. Equilibrium values for the tumor cell fractions along SL as a function of b. The fraction of uninfected tumor cells X is always greater than the fraction of dead tumor cells N .

Fig. 7. Equilibrium values for the tumor cell fractions along F c ¼ 0:5 as a function of b. The fraction of dead tumor cells N is always above fraction of uninfected tumor cells X .

increases the risk of metastasis. The opposite occurs for F c ¼ 0:5 in Fig. 7, where N is significantly greater than X . Based on Figs. 6 and 7, the minimum values allowed for b are 5.77 and 7.35 for SL and F c ¼ 0:5, respectively. The upper limits for b are provided by the condition of P Z 0 in Fig. 2, which indicates b less than 9.77 and 12.76 for SL and F c ¼ 0:5, respectively. Such limits define the evaluation range for Ac in order to establish local stability conditions. Figs. 8 and 9 illustrate the location of the real part of the Ac eigenvalues as a function of b along SL and F c ¼ 0:5, respectively. Note that all real parts remain negative, which indicates that all equilibrium points along SL and F c ¼ 0:5 are locally stable. The existence of a stable conjugate pair of eigenvalues in Figs. 8 and 9 indicates oscillations in dynamic responses, as it is visualized in Fig. 3. Controllability properties were also verified for the pair ðAc ,Bc Þ at SL and F c ¼ 0:5, which indicates that P can

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

927

Z and the tumor cells density y: Z¼

Z

ð21Þ

y

where y is assumed constant [16] during tumor treatment. This assumption allows keeping the same set of model parameters used in the convection case for the diffusion model. The immunization cells, as well as the virus cells, migrate by diffusion along the tumor radial direction. Therefore Eq. (5), which describes the convective migration of the immune cells inside the tumor, is substituted by the following equation: ! @Z r R_ 2Dz @Z Dz @2 Z  þ 2  ¼ sYZwZ 2 PðtÞZ ð22Þ @t R R r @r R2 @r 2 where Dz represents the diffusion coefficient of the immune system cells. This coefficient is about two orders of magnitude less than the virus diffusion coefficient Dv. Symmetry and boundary conditions associated to Eq. (22) are similar to the ones applied for Eq. (6). An important adjustment in the application of these equations to the diffusion model is the substitution of Fc by Fd, where Fd ¼ lXmN Fig. 8. Real part of Ac eigenvalues at the Static Line (SL) as a function of b. All points along SL have negative real part eigenvalues. A pair of complex conjugate eigenvalues is always present.

ð23Þ

The assumptions made above are consistent with the work developed by Tao [18], in which the immune cell migration is based on a diffusion mechanism. In a similar way to the convection model, the equilibrium expressions for the immunization, virus density and tumor cell fractions are obtained by setting the right-hand side of all PDEs to zero: V ¼ ðlF d Þ=b

ð24Þ

Z¼

sV gP d wdk0 sV

ð25Þ

Y¼

P þ wZ s

ð26Þ

X¼

Y ðkZ þ d þ F d Þ

bV

N ¼ 1X Y

ð27Þ ð28Þ

where F d ¼ lX mN

Fig. 9. Real part of Ac eigenvalues for F c ¼ 0:5 as a function of b. All points have negative real part eigenvalues, and a pair of complex conjugate eigenvalues is always present.

Applying similar reasoning as for the convection model, the tumor radius will reach steady state when F d ¼ 0. Substitution of Eqs. (29), (27), (26), (28) in this last identity provides a secondorder polynomial of the form: ad P 2 þ bd P þ cd ¼ 0

effectively be used as the manipulated variable to maintain specific treatment conditions around equilibrium points.

4. Diffusion model The immunization cells are not considered part of the tumor tissue in the diffusion model. Therefore, only the infected (Y), uninfected (X) and dead tumor cell fractions (N) are taken into account in the overall tumor mass balance: X þY þN ¼ 1

ð20Þ

This assumption is consistent with the notion of having the immune system cells penetrating the tumor tissue. Therefore, Z does not represent a fraction with an upper bound limit of one in the diffusion model, but a ratio between the immune cells density

ð29Þ

ð30Þ

where the general expression for the quadratic polynomial coefficients is provided in Appendix B.1. Substitution of the parameters in Table 1 makes the polynomial coefficients only dependent on b. The Static Line (SL) for the diffusion case is obtained by calculating the solution of the quadratic polynomial for different values of b. Notice the similarity between the convection and diffusion models when comparing Fig. 2 with Fig. 10. The Static Line as well as the contour lines in the diffusion model tend to be steeper than in the convective model. This difference makes the control of tumor growth more difficult in the diffusion case when using P as manipulated input. Nevertheless, larger values of b are required in the convective case to shrink the tumor when compare to the diffusion case. Points fA,B,Cg in Fig. 10 share the same drug delivery dose (P ¼4) but demonstrate different types of treatment responses. Point A is to the left of SL (where Fd 4 0) and lacks an infection

928

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Fig. 10. Static Line (SL) for tumor treatment in the diffusion model. Points fA,B,Cg illustrate cases in which the tumor size growths, stabilizes or shrinks in time, respectively. Their respective b values are f5:418,6:645,8:392g.

Fig. 12. Tumor cell equilibrium fractions for the diffusion model along SL. The fraction of uninfected tumor cells X is always greater than the fraction of dead tumor cells N .

Fig. 13. Tumor cell equilibrium fractions for the diffusion case along F d ¼ 0:5. The fraction of dead tumor cells N is always above the fraction of uninfected tumor cells X . Fig. 11. Dynamic responses for the tumor radius at the conditions given by points fA,B,Cg in Fig. 10. In case B, the steady state radius is 1.29 mm.

rate level capable to prevent tumor enlargement. On the contrary, point C is on the contour line F d ¼ 0:5 at the right of SL. A negative F d makes R_ o0 after few hours of treatment due to the fast proliferation of the virus infection. Finally, points at SL (like point B) reach a steady radius after few hours of treatment. Fig. 11 illustrates the dynamic response of the tumor radius for the points highlighted in Fig. 10. Larger b increases the chances to avoid tumor enlargement for a given drug delivery amount P. The equilibrium cell fraction profiles along SL (F d ¼ 0) and along F d ¼ 0:5 are given in Figs. 12 and 13, respectively. Note that X is drastically reduced compared to N when the treatment point moves to the right of SL. This effect results in the same conclusion as for the convective model: the more negative is F d the less is the risk of metastasis.

The limits required to obtain feasible equilibrium fractions and positive drug delivery are considered to study stability conditions in the vicinity of SL and contour line F d ¼ 0:5 of Fig. 10. These limits indicate that the infection rate parameter b should be contained in [5.77, 7.45] and [7.30, 9.42] for the SL and 0.5 contour line, respectively. Appendix B.2 provides the expressions for the system matrices Ad and Bd to determine the dynamic stability and controllability of the linearized model around equilibrium conditions. Figs. 14 and 15 show that all eigenvalues (real part) of matrix Ad are negative, which indicates that all feasible responses are locally stable. These figures also demonstrate the presence of a conjugate pair of eigenvalues that provide response oscillations, like the ones shown in Fig. 10. Finally, local controllability conditions were satisfied for the pair ðAd ,Bd Þ at all feasible b values along SL and F d ¼ 0:5.

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Fig. 14. Real part of Ad eigenvalues at the Static Line (SL) as a function of b. All points along SL have negative real part eigenvalues. A pair of complex conjugate eigenvalues is always present.

929

with respect to the diffusion model is the application of CPA and generation of immune cells at the blood stream level. Therefore, the expression that corresponds for the migration of immune cells Z is different for all three models. Although the immune cell migration in the irrigation model is also based on a diffusion expression, it partitions the spherical tumor into shells. Such shells relocate the generation of immune system cells to the blood circulation level. Therefore, the generation term for immune cells occurs not inside the tumor, as it was considered in the convection and diffusion models, but at the blood level compartment. Fig. 16 illustrates the notion of an irrigated tumor divided into spherical shells. Inner tumor shell compartments have both faces exposed to irrigation layers. The core inner portion of the tumor is equivalent to a spheroid, while the most external tumor shell is exposed to a immune cell surface density, ZR. The location of the tumor determines the value of the boundary value ZR. The immune system cells migrate from the blood irrigation layers to the inner tumor tissue by diffusion. The irrigation layers are a subset of the actual nodes used for the calculation of radial profiles. These irrigation layers or nodes are considered in contact with the blood stream. The partial differential equation for the immune cell density at the irrigation nodes is given by: ! @Z r R_ 2Dz @Z Dz @2 Z ðZZe Þ  þ 2  ¼ ð31Þ @t R tt R r @r R2 @r 2 where Ze is the immunization density at the blood stream. The time constant tt is in the order of 1 h and determines the rate at which the irrigation node density Z is impacted by the blood stream density, Ze. A diffusion expression, like the one obtained for the diffusion model, governs the migration of the immunization cells to the non-irrigated nodes: ! @Z r R_ 2Dz @Z Dz @2 Z  þ 2  ¼0 ð32Þ @t R R r @r R2 @r 2 This type of model considers no creation or suppression of the immune system inside the tumor. Such effects take place in the

Fig. 15. Real part of Ad eigenvalues for F d ¼ 0:5 as a function of b. All points have negative real part eigenvalues, and a pair of complex conjugate eigenvalues is always present.

5. Irrigation model The immune system response consists mainly of four phases: recognition, amplification of defense, attack and suppression. Some of these phases take place outside the tumor boundary. Body organs, like the bone marrow and the thymus are the main producers of immunization cells. Other organs like the spleen produces large amounts of antibodies to fight against antigens, and lymph nodes filter the antigens to evoke a full-fledged immune response. Therefore, it is natural to consider that the dominant dynamics of the immune system response, as well as the application of CPA are connected at the blood level compartment to the irrigated tumor. The convection, diffusion and irrigation models share the same equations for X, Y and V. The irrigation model is introduced in this section as a variant of the diffusion model. Its main differentiator

Fig. 16. Irrigation model for virotherapy. The immune system cells circulate in the blood stream and reach the tumor internals through blood irrigation capillarity. Irrigation layers are placed in predefined node locations, and immune cells migrate from irrigation layers to the tumor inner tissue by diffusion.

930

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

innate blood stream compartment outside the tumor confinement, and directly impacts the immune cell density rate of change at the blood stream level:

te

dZ e ¼ ½sY^ wðZe Zn ÞPðtÞZe dt

ð33Þ

where Zn represents the nominal immune cell density at the blood stream compartment in the absence of tumor cells and immune suppression drugs. Such a parameter indicates a minimum immunization cell density at the blood stream level during normal conditions. The time constant te is on the order of 50 h, typical for the time required by CPA to impact the immune system. The average infected cell fraction Y^ is used in the place of the local fraction Y to account for an overall effect of the infected cells on the immunization system generation. The equilibrium expressions based on the irrigation model for ZðRÞ ¼ Ze are given by Y¼

V ðwgk0 P þk0 wZ n Þ wdk0 sV

Z ¼ Ze Ze ¼

ð34Þ

ð35Þ

sV gP d þ wZ n d wdk0 sV

ð36Þ

where X , V , F d and N have the same expressions as for the diffusion model. The quadratic polynomial form for the irrigation model under steady tumor size is given by: ai P 2 þbi P þ ci ¼ 0

ð37Þ

enlargement (F d ¼ 0:5) to tumor reduction (F d ¼ 0:5) by increasing b and keeping P ¼4. Point B is at the SL and represents the response for b ¼ 6:88, which makes the tumor reach a constant radius R a0 after few hours of treatment. These results are similar to the diffusion model when compared to Fig. 10. The assumption of a nominal density Zn 4 0 increases the requirements for larger b in the irrigation model when compared to the diffusion model. This shift in b is verified by the contour lines in Fig. 17, when compared to the ones in Fig. 10. The dynamic responses for the tumor radius at points fA,B,Cg in Fig. 17 are illustrated in Fig. 18. Although these points share the same immune suppression drug delivery, their b values of f5:61,6:88,8:694g make the responses significantly different. These results also show slow responses compared to the diffusion model responses in Fig. 11. The addition of the blood stream compartment in the irrigation model has increased the time required for the CPA dose to affect the immune system, which will eventually improve the virus replication. A delay term was introduced in the diffusion model presented by Tao [18] to account for the lapse of time that takes for P to have an effect in the immune system, Z. The irrigation model takes into account such an effect by bringing the generation and attenuation of the immune system to the blood stream level. Stability analysis for the cell fraction dynamics is achieved by linearizing the irrigation model around feasible equilibrium conditions. Table 2 shows b ranges to reach equilibrium conditions along SL and F d ¼ 0:5. Appendix C.2 provides the linearized system matrices Ai and Bi to test for stability and controllability of the irrigation model. Figs. 19 and 20 show negative real part eigenvalues for Ai along the feasible range of b. Such a result

where the expressions for the polynomial coefficients are provided in Appendix C.1. The Static Line (SL) in Fig. 17 represents the solution of the quadratic polynomial in Eq. (37) for different b. This solution represents the constant tumor size condition for the irrigation model. Points to the right of SL (F d o 0) define treatment conditions in which the tumor shrinks in time, while the region to the left of SL (F d 40) represents treatment conditions for tumor enlargement. Points fA,B,Cg show the progress from tumor

Fig. 18. Tumor radius dynamic responses at points fA,B,Cg of Fig. 17.

Table 2 Irrigation model b limits. Cell fraction equilibrium is not achieved outside these limits. Fig. 17. Static Line (SL) for the irrigation model as a function of the infection rate parameter, b. Point A is at the left of SL, which indicates tumor enlargement. The contrary occurs at the right of SL (point C), where the tumor shrinks after some hours of treatment. Points at the SL (like point B) prevents tumor enlargement by keeping the tumor size constant after days of treatment.

Fd

Min

Max

 0.5 0

7.30 5.77

9.75 7.70

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Fig. 19. Eigenvalues (real part) for the linearized model along the contour line F d ¼ 0:5 in Fig. 17.

Fig. 20. Eigenvalues (real part) for the linearized model around points along the SL in Fig. 17.

indicates stable dynamic responses along SL and F d ¼ 0:5. A complex pair of eigenvalues is also present in the irrigation model, which makes the response oscillate as it shown in Fig. 18. Stable dynamic responses indicate that Eqs. (34) to (36) provide stable cell equilibrium fractions and immune system density. Figs. 21 and 22 illustrate such steady state values for SL and F d ¼ 0:5, respectively. Notice that the fraction of dead cells N dominates the cell fraction distribution for F d ¼ 0:5. These results match the ones obtained for the convection and diffusion models. Controllability conditions were verified for the pair ðAi ,Bi Þ along SL and F d ¼ 0:5. The immune cell density in the blood stream (Ze) can be measured and used to estimate the remaining states in the irrigation model. Therefore, the system matrix Ci has the form [0 0 0 0 1] and indicates that only Ze is measured. The pair ðAi ,Ci Þ are proved to be observable along (see Appendix B.2). This result indicates the potential use of Ze measurements to estimate the tumor cells fractions toward an optimal treatment.

931

Fig. 21. Tumor cell steady state fractions for the irrigation model along SL.

Fig. 22. Tumor cell steady state fractions for the irrigation model along F d ¼ 0:5.

Next section compares the proposed immune cell migration models based on the virotherapy treatment effect on the cancer growth.

6. Model comparison for cancer treatment The different model assumptions may impact the predicted success of virotherapy. Therefore, it is important to compare them in order to determine which assumption gives the most beneficial result based on the following treatment aspects:

 The relative location of the Static Lines, which determines the amount of CPA required to avoid tumor growth for a given b.

 The dynamic response time, which establishes the response time for treatment progress.

 The reduction of metastasis risk, which determines the chances of cancer propagation in other areas of the host body.

932

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

The comparison will determine predictions for tumor treatment. Furthermore, the fact that all three models share equivalent model parameters helps to determine the effect of the different immune system cell migration in the effectiveness of cancer treatment. 6.1. CPA dose requirements A way to compare the amount of immune suppressive drug requirement for the different models is by plotting the static lines on the same graph. Fig. 23 illustrates how the SL of the convection model tends to require larger doses of CPA than for the diffusion and irrigation models. These last two models have a very similar SL profile shift by a constant value Db ¼ 0:22. Such a small difference is due to the nominal immune suppressive drug density in the blood stream, Zn ¼0.06. Large CPA dose requirements for the convection model indicates that the immune suppression drug is more effective when the immune system cells migrate by a diffusion mechanism. The immune cell diffusivity Dz was made two orders of magnitude lower than the virus diffusivity coefficient, Dv. Such an assumption is in accordance with Tao [18], but might have reduced drastically the effectiveness of the immune system for the diffusion model. In terms of the shape of the PðbÞ function that defines the SL plots, the diffusion and irrigation lines have steeper profiles than the convection one. Steep SL profiles indicate that larger changes in P are required for a given Db. Such a characteristic makes the cancer more difficult to control under the diffusion and irrigation models than using the convection model. 6.2. Treatment response time Treatment response time is related to the effectiveness of virotherapy in reducing the tumor size. For that reason, a treatment point to the right of SL is considered for all three type of models. The treatment point under consideration is given by b ¼ 7 and P¼8 as it lies to the right of the SL for all three models. Such a condition allows to compare radial size reduction for all three models using the same infection rate and drug delivery. Fig. 24 shows the phase diagram for all three models with time stamp locations for 100, 500, 1000 and 1500 h of treatment. The colored dots used for each model tend to be close to each other at every

Fig. 23. Comparison of the convection, diffusion and irrigation Static Lines.

Fig. 24. Phase diagram comparison for the convection, diffusion and irrigation models with time stamp locations of 100, 500, 1000 and 1500 h of treatment. Fractions near to equilibrium conditions are achieved at around 750 h of treatment.

time stamp considered in the figure. Nevertheless, differences can be appreciated, and the diffusion model shows the most favorable progress in tumor reduction size as the green dot looks ahead of the other two models at every stamp time considered. It is important to mention that even though the equilibrium region seems to represent a small portion of the response in the phase diagram, it becomes more than half of the virotherapy treatment in time domain because the last 750 h of therapy occur inside the equilibrium region. 6.3. Risk of metastasis The risk of having metastasis is proportional to the fraction of uninfected tumor cells, X. Therefore, it is important to keep this fraction low during treatment and reduce oscillations that could spark metastasis. The equilibrium fractions provide the final values for the fraction X, while phase diagrams are useful to determine the magnitude of the response oscillations that are detrimental for cancer propagation. Fig. 25 compares X for the different models along SL. The equilibrium fractions of uninfected cells seems to be minimum for the convection model along SL. Irrigation and diffusion models share the same equilibrium fractions because the equations that describe their equilibrium conditions are similar. Furthermore, the larger b, the larger X for the irrigation and diffusion models. The contrary happens to the convection model, where large b values provide lower X . Such a behavior tends to support the use of the convection model for more favorable results to prevent metastasis. A second way to determine the effectiveness of the virotherapy in reducing metastasis is by looking at the oscillations in the phase diagram defined by X^ ðtÞ vs. R(t). Fig. 26 shows that the tumor radius under the convection model has the slowest initial response to the treatment because it goes well above the tumor initial radius of 2 mm. Nevertheless, the treatment is very effective at reducing X^ ðtÞ as soon as the virotherapy starts for all three models. The largest oscillation occurs for the irrigation model close to the 250 h of treatment, where X^ ðtÞ reaches values above 0.65. Therefore, the irrigation model tends to predict larger chances of metastasis that the other two models.

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Fig. 25. Uninfected equilibrium fraction for the convection, diffusion and equilibrium models along SL.

933

virotherapy treatment estimation and control, where the drug delivery takes the place of the manipulated variable and the immune system cell density represents the measured variable. Stability, controllability and observability properties were here verified in a linearized model around equilibrium points. The irrigation model equations were based on the diffusion model [18], with the difference that an external blood stream compartment was created. This blood compartment is in charge of the generation and consumption of immune cells, and interacts with the tumor cells by means of irrigation layers. The irrigation model was also analyzed in terms of its dynamics and equilibrium aspects. Its characteristics resembles the diffusion model because the immune cells migrates from the irrigation layers to internal tumor nodes following a diffusion mechanism. Convection [17], diffusion and irrigation models were compared to determine the effect of virotherapy treatment under the different model assumptions. The convection model seems to require the largest amount of immune suppression drug to prevent tumor growth. However, it provides the lowest uninfected tumor cell equilibrium fractions, which is beneficial to avoid cancer metastasis. The diffusion model provides the most prompt response to the drug therapy, but at the same time it requires larger drug dose changes to adjust and bring the tumor system to desirable treatment points. Nevertheless, there is no major discrepancy among the models, as these adjust accordingly to experimental data collected from animal labs. The controllability and observability properties of the linearized irrigation model are promising for advanced control applications. Recent publications have shown major cancer treatment improvements when using model based control for radiotherapy [20] and chemotherapy [21]. The irrigation model accounts for possible manipulated and measurement variables taken at the blood stream level, which could allow continuous monitoring and control action.

Acknowledgments The authors acknowledge the suggestions made by the reviewers and the guidance given by Dr. Marisol Fernandez from the infectious disease staff at Dell Children’s Hospital in Austin, Texas.

Appendix A. Convection model Fig. 26. X^ ðtÞ vs. R(t) phase diagram comparison for P ¼ 8 and b ¼ 7. The irrigation model predicts large oscillations in the average uninfected cell fractions.

7. Conclusions and future work An analysis and comparison of the spatial cancer models provided in the literature have given insight into the equilibrium aspects of virotherapy treatment. Equilibrium fractions have been obtained for a constant immune suppression drug delivery. A static line is defined for CPA drug delivery to determine its minimum dose to prevent the tumor enlargement after few hours of virotherapy treatment. Higher doses than the ones specified by the static line will make the cancer shrink, while doses below the ones given by the static line make the tumor growth. The models proposed in the literature consider the confined tumor as the place were the immune system cells are generate and consumed. An irrigation model is introduced in this paper to account for the effect of the innate blood stream compartment. This model permits the application of the immune suppression drug as well as the measurement of the immune cells density at the blood stream level. Such aspects provide a realistic view of the

A.1. Static line The Static Line for the convection model indicates that R_ ¼ 0 after some hours of virotherapy treatment. Eq. (19) indicates that such a condition is represented by a quadratic expression in P, where the polynomial coefficients are given by ac ¼ kk0 bdðl þ mÞ

ðA:1Þ

bc ¼ kgðwbd þ k0 slÞðl þ mÞ þ ðwbd þ k0 slÞðbdm þ k0 ðlm þ dðl þ mÞÞÞ

ðA:2Þ cc ¼ ðk0 s þ sgÞlmðwbdk0 slÞ þkswg2 lðl þ mÞ wðwbd þ k0 slÞðbdm þ gðlm þ dðl þ mÞÞÞ

ðA:3Þ

The substitution of the parameters given in Table 1 provides a relation of the form PðbÞ ac ðbÞ ¼ 0:00366b

ðA:4Þ

bc ðbÞ ¼ 0:105ðb0:1754Þðb þ 3:2282Þ

ðA:5Þ

cc ðbÞ ¼ 2:1ðb9:771Þðb0:503Þ

ðA:6Þ

934

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

Note that P becomes a function of b through the calculation of the quadratic polynomial coefficients. Only one of the quadratic roots provide a feasible value for P. A.2. Linearized model A linearized model is developed for the convection model around the state vector x at equilibrium conditions: x ¼ ½X

Y

T

Z

V

ðA:7Þ

Note that only four states are considered because the summation of all tumor cell fractions is one. The linearized PDE right-hand ~ where x~ represents the state vector in side has the form Ac x~ þ Bc P, deviation variables: x~ ¼ ½X~

Y~

Z~

V~ T

ðA:8Þ

while P~ are the variations of P(t) from the equilibrium condition used for linearization. The elements of the system matrix Ac are

elements of the system matrix Ad are ð1,1Þ : F d V b þ lX ðl þ mÞ

ð2,1Þ : V bY ðl þ mÞ

ð1,2Þ : X m

ð2,2Þ : F d kZ dY m

ð1,3Þ : 0

ð2,3Þ : kY

ð1,4Þ : X b

ð2,4Þ : X b

ð3,1Þ : 0

ð4,1Þ : 0

ð3,2Þ : Z s ð3,3Þ : Z w

ð4,2Þ : d ð4,3Þ : k0 V

ð3,4Þ : 0

ð4,4Þ : k0 Z g

where the pair (i,j) indicates the matrix element located at the ith row and jth column. The system matrix Bd is given by BTd ¼ ½0

0

Z

0

ðB:8Þ

Local controllability conditions were verified for the pair ðAd ,Bd Þ along the feasible range of cell equilibrium fractions.

ð1,1Þ : F cV b þ lX ðl þ mÞ ð1,2Þ : X ðsZ þ mÞ

ð2,1Þ : V bY ðl þ mÞ ð2,2Þ : F c kZ sY Z dY m

Appendix C. Irrigation model

ð1,3Þ : X ðsY þ 2wZ m þ PÞ ð1,4Þ : X b

ð2,3Þ : Y ðk þ PsY þ 2wZ mÞ ð2,4Þ : X b

C.1. Static line

ð3,1Þ : Z ðlmÞ ð3,2Þ : Z ðssZ mÞ

ð4,1Þ : 0 ð4,2Þ : d

ð3,3Þ : Z ðPwsY þ2wZ mÞ ð3,4Þ : 0

ð4,3Þ : k0 V ð4,4Þ : k0 Z g

ðA:9Þ

where the pair (i,j) indicates the matrix element located at the ith row and jth column. The system matrix Bc is given by Bc ¼ Z ½X

Y

ðZ 1Þ

0

T

ðA:10Þ

Controllability conditions were verified for the pair ðAc ,Bc Þ along the feasible range of cell equilibrium fractions.

ðB:7Þ

The Static Line for the irrigation model have the same quadratic form than for the convection and diffusion models. The quadratic coefficients are given by the following expressions: ai ¼ kk0 bdðl þ mÞ

ðC:1Þ

bi ¼ kð2k0 wZ n bd þwbgd þ k0 sglÞðl þ mÞ þk0 ðwbd þ k0 slÞðlm þ dðl þ mÞÞ

ðC:2Þ

2

ci ¼ k20 s2 l m þw2 bdðbdm þ gðlm þ dðl þ mÞÞÞ þswlðkg2 ðl þ mÞk0 ð2bdm þ gðlm þ dðl þ mÞÞÞÞ

ðC:3Þ

where the substitution of the parameters in Table 1 provides coefficients that are solely dependent of b: Appendix B. Diffusion model

ai ðbÞ ¼ 0:005143b

ðC:4Þ

B.1. Static line

bi ðbÞ ¼ 0:6102ðb0:137Þ

ðC:5Þ

The Static Line for the diffusion model is given by the solution of the quadratic Eq. (30), where the polynomial coefficients are given by

ci ðbÞ ¼ 2:95ðb7:7014Þðb0:3181Þ

ðC:6Þ

ad ¼ kk0 bdðl þ mÞ

ðB:1Þ

C.2. Linearized model

bd ¼ kgðwbd þk0 slÞðl þ mÞ þk0 ðwbd þ k0 slÞðlm þ dðl þ mÞÞ

ðB:2Þ

The elements of the system matrix Ai are ð1,1Þ : F d V b þ lX ðl þ mÞ

ð2,1Þ : V bY ðl þ mÞ

ðB:3Þ

ð1,2Þ : X m ð1,3Þ : 0

ð2,2Þ : F d kZ dY m ð2,3Þ : kY

The substitution of the parameters given in Table 1 provides a relation of the form PðbÞ

ð1,4Þ : X b

ð2,4Þ : X b

ð3,1Þ : 0

ð4,1Þ : 0

ad ðbÞ ¼ 0:00366b

ðB:4Þ

ð3,2Þ : 0

ð4,2Þ : d

bd ðbÞ ¼ 0:425536ðb0:137Þ

ðB:5Þ

ð3,3Þ :  t1t

ð4,3Þ : k0 V

ð3,4Þ : 0

ð4,4Þ : k0 Z g

cd ðbÞ ¼ 2:1ðb7:44952Þðb0:324292Þ

ðB:6Þ

ð5,1Þ : 0

ð5,2Þ : sZ te

Note the resemblance of this quadratic from with the one obtained for the convection model. In particular, ad ¼ac, and the other coefficients conserve the same order in b.

ð5,3Þ : 0

cd ¼ swlðkg2 ðl þ mÞk0 ð2bdm þ gðlm þ dðl þ mÞÞÞÞ 2

þ w bdðbdm þ gðlm þ dðl þ m

ÞÞÞk20

2 2

s l m

B.2. Linearized model Similar to the convective model, a linearized model is developed for the diffusion model around equilibrium conditions. The

ðC:7Þ

ð5,4Þ : 0 wZ

ð5,5Þ :  te ð2,5Þ : 0

ð1,5Þ : 0 ð3,5Þ : t1t

ð4,5Þ : 0 where the pair (i,j) indicates the matrix element located at the ith row and jth column. Note that five states are required to describe the dynamic response of this model because of the inclusion of

R. Dunia, T.F. Edgar / Computers in Biology and Medicine 41 (2011) 922–935

the external immune cell density Ze. The system matrix Bi is given by BTi ¼ Z ½0

0

0

0

t1e 

ðC:8Þ

Local controllability conditions were verified for the pair ðAi ,Bi Þ along the feasible range of cell equilibrium fractions. Furthermore, with the assumption that the external immune cell density Ze is a measurable quantity, the system matrix Ci is given by Ci ¼ ½0

0

0

0

1

ðC:9Þ

and the pair ðAi ,Ci Þ is observable. This property indicates that the continuous measurement of Ze can be used to estimate the remaining set of states during virotherapy. Such an estimation is subject to the local linearized model accuracy. References [1] S.J. Ries, C.H. Brandts, Oncolytic viruses for the treatment of cancer: current strategies and clinical trials, Drug Discovery Today 9 (2004). [2] M. Vaha-Koskela, J.E. Heikkila, A.E. Hinkkanen, Oncolytic viruses in cancer therapy, Cancer Letters 254 (2007) 178–216. [3] L.M. Wein, J.T. Wu, D.H. Kirn, Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: implications for virus design and delivery, Cancer Research 63 (2003) 1317–1324. [4] F. Khuri, J. Nemunaitis, I. Ganly, J. Arseneau, I. Tannock, L. Romel, M. Gore, J. Ironside, R. MacDougall, C. Heise, B. Randlev, A. Gillenwater, P. Bruso, S. Kaye, W. Hong, D. Kirn, A controlled trial of intratumoral onyx-015, a selectively-replicating adenovirus, in combination with cisplatin and 5-fluorouracil in patients with recurrent head and neck cancer, Nature Medicine 6 (2000). [5] J. Markert, M. Medlock, S. Rabkin, G. Gillespie, T. Todo, W. Hunter, C. Palmer, F. Feigenbaum, C. Tornatore, F. Tufaro, R. Martuza, Conditionally replicating herpes simplex virus mutant, g207 for the treatment of malignant glioma: results of a phase i trial, Gene Therapy 7 (2000). [6] L. Csatary, G. Gosztonyi, J. Szeberenyi, Z. Fabian, V. Liszka, B. Bodey, C. Csatary, Mth-68/h oncolytic viral treatment in human high-grade gliomas, Journal of Neuro-Oncology 67 (2004).

935

[7] A. Friedman, J.P. Tian, G. Fulci, E.A. Chiocca, J. Wang, Glioma virotherapy: effects of innate immune suppression and increase viral replication capacity, Cancer Research 66 (2006) 2314–2319. [8] H. Wakimoto, G. Fulci, E. Tyminski, E.A. Chiocca, Altered expression of antiviral cytokine mrnas associated with cyclophosphamide’s enhancement of viral oncolysis, Gene Therapy 11 (2004) 214–223. [9] H. Vladar, J. Gonzalez, Dynamic response of cancer under the influence of immunological activity and therapy, Journal of Theoretical Biology 227 (2004) 335–348. [10] D. Wodarz, Gene therapy for killing p53-negative cancer cells: use of replicating versus nonreplicating agents, Human Gene Therapy 14 (2003) 153–159. [11] S.C. Ferreira, M.L. Martins, M.J. Vilela, Fighting cancer with viruses, Physica 345 (2004) 591–602. [12] D. Wodarz, Viruses as antitumor weapons: defining conditions for tumor remission, Cancer Research 61 (2001) 3501–3507. [13] Z. Bajzer, T. Carr, K. Josic, S. Russell, D. Dingli, Modeling of cancer virotherapy with recombinant measles viruses, Journal of Theoretical Biology 252 (2008) 109–122. [14] J. Wu, H.M. Byrne, D.H. Kirn, L.M. Wein, Modeling and analysis of a virus that replicates selectively in tumor cells, Bulletin of Mathematical Biology 63 (2001) 731–768. [15] J. Wu, D.H. Kirn, L.M. Wein, Analysis of a three-way race between tumor growth, a replication-competent virus and an immune response, Bulletin of Mathematical Biology 66 (2004) 605–625. [16] A. Friedman, Y. Tao, Analysis of a model of a virus that replicates selectively in tumor cells, Journal of Mathematical Biology 47 (2003) 391–423. [17] J. Wang, J. Tian, Numerical study for a model of tumor virotherapy, Applied Mathematics and Computation 196 (2008) 448–457. [18] Y. Tao, Q. Guo, The competitive dynamics between tumor cells, a replicationcompetent virus and an immune response, Journal of Mathematical Biology 51 (2005) 37–74. [19] D. Dingli, M. Cascino, K. Josic, S. Russell, Z. Bajser, Mathematical modeling of cancer radiovirotherapy, Mathematical Biosciences 199 (2006) 55–78. [20] L.M. Wein, J.E. Cohen, J.T. Wu, Dynamic optimization of a linear-quadratic model with incomplete repair and volume-dependent sensitivity and repopulation, International Journal of Radiation Oncology Biology Physics 47 (2000) 1073–1083. [21] J.A. Florian, J. Eisemanb, R.S. Parker, Nonlinear model predictive control for dosing daily anticancer agents using a novel saturating-rate cell-cycle model, Computers in Biology and Medicine 38 (2008) 339–347.