Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Wiley Interdiscip Rev Syst Biol Med. Author manuscript; available in PMC 2017 May 1.
Published in final edited form as:
PMCID: PMC4896310

Computational modeling approaches to the dynamics of oncolytic viruses


Replicating oncolytic viruses represent a promising treatment approach against cancer, specifically targeting the tumor cells. Significant progress has been made through experimental and clinical studies. Besides these approaches, however, mathematical models can be useful when analyzing the dynamics of virus spread through tumors, because the interactions between a growing tumor and a replicating virus are complex and nonlinear, making them difficult to understand by experimentation alone. Mathematical models have provided significant biological insight into the field of virus dynamics, and similar approaches can be adopted to study oncolytic viruses. The review discusses this approach and highlights some of the challenges that need to be overcome in order to build mathematical and computation models that are clinically predictive.


Modern treatment approaches against cancer are moving beyond traditional chemotherapies and aim to specifically target cancer cells while sparing healthy cells. One such approach that has shown great promise is the use of targeted therapies or small molecule inhibitors1. These drugs specifically target molecular pathways that are responsible for the aberrant behavior of cancer cells. They can abolish the ability of cancer cells to divide and can also cause cell death, among other effects. They have shown great promise in leukemias, and are emerging to be important in the treatment of other cancers as well. Another treatment approach that also specifically targets cancer cells, but that is less well established, is the use of oncolytic viruses26. These are viruses that have been specifically engineered to infect cancer cells, while sparing the healthy cells. Common mutations in cancer cells can allow virus replication, while the absence of those mutations in healthy cells tend to abort the replication cycle of the virus. The virus population multiplies in the cancer cells, offspring virus is released, and these offspring viruses go on to infect further cancer cells. In this way, the virus population is supposed to spread through the tumor cell population, kill the tumor cells, and drive the disease into remission. While this is theoretically a desirable idea, the treatment of cancers with oncolytic viruses has encountered a variety of difficulties. For example, anti-viral immune responses can have a negative impact on virus spread, and so can barriers within the tumor, such as the presence of extracellular matrix. Nevertheless, success has been observed, and oncolytic virus therapy is gaining in clinical relevance. In China, the adenovirus H101 (Shanghai Sunway Biotech, Shanghai, China) was approved for the treatment of head and neck cancer in combination with chemotherapy. In the United States, the herpes-based virus T-VEC (Amgen Inc)7 has recently been approved by the FDA. While anti-viral immune response can hinder treatment success, oncolytic viruses are also thought to induce immune responses against tumor antigens, presumably by enhancing presentation of tumor antigens. It is possible that this is a critical mechanism underlying the success of oncolytic virus therapy8.

Experiments and clinical studies have been driving progress in this field. In addition to these approaches, however, mathematical modeling has emerged as a capable tool that can help us understand why certain treatment regimens might not lead to success, and to suggest mechanisms by which the outcomes can be improved9. The interactions between a growing cancer cell population and a virus population, and possibly a variety of immune responses, are highly complex and non-linear. This makes it impossible to use verbal or graphical reasoning to predict outcomes. Instead, mathematical models can capture crucial biological assumptions in equations, and follow them to their precise logical conclusions. Mathematical models can be used to interpret biological data, to make specific, testable predictions, and to measure biologically important parameters. Mathematical models have a long history in the field of virus dynamics in general10, and equivalent approaches have been used to study the dynamics of oncolytic viruses. As will be discussed throughout this review, model properties and results depend on the exact assumptions and formulations upon which the models are built. This in turn depends on our biological/experimental understanding of the virus, the cancer cells, and the virus-cell interactions, which is of course still incomplete.

Much work has been performed in this field of research. Rather than providing a comprehensive review of this growing area of research, this article will summarize the basic approach and highlight some important challenges and difficulties that need to be investigated further in order to construct models that are truly predictive and are useful for clinical research.

Basic model of virus dynamics, and its application to oncolytic virus dynamics

A basic framework to study virus dynamcis in the context of a variety of human pathogens11 takes into account three populations: uninfected cells, x, infected cells, y, and free virus, v. The model is given by ordinary differential equations, that describe the average behavior of these populations over time. The equations are given as follows.


Uninfected cells are produced with a rate λ, die with a rate δ, and become infected by virus with a rate β. Infected cells produce offspring virus with a rate k, and die with a rate a. Free offspring virus decays with a rate u. This model, and variations thereof, have been used to study the dynamics of human pathogens, such as human immunodeficiency virus (HIV), and Hepatitis B (HBV) and Hepatitis C (HCV) viruses, as well as other viral infections10,12,13. A common simplification with such a model is to assume that the turnover of the free virus population is much faster than that of the infected cells. In this case, it can be assumed that free viruses are in a quasi-steady state, and the model can be re-written in the following, more compact form.


The same basic framework can be used to study oncolytic viruses. While tissue target cells are produced with a constant rate in model (2), this assumption does not apply to tumor cells. Tumor cell population dynamics are governed by cell division. Tumor growth is initially exponential, but can follow a variety of laws in the long term14. As an example, logistic growth of the tumor cell population can be assumed. In this case, the basic model is given by the following pair of ordinary differential equations.


The division rate of tumor cells is denoted by r, and the parameter ω described the carrying capacity of the cells, i.e. the maximum number of tumor cells that can be sustained (e.g. due to nutrient or spatial limitations). This model has been analyzed in some detail15,16.

An important concept with virus dynamics models is the basic reproductive ratio of the virus. This describes the average number of newly infected cells generated by a single infected cell when placed into an uninfected tumor cell population. It is denoted by R0 and decides whether the virus will establish a successful infection or not. If R0 > 1, then one cell on average gives rise to more than one newly infected cell, and a successful infection is established. In the opposite case, the infection goes extinct. For the above model, the basic reproductive ratio can be calculated and is given by R0=βω(r-δ)/ra. This expression can be obtained by standard stability analysis of this system of ODEs.

Basic insights

This type of model can be used to study the spread of an oncolytic virus through a tumor cell population, starting with a relatively low inoculum, which likely corresponds to relevant clinical settings. Because the model is based on ordinary differential equations, note that perfect mixing of cells and viruses is assumed, which can be unrealistic for solid tumors. Nevertheless, such approaches can be insightful, and a spatially explicit modeling approach is reviewed below. In the current framework, if the basic reproductive ratio of the virus (defined above) is greater than one, the virus will spread through the tumor in the model, and in the long term the populations will converge to a stable equilibrium. In the current modeling framework (based on ordinary differential equations), the virus cannot drive the tumor cell population extinct (i.e. to a population size of zero), because the model is deterministic. Instead, a very efficient virus can only reduce the target cell population to very low levels. If the number of tumor cells fall below a threshold value in the model (i.e. below an average population size of one cell), then extinction is the likely outcome in a stochastic setting. Equilibrium tumor load in the model is given by x(1)+y(1)=ω(a+r)r+βω. An important result is that a lower virus-induced death rate of infected cells (small a) results in lower equilibrium tumor load15,16 (Figure 1). If the rate of virus-induced tumor cell killing is too high, the outcome is persistence of the tumor in the face of ongoing viral replication (Figure 1). The reason is as follows. Low viral cytopathicity increases virus load. Higher virus load results in more infection and in a greater decline of the uninfected cancer cells. Higher viral cytopathicity results in lower virus load. Low virus load results in less infection and in less reduction of the uninfected tumor cells. Note that this argument assumes that an increase in viral cytopathicity is not accompanied by a higher rate of virus replication. If a higher viral cytopathicity correlates with a faster viral replication rate, this result holds as long as the increase in cytopathicty is stronger than the increase in the replication rate of the virus. In addition, equilibrium tumor load is reduced by other parameters, most notably by a high replication rate of the virus (high value of β’′), and a slow growth rate of the tumor (low value of r). As expected, the replication rate of the virus (value of β) has to lie above a threshold for tumor load to become low enough that tumor eradication is feasible in a stochastic setting15,16. If the replication rate of the virus is too low, the virus does not spread sufficiently through the population of tumor cells. In summary, if the virus is replicating, the best strategy is to (i) use a weakly cytopathic virus, (ii) use a fast replicating virus, and (iii) reduce the growth rate of the tumor by alternative therapeutic means.

Figure 1
Simulation of therapy using a replicating virus, model (1). The virus is administered once, as indicated by the arrow. Shading indicates the phase of the dynamics following administration of the virus. (i) Use of a weakly cytopathic virus results in sustained ...

Immune responses have been incorporated into this basic model, and again model properties were explored in relatively simple settings15. Two types of immune responses were included: (i) virus-specific immune responses, and (ii) tumor-specific immune responses, particularly concentrating on cytotxic T lymphocytes (CTL). Anti-viral CTL can both promote and counter treatment success. They promote treatment success because they kill infected tumor cells. They counter treatment success because they limit virus spread and can potentially clear the virus from the tumor. A basic model suggested that there is an optimal strength of an anti-viral CTL response that maximizes the chances of tumor control15. Regarding tumor specific CTL responses, it is thought that they can be stimulated or induced by an oncolytic virus because virus-induced cell death triggers the availability of tumor antigen to the immune system. These dynamics were also included in the model, focusing on conditions that maximize this effect. Interestingly, if different immune responses arise that both fight a common enemy, competition between immune responses for antigen can occur17. In such a situation, mathematical models have shown that the stronger response that requires less stimulation can outcompete and exclude a weaker response that requires higher levels of antigen for stimulation. Such competition dynamics can also occur between tumor-specific and virus-specific immune response if both responses require the presence of the virus for stimulation. This can lead to intricate interactions and complex outcomes, which have been explored by modeling15.

These are a few insights that have emerged from the simplest possible models of oncolytic virus spread and can be useful to keep in mind when interpreting experimental data. Many additional biological complexities can be incorporated into this framework, and depending on the exact assumptions, some of the conclusions may or may not change. Different mathematical models of oncolytic virus spread, based on the framework described here, have been explored and applied to specific biological systems, with examples given in references1831,26,3237. Rather than reviewing all these insights here, the emphasis of the current review is different. It discusses some important difficulties that are inherent in such modeling approaches, and these difficulties currently pose a significant barrier to further progress, especially when it comes to fitting mathematical models to experimental or clinical data, estimating parameters, and making robust predictions. I hope that this review will contribute to the increased awareness of these difficulties and motivate further work to overcome them.

Model uncertainties

Many terms in the mathematical models are of uncertain nature. For example, the infection process of cells can be described mathematically in a variety of ways, and it is not known which description is the “correct one”. Because many mathematical models assume that the rate of infection is simply proportional to the number of target cells and the amount of virus, it was tested whether this can describe virus infection experiments performed in vitro38. This work was not done with oncolytic viruses specifically, but utilized human immunodeficiency virus infection of jurkat cells as a model system. It is nevertheless highly relevant because the aim was to investigate how to mathematically describe the rate of infection in an in vitro setting. Single round replication experiments were performed, where different amounts of virus were used to infect the target cell population, and the number of resulting productively infected cells were recorded. The virus was labeled with fluorescent reporters, and the number of productively infected cells was determined by measuring the number of fluorescent cells. This was repeated, using different number of target cells in the culture. The basic virus dynamics model was fit to the data with the lowest target cell numbers, and the best fit provided parameter estimates. The parameterized model was then used to predict the data for the experiments with larger target cells numbers in order to determine whether the model could make the correct prediction or not. It turned out that the standard model with the infection term βxy could not correctly predict the data. While it could be fit to the data for a given number of target cells, the prediction for the experiments with the other target cell numbers was incorrect. The model consistently over-predicted numbers of productively infected cells observed in the experiments (Figure 2). This indicated that the rate of infection is not proportional to the number of target cells available for virus infection.

Figure 2
Infection terms and model predictions. Single-round virus infection experiments were performed. Different doses of the virus were used to infect the target cell population, and the number of productively infected cells was determined by fluorescence after ...

These data suggest that the rate of infection is a saturating function of the number of target cells available to the virus. This can be described phenomenologically by expressing the infection term as (1+epsilon)βxy/(x+epsilon). For large values of epsilon, the infection term converges to the standard model, βxy. For small values of epsilon, the infection term converges to βy. When such a model was fitted to the data for a given target cell number, it could correctly predict the infection levels for the other target cell numbers, by adjusting the value of x in the model (Figure 2). Therefore, this is a more accurate term to describe infection than the simplest term that is used in many mathematical models of virus dynamics. Because this term is phenomenological in nature, however, it is not based on an underlying biological mechanism. It was hypothesized that a reason for a saturated rate of infection could be the following. Standard models assume that during a time step a virus only has a chance to infect one target cell. In reality, however, it is possible that a virus has simultaneously access to multiple target cells in the neighborhood, especially if a certain degree of mixing occurs in the system. This can be described as follows38. Suppose at each time-step, a number of viruses are randomly applied to a set of cells randomly distributed on a grid. At each attempt, a virus picks a randomly chosen spot, and if the spot contains a cell (infected or uninfected), with probability γI (where i is the serial number of the attempt), the cell becomes infected, and otherwise the next attempt follows. Set 1=γ1γ2γ3γm. A particular (exponential) example of a model for the decay of infection probabilities is given by


where parameter κ measures the rate of decay. Large values of κ correspond to effectively having only very few attempts of infection per virus particle per update; small values of κ correspond to a large number of attempts with similar probabilities of infection. The infection term for such a model has been derived, and under biologically realistic regimes, it can be given by (1+α)βv(1 – e−αS )/α. The larger the value of α, the more cells a virus can access. For small values of α, the infection term converges to the standard expression βxy, and for large values of α, it converges to βy, as in the model above. Figure 2 shows that this model can also accurately predict the experimental data.

Note, that in the above analysis, the properties of different infection terms was compared, each term expressing basic spread dynamics in mathematically different ways. All terms describe a relatively simple biological scenario, most likely corresponding to virus infection of a homogeneous cell culture. In more complex settings and in vivo, other biological processes likely influence the laws of virus infection, such as receptor densities or blood flow. Such more complicated scenarios can be considered once we have obtained a very solid understanding of the more basic dynamics.

With this example, it was shown that modifications of the standard infection model are required to accurately describe the experimental data. At the same time, however, it is clear that different models can describe the experimental data equally well. It is probable that there are further infection models not explored here that can also accurately describe the experiments. Importantly, each of these models has different properties and can lead to different dynamics. This created problems with the robustness and validity of any modeling results, even if a particular model fits a particular set of experimental data. The next section discusses an axiomatic modeling approach, which can be more appropriate in such situations.

Axiomatic modeling

As discussed in the last section, the “biologically correct” term to describe infection is currently not known. Other processes assumed in the model are also described by arbitrary mathematical expressions, such as the term describing cell growth, given by the logistic equation. While this is assumed in many models of tumor growth, the laws underlying tumor growth are likely to be more complex, and this could impact the properties of the virus dynamics model. To address this issue, it might be useful to avoid concentrating on a particular model, and take a more general approach. Processes such as infection and cell growth were formulated in terms of general functions and results were thus independent from biologically uncertain and arbitrary mathematical formulations29. Through specific restrictions about biological assumptions, a class of mathematical models was analyzed that aim to describe viral spread through a tumor in different settings. Details of the mathematical analysis are given in reference24.

It was found that based on the infection term, models can be divided into two categories with fundamentally different behavior29. In one group, virus growth is relatively fast, which could be characteristic of a relatively well mixed system. In these models, there is a clear viral replication rate threshold beyond which the number of cancer cells drops to levels of the order of one or less, corresponding to extinction in practical terms. Under this parameter region, this is the only outcome in this class of model. In the other category, virus growth is relatively slow, which could be characteristic of more intricate tumor cell arrangements. In this class of model, the larger the number of infected cells, the smaller the proportion of cells that can pass on the virus. In this scenario, virus therapy is more difficult. If tumor growth saturates only at relatively large sizes or does not saturate, then even in the parameter regions where the dynamics can converge to a tumor control or eradication outcome, there can be the possibility that the cancer can outrun the virus if the number of cancer cells lies above a threshold at the start of virus therapy. This is because of the existence of the saddle node equilibrium which ensures dependence of the outcome on initial conditions. This might be problematic in clinical settings, because there is only a relatively small window between the size at which the tumor becomes detectable (about 1010 cells) and the size at which it can induce mortality (around 1013 cells). Tumor growth saturation at lower levels introduces a parameter region in which only the tumor control outcome is possible. A further reduction in the number of tumor cells at which growth saturation occurs can abolish the existence of the saddle node equilibrium altogether. In this case, the only outcome is tumor control. This result makes intuitive sense: earlier saturation of tumor growth slows down the cancer, and makes it easier for the virus to gain the upper hand. It also means that if the tumor is found early, it might be possible to slow down tumor growth by means of more conventional drug therapy, enabling the virus to control the cancer and to prevent runaway growth. There is indication in clinical data that a combination of chemotherapy and oncolytic virus therapy leads to better results than either approach alone39.

Another important finding of this study is that the basic results regarding the outcome of oncolytic virus therapy do not depend on the particular tumor growth terms used in the model29. The exact kinetics of tumor growth are still poorly understood and a source of uncertainty40. While there are minor differences in model behavior when comparing different tumor growth laws, the properties of the tumor control equilibrium are largely independent from the exact way in which tumor growth is modeled29.

This analysis showed that models can be grouped into classes that are characterized by similar properties, and this helps understand better how certain results depend on the exact form of the model and whether given results are robust or not. Significant difficulties, however, still remain when mathematical models need to be fitted to experimental or clinical data. This has to rely on specific mathematical formulations, and it has been shown that even if different specific models within the same class are fit to the same data, the parameter estimates and long-term projections can be quite different29. As an example, models of both classes have been fitted to data of oncolytic virus infection in mice. In particular, A549 human lung cancer nude xenografts were infected with the adenovirus Ad309, and the tumor volume over time was recorded in the presence of the infection41. Figure 3a shows that the data can be fit relatively well with two specific models that belong to the group of “slow infection” terms, and one model that belongs to the group of “fast infection” terms. These are arbitrarily chosen models defined in reference41. Figure 3b shows the same fit, with the dynamics continued for a longer time frame that goes beyond the duration of the experiments. As is obvious from the graphs, the three particular models give rise to very different long term dynamics. This applies even to models belonging to the same class, one predicting relatively low, and the other relatively high tumor loads as time progresses. Therefore, axiomatic modeling facilitates our understanding of model properties and the robustness of these properties. When it comes to fitting experimental or clinical data, however, specific models need to be considered and the best strategy for model selection is to collect data that are as detailed as possible and continue for as long a duration as the experimental setup permits.

Figure 3
Fitting mathematic models to data. 549 human lung cancer nude xenografts were infected with the adenovirus Ad309, and tumor volume was recorded over time41. Different mathematical models were fitted to these data, including two models from the “slow” ...

Spatial modeling

So far, the review discussed mathematical modeling approaches that are based on ordinary differential equations. Such models assume that cells and viruses mix well within the tumor and that no spatial structure exits (perfect mixing or mass action). Most tumors, however, have intricate spatial structure and this can influence the dynamics of virus spread. Different spatial modeling approaches have been used that took into account a relatively high level of biological complexity, e.g. references26,3234. Again, issues of robustness can arise since there are biological uncertainties, and the same assumptions about spatial dynamics can be formulated computationally in different ways and can give rise to different properties. In the context of the ordinary differential equations, the discussion focused on the infection term and the uncertainty connected with its formulation. Similarly, a crucial question is how to accurately describe spatially restricted virus replication and spread with computational models. This is best investigated if a simple scenario is considered that is stripped of most biological complexity present in tumors. Once the laws of spatial virus spread have been understood in the simplest setting, this can form the basis for the gradual introduction of further biological complexities. Such a simple setting has been investigated experimentally, and models have been built in order to describe those dynamics. Experimentally, an adenovirus was used to infect a monolayer of 293 cells (a cell type) with agar layover42. This setting prevented the virus from mixing in the culture. Instead, an infected cell could pass its virus only to the nearest neighboring target cells. The virus was labeled with green fluorescent protein, allowing the spatial patterns of virus spread to be investigated in addition to the time evolution of the infection. In these experiments, two patterns of virus spread were observed42. The first pattern was called “robust infection”. The virus population was observed to spread in the form of an expanding hollow ring, where infected cells on the surface of the ring spread the infection outward, while cells in the middle died and gave rise to a hollow core. This represents typical plaque growth. The second pattern was called “diffuse growth”. In such cases, the virus population grew more slowly in a spatial configuration in which infected and uninfected cells were interspersed among each other. No plaque pattern was observed in such cases. The question arose what those different patterns mean in terms of the outcome of the virus infection in this system.

In order to address this question, an agent-based computational model was considered that captured spatial infection dynamics in two dimensions, with virus spread limited to the eight nearest neighbors of the source cell. This model was fitted to the experimental data42, and this is shown in Figure 4. The model could describe the time evolution of cells for the two growth patterns well, and the difference was explained by differences in virus spread parameters. The diffuse pattern arose for parameters that described slower overall virus spread, i.e. a slower rate of virus replication and a higher rate of infected cell death. The model was only fitted to the numbers of cells over time, not to the spatial configurations, since this provided too many technical challenges. Nevertheless, the spatial patterns that arose from those fits corresponded to the experimentally observed patterns on a qualitative level (Figure 4). This model was analyzed further, predicting the long-term outcomes of these two patterns. It was shown that the diffuse infection spread patterns typically resulted in the failure of the virus to eliminate the target cell population. The robust infection pattern, on the other hand, could result either in the virus-mediated elimination of the target cells, or in the long-term persistence of both tumor cells and the virus depending on the exact long-term spatial pattern that developed. Further, longer term experiments will be necessary to investigate this further.

Figure 4
Experimental outcomes of spatial virus spread and model fitting. (a) A “ring structure” can develop where the virus population expands as a growing ring. The model can accurately describe the time evolution of cells and qualitatively reproduces ...

Follow-up work, however, revealed that even this simplest spatial in vitro setting, the dynamics were driven by considerable complexity. According to the model, the difference between the robust and the diffuse growth pattern arose from differences in parameters. Experiments, however, revealed that the two growth patterns could be observed with the same virus and the same target cells. Indeed, the two growth patterns were observed in the same experimental dish when multiple infection foci were seeded. Further modeling and experimental work revealed that this could be explained by a stochastic race between the initial spread of the virus population, and the initial spread of anti-viral states induced by interferon responses. This is subject to current investigation. The central message of this work, however, is that even the simplest spatial virus dynamics scenario is characterized by a high degree of complexity, and that a thorough understanding of the principles underlying simple spatial virus spread will be essential to construct more robust models that contain additional biological complexity. It will also be important to extend these result into a 3D setting, both experimentally (e.g. in a matrigel setting) and computationally. The key dynamics described here are likely to hold if infection is still restricted to the nearest neighboring cells. The difference is that in a 3D setting, a cell will have more immediate neighbors, thus relaxing the degree of spatial restriction to a certain extent.

Three-dimensional computational models of virus spread have been considered before in different contexts. The most closely related study examined oncolytic virus spread in myeloma animal models43. Data from murine studies clearly showed virus spread to be spatially restricted, thus supporting the type of nearest neighbor dynamics discussed above. These biological data were used to construct a computational model of spatial virus spread in vivo, which gave rise to a variety of interesting insights. The model calculated the survival probability of any cell in the tumor. Importantly, the model suggested that relatively small changes in the initial density of infected cells as well as in the radius of the infection centers could have a major impact on the outcome. This clearly has practical implications for optimizing virus delivery. A strong sensitivity of outcome on the initial virus density has also been demonstrated in vitro, using experiments where cell cultures were infected with low infection multiplicities44.

Next level of complexity: Heterogeneity

A further issue that can arise, especially in a spatially restricted scenario, is that there is significant heterogeneity among parameters throughout the tumor. On some level, it can be argued that the model parameters describe average rates across the whole tumor cell population. However, it is possible that different spatial locations of the tumor are characterized by genetic differences in the tumor cell population. This can lead to locally different rates of tumor cell division, cell death rates, and viral replication kinetics, and can in principle have significant impacts on the outcome of the dynamics. The effect of such heterogeneity on the dynamics between a growing tumor and a replicating oncolytic virus has been studied with a mathematical model23. This work showed that several phenomena, such as tumor dormancy, occur in the context of heterogeneity, but did not occur in homogeneous models. This topic, however, remains under-explored, partly due to the fact that no systematic data exist that parses the effect of such heterogeneity. A systematic attempt to study this experimentally could involve cell cultures that contain different tumor cell lines in different parts of the culture that differ in susceptibility to infection. Related to this issue is the presence of physical barriers within the tumors, such as stroma cells, that are not susceptible to virus infection2. Exploring heterogeneity with a combination of experiments and mathematical models perhaps represents the next layer of complexity that can be incorporated into the frameworks reviewed here, once these very basic dynamics are understood in sufficient detail. Heterogeneity can be explored in the context of in vitro settings, which allows for systematic model testing and validation, and can be investigated further in mouse xenografts.


In this paper, different mathematical modeling approaches to describe the dynamics of oncolytic virus spread were discussed. Ordinary differential equations have been an important tool to gain biologically relevant basic insights, and to analyze particular treatment scenarios. The most basic model was discussed and some insights highlighted. It has become clear, however, that several model components remain uncertain and that the same biological assumptions can be described mathematically in different ways. This leads to the development of different models that all describe the same biology, but that differ in properties and predictions, thus casting doubt on issues related to model robustness. Axiomatic modeling approaches can address some of these difficulties, and it was discussed how “classes” of models with common properties can be identified and analyzed. Another problem related to the use of ordinary differential equations is that they assume well-mixed populations, which rarely applies to cancer. Hence, spatially explicit models were discussed, again concentrating on the simplest possible setting. Even in this simplest setting, however, the dynamics turned out to be highly complex and counter-intuitive, and an understanding of these dynamics is vital before further biological complexities can be modeled.

The fact that different model types, and indeed different mathematical formulations of the same model type can lead to different properties and to differing long-term predictions poses a dilemma for this area of research and requires thorough analysis of model robustness before researchers can come to solid biological conclusions. If the aim of mathematical modeling is to qualitatively interpret biological data and gain insights into possible outcomes under a defined set of biological assumptions, then the axiomatic modeling approach discussed here provides a good tool. This is less suitable if models are considered that are explored primarily by numerical simulations, such as spatially explicit agent-based models. In this case, it is important not to base the conclusion on a single formulation of this model, but to consider alternative computational implementations of the same biological assumptions in order to examine whether the results remain robust. If the goal of the mathematical modeling is to accurately describe and predict experimental data, it is equally essential to check for robustness. As discussed in this review, different model formulations can fit a given set of experimental data equally well, but can lead to different dynamics when extrapolating the dynamics beyond the time frame of the available data. Rather than using one particular model, different variations should be fit to the data and long term predictions should be compared. An important aspect in model validation would be to use the mathematical model that was parameterized by data fitting in order to make specific predictions under slightly different parameters, and then to perform corresponding experiments to investigate whether the model prediction is accurate. If the data can be predicted in this way without fitting the model to the new experimental data set, we can have higher confidence that the mathematical model indeed has predictive power. The finding that a single model formulation can be successfully fit to data is by itself not conclusive and does not mean that the model is predictive.


This work was funded in part by 1U01CA187956.


1. Sawyers C. Targeted cancer therapy. Nature. 2004;432:294–297. [PubMed]
2. Bell JC, Lichty B, Stojdl D. Getting oncolytic virus therapies off the ground. Cancer Cell. 2003;4:7–11. [PubMed]
3. Parato KA, Senger D, Forsyth PA, Bell JC. Recent progress in the battle between oncolytic viruses and tumours. Nature reviews Cancer. 2005;5:965–976. [PubMed]
4. Bell JC, Garson KA, Lichty BD, Stojdl DF. Oncolytic viruses: programmable tumour hunters. Current gene therapy. 2002;2:243–254. [PubMed]
5. McCormick F. Cancer-specific viruses and the development of ONYX-015. Cancer Biol Ther. 2003;2:S157–160. [PubMed]
6. McCormick F. Future prospects for oncolytic therapy. Oncogene. 2005;24:7817–7819. [PubMed]
7. Kaufman HL, Bines SD. OPTIM trial: a Phase III trial of an oncolytic herpes virus encoding GM-CSF for unresectable stage III or IV melanoma. Future oncology. 2010;6:941–949. [PubMed]
8. Lichty BD, Breitbach CJ, Stojdl DF, Bell JC. Going viral with cancer immunotherapy. Nature reviews Cancer. 2014;14:559–567. [PubMed]
9. Komarova N, Wodarz D. Modeling and simulation in science, engineering and technology. Springer; 2014.
10. Nowak MA, May RM. Virus dynamics Mathematical principles of immunology and virology. Oxford University Press; 2000.
11. Nowak MA, Bangham CR. Population dynamics of immune responses to persistent viruses. Science. 1996;272:74–79. [PubMed]
12. Perelson AS. Modelling viral and immune system dynamics. Nature Rev Immunol. 2002;2:28–36. [PubMed]
13. Perelson AS, Ribeiro RM. Modeling the within-host dynamics of HIV infection. BMC Biol. 2013;11:96. [PMC free article] [PubMed]
14. Rodriguez-Brenes IA, Komarova NL, Wodarz D. Tumor growth dynamics: insights into evolutionary processes. Trends in ecology & evolution. 2013;28:597–604. [PubMed]
15. Wodarz D. Viruses as antitumor weapons: defining conditions for tumor remission. Cancer Res. 2001;61:3501–3507. [PubMed]
16. Wodarz D. Gene Therapy for Killing p53-Negative Cancer Cells: Use of Replicating Versus Nonreplicating Agents. Hum Gene Ther. 2003;14:153–159. [PubMed]
17. Nowak MA, et al. Antigenic Oscillations and Shifting Immunodominance in Hiv-1 Infections. Nature. 1995;375:606–611. [PubMed]
18. Bajzer Z, Carr T, Josic K, Russell SJ, Dingli D. Modeling of cancer virotherapy with recombinant measles viruses. Journal of theoretical biology. 2008;252:109–122. [PubMed]
19. Biesecker M, Kimn JH, Lu H, Dingli D, Bajzer Z. Optimization of virotherapy for cancer. Bull Math Biol. 2010;72:469–489. [PubMed]
20. Dingli D, Cascino MD, Josic K, Russell SJ, Bajzer Z. Mathematical modeling of cancer radiovirotherapy. Mathematical biosciences. 2006;199:55–78. [PubMed]
21. Dingli D, et al. Dynamics of multiple myeloma tumor therapy with a recombinant measles virus. Cancer Gene Ther. 2009;16:873–882. [PMC free article] [PubMed]
22. Friedman A, Tian JP, Fulci G, Chiocca EA, Wang J. Glioma virotherapy: effects of innate immune suppression and increased viral replication capacity. Cancer Res. 2006;66:2314–2319. [PubMed]
23. Karev GP, Novozhilov AS, Koonin EV. Mathematical modeling of tumor therapy with oncolytic viruses: effects of parametric heterogeneity on cell dynamics. Biology direct. 2006;1:30. [PMC free article] [PubMed]
24. Komarova NL, Wodarz D. ODE models for oncolytic virus dynamics. Journal of theoretical biology. 2010;263:530–543. [PMC free article] [PubMed]
25. 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;1:6. [PMC free article] [PubMed]
26. Wein LM, Wu JT, Kirn DH. Validation and analysis of a mathematical model of a replication-competent oncolytic virus for cancer treatment: implications for virus design and delivery. Cancer Res. 2003;63:1317–1324. [PubMed]
27. Wodarz D. Computational approaches to study oncolytic virus therapy: insights and challenges. Gene therapy and molecular biology. 2004;8:137–146.
28. Wodarz D. Use of oncolytic viruses for the eradication of drug-resistant cancer cells. Journal of the Royal Society, Interface/the Royal Society. 2009;6:179–186. [PMC free article] [PubMed]
29. Wodarz D, Komarova N. Towards predictive computational models of oncolytic virus therapy: basis for experimental validation and model selection. PloS one. 2009;4:e4271. [PMC free article] [PubMed]
30. Bagheri N, Shiina M, Lauffenburger DA, Korn WM. A dynamical systems model for combinatorial cancer therapy enhances oncolytic adenovirus efficacy by MEK-inhibition. Plos Comput Biol. 2011;7:e1001085. [PMC free article] [PubMed]
31. Zurakowski R, Wodarz D. Model-driven approaches for in vitro combination therapy using ONYX-015 replicating oncolytic adenovirus. Journal of theoretical biology. 2007;245:1–8. [PMC free article] [PubMed]
32. Mok W, Stylianopoulos T, Boucher Y, Jain RK. Mathematical modeling of herpes simplex virus distribution in solid tumors: implications for cancer gene therapy. Clin Cancer Res. 2009;15:2352–2360. [PMC free article] [PubMed]
33. Paiva LR, Binny C, Ferreira SC, Jr, Martins ML. A multiscale mathematical model for oncolytic virotherapy. Cancer Res. 2009;69:1205–1211. [PubMed]
34. Reis CL, Pacheco JM, Ennis MK, Dingli D. In silico evolutionary dynamics of tumour virotherapy. Integr Biol (Camb) 2010;2:41–45. [PubMed]
35. Crivelli JJ, Foldes J, Kim PS, Wares JR. A mathematical model for cell cycle-specific cancer virotherapy. J Biol Dynam. 2012;6:104–120. [PubMed]
36. Kim PS, Crivelli JJ, Choi IK, Yun CO, Wares JR. Quantitative Impact of Immunomodulation Versus Oncolysis with Cytokine-Expressing Virus Therapeutics. Math Biosci Eng. 2015;12:841–858. [PubMed]
37. Camara BI, Mokrani H, Afenya E. Mathematical Modeling of Glioma Therapy Using Oncolytic Viruses. Math Biosci Eng. 2013;10:565–578. [PubMed]
38. Wodarz D, Chan CN, Trinite B, Komarova NL, Levy DN. On the laws of virus spread through cell populations. Journal of virology. 2014;88:13240–13248. [PMC free article] [PubMed]
39. Chahlavi A, Todo T, Martuza RL, Rabkin SD. Replication-competent herpes simplex virus vector G207 and cisplatin combination therapy for head and neck squamous cell carcinoma. Neoplasia. 1999;1:162–169. [PMC free article] [PubMed]
40. Rodriguez-Brenes IA, Komarova NL, Wodarz D. Evolutionary dynamics of feedback escape and the development of stem-cell-driven cancers. Proc Natl Acad Sci U S A. 2011;108:18983–18988. [PubMed]
41. Harrison D, et al. 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:1323–1332. [PubMed]
42. Wodarz D, et al. Complex spatial dynamics of oncolytic viruses in vitro: mathematical and experimental approaches. Plos Comput Biol. 2012;8:e1002547. [PMC free article] [PubMed]
43. Bailey K, et al. Mathematical model for radial expansion and conflation of intratumoral infectious centers predicts curative oncolytic virotherapy parameters. PloS one. 2013;8:e73759. [PMC free article] [PubMed]
44. Hofacre A, Wodarz D, Komarova NL, Fan H. Early infection and spread of a conditionally replicating adenovirus under conditions of plaque formation. Virology. 2012;423:89–96. [PubMed]