Abstract
One of the most promising strategies to treat cancer is attacking it with viruses. Oncolytic viruses can kill tumor cells specifically or induce anticancer immune response. A multiscale model for virotherapy of cancer is investigated through simulations. It was found that, for intratumoral virus administration, a solid tumor can be completely eradicated or keep growing after a transient remission. Furthermore, the model reveals undamped oscillatory dynamics of tumor cells and virus populations, which demands new in vivo and in vitro quantitative experiments aiming to detect this oscillatory response. The conditions for which each one of the different tumor responses dominates, as well as the occurrence probabilities for the other nondominant therapeutic outcomes, were determined. From a clinical point of view, our findings indicate that a successful, single agent virotherapy requires a strong inhibition of the host immune response and the use of potent virus species with a high intratumoral mobility. Moreover, due to the discrete and stochastic nature of cells and their responses, an optimal range for viral cytotoxicity is predicted because the virotherapy fails if the oncolytic virus demands either a too short or a very large time to kill the tumor cell. This result suggests that the search for viruses able to destroy tumor cells very fast does not necessarily lead to a more effective control of tumor growth. [Cancer Res 2009;69(3):1205–11]
Introduction
One of the most promising strategies to treat cancer is attacking it with viruses. Oncolytic viruses are biological nanomachines that can be programmed to specifically target, replicate in, and ultimately kill cancer cells. The new viruses released from the lysed cells can subsequently infect adjacent or distant tumor cells and trigger multiple cycles of infection. All the recombinant DNA technology can be readily used to enhance the selectivity of oncolytic viruses for tumor cells, thereby limiting the infection of normal cells and minimizing therapeutic side effects (1). Several approaches to selectivity have been explored; for example, specifically retargeting viruses to cancer cells (2); linking virus replication to cancer-specific transcription factors such as telomerase (3) or prostate-specific antigen (4); and deleting viral genes for advancing cell cycle and/or preventing apoptosis, typically unnecessary for infecting cancer cells (1). An infected tumor cell is eventually killed by the oncolytic virus through several routes such as lysis, apoptosis, immune-modulated cell killing, or other forms of programmed cell death (5).
Over the past decade, clinical trials have been conducted with at least 14 distinct viruses from six different viral strains (6). Although the arsenal of oncolytic viruses expands continuously, to date only two viruses (Onyx-015 and H101) have entered randomized phase III clinical testing (1, 6). The data from these clinical trials have established the safety, effectiveness, and lack of cross-resistance to currently available treatments of virotherapy. Hence, in the near future, this approach has the potential to supplant many of the current standard anticancer therapies, mainly if combined with existing treatments. However, several hurdles such as the host immune response, systemic distribution, and intratumoral spread must be overcome. Because our clinical experience is limited, mathematical models may provide relevant insights about the nonlinear dynamics ruling the interplay between normal and cancer cells and their viral parasites. A quantitative understanding of this complex network of interactions will certainly accelerate the further development necessary to turn virotherapy into a fully functional and effective therapeutics.
In this article, a virotherapy model recently proposed (7) is investigated. It considers normal cells, uninfected and infected tumor cells, necrotic cells, and free viruses. At the tissue level, reaction-diffusion equations describe the dynamics of nutrients and viruses, both small particles in comparison with cells. In these reaction-diffusion equations, cells act as sinks and sources of nutrients and viruses depending on their internal states. In turn, cell dynamics is a stochastic process controlled by the local concentration of nutrients and free viruses determined at the tissue level. Thus, the model involves partial differential equations and probabilistic cellular automata rules to describe the multiscale dynamics of tumor growth (8, 9). It is worth mentioning that this model was designed for the avascular stages of cancer growth. However, the effects of therapies even at these stages are crucial for the oncologic practice because usually clinically detected tumors have often already generated micrometastases that are growing avascularly in other tissues of the organism. Thus, for instance, of the hundreds of thousand of people who develop colorectal cancer every year, as many as 50% of them have synchronous or will have later hepatic metastases (10). Furthermore, the current available treatments offer essentially no hope for cure to those patients with unresectable liver metastases.
Here, the aim is to perform an extensive analysis of the system behavior in a subspace of the model parameter space. Our results reveal that, for intratumoral virus administration, the solid tumor can be completely eradicated or keep growing after a transient remission depending on the parameter values. The regions in parameter space associated to each one of such tumor responses were determined. Finally, our findings are compared with the reported results from clinical trials and the predictions of two other mathematical models based on systems of ordinary (11) or partial (12) differential equations. Differently from them, the present model assumes that each cell is an individual agent and uses a sigmoidal cell infection rate (instead of a linear one).
Materials and Methods
The mathematical model represents the tissue as a square lattice fed by a single capillary vessel at the top of the lattice. As illustrated in Fig. 1, four different cell types are considered: normal cells (σn), dead cells (σd), and uninfected (σc) and infected (σv) tumor cells. In contrast to the normal and dead cells, one or more uninfected or infected cancer cells can pile up in a given site, reflecting the fact that the division of tumor cells is not constrained by contact inhibition. Normal cells cannot pile up because they are not dividing. Piled, infected cancer cells result from infection of uninfected cancer cells in a pile. The nutrients diffuse from the capillary vessel through the tissue. Each cancer cell, randomly selected with equal probability, can carry out one of three actions: replication, death, or becoming infected if yet uninfected. In contrast, an infected cancer cell does not divide because its slaved cellular machinery is focused on virus replication. It is assumed that infected cancer cells sustain their metabolism until lysis. In addition, infected cancer cells mainly die by lysis, and the death due to insufficient nutrient supply is neglected. The lysis of each infected cancer cell releases v0 viruses to the extracellular medium. Because virions are present throughout the infected cell, a random fraction of these viruses remains on the site of the lysed cell, and the rest is uniformly distributed among its four nearest neighbors at the time of lysis. The infection of a normal cell by the oncolytic virus is neglected. Thus, the model assumes perfect viral selectivity for cancer cells. At last, the clearance rate (γv) embodies the complex innate and adaptive immune responses to a virus. Such response involves the synthesis of antiviral cytokines, activation/selection of immune cells, and production of antiviral antibodies (13). The free viruses diffuse through the tissue. The mathematical details of the rules and partial differential equations controlling cellular actions, nutrient diffusion, and virus spreading are described in Appendix.
The model considers normal cells, tumor cells, and oncolytic viruses. Uninfected cancer cells can replicate by mitosis or die due to starvation (lack of nutrients) or be infected by viruses. Such cell actions are stochastic and governed by probabilities Pdiv, Pdel, and Pinf dependent on the nutrient and virus concentrations per cancer cell. Analogously, normal and infected cancer cells die due to starvation and lysis with probabilities Pdel and Plysis, respectively. The functional forms of such cell action probabilities are stated in Appendix. After mitosis, the new tumor cell piles up (loss of contact inhibition) at the site of its parent cell if the last is inside the tumor. Otherwise, if the dividing cell is on the tumor border, its daughter cell will randomly replace one of their normal or necrotic nearest neighbor cells. Nutrients and viruses are modeled as continuous concentrations that evolve in space (tissue) and time according to reaction-diffusion equations described in Appendix. Essentially, nutrients diffuse at a rate D from the capillary vessel through the tissue and are taken up by normal and tumor cells at rates γ and λγ, respectively. In addition, free viruses diffuse at a rate Dv, are cleared at a rate γv, and have as source infected cancer cells, which release v0 new viruses after their lysis. The virus clearance rate represents the simplest way to model viral loss due mainly to the innate and adaptive immune responses mediated by antibodies, CD8 T cells, IFNs, and other cytokines or inherent viral instability. It can also take into account other mechanisms such as the virus outflow from the tumor driven by the gradient of interstitial fluid pressure and eventual virus trapping by the extracellular matrix. γv is always a constant, but several distinct values were simulated.
The model considers normal cells, tumor cells, and oncolytic viruses. Uninfected cancer cells can replicate by mitosis or die due to starvation (lack of nutrients) or be infected by viruses. Such cell actions are stochastic and governed by probabilities Pdiv, Pdel, and Pinf dependent on the nutrient and virus concentrations per cancer cell. Analogously, normal and infected cancer cells die due to starvation and lysis with probabilities Pdel and Plysis, respectively. The functional forms of such cell action probabilities are stated in Appendix. After mitosis, the new tumor cell piles up (loss of contact inhibition) at the site of its parent cell if the last is inside the tumor. Otherwise, if the dividing cell is on the tumor border, its daughter cell will randomly replace one of their normal or necrotic nearest neighbor cells. Nutrients and viruses are modeled as continuous concentrations that evolve in space (tissue) and time according to reaction-diffusion equations described in Appendix. Essentially, nutrients diffuse at a rate D from the capillary vessel through the tissue and are taken up by normal and tumor cells at rates γ and λγ, respectively. In addition, free viruses diffuse at a rate Dv, are cleared at a rate γv, and have as source infected cancer cells, which release v0 new viruses after their lysis. The virus clearance rate represents the simplest way to model viral loss due mainly to the innate and adaptive immune responses mediated by antibodies, CD8 T cells, IFNs, and other cytokines or inherent viral instability. It can also take into account other mechanisms such as the virus outflow from the tumor driven by the gradient of interstitial fluid pressure and eventual virus trapping by the extracellular matrix. γv is always a constant, but several distinct values were simulated.
Figure 1 also illustrates some of the multiple processes involved in the model. At the tissue level (macroscopic scale), the nutrient molecules and viruses are both very small particles in comparison with cells. Indeed, virus species tested in cancer clinical trials form particles with sizes that vary from 18 to 30 nm (picornaviruses) to 300 to 400 nm (vaccinia viruses; refs. 1, 6), whereas a cell is typically 10 μm in diameter. These molecules and viruses disperse in the tissue like Brownian particles. Thereby, reaction-diffusion equations are used to describe their dynamics. In these equations, normal and cancer cells are sinks of nutrients; cancer cells become sinks of viruses on infection; and infected tumor cells act as sources of viruses at lysis. Thus, the evolution in time of cells determines the macroscopic distribution of nutrients and free viruses. In turn, cell division, death, infection, and lysis, occurring at the cellular or mesoscopic scale, are stochastic processes controlled by the local concentration of nutrients and viruses determined at the macroscopic scale. Such cell actions occur at a time scale (∼10 h) much greater than that associated to nutrient and virus diffusions (∼10 s). Hence, cell dynamics evolves in a stationary state of these concentration fields.
Summarizing, both mesoscopic and macroscopic scales are inexorably interwoven, although described in terms of distinct physical models: partial differential equations (macroscopic level) and probabilistic cellular automata rules (mesoscopic level) coupled in a single model through division, death, infection, lysis, uptake, release and adsorption rates, and diffusion coefficients.
Results
We are interested in the effects of virotherapy on solid tumors commonly used in experimental assays. Figure 2 shows the total number of uninfected cancer cells and corresponding virus concentrations as functions of time for different virus diffusion and clearance rates. The delay between virus and uninfected tumor cell dynamics is due to the period for cell lysis (Tl). Depending on the parameter values (Table 1), the tumor can be completely eradicated few time steps following a single intratumoral virus administration, undergo damped oscillations leading to eradication at a later time, or sustain an overall tumor growth. In this case, after a transient shrinkage of the neoplastic mass, the growth can be oscillatory and at a reduced rate or monotonous (nonoscillatory) and even at an increased rate. Without treatment, the progress in time of cancer cell populations for all the simulated patterns is Gompertzian (9, 14). The curves in Fig. 2C and D correspond to failed therapies because the cancer keeps growing with time.
Temporal evolution of the number of uninfected cancer cells (Nc) and total free virus concentration with fixed virus infectivity 𝛉inf = 0.03 and clearance rate γv = 0.01. A single intratumoral virus load is administered at the time indicated by the arrow. The following virus parameters were used: Dv = 0.7 and Tl = 4, rapid eradication (A); Dv = 0.8 and Tl = 2, slow oscillatory eradication (B); Dv = 0.2 and Tl = 2, oscillatory growth at a reduced rate (C); and Dv = 0.9 and Tl = 32, exponential growth (D). These curves represent the most probable behavior observed for the corresponding sets of model parameters. Each one refers to a single simulation. In all of them, the viruses inside infected cells are not counted for determining free virus concentration in the system.
Temporal evolution of the number of uninfected cancer cells (Nc) and total free virus concentration with fixed virus infectivity 𝛉inf = 0.03 and clearance rate γv = 0.01. A single intratumoral virus load is administered at the time indicated by the arrow. The following virus parameters were used: Dv = 0.7 and Tl = 4, rapid eradication (A); Dv = 0.8 and Tl = 2, slow oscillatory eradication (B); Dv = 0.2 and Tl = 2, oscillatory growth at a reduced rate (C); and Dv = 0.9 and Tl = 32, exponential growth (D). These curves represent the most probable behavior observed for the corresponding sets of model parameters. Each one refers to a single simulation. In all of them, the viruses inside infected cells are not counted for determining free virus concentration in the system.
Parameter values used for a compact morphology (𝛉div = 0.3, 𝛉del = 0.03 without cell migration)
Parameter . | Description . | Value . |
---|---|---|
D | Glucose diffusion constant | 3.6 × 10−1 mm2 h−1 |
γ | Nutrient uptake rate of normal cells | 10−6–5 × 10−4 s−1 |
λ | γc/γ | 10 |
Dv | Virus diffusion coefficient | 10−2 mm2 h−1 |
γv | Virus clearance rate | 2.5 × 10−2 h−1 |
Tl | Characteristic time for cell lysis | 20 h |
𝛉inf | 0.911 × MOI, assuming Pinf = 70% just after intratumoral injection | MOI ϵ [0.01,0.1] |
Δ | Linear cell size (lattice constant) | 10 μm |
Parameter . | Description . | Value . |
---|---|---|
D | Glucose diffusion constant | 3.6 × 10−1 mm2 h−1 |
γ | Nutrient uptake rate of normal cells | 10−6–5 × 10−4 s−1 |
λ | γc/γ | 10 |
Dv | Virus diffusion coefficient | 10−2 mm2 h−1 |
γv | Virus clearance rate | 2.5 × 10−2 h−1 |
Tl | Characteristic time for cell lysis | 20 h |
𝛉inf | 0.911 × MOI, assuming Pinf = 70% just after intratumoral injection | MOI ϵ [0.01,0.1] |
Δ | Linear cell size (lattice constant) | 10 μm |
Abbreviation: MOI, multiplicity of infection.
A striking and novel result shown in Fig. 2 is the oscillatory behavior of both cancer cells and viruses seen at higher viral diffusivity and clearance rates. It contrasts to the long-term, spatially uniform profiles for cells and free viruses obtained by Wu and colleagues (15). Instead of only damped oscillations, the total numbers of viruses and uninfected and infected cancer cells can also undergo stable aperiodic oscillations superimposed on a steady growth.
In Fig. 3, snapshots of the simulated spatial patterns of uninfected tumor cells are shown. The parameter values used are those for the curve C in Fig. 2. Successive waves of virus infection and therefore cell death by lysis are neatly evident. However, these rounds of infection are not strong enough to eradicate the tumor. Videos S1 and S2 for unsuccessful and successful treatments, respectively, are provided as supplementary information.
Snapshots of cancer cell patterns at intervals of 3 time steps after virus injection. The first pattern corresponds to the solid tumor immediately before the treatment. The parameters are the same used for the oscillatory growth at a reduced rate shown in Fig. 2C. Such spatiotemporal evolution can hardly be observed experimentally either in vitro or in vivo. Indeed, histologic samples require killing of the cells, and noninvasive methods (PET/fluorescently labeled virus) often provide fuzzy pictures characterized by a low spatial accuracy.
Snapshots of cancer cell patterns at intervals of 3 time steps after virus injection. The first pattern corresponds to the solid tumor immediately before the treatment. The parameters are the same used for the oscillatory growth at a reduced rate shown in Fig. 2C. Such spatiotemporal evolution can hardly be observed experimentally either in vitro or in vivo. Indeed, histologic samples require killing of the cells, and noninvasive methods (PET/fluorescently labeled virus) often provide fuzzy pictures characterized by a low spatial accuracy.
A fundamental feature uncovered by the simulations of the model is its great sensitivity to the local, stochastic fluctuations of uninfected cell population and free virus concentration. Thus, for instance, although the complete eradication of the solid tumor is the most frequent therapeutic outcome, the other two types of tumor response are also observed for γv = 0.01, Dv = 0.7, and Tl = 4. Hence, it becomes imperative to find the dominant tumor response in every region of the model parameter space. For a given tumor cell type (𝛉inf and Tl fixed), the dominant response was determined for each pair of parameters Dv and γv controlling free virus spreading. Whatsoever the values of Dv and Tl, the therapy always fails if the virus clearance rate is large (γv ≥ 0.03). However, a small decrease in γv with all the other parameters (Dv, Tl, and 𝛉inf) fixed can lead to complete tumor eradication (see Supplementary Fig. S1). Thereby, virus clearance rate is a key factor determining therapeutic outcome. In turn, in a tissue characterized by a moderate clearance rate (0.01 ≤ γv < 0.03), if 𝛉inf ≥ 0.1 (i.e., the virus has a low efficacy to invade uninfected cancer cells), the tumor mass steadily grows and, again, the treatment fails. In Fig. 4, the boundaries corresponding to the various therapeutic outcomes of an aggressive virotherapy are shown.
Boundaries in parameter space associated to different tumor responses to an aggressive virotherapy (γv = 0.01, low virus clearance rate; 𝛉inf = 0.03, high virus infectivity). In regions labeled A, the cancer is eradicated (absolutely successful treatment); in C, the tumor growth progress at a reduced rate and oscillating, whereas in D the growth is monotonic at a rate that can even be greater than that of untreated tumors. In both cases the therapy fails. This “phase diagram” was drawn from simulations involving 200 independent samples (runs) for 50 points in parameter space arranged as a square lattice. Each simulation lasts for 500 time steps or ∼3 mo, because a time step is estimated as ∼4 h.
Boundaries in parameter space associated to different tumor responses to an aggressive virotherapy (γv = 0.01, low virus clearance rate; 𝛉inf = 0.03, high virus infectivity). In regions labeled A, the cancer is eradicated (absolutely successful treatment); in C, the tumor growth progress at a reduced rate and oscillating, whereas in D the growth is monotonic at a rate that can even be greater than that of untreated tumors. In both cases the therapy fails. This “phase diagram” was drawn from simulations involving 200 independent samples (runs) for 50 points in parameter space arranged as a square lattice. Each simulation lasts for 500 time steps or ∼3 mo, because a time step is estimated as ∼4 h.
Finally, in addition to determining the most probable prognosis after virotherapy, it is relevant to know the chances for the occurrence of the other tumor responses in every region in the parameter space. This is especially true for the probability of therapeutic success. For the aggressive treatment considered here, those probabilities are shown in Fig. 5.
Probability for tumor eradication (A), monotonic growth (B), and oscillatory growth (C) as a function of the parameters Tl and Dv. Both 𝛉inf = 0.03 and γv = 0.01 remain fixed. The probabilities were evaluated from samples including 200 independent, 500 time-step runs for each pair (Tl, Dv) considered.
Probability for tumor eradication (A), monotonic growth (B), and oscillatory growth (C) as a function of the parameters Tl and Dv. Both 𝛉inf = 0.03 and γv = 0.01 remain fixed. The probabilities were evaluated from samples including 200 independent, 500 time-step runs for each pair (Tl, Dv) considered.
Discussion
Computer simulations of our multiscale model revealed that, after intratumoral virus administration, a solid tumor can completely be eradicated or keep growing after a transient remission. Furthermore, depending on the parameter values, the tumor cell and virus populations exhibit undamped oscillatory dynamics around increasing, nonstationary values (Fig. 2). The most likely therapeutic responses are determined by the oncolytic activity (𝛉inf and Tl) and spreading properties (Dv and γv) of the virus. These are the key parameters that control treatment success. For example, independently of the viral diffusivity, infectivity, and time for cell lysis, a high clearence rate (γv ≥ 0.03) leads to a neoplastic mass and an invaded area that grow monotonously and even faster than without treatment. Because the immune response directed at free virus and viral antigens expressed on the surface of infected cancer cells enhances γv, this result indicates the necessity of a strong ablation of the patient immune system. Virotherapy also fails when the tumor growth is oscillatory and progresses at a small rate. This is the most frequent outcome of varied treatments. For instance, the use of an oncolytic virus with a low infectivity (𝛉inf ≥ 0.1) in a tissue characterized by a moderate virus clearance rate (0.01 ≤ γv ≤ 0.03). Here, the failure is caused by the inefficiency of the viruses at infecting tumor cells, not by the strength of the host immune response. The same outcome can be obtained in such a tissue for oncolytic viral agents with greater infectivity (0.01 ≤ 𝛉inf < 0.03) but a small diffusivity (Dv ≤ 0.2; Fig. 3). Now, instead of the host immune response or virus inefficacy to infect tumor cells, the therapeutic failure is due to infection waves expanding slower than the neoplastic mass. Moreover, because a fraction of the tumor cells are killed by the virus, the nutrient supply available for the uninfected cells becomes larger. In response, they divide more frequently, spread, and invade a larger area of the tissue.
With regard to a successful therapy, our results indicate that the key factor is really the virus clearance rate. Indeed, in a host tissue with small virus clearance rate (0.005 ≤ γv < 0.01), even mild oncolytic viruses (0.03 ≤ 𝛉inf < 0.1) with a rapid virus spreading (Dv ≥ 0.3) can eradicate the tumor. In agreement, in vivo experiments reveal that the potency of oncolytic viruses is markedly reduced when moved to an immunocompetent model (16). It is also important to note that, in an immunocompentent host, the evolving immune response to the virus will markedly change the effective viral clearance rate in a complex way that is currently difficult to predict (17).
The time for lysis cannot be either too short (Tl < 2 time steps) or very long (Tl > 16 time steps). Indeed, due to the nonhomogeneous virus dynamics, new infection waves cannot reach the tumor border because either they are infrequent and thus start far from this border or the large “voids” generated inside the tumor by rapid cell lysis cannot provide further fuel for such waves. Hence, from a clinical point of view, our findings indicate that a successful, single agent virotherapy requires a strong inhibition of the host immune response and the use of potent virus species with a high intratumoral mobility. A rapid immune-mediated clearance can be avoided through maintenance or incorporation of virus-encoded immune response modifier genes or nongenetic approaches such as coating the viruses with polymers (1, 6). Possibly, viral genes associated with immune evasion could be used to reduce the viral clearance rate (18). In contrast, deliberate up-regulation of the immune response to target cancer cells is also being considered (19). Approaches to increasing oncolytic efficacy include arming the viruses with directly lethal (20) or prodrug-converting (21) enzymes. Enhanced intratumoral diffusivity is possible if, for instance, the virus is engineered to incorporate a gene encoding for a matrix metalloprotease. This enzyme is secreted by the infected cancer cell before its lysis and breaks down the surrounding extracellular matrix, facilitating the spreading of the viruses released after cell death (22).
The above scenario about the most frequent therapeutic outcomes is a simplified one because their corresponding regions in the model parameter space are irregular and intertwined (Fig. 4). Furthermore, the model is sensitive to the local, stochastic fluctuations of uninfected cell population and free virus concentration. As a result of such sensitivity, the history of cancer progression and therapeutic outcome can drastically vary even for a given set of fixed parameters.
A relevant issue is how our results compare with those obtained from clinical trials (1, 6). Even with intratumoral administration, wild or engineered adenoviruses, herpes viruses, and reoviruses have characteristically induced only transient tumor remissions (6). Thus, for instance, single intratumoral injection of hrR3 virus into HT29 flank tumors in nude mice has only reduced the rate of monotonic growth of this human colon carcinoma (10), a response predicted in Fig. 2D. In addition, transient antitumoral responses were observed after single intratumoral administration of Reolysin (reovirus) and dl1520 (adenovirus ONYX-015) into superficial (23) and head and neck (24–26) cancers. Because intratumoral spread seems to be inefficient with adenovirus (ref. 6; small diffusivity), Figs. 4 and 5 can explain why virotherapy using dl1520 fails. Unfortunately, in these phase I and phase II clinical trials with 11% and 14% of positive responses, respectively, it was not reported if the tumor growth progress either in a monotonic or oscillatory way. A neat indication of such oscillatory tumor growth is observed in an ovarian cancer xenograft model attacked with MV-CEA (an engineered measles virus) administered intratumorally (27). By contrast, vaccinia viruses have induced durable complete responses in patients after intratumoral or even i.v. injection. Specifically, the BP (28) and the engineered JX-594 (29) vaccinia viruses eradicated melanomas in 4 (13%) and 2 (29%) patients, respectively. These successful therapeutic outcomes, corresponding to Fig. 2A, or B, have frequencies consistent with eradication probabilities smaller than 35% obtained from the model simulations for Dv ≤ 0.6 and Tl > 5 h (Fig. 5A). Again, a long-term oscillatory eradication was not evidenced in reported clinical trails. Accurate monitoring of tumor burden is often difficult due to irregularly shaped and inaccessible tumors. This, combined with the contribution of dead cells to the tumor mass, would make detection of oscillating responses very difficult. New noninvasive imaging techniques such as high-sensitivity fluorescence imaging and positron emission tomography (PET) scanning (30) could facilitate quantitative monitoring of both living tumor mass and viral load in vivo. This approach could be used to test the predictions of our model.
Very recent results obtained by Nagano and colleagues (31) reveal the spatial distribution of herpes simplex virus inside tumors developed in mice. These authors show that the induction of cancer cell apoptosis before the virus administration enhances intratumoral virus spreading and antitumor efficacy by generating voids and channels inside the tumor. With the main barriers against virus spreading disrupted, intratumoral diffusion is likely to be closer to the regimen assumed in our model (i.e., viruses can always diffuse through the neoplastic mass). Indeed, the spatial patterns of virus distribution shown in Fig. 2A and B of Nagano's experiment are very similar to the free virus simulated patterns exhibited in our Supplementary videos S1 and S2.
Although the present model does not consider the natural barriers to virus spreading (fibrillar collagen in the extracellular matrix, narrow space between tumor cells, and contrary interstitial pressure gradient) or intratumoral heterogeneity in the virus diffusivity (as in necrotic areas, for instance), it can predict that virotherapy has a great chance to fail if virus diffusivity is very low, in complete agreement with the findings of Nagano and colleagues (31). Furthermore, our results indicate that virus diffusivity is only part of the puzzle. Larger diffusivity is desirable, but a low virus clearance rate and an adequate lysis time are also necessary to sustain successive waves of infection required for tumor eradication.
Finally, our results are compared with those provided by two other mathematical models. Wein and colleagues (12) modeled the dynamics of oncolytic viruses through a deterministic set of partial differential equations. They found that, in the case of a uniform virus injection, the tumor mass either grows exponentially or is eradicated if the waves of infection can or cannot outpace the tumor proliferation. Oscillations in time of the virus population were explicitly shown. Furthermore, because the wave speed of the infection is inversely proportional to Tl (the time for cell lysis), these authors suggest that the most efficacious viruses will be those that kill an infected tumor cell faster. In contrast, our results indicate that there is an optimal range for Tl. Indeed, if the time for cell lysis is very short, the waves of infection cannot control the neoplastic growth. A few isolated cancer cells, distant from the large voids produced by rapid cell lysis and plenty of free viruses looking for targets, provide the opportunity to tumor escape. This fundamental and novel prediction seems to emerge from the stochastic nature of our model in which the entry of a virus in a cell, the lysis of an infected cell, the replication or death of a uninfected cell, etc., are all random, individual events. Cell discreteness and stochasticity of cellular responses deeply affect the strength of the fluctuations, particularly at tumor regions in which the number of cancer cells is very low. Surprisingly, the deterministic, ordinary differential equation model proposed by Wodarz (11) also predicts an optimal viral cytotoxicity. In this model, the infection of a cancer cell occurs only through contact with an infected cell, which, in turn, is killed by the virus at a rate a, the viral cytotoxicity. According this model, if there are uninfected cancer cells, an increase in the viral cytotoxicity increases the tumor mass because infected cells are eliminated before the virus had a chance to significantly spread. The reciprocal of a should be interpreted as the time for cell lysis Tl. However, this model can only exhibit damped oscillations in tumor cell populations.
In summary, we investigate a multiscale model for virotherapy of cancer in which the cells are discrete agents. Our predictions are in qualitative agreement with results from clinical reports. Through extensive simulations, it was found that, for a single intratumoral virus administration, a solid tumor can completely be eradicated or keep growing after a transient remission. Furthermore, the model reveals undamped oscillatory dynamics of tumor cells and virus populations, which demands for new in vivo and in vitro quantitative experiments aiming to detect this oscillatory response. The conditions (regions in the model parameter space) for domination of each one of the different tumor responses, as well as the occurrence probabilities for the other nondominant therapeutic outcomes, were determined. From a clinical point of view, our findings indicate that a successful, single agent virotherapy requires a strong inhibition of the host immune response and the use of potent virus species with intratumoral high mobility. Moreover, due to the discrete and stochastic nature of cells and their responses, an optimal range for viral cytotoxicity is predicted because the virotherapy fails if the oncolytic virus demands either a too short or a very large time for killing the tumor cell. This finding suggests that the virus that kills cancer cells most rapidly is not necessarily the more effective agent to eradicate the tumor. The implications of such a result for the design of new replication-competent viruses are clear.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).
Acknowledgments
Grant support: Brazilian Agencies Coordenação de Aperfeiçoamento de Pessoal de Nível Superior, Conselho Nacional de Desenvolvimento Científico e Tecnológico, and Fundação de Amparo à Pesquisa do Estado de Minas Gerais.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
We thank M.L. Munford (Department of Physics, UFV) for preparing Fig. 1.