- Open Access
Mathematical modeling of tumor therapy with oncolytic viruses: effects of parametric heterogeneity on cell dynamics
© Karev et al; licensee BioMed Central Ltd. 2006
- Received: 28 September 2006
- Accepted: 03 October 2006
- Published: 03 October 2006
One of the mechanisms that ensure cancer robustness is tumor heterogeneity, and its effects on tumor cells dynamics have to be taken into account when studying cancer progression. There is no unifying theoretical framework in mathematical modeling of carcinogenesis that would account for parametric heterogeneity.
Here we formulate a modeling approach that naturally takes stock of inherent cancer cell heterogeneity and illustrate it with a model of interaction between a tumor and an oncolytic virus. We show that several phenomena that are absent in homogeneous models, such as cancer recurrence, tumor dormancy, and others, appear in heterogeneous setting. We also demonstrate that, within the applied modeling framework, to overcome the adverse effect of tumor cell heterogeneity on the outcome of cancer treatment, a heterogeneous population of an oncolytic virus must be used. Heterogeneity in parameters of the model, such as tumor cell susceptibility to virus infection and the ability of an oncolytic virus to infect tumor cells, can lead to complex, irregular evolution of the tumor. Thus, quasi-chaotic behavior of the tumor-virus system can be caused not only by random perturbations but also by the heterogeneity of the tumor and the virus.
The modeling approach described here reveals the importance of tumor cell and virus heterogeneity for the outcome of cancer therapy. It should be straightforward to apply these techniques to mathematical modeling of other types of anticancer therapy.
Leonid Hanin (nominated by Arcady Mushegian), Natalia Komarova (nominated by Orly Alter), and David Krakauer.
- Initial Distribution
- Tumor Heterogeneity
- Homogeneous Model
- Uninfected Cell
- Heterogeneous Model
Reviewed by Leonid Hanin (nominated by Arcady Mushegian), Natalia Komarova (nominated by Orly Alter), and David Krakauer.
For the full reviews, please go to the Reviewers' comments section.
Cancers are extremely complex systems that possess many features of robustness, i.e., they are systems that tend to maintain stable functioning despite various perturbations [1, 2]. The main two mechanisms that enable cancer robustness are functional redundancy, that comes from inherent heterogeneity of tumor cells, and feedback-control systems that facilitate survival of a tumor under adverse conditions, e.g., caused by anticancer drugs [3–6]. Many approaches to anticancer treatment have had limited success due to cancer robustness. An important challenge is to identify fragilities of cancers as robust biological systems and develop treatment strategies that take advantage of these weak links. Thus, a better understanding of the mechanisms that yield cancer robustness, particularly, functional redundancy, is vital.
Redundancy in cancers occurs at two levels. First, multiple copies of identical cells increase the likelihood of restoration of a tumor after treatment. Second, functional redundancy can be mediated by functionally equivalent but heterogeneous components (known as heterogeneous redundancy). Heterogeneous redundancy is thought to be central to tumor robustness. Although heterogeneity is a well-recognized characteristic of tumors that can lead to drug resistance, the current body of experimental and clinical data that relate to dynamic changes in intratumoral heterogeneity during progression, as well as to responses to various therapeutic strategies, is insufficient .
The term 'tumor heterogeneity' means the existence of distinct subpopulations of tumor cells with specific characteristics within a single neoplasm [7, 8]; in particular, surviving subpopulations of cells with metastatic potential after anticancer therapy can lead to tumor recurrence [9–11]. Tumor heterogeneity is a well documented phenomenon  as demonstrated by extensive cytogenetic analysis [12, 13]. Furthermore, many tumors are characterized by heterogeneous, often non-random, spatio-temporal distribution of genetically heterogeneous tumor cells . Distinct subpopulations of cells have been isolated from experimental and human neoplasms of every major histological type and location. These subpopulations are heterogeneous for many characteristics, such as morphology, growth rate, metastatic potential, karyotype, antigenicity, immunogenicity, biochemical properties, sensitivity to chemotherapeutic agents and radiation, etc. Genetic heterogeneity is, arguably, a major cause of acquired drug resistance of tumors [14, 15].
The dynamics of host-tumor system, which entails co-evolution under the selective pressure that is imposed by host environments, including antitumor drugs, is highly complex and nonlinear. Thus, to precisely define the conditions for successful therapy, mathematical models are needed. Extensive efforts have been dedicated over many years to mathematical modeling of cancer development and anticancer therapy. Stochastic models that take into account random mutations and cell proliferation proved to be useful in the context of epidemiology and statistical data , and for modeling cancer initiation and progression in terms of somatic evolution . Deterministic models of tumor growth have proved valuable as well. Many of these have addressed avascular and vascular tumor growth taking advantage of methods borrowed from physics  but some use models from population biology to treat a tumor as a dynamic society of interacting cells [19–21]. A variety of mathematical approaches contribute to modeling cancer progression from different standpoints and take stock of various factors affecting tumor growth [ and references therein, ]. Combined analysis of tumor growth and anticancer therapies within the framework of mathematical modeling also produced a number of significant results [24, 25].
Inasmuch as tumor heterogeneity is one of the crucial factors in determining possible outcomes of anticancer treatment, it has to be incorporated and investigated in mathematical models. Although heterogeneity of tumor cell populations has been widely studied, there is no general approach that would provide a comprehensive description of intrinsic heterogeneity of tumors and would be amenable to qualitative and quantitative mathematical analysis. The existing mathematical models of anticancer therapy either do not include the effects of heterogeneity [e.g., [26, 27]], or consider finite, usually, small, number of subpopulations, e.g., sensitive and resistant cells or proliferating and quiescent cells [e.g., [28, 29]]. Other models take into account only spatial heterogeneity and do not address the important subject of genetic heterogeneity [e.g., [30–32]], or use individual-based approach and simulation techniques and thus are, practically, not amenable to theoretical mathematical analysis [e.g., [33–35]].
The main goal of the present paper is to introduce a novel approach to model heterogeneity into the field of cancer modeling. We show that heterogeneous models, although still oversimplified, reflect qualitatively new phenomena, some of which are observed in experiments and clinical trials. We present an effective method to analyze these models and examine the predictions that such heterogeneous models can yield.
We illustrate our approach by addressing a complex process that involves both virus-cell interaction and tumor growth, namely, the interaction of the so-called oncolytic viruses with tumors. Oncolytic viruses are viruses that specifically infect and kill cancer cells but do not affect normal cells [36–39]. Many types of oncolytic viruses have been studied as candidate therapeutic agents, including adenoviruses, herpesviruses, poxviruses, reoviruses, paramyxoviruses, and retroviruses [37, 39]. Probably, the best-characterized oncolytic virus, that has drawn much attention, is ONYX-015, an attenuated adenovirus that selectively infects tumor cells with a defect in the p53 gene [38, 40]. This virus has been shown to possess substantial antitumor activity and has proven relatively effective at reducing or eliminating tumors in clinical trials [41–43]. Although safety and efficacy remain major 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 or are about to enter clinical trials [37, 44–47]. Recently, synergistic use of immune therapy and oncolytic viruses has shown particular promise in cancer treatment .
The oncolytic effect can result from at least three distinct modes of virus-host interaction [37, 39]. The first mode involves repeated cycles of viral replication in the tumor cells leading to cell death and, consequently, to tumor reduction and, potentially, elimination. The second mode involves low-level virus reproduction that, however, results in the production of a cytotoxic protein that causes cell damage. The third mode consists in induction of antitumor immunity by virus infection of cancer cells. Cancer cells possess weak antigens for host immune sensitization. Virus infection causes inflammation and lymphocyte penetration into the tumor, with the virus antigens eliciting increased sensitivity to tumor necrosis factor-mediated killing.
Although the indirect modes of virus cancer therapy based on production of cytotoxic proteins or antitumor immunity might be promising, direct lysis of tumor cells by an oncolytic virus is the current mainstream strategy. Experiments on human tumor xenografts in nude mice have shown that the effect of oncolytic virus infection on tumors can range from no apparent effect, to reduction and stabilization of the tumor load (i.e., the overall size of a tumor), to elimination of the tumor . Complete regression of tumors has been reported also in some patients treated with oncolytic viruses as part of clinical trials . In a previous study , we presented a conceptual mathematical model of tumor cells-virus interaction which, depending on system parameter values, exhibits various behaviors including deterministic elimination of the cancer cells. Here, we further examine this model in conjunction with different possible heterogeneities of tumor cells, such as different susceptibility to infection, different death rates of infected cells, and others.
Homogeneous mathematical models (phase-parameter portrait)
In this section, we briefly present the main results of the analysis of a simple conceptual model of virus-tumor interaction ; these will be important for the further analysis described in the present work. The model allows for two populations of cells: uninfected tumor cells and infected tumor cells and is an extension of the previous work of Wodarz  to include an alternative non-linear functional response for infection rate.
The model, which considers two types of cells growing in the logistic fashion, has the following form:
where X is the size of the uninfected cell population; Y is the size of the infected cell population; r1 and r2 are the maximum per capita growth rates of uninfected and infected cells, respectively; K is the carrying capacity; b is the transmission coefficient (this parameter may also include the replication rate of the virus); and a is the rate of infected cell killing by the virus (cytotoxicity). All the parameters of the model are supposed to be non-negative. Model (1) is subject to initial conditions X(0) = X0 > 0 and Y(0) = Y0 > 0. The concentration of viral particles is not explicitly included; it is assumed that virus abundance is proportional to infected cell abundance .
Rescaling model (1) by letting X* (t*) = X(t)/K, Y* (t*) = Y (t)/K, t* = r1t leads to the system
where β = b/r1, γ = r2/r1, and δ = a/r1. In the following analysis, we suppress the asterisks to simplify the notations, but it should be clear that we use non-dimensional parameters and scaled sizes of cell populations, so that X + Y ≤ 1 for any t.
The model (2) exhibits all possible outcomes of oncolytic virus infection, i.e., no effect on the tumor (domains I and II in Fig. 1), stabilization or reduction of the tumor load (domains IV and V), and complete elimination of the tumor (domain VIII). Moreover there are two domains (domains III and VII) where the final outcome crucially depends on the initial conditions and can result either in failure of virus therapy or in stabilization (domain III) and elimination (domain VII) of the tumor.
The model (1) consists of two coupled, deterministic differential equations allowing for cell reproduction and death, and cell infection. This model is one of the conceptual mathematical models of tumor growth that treat a tumor as a dynamic society of interacting cells (e.g., [19, 20, 22]). Most population models suppose that all individuals have identical attributes, in particular, identical rates of growth, death and birth. This assumption simplifies computation, albeit at the cost of realism.
There are different ways to model heterogeneity. If it is assumed that heterogeneity of population is mediated by a structured variable (such as explicit space or age), we have to deal with partial differential equations (PDEs) (for application of space explicit models to tumor-virus interaction, see  and references therein; cell age is modeled, e.g., in ). By contrast, here, we explore heterogeneity of parameters (such as birth and death rates, and infection rate), i.e., we consider parameters as inherent and invariable properties of individuals, whereas parameter values can vary between individuals; any changes of mean, variance and other characteristics of the parameter distribution with time are caused only by variation of the population structure due to the system dynamics; the distribution of parameter values changes with time because different subpopulations evolve with different, non-constant rates. This type of heterogeneity has been designated 'parametric' by Dushoff .
The most common method to take account of parametric heterogeneity is to divide a population into groups [28, 57, 58]. An important disadvantage of this approach is that heterogeneity within a group cannot be incorporated and the number of groups usually should be small because, otherwise, the dimension of the model is too large for qualitative analysis. Another approach is to consider the population as having a continuous distribution of parameters [59–62]. In the latter case, one usually has to deal with infinite-dimensional dynamical systems but, in a wide range of specific cases, there is an efficient analytical theory that allows one to reduce an initial distributed system to a system of ordinary differential equations (ODEs) [61, 62] (see also Mathematical Appendix (hereinafter MA) [see Additional file 1]). Here, we adopted this approach and emphasize that the applied reduction permits us to follow not only the total sizes of the populations, but also the time-dependent dynamics of the distributions, i.e., no information is lost when we switch from infinite-dimensional dynamical system to system of ODEs. It should be further emphasized that this reduction is not an approximation but results in an equivalent model, in contrast to the approaches commonly used to solve similar systems numerically [59, 63, 64]. Moreover, this approach can be used for arbitrary distributions, both continuous and discrete, with the mathematical formalism remaining the same.
The starting point of our analysis is the model (2) which is a non-dimensional version of the initial model (1). As a matter of fact, by assuming that one or more parameters are distributed through the populations of cells, we have to deal with the initial system (1) with dimensional parameters. However, the two parameters that are of particular interest are transmission coefficient b and cytotoxicity a, and in the following we assume that these parameters (or either of them) are distributed, whereas the net growth rates and carrying capacity are supposed to be the same for all tumor cells. Since non-dimensional parameters β and δ differ from b and a only by linear scaling, the system with distributed non-dimensional parameters β and δ can be analyzed without loss of generality.
Let us assume that transmission coefficient β is distributed through the population of uninfected tumor cells. Denote B the set of possible values of β and x(t, β) the density of uninfected cells that have a given value of β at moment t (i.e., for any subset ⊆ B, the part of the population with trait values belonging to is given by , and the total size of the population, X (t), is given by ).
The number of uninfected cells with a given value of β that become infected per time unit is given by ρ x(t, β), where ρ = ρ(β, X, Y) is the rate of infection. We assume that the rate of infection depends only on the value of β and the total sizes of the infected and uninfected cell subpopulations. This assumption reflects the fact that different cells can have different susceptibilities to virus infection, i.e., some of them can be infected with a relatively high probability whereas others are infected with a low or even arbitrary, close to zero probability; the latter case corresponds to the situation when a particular subpopulation of tumor cells can hardly be infected by the virus.
Using the nonlinear transmission function in (2) we obtain that the change in subpopulation of uninfected cells with parameter value β is given by β x(t, β)Y(t)/(X(t) + Y(t)). We also need to specify the law of growth of tumor cells with a fixed value of β. Let us assume that uninfected cells with a particular value of β beget daughter cells with the same parameter value. For the net growth rate, we use the logistic law that depends on the total population sizes, X and Y. We do not consider heterogeneity in the infected cell population, and, hence, the total amount of newly infected cells per time unit is given by Y(t) /X(t) + Y(t)), which should be taken into account when writing down the equation for the change in the infected cell population. By specifying initial conditions for x(t, β) and Y(t), we obtain the model in the form
where the following notation was used: E β (t) = ()/X(t).
The initial conditions are
x(0, β) = x0 (β) = X0p(0, β), Y(0) = Y0. (4)
Here x0(β) is the initial distribution of β in the population of uninfected cells, X0 is the total number of uninfected cells at the initial moment, and p(0, β) is the probability density function (pdf) of the initial distribution of. β. Note that E β (t) can be rewritten in the form E β (t) = , where p(t, β) is the pdf of the distribution of β at time t, and, thus, E β (t) is the mean value of β at moment t.
Integrating the first equation in (3) over β we obtain the system
where we suppressed the dependence of phase variables on t, and the initial conditions are X(0) = X0, Y(0) = Y0. In order to solve this system, we only need to know the explicit expression for E β (t). If we use the usual notation
for the moment generation function (mgf) of the initial pdf, it can be shown ([61, 62], see also MA) that the current mean value of β can be calculated using the mgf of the initial distribution of β ; the exact expression is
where q(t) is an auxiliary variable that satisfies the equation
and M β (λ) is given as far as we specify the initial pdf p(0, β).
In other words, if we want to keep track of the total population sizes in the distributed model (3)–(4), then this model implies the system of ODEs (5)–(6) that can be analyzed in the usual way; moreover, it is very easy to solve this system numerically (there are methods to obtain the solution of (3)–(4), e.g., , but these methods are usually computationally intensive and not so effective as numerical methods for ODEs). It is worth noting that the current population density x(t, β) also can be computed (see Theorem 1 in MA).
Note that system (5) differs from (2) in that the fixed value of parameter β is replaced by its current mean value E β (t) in the distributed model. This mean value, obviously, depends on the system dynamics and the population sizes (see (6) and above), and its calculation at any time moment gives us solution of problem (5).
It can be shown that, for model (3)–(4), the mean of the parameter distribution decreases monotonically with time (dE β (t)/dt ≤ 0 for any t, see Corollary 2 in MA) from the initial value E β (0) to the final value η, which shows the direction of selection: the cells with lower parameter values are selected for. More importantly, however, the simultaneous knowledge of the time dependence of E β (t) and the parametric portrait of the homogeneous model (2) (Fig. 1) allows one to identify transient behavior of the solutions of (5) (and, accordingly, of model (3)), i.e., the behavior of solutions prior to the time moment when E β (t) becomes constant. This behavior can be of particular interest because this is one of the features that distinguishes model (3) from (2) and allows us to explore in detail the temporal dynamics of inhomogeneous biological systems.
Due to the non-negativity of β, we should consider only distributions with support on half-axis β ≥ 0. To illustrate the possible dynamical behavior of the cell populations, we need the initial distribution of β, which is unknown in practice. In the following, we assume that the initial pdf of β is the gamma distribution on [η, ∞) with positive parameters k, s and η ≥ 0. The choice of the gamma distribution can be justified by the fact that it can well approximate almost any unimodal distribution concentrated on the positive half-line. Other possible initial distributions include uniform, log-normal, Pareto, and many others; any of them can be incorporated in the model as far as there exist mgfs for the given q(t) (we note, however, that, in some particular cases, mgf does not exist for some values of q(t), see MA).
For the case of gamma distributed β at the initial time moment, it can be proved (see MA, Example 1) that β is gamma-distributed at any time moment with E β (t) = η + k /(s-q(t)). Typically, when Y(t) > const > 0, then q(t) → -∞ for t → ∞ (see (6)) and hence E β (t) → η. Thus, the knowledge of the support of the initial distribution and the direction of selection allows us, in a number of cases, to predict the final state of the system.
In brief, we use the following approach to perform and present numerical simulations. First, we need to specify the initial conditions and the initial distribution of β, i.e., assign values of η, k, s. Here, for the sake of transparency, we have chosen to set the values of the mean E β = and variance = - E β (0)2 of the initial distribution, as well as the bounds of the distribution if they are different from zero or infinity (for the gamma distribution, the bound that has to be specified is its left boundary, η). In the analyzed case (when we work with the distribution that has two free parameters, k and s), this is sufficient to unambiguously determine the values of the parameters (obviously, however, this is not the case in the general situation).
Loosely speaking, during the simulation, the parametric point with coordinates (γ, δ, E β (t)) travels in the parameter space of model (2); the movement occurs along the line connecting points (γ, δ, E β (0)) and (γ, δ, η) in the parameter portrait of (2) (see Fig. 1 for the cut of the parameter space for fixed γ). The speed of movement is determined, together with the initial conditions X0 and Y0, by the parameters of the initial distribution, or, equivalently, by the initial mean and variance of the gamma distribution with the given left boundary.
Fig. 2 shows the phenomenon of tumor recurrence after a relatively long period of small tumor load, a phenomenon that is not seen in the original homogeneous model (2). Clearly, reappearance of the tumor is a probable outcome of tumor-specific virus treatment if there are cancer cells that are inaccessible to the virus (or are accessible at an extremely slow rate compared to the characteristic time of virus propagation): after the virus kills off the susceptible tumor cells and is cleared (given that it cannot infect healthy tissues), these resistant cancer cells can develop into a new tumor. It should be emphasized that, in the simulations in Fig. 2, the final outcome is tumor recurrence, i.e., the tumor dynamics is strongly affected by cell heterogeneity although all cells have a positive probability to be infected (η > 0). This example shows that heterogeneous model (3), together with dynamical regimes inherited from model (2), possesses new regimes, and, thus, is more general.
Distributed susceptibility and distributed cytotoxicity
Together with differential susceptibility considered in the previous section, parameter δ (the rate of killing of infected cells) also may be assumed to take different values in the infected cell population. The population of infected cells can be heterogeneous on its own such that some cells are more likely to die faster than others when infected by the virus. From different perspective, we can attribute this heterogeneity to the virus, i.e., the virus population is assumed to consist of strains that differ in their ability to kill cells. In previous work [51, 52], it has been shown that, in some parameter domains, there is an optimal level of cytotoxicity that is necessary to maximally reduce the tumor size if other parameter values are fixed. In addition, it is sometimes desirable to kill as many cells as possible during a short period of time; such situation would favor using different viruses for therapeutic purposes. Moreover, the latter assumption on virus heterogeneity can help one to interpret the model equations which, assuming that δ takes values from set Δ, and y(t, δ) is the density of infected cells having parameter value δ, take the form:
where E β (t) is defined as above, and the initial conditions are
x(0, β) = x0(β) = X0 p1 (0, β), y(0, δ) = y0(δ) = Y0 p2 (0, δ). (8)
Note that in (7) both cell populations are closed under reproduction; that is, e.g., the uninfected cell that has parameter value β* can produce daughter cells only with the same parameter value (the same holds for the infected cell population), moreover, if an uninfected cell was infected by the virus from an infected cell with a parameter value δ*, it falls into the same class. If we assume that δ is an attribute of the cell population, this assumption is difficult to justify, whereas if δ is an attribute of the virus, the assumption follows from obvious considerations.
System (7)–(8) implies the following system of ODEs (see MA, Corollary 1):
where the mean parameter values are
for the given mgf of p1 (0, β) and p2(0, δ), and the auxiliary variables can be found from
To analyze system (9)–(10), one has to specify initial conditions and initial distributions for parameters β and δ. As before, to obtain some insight into the possible dynamical behavior of (9), we can consider the parameter portrait (Fig. 1) and time-dependent changes of the mean parameter values. In the case under consideration E β (t), as before, moves from top to bottom in Fig. 1, whereas E δ (t) goes from right to left. Note that the speed of movement in the δ -direction is determined only by the properties of the mgf, because q2 (t) = -t from (10) and q2(t) does not depend on the population sizes.
In Fig. 6, the initial and final parameter values belong to the most favorable domain VIII (Fig. 1), in which complete eradication of the tumor cells occurs. On its way to extinction, however, the tumor load behaves in an irregular, complex way. This example shows the possibility that complex and erratic observational data can be explained not only by random effects and noise but also by the innate heterogeneity of the cell and virus populations.
In general, both distributed susceptibility of uninfected cells and distributed cytotoxicity of the virus do not favor the tumor treatment and decrease chances for effective cure by increasing tumor robustness; one of the possible ways to overcome this problem (within the framework of the considered models) is suggested in the next section.
Distributed susceptibility and distributed virulence
The inherent heterogeneity of tumor cells implies that cells differ in their susceptibility to virus infection, and this situation was modeled in two previous sections. In general, it can be assumed that the virus population is also heterogeneous in the sense that it contains viruses with different ability to infect tumor cells, which we denote virulence. Inasmuch as model (1) does not explicitly model the virus dynamics, we can incorporate the effect of virus heterogeneity into the infected cell population. Assuming the susceptibility of uninfected cells β1 and the virus virulence β2 to be independent, we can consider the transmission coefficient β as being proportional to the product of β1 and β2 (moreover, the proportionality constant can be incorporated in one of the parameters, thus we have β = β1β2).
As before, denote B1 the set of possible values of β1 and x(t, β1) the density of uninfected cells having the given value of β1 at the moment t; similarly, denote B2 the set of possible values of β2 and y(t, β2) the density of infected cells with the given value of β2 at the time moment t. The number of cells with susceptibility β1 newly infected by the virus produced from previously infected cells with virulence β2 is given by β1β2x(t, β1) y (t, β2)/(X(t) + Y(t)). The total number of cells that leave the uninfected population and have susceptibility β1 is β1x(t, β1) /(X(t) + Y(t)); and the total change in the infected cell population with parameter value β2 due to the infection process is given by β2 y(t, β2) /(X(t) + Y(t)). Combining the above assumptions, we obtain the model
where the notations (t) = ()/X(t), (t) = ()/Y(t) were used for the current mean values of the corresponding parameters, and the initial conditions are
x(0, β1) = x0(β1) = X0 p1(0, β1), y(0, β2) = y0(β2) = Y0 p2(0, β2), (12)
The model (11)–(12) implies the system of ODEs (Corollary 1, MA)
where the mean parameter value is E(t) = (t) (t), and
are expressed with the mgfs of p1 (0, β1), p2 (0, β2). The auxiliary variables can be found from the equations
As before, we have to specify the initial pdfs for β1 and β2 Now we have dq2 (t)/dt ≥ 0, and it can be shown that d (t)/dt ≥ 0 (Corollary 2, MA), that is, formally, the movement in the β2-direction goes from smaller to greater mean values. In order not to deal with an infinitely large coefficient, it is reasonable to assume that the initial distribution of β2 is determined on a closed interval [c1, c2]. In this case, we can choose as the initial distribution the beta distribution on [c1, c2].
In Fig. 7, the results of numerical simulation of system (13)–(14) are shown together with functions (t), (t) and E(t) = (t) (t). We start in domain VIII (elliptic sector in Fig. 1), and the asymptotic state is in domain VII where two opposite outcomes are possible, namely, complete tumor eradication or logistic tumor growth. The initial conditions and parameter values are the same for both cases; the two cases differ only in the initial variance of the β2 distribution. It can be seen from Fig. 7 that even a small difference in the variance of β2 may yield dramatically different results: in the first case, the tumor is cured, whereas, in the second case, virus therapy fails. This example emphasizes that, to predict the outcome of oncolytic virus therapy, it is necessary to know not only the initial sizes of the cell populations but also the degree of heterogeneity of all the parameters under consideration.
We presented a general framework within which to construct and analyze mathematical models of anticancer treatment, with a special emphasis on tumor heterogeneity. Conceptually, this approach is connected to the theoretical considerations of Kitano on cancer robustness [1, 2]. At this initial step of analysis, it appears most important that the mathematical models reproduce important behaviors of tumors at the qualitative level. Using the previously developed model of tumor cell-oncolytic virus interaction , we show that qualitatively distinct behaviors that are absent in the homogeneous model appear as natural consequences of tumor and oncolytic virus heterogeneity. For example, our model accounts for cancer recurrence, and the time until reappearance of the tumor crucially depends on the heterogeneity of the transmission coefficient (Fig. 2). It should be emphasized that the effects of the tumor cell heterogeneity are not limited to trivial resistance of a subpopulation of cells to the virus, i.e., recurrence might ensue even when all cells have a non-zero probability to be infected. Another phenomenon missing in the homogeneous model is constant tumor size during a particular period of time (Fig. 4). Heterogeneity in parameters can lead to complex, irregular evolution of the tumor (Fig. 5 and 6). Thus, interpretation of the results of anticancer therapy should take into account the possibility that irregular, quasi-chaotic behavior can be caused not only by random fluctuations but also by the heterogeneity of the tumor and the virus (Fig. 6). All the simulations presented here reveal the effect of the level of heterogeneity on tumor dynamics. The most obvious case is shown in Fig. 7 where different initial variances of parameter distribution lead either to tumor eradication or to logistic tumor growth, i.e., failure of virus therapy.
Analysis of the models proves that tumor heterogeneity increases cancer robustness, in agreement with the theoretical considerations of Kitano . The results presented here further show that, to counter this adverse effect of tumor heterogeneity, it should be possible to employ a heterogeneous population of an oncolytic virus (see model (11) and Figure 7). The interaction of the two non-homogeneous populations, the tumor and the oncolytic virus, may result in complete elimination of the tumor.
We applied a previously developed, general mathematical technique [61, 62] to investigate non-homogeneous models of complex tumor cells-virus interaction. The advantages of this modeling approach are as follows. This approach can account for different types of parametric heterogeneity of the analyzed populations. To infer the consequences of heterogeneity, we use well-known mathematical tools, such as bifurcation analysis, which identifies points of qualitative change in the system dynamics. The theory of heterogeneous populations [61, 62, 65] allows one to reduce the models to systems of ODEs that, in many cases, can be explored analytically or, if this is not possible, can be solved numerically with high precision.
Cancer is an evolving system, and the main evolutionary forces are selection, mutation and random drift. It should be noticed that our approach explicitly examines only selection. The techniques is applied to deterministic systems, and it is our belief that such system are important for modeling purposes, although analysis of extinction phenomena, such as tumor eradication, may require stochastic factors to be considered. The principal obstacle to the validation of these models against empirical data and their use for prediction purposes is the necessity to know the initial distribution of the model parameters. It should be emphasized that knowledge of mean, variance, and any finite number of the moments of the distribution is not sufficient to describe the evolution of the system over indefinite time . The entire distributions have to be known. However, there is a way to bypass this problem. First, knowledge of several first moments is sufficient to estimate tumor evolution over short time spans. Second, to model the evolutionary process, we need to know how the mean parameter values behave with time. Using the theory of heterogeneous populations, we can infer generic properties of these functions and identify them from empirical data .
To conclude, the approach developed here allows one to take stock of such a complex aspect of cancer as tumor heterogeneity and apply effective analytical techniques to the analysis of heterogeneous models of tumor evolution. Although, in this study, we analyzed the specific case of a tumor interaction with an oncolytic virus, it is our hope that these techniques will prove useful in other systems that include interaction between tumor cells and anticancer agents.
Reviewer's report 1
Leonid Hanin, Department of Mathematics, Idaho State University (nominated by Arcady Mushegian)
This work was supported, in part, by the Intramural Research Program of the National Library of Medicine at National Institutes of Health/DHHS.
- Kitano H: Cancer robustness: tumour tactics. Nature 2003,426(6963):125. 10.1038/426125aPubMedView ArticleGoogle Scholar
- Kitano H: Cancer as a robust system: implications for anticancer therapy. Nat Rev Cancer 2004,4(3):227-235. 10.1038/nrc1300PubMedView ArticleGoogle Scholar
- Alexandrova R: Tumour heterogeneity. Exper Path and Paras 2001,4(6):57-67.Google Scholar
- Gonzalez-Garcia I, Sole RV, Costa J: Metapopulation dynamics and spatial heterogeneity in cancer. Proc Natl Acad Sci U S A 2002,99(20):13085-13089. 10.1073/pnas.202139299PubMedPubMed CentralView ArticleGoogle Scholar
- Maley CC, Galipeau PC, Finley JC, Wongsurawat VJ, Li X, Sanchez CA, Paulson TG, Blount PL, Risques RA, Rabinovitch PS, Reid BJ: Genetic clonal diversity predicts progression to esophageal adenocarcinoma. Nat Genet 2006,38(4):468-473. 10.1038/ng1768PubMedView ArticleGoogle Scholar
- Maley CC, Reid BJ: Natural selection in neoplastic progression of Barrett's esophagus. Semin Cancer Biol 2005,15(6):474-483. 10.1016/j.semcancer.2005.06.004PubMedView ArticleGoogle Scholar
- Heppner GH: Tumor heterogeneity. Cancer Res 1984,44(6):2259-2265.PubMedGoogle Scholar
- Heppner GH, Miller FR: The cellular basis of tumor progression. Int Rev Cytol 1998, 177: 1-56.PubMedView ArticleGoogle Scholar
- Harris JF, Chambers AF, Ling V, Hill RP: Dynamic heterogeneity: characterization of two cell lines derived from experimental lung metastases of mouse KHT fibrosarcoma. Invasion Metastasis 1987,7(4):217-229.PubMedGoogle Scholar
- Hill RP, Chambers AF, Ling V, Harris JF: Dynamic heterogeneity: rapid generation of metastatic variants in mouse B16 melanoma cells. Science 1984,224(4652):998-1001. 10.1126/science.6719130PubMedView ArticleGoogle Scholar
- Nowell PC: Mechanisms of tumor progression. Cancer Res 1986,46(5):2203-2207.PubMedGoogle Scholar
- Gorunova L, Dawiskiba S, Andren-Sandberg A, Hoglund M, Johansson B: Extensive cytogenetic heterogeneity in a benign retroperitoneal schwannoma. Cancer Genet Cytogenet 2001,127(2):148-154. 10.1016/S0165-4608(00)00440-4PubMedView ArticleGoogle Scholar
- Fujii H, Yoshida M, Gong ZX, Matsumoto T, Hamano Y, Fukunaga M, Hruban RH, Gabrielson E, Shirai T: Frequent genetic heterogeneity in the clonal evolution of gynecological carcinosarcoma and its influence on phenotypic diversity. Cancer Res 2000,60(1):114-120.PubMedGoogle Scholar
- Goldie JH, Coldman AJ: A mathematic model for relating the drug sensitivity of tumors to their spontaneous mutation rate. Cancer Treat Rep 1979,63(11-12):1727-1733.PubMedGoogle Scholar
- Skipper HE, Schabel FMJ: Tumor stem cell heterogeneity: implications with respect to classification of cancers by chemotherapeutic effect. Cancer Treat Rep 1984,68(1):43-61.PubMedGoogle Scholar
- Moolgavkar SH, Knudson AGJ: Mutation and cancer: a model for human carcinogenesis. J Natl Cancer Inst 1981,66(6):1037-1052.PubMedGoogle Scholar
- Michor F, Iwasa Y, Nowak MA: Dynamics of cancer progression. Nat Rev Cancer 2004,4(3):197-205. 10.1038/nrc1295PubMedView ArticleGoogle Scholar
- Preziosi L: Cancer Modeling and Simulation. Boca Raton, FL , CRC; 2003.View ArticleGoogle Scholar
- Gatenby RA: Models of Tumor-Host Interaction as Competing Populations: Implications for Tumor Biology and Treatment. Journal of Theoretical Biology 1995,176(4):447-455. 10.1006/jtbi.1995.0212PubMedView ArticleGoogle Scholar
- Gatenby RA, Vincent TL: An evolutionary model of carcinogenesis. Cancer Res 2003,63(19):6212-6220.PubMedGoogle Scholar
- Gatenby RA, Vincent TL: Application of quantitative models from population biology and evolutionary game theory to tumor therapeutic strategies. Mol Cancer Ther 2003,2(9):919-927.PubMedGoogle Scholar
- Wodarz D, Komarova N: Computational biology of cancer. World Scientific; 2005.View ArticleGoogle Scholar
- Komarova NL: Mathematical modeling of tumorigenesis: mission possible. Curr Opin Oncol 2005,17(1):39-43. 10.1097/01.cco.0000143681.37692.32PubMedView ArticleGoogle Scholar
- Goldie JH, Coldman AJ: Drug Resistance in Cancer: Mechanisms and Models. Cambridge University Press; 1998.View ArticleGoogle Scholar
- Martin R, Teo KL: Optimal control of drug administration in cancer chemotherapy. World Scientific Publishing; 1993.Google Scholar
- de Pillis LG, Gu W, Radunskaya AE: Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. Journal of Theoretical Biology 2006,238(4):841-862. 10.1016/j.jtbi.2005.06.037PubMedView ArticleGoogle Scholar
- Liu W, Freedman HI: A mathematical model of vascular tumor treatment by chemotherapy. Mathematical and Computer Modelling 2005,42(9-10):1089-1112. 10.1016/j.mcm.2004.09.008View ArticleGoogle Scholar
- Birkhead BG, Rankin EM, Gallivan S, Dones L, Rubens RD: A mathematical model of the development of drug resistant to cancer chemotherapy. European Journal of Cancer and Clinical Oncology 1987,23(9):1421-1427. 10.1016/0277-5379(87)90133-7PubMedView ArticleGoogle Scholar
- Lakmeche A, Arino O: Nonlinear mathematical model of pulsed-therapy of heterogeneous tumors. Nonlinear Analysis: Real World Applications 2001,2(4):455-465. 10.1016/S1468-1218(01)00003-7View ArticleGoogle Scholar
- Norris ES, King JR, Byrne HM: Modelling the response of spatially structured tumours to chemotherapy: Drug kinetics. Mathematical and Computer Modelling In Press, Corrected Proof:Google Scholar
- Stephanou A, McDougall SR, Anderson ARA, Chaplain MAJ: Mathematical modelling of flow in 2D and 3D vascular networks: Applications to anti-angiogenic and chemotherapeutic drug strategies. Mathematical and Computer Modelling 2005,41(10):1137-1156. 10.1016/j.mcm.2005.05.008View ArticleGoogle Scholar
- Chaplain MAJ: Avascular growth, angiogenesis and vascular growth in solid tumours: The mathematical modelling of the stages of tumour development. Mathematical and Computer Modelling 1996,23(6):47-87. 10.1016/0895-7177(96)00019-2View ArticleGoogle Scholar
- Lankelma J, Fernandez Luque R, Dekker H, Pinedo HM: Simulation model of doxorubicin activity in islets of human breast cancer cells. Biochim Biophys Acta 2003,1622(3):169-178.PubMedView ArticleGoogle Scholar
- Drasdo D, Hohme S: Individual-based approaches to birth and death in avascu1ar tumors. Mathematical and Computer Modelling 2003,37(11):1163-1175. 10.1016/S0895-7177(03)00128-6View ArticleGoogle Scholar
- Quaranta V, Weaver AM, Cummings PT, Anderson ARA: Mathematical modeling of cancer: The future of prognosis and treatment. Clinica Chimica Acta 2005,357(2):173-179. 10.1016/j.cccn.2005.03.023View ArticleGoogle Scholar
- Kirn DH, McCormick F: Replicating viruses as selective cancer therapeutics. Mol Med Today 1996,2(12):519-527. 10.1016/S1357-4310(97)81456-6PubMedView ArticleGoogle Scholar
- Parato KA, Senger D, Forsyth PA, Bell JC: Recent progress in the battle between oncolytic viruses and tumours. Nat Rev Cancer 2005,5(12):965-976. 10.1038/nrc1750PubMedView ArticleGoogle Scholar
- McCormick F: Cancer-specific viruses and the development of ONYX-015. Cancer Biol Ther 2003,2(4 Suppl 1):S157-60.PubMedGoogle Scholar
- Kasuya H, Takeda S, Nomoto S, Nakao A: The potential of oncolytic virus therapy for pancreatic cancer. Cancer Gene Ther 2005,12(9):725-736. 10.1038/sj.cgt.7700830PubMedView ArticleGoogle Scholar
- Kaplan JM: Adenovirus-based cancer gene therapy. Curr Gene Ther 2005,5(6):595-605. 10.2174/156652305774964677PubMedView ArticleGoogle Scholar
- Kirn D, Hermiston T, McCormick F: ONYX-015: clinical data are encouraging. Nat Med 1998,4(12):1341-1342. 10.1038/3902PubMedView ArticleGoogle Scholar
- Khuri FR, Nemunaitis J, Ganly I, Arseneau J, Tannock IF, Romel L, Gore M, Ironside J, MacDougall RH, Heise C, Randlev B, Gillenwater AM, Bruso P, Kaye SB, Hong WK, Kirn DH: 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. Nat Med 2000,6(8):879-885. 10.1038/78638PubMedView ArticleGoogle Scholar
- Nemunaitis J, Khuri F, Ganly I, Arseneau J, Posner M, Vokes E, Kuhn J, McCarty T, Landers S, Blackburn A, Romel L, Randlev B, Kaye S, Kirn D: Phase II trial of intratumoral administration of ONYX-015, a replication-selective adenovirus, in patients with refractory head and neck cancer. J Clin Oncol 2001,19(2):289-298.PubMedGoogle Scholar
- Shah AC, Benos D, Gillespie GY, Markert JM: Oncolytic viruses: clinical applications as vectors for the treatment of malignant gliomas. J Neurooncol 2003,65(3):203-226. 10.1023/B:NEON.0000003651.97832.6cPubMedView ArticleGoogle Scholar
- Kaufman HL, Deraffele G, Mitcham J, Moroziewicz D, Cohen SM, Hurst-Wicker KS, Cheung K, Lee DS, Divito J, Voulo M, Donovan J, Dolan K, Manson K, Panicali D, Wang E, Horig H, Marincola FM: Targeting the local tumor microenvironment with vaccinia virus expressing B7.1 for the treatment of melanoma. J Clin Invest 2005,115(7):1903-1912. 10.1172/JCI24624PubMedPubMed CentralView ArticleGoogle Scholar
- Reid T, Warren R, Kirn D: Intravascular adenoviral agents in cancer patients: lessons from clinical trials. Cancer Gene Ther 2002,9(12):979-986. 10.1038/sj.cgt.7700539PubMedView ArticleGoogle Scholar
- Shen Y, Nemunaitis J: Herpes simplex virus 1 (HSV-1) for cancer treatment. Cancer Gene Ther 2006.Google Scholar
- Thorne SH, Negrin RS, Contag CH: Synergistic antitumor effects of immune cell-viral biotherapy. Science 2006,311(5768):1780-1784. 10.1126/science.1121411PubMedView ArticleGoogle Scholar
- Harrison D, Sauthoff H, Heitner S, Jagirdar J, Rom WN, Hay JG: Wild-type adenovirus decreases tumor xenograft growth, but despite viral persistence complete tumor responses are rarely achieved--deletion of the viral E1b-19-kD gene increases the viral oncolytic effect. Hum Gene Ther 2001,12(10):1323-1332. 10.1089/104303401750270977PubMedView ArticleGoogle Scholar
- Lorence RM, Pecora AL, Major PP, Hotte SJ, Laurie SA, Roberts MS, Groene WS, Bamat MK: Overview of phase I studies of intravenous administration of PV701, an oncolytic virus. Curr Opin Mol Ther 2003,5(6):618-624.PubMedGoogle Scholar
- Novozhilov AS, Berezovskaya FS, Koonin EV, Karev GP: Mathematical modeling of tumor therapy with oncolytic viruses: Regimes with complete tumor elimination within the framework of deterministic models. Biology Direct 2006.Google Scholar
- Wodarz D: Viruses as antitumor weapons: defining conditions for tumor remission. Cancer Res 2001,61(8):3501-3507.PubMedGoogle Scholar
- Nowak MA, May RM: Virus Dynamics: mathematical principles of immunology and virology. New-York , Oxford; 2000.Google Scholar
- Tao Y, Guo Q: The competitive dynamics between tumor cells, a replication-competent virus and an immune response. J Math Biol 2005,51(1):37-74. 10.1007/s00285-004-0310-6PubMedView ArticleGoogle Scholar
- Dyson J, Villella-Bressan R, Webb GF: Asynchronous exponential growth in an age structured population of proliferating and quiescent cells. Mathematical Biosciences 2002, 177-178: 73-83. 10.1016/S0025-5564(01)00097-9PubMedView ArticleGoogle Scholar
- Dushoff J: Host Heterogeneity and Disease Endemicity: A Moment-Based Approach. Theoretical Population Biology 1999,56(3):325-335. 10.1006/tpbi.1999.1428PubMedView ArticleGoogle Scholar
- Hsu Schmitz SF: Effects of genetic heterogeneity on HIV transmission in homosexual populations. In Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods and Theory. Edited by: Castillo-Chavez C, Blower S, van den Driessche P, Kirschner D, Yakubu AA. Springer-Verlag; 2002:245-260.View ArticleGoogle Scholar
- Dushoff J, Levin S: The effects of population heterogeneity on disease invasion. Mathematical Biosciences 1995,128(1-2):25-40. 10.1016/0025-5564(94)00065-8PubMedView ArticleGoogle Scholar
- Ackleh AS: Estimation of rate distributions in generalized Kolmogorov community models. Nonlinear Analysis 1998,33(7):729-745. 10.1016/S0362-546X(97)00665-2View ArticleGoogle Scholar
- Diekmann O, Heesterbeek H, Metz H: The Legacy of Kermack and McKendrick. In Epidemic Models: Their Structure and Relation to Data Edited by: Mollison D. 1995, 95-115.Google Scholar
- Karev G: The effects of non-homogeneity in population dynamics. Doklady Matematica 2000, 62: 141-145.Google Scholar
- Karev G: Dynamics of heterogeneous populations and communities and evolution of distributions. Dis and Cont Dyn Sys 2005, Suppl: 487-496.Google Scholar
- Veliov VM: Newton's method for problems of optimal control of heterogeneous systems. Optimization Methods & Software 2003,18(6):689-703. 10.1080/10556780310001639753View ArticleGoogle Scholar
- Ackleh AS, Marshall DF, Heatherly HE, Fitzpatrick BG: Survival of the fittest in a generalized logistic model. Mathematical Models & Methods in Applied Sciences 1999,9(9):1379-1391. 10.1142/S0218202599000610View ArticleGoogle Scholar
- Novozhilov AS: Analysis of a generalized population predator-prey model with a parameter distributed normally over the individuals in the predator population. Journal of Computer and Systems Sciences International 2004,43(3):378-382.Google Scholar
- Veliov VM: On the effect of population heterogeneity on dynamics of epidemic diseases. J of Math Bio 2005, 51: 123-143. 10.1007/s00285-004-0288-0View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.