Fighting Cancer with Mathematics and Viruses.

After decades of research, oncolytic virotherapy has recently advanced to clinical application, and currently a multitude of novel agents and combination treatments are being evaluated for cancer therapy. Oncolytic agents preferentially replicate in tumor cells, inducing tumor cell lysis and complex antitumor effects, such as innate and adaptive immune responses and the destruction of tumor vasculature. With the availability of different vector platforms and the potential of both genetic engineering and combination regimens to enhance particular aspects of safety and efficacy, the identification of optimal treatments for patient subpopulations or even individual patients becomes a top priority. Mathematical modeling can provide support in this arena by making use of experimental and clinical data to generate hypotheses about the mechanisms underlying complex biology and, ultimately, predict optimal treatment protocols. Increasingly complex models can be applied to account for therapeutically relevant parameters such as components of the immune system. In this review, we describe current developments in oncolytic virotherapy and mathematical modeling to discuss the benefit of integrating different modeling approaches into biological and clinical experimentation. Conclusively, we propose a mutual combination of these research fields to increase the value of the preclinical development and the therapeutic efficacy of the resulting treatments.

1. Introduction

In vitro and in vivo experiments are model systems of the real world, and given the understanding of the limitations of such approaches, study results are to be interpreted as simplifications that may not hold truth in the real world, which for cancer biology is ultimately the patient. Whilst there is no substitute for wet lab experiments, for many biological problems it is not only unfeasible but actually impossible to test all possible experimental conditions, population heterogeneities or drug combinations for an infinite amount of time. Mathematical models may synergize with experiments in helping analyze experimental data to iteratively inform subsequent experiments with the highest likelihood of advancing our knowledge [1]. Mathematics is a formal language that can be used to describe a complex, non-linear system to, amongst others, evaluate the mechanistic underpinnings of dynamic behaviors and predict long-term behavior in silico. If fitted to experimental data, mathematical models could identify appropriate parameters to estimate contributions of different mechanisms to simulated emergent dynamics. From that, sensitivity analysis may be deployed to investigate which mechanisms truly drive experimental outcome to identify promising treatment targets for example. For cancer biology, it may be cell migration rather than cell proliferation that drives an aggressive phenotype [2]. The construction of a mathematical model with given assumptions may reveal useful information such as unknown rates that may be found during a parameter search. For example, if the proliferation rate is known for a cancer cell line, the death rate may be approximated by a tumor growth model by modeling a net growth rate parameter. The fitting of experimental data to a math model may either confirm negligible impact or the necessity for a death/exhaustion rate to be part of the model. Once the fitting is done, a model may then help interpret observations. Here, conditions of the math model may help validate the model where confounding or confirming characteristics coupled with the experimental data help to refine or affirm the model, respectively. Finally, a mathematical model may predict outcomes, albeit a complicated task. At this stage, predicted outcomes might be rejected due to impossible or improbable biological conditions necessary for the outcomes (much like an imaginary solution is not considered for an algebraic equation when used to solve a word problem that necessitates a real solution). Elaborate discussions on parameter fitting and biological implications on their virotherapy models can be found elsewhere [3,4]. It is of importance to understand that a mathematical model, often informed by limited experimental data, can never prove a biological hypothesis to be right—but it may prove a hypothesis false or support its plausibility.

A myriad of work has been done on mathematical modeling in immunology and virology. Alan Perelson and colleagues published a series of articles on optimal strategies in immunology [5,6,7]. Such models helped identify the dynamics and time course of HIV infection [8,9,10] and response to treatment [11]. Mathematical principles of immunology and virology are summarized in an excellent textbook [12]. In this review, we aim to describe and discuss the contribution that integrated mathematical modeling could make to the advancement of oncolytic virotherapy (OVT)—with emphasis on the immunological aspects of OVT.

More than 100 years passed between the first reported oncolytic effect of a virus [13] and the approval of virotherapy for cancer, first in China (H101 by Shanghai Sunway Biotech, 2005) then later by the U.S. Food and Drug Administration and European Medicines Agency (T-VEC by Amgen, 2015 and 2016, respectively). The first reported cancer remission in the context of viral infection was described in 1904 for a woman with myelogenous leukemia after being infected with influenza [13] (almost three decades before influenza was found to be a viral infection). Organized efforts to unlock the potential anticancer effects of viruses have been under way since the 1940–1950s [14,15,16,17]. In a 1949 clinical trial of Hepatitis B virus applied to Hodgkin’s lymphoma, Hoster, Zanes, and Von Hamm noticed that 7 of 22 patients improved in the clinical aspect of the disease; a reduction of tumor volume was observed in 4 of 22 patients [14]. Therapeutic effects were found to rely on natural selective replication of viruses in tumor cells. Moreover, mechanisms of action are being investigated to further develop OVT. Modes of effective OVT may include any individual or combination of cellular lysis, apoptosis, and innate or adaptive immune responses [18]. Various genetically engineered oncolytic viruses (OVs) are available, and the diversity of OVT combination regimens has grown rapidly in recent years, especially with emphasis on personalized cancer therapy. In the hope of guiding and optimizing these more complex therapies, mathematical models are deployed to understand the key mechanisms underlying the complex biological interactions. Having quantitative models in place that can simulate effective OVT could greatly reduce time and effort in the search for optimal OVT for subpopulations of cancer patients.

2. Mathematical Modeling of Tumor Growth

Oncolytic efficacy depends on tumor growth dynamics, of which mathematical modeling has a long history [19,20]. Differential equations are commonly used to describe the mechanisms that govern change in tumor cell number, or tumor volume: where the variable c represents the number of cells, dc denotes the change (difference, d) in cell number during a specific time interval, dt. The rate of change r describes the net change of the cell population. In the simplest case, r could be a constant number (such as r = 0.1), which describes the rate at which the number of existing cells increases. This will lead to exponential growth as in each time interval dt the number of cells increases by r × c. Intuitively, if r is a negative number, the population will decay. Experimental measurements such as cell numbers in vitro or volumes in vivo may help identify the value of such parameter r. Of note, r itself could be a function that describes in more detail how the number of cells could change, for example r = f − g, where f could describe the rate at which cells divide per unit time, and g the corresponding rate of cell death. The results of net growth is identical to only using the population net growth rate r, but the separation of cell division and cell death events may become very important in virology where viral dynamics are strongly influenced by those processes. Different growth laws have been developed and successfully fit to experimental and clinical tumor volume data. These growth laws can be as simple as above-discussed exponential growth. More complicated models may include increasingly complex biology. Initial exponential growth at low densities, when most cells have access to ample resources, decelerates when cells at the core of the tumor become growth-arrested. This is largely due to limited space and exhausted intratumoral nutrient supply as resources are consumed by cells closer to the tumor surface [21,22,23]. This established the notion of a tumor carrying capacity (K) as the maximum cell number that can be supported by a given environment [24]. The rate at which tumor growth saturates as the cell number approaches its carrying capacity can be shaped differently including linear (logistic growth; f = a × (1 − c/K)) and logarithmical (Gompertzian; f = a × ln(K/c)) functions, etc. A tumor carrying capacity may evolve with changing oxygen and nutrient supply through tissue vascularization, removal of metabolic waste products, and evasion of immune surveillance [25,26]. Then, carrying capacity itself will become a variable, i.e., dK/dt, whose rate of change can be described with a differential equation [25].

The most likely growth law for a specific tumor can be obtained by fitting the different models to experimental/clinical data and comparing the regression results [27]. The carrying capacity of the tumor may be more important as cancer cells multiply and the tumor grows; consequently, time constraints of the experiment or mathematical modeling may influence the decision to incorporate a tumor carrying capacity. Ordinary differential equation (ODE) systems are typically used for relatively simple models of growth dynamics. In a recent study, Murphy and colleagues showed that different tumor growth models fit retrospective experimental and clinical data equally well, but forward predictions in time may significantly vary [28]. There is a need to identify the most applicable growth dynamics [29], with explicit consideration of available data and the number of undetermined mathematical model parameters and their identifiability [30,31,32]. The Akaike information criterion (AIC) is often utilized to correlate model complexity with fit to data. This penalizes models with too many degrees of freedom and only marginal improvements in data fit [33].

While most examples in this review will be ODE systems, mathematical modeling of tumor growth and oncolytic virotherapy may also be described by partial differential equations (PDE) or agent-based model (ABM) systems. While ODE systems describe quantities over time, PDE systems describe quantities over time and relative dimension in space. A spatially explicit PDE system would become necessary when the biological effects are predominantly space-dependent (such as diffusion of cancer drugs or viral penetration into a tumor) and cannot be averaged in an ODE model (see de Pillis et al. [34]). Agent-based models (or single cell models; each cell is an individual “agent”) may be used to simulate individual cell behaviors and cell–cell or cell–environment interactions that propagate to produce emerging population-level dynamics. From such an approach, one can exclude or identify likely cellular mechanisms that underlie observed complex system-level behavior. See Wang et al. [35] for a recent review on ABM models of tumor cells and Wodarz [36] for a review on overall approaches to OV dynamics.

3. Viral Life Cycle

Viral entry into a host cell, replication within the cell and finally release of progeny particles is often referred to as the “life cycle” of a virus. Cellular requirements for the completion of the viral life cycle compare to the hallmarks of cancer [37]. Both processes benefit from pro-mitogenic, antiapoptotic and metabolic alterations promoting cell survival, proliferation and protein biosynthesis. Inflammation provides further stimuli. Induction of angiogenesis improves supply of nutrients and oxygen and furthermore allows spread of both viruses and tumor cells. The disruption of innate antiviral pathways, namely of the interferon (IFN) response, by mutations or virulence factors, is another important common mechanism of both viral infection and tumorigenesis. Viral infection and malignant transformation therefore share important signaling pathways and core elements required for successful progression (reviewed in detail in [38]).

Consequently, the first known oncogenes, amongst them the protein kinase-encoding genes in Rous sarcoma virus (v-Src [39,40]) and Abelson murine leukemia virus (v-Abl [41]), were found to be acquired from viruses. There are several ways in which viral infection can support malignant transformation. Both expression of oncogenic viral proteins, e.g., E6 and E7 in human papillomavirus (HPV)-associated cancers [42,43], and insertional mutagenesis—as observed in a gene therapy trial with retroviral vectors [44]—can provide a survival advantage to infected cells by promoting genetic instability, activating mitogenic signaling pathways and inhibiting apoptosis. Persistent inflammation as a result of viral infection also contributes to tumorigenesis in several cancer types such as hepatocellular carcinoma induced by chronic infection with hepatitis B or C viruses (reviewed in [45]).

On the other hand, anecdotic observations of cancer remissions associated with viral infections have been reported historically [13,46], indicating a potential for using these “culprits” as a strategy to cure cancer. This approach is also linked to the close resemblance of mechanisms required for viral replication and malignant transformation. Facilitation of productive viral infection by changes in cellular signaling, metabolism and innate immune responses in the course of malignant transformation has been referred to as “phenotypic complementation” [38]. Providing cellular requirements for productive infection can thereby enable viruses that are not able to replicate within healthy tissue to selectively infect and kill tumor cells.

Such viruses, which can either have a natural tropism towards cancer cells or be genetically engineered to enhance tumor-specific replication, are termed OVs and can be found throughout different virus classes [47]. In contrast to oncogenic viruses that cause latent infections allowing host cells to survive and accumulate mutations, OVs usually have a lytic replication cycle leading to the death of infected cells. Oncolytic virotherapy is a promising approach to treat cancer by making use of these agents, relying on a variety of mechanisms of action differing from those of conventional treatment options such as surgery, chemotherapy and radiotherapy.

A unique feature of OVT is the amplification of the agent within the tumor, increasing the therapeutic potential of the initially applied dose [38]. Killing of infected cells typically occurs by either extensive budding of viral progeny or expression of viral proteins on the cell surface and subsequent fusion with neighboring cells, both finally resulting in bursting of the host cell, i.e., oncolysis [48]. The unique mechanisms of OV cytotoxicity generally do not completely rely on intrinsic cell death programs, thus providing potential to overcome typical forms of treatment resistance observed in chemotherapy and radiotherapy [38]. Both viral amplification and the individual modes of oncolysis are important factors that require consideration in the mathematical modeling of OV infection.

The sequential steps in the reproductive cycle of a virus (Figure 1) can differ substantially between individual viruses and influence their rate of cell killing and spread, accordingly affecting their oncolytic potential. For a productive infection, one or more virus particles must enter the host cell. This process differs between classes of viruses, especially regarding the presence or absence of a viral membrane envelope. Attachment to cells is mediated via viral surface proteins targeting molecules accessible on cell membranes. Adaptation to cell entry receptors contributes to viral tropism. The host cell range can be broad, as in the case of vesicular stomatitis virus (VSV), which binds to low-density lipoprotein receptor (LDL-R) via the vesicular stomatitis virus G glycoprotein (VSV-G) [49], or more limited: wild type measles virus requires interaction of the hemagglutinin (H) protein with signaling lymphocyte activation molecule (signaling lymphocytic activation molecule (SLAM), cluster of differentiation (CD) 150) [50] or the adherens junction protein Nectin-4 [51] for entry and is restricted to immune and epithelial cells, which express these molecules, respectively. For some viruses, the host cell specificity can be modified by genetic engineering. Serial passaging of wild type measles virus has resulted in the generation of live attenuated vaccine strains [52]. These have adapted to usage of the complement-regulatory protein CD46 [53,54,55], which is frequently overexpressed on tumors, in addition to its natural tropism. Point mutations in the H gene can abrogate binding to these receptors, and introduction of transgenes encoding single chain antibodies or receptor ligands can be applied to retarget measles virus [56,57,58]. Pseudotyping with attachment proteins of different virus families is another tool for retargeting of viruses which can be applied e.g., to broaden the host cell range of lentiviral vectors via VSV-G (reviewed in [59]) or measles virus glycoproteins [60].

Figure 1

Schematic representation of the reproductive cycle of viruses. After successful attachment to a cell entry receptor (A), entry via membrane fusion or endocytosis (B), and uncoating of nucleic acids (C), viral gene expression (D1) and genome amplification (D2) are initiated. These processes can be complex and may include reverse transcription of the genome, shuttling to the nucleus, and further processing and modification of generated nucleic acids and proteins. Assembly (E) and subsequent release (F) of viral progeny complete the viral “life cycle”. Each of these steps contributes to (tumor) cell specificity of a particular virus as well as its replicative and cytolytic potential and, depending on the scientific question, may need to be considered in mathematical modeling.

Upon attachment, the viral particle passes the host cell membrane in a process termed penetration, involving different mechanisms depending on the virus type. Measles virus and most other members of the family Paramyxoviridae require interaction of two distinct surface proteins for entry. Binding of the H protein to a cell entry receptor induces a change in the structural conformation of the F (fusion) protein, leading to an approximation of viral particle and host cell. This finally results in membrane fusion, allowing entry of the viral particle into the cell (Plattet et al. [61], and Lamb and Parks [62]). Endocytosis provides an alternative route of cell entry. Acidification of maturing endosomes can induce conformation changes in viral attachment proteins to promote membrane approximation and fusion [63]. An example for this mechanism is the H protein of influenza viruses, which are also studied as potential oncolytic vectors [64,65,66]. Non-enveloped viruses omit the need for membrane fusion and rely on endocytic pathways for entry (reviewed in [67]).

In line with the heterogeneity of genome sizes among OV, ranging from approximately 5 kb of oncolytic parvovirus [68] to 300 kb of oncolytic vaccinia viruses [69], the number of genes and complexity of their regulation also differs greatly. For more detailed information on genome organization and replication of particular viruses, please refer to comprehensive reviews, e.g., for herpes simplex virus [70], poxvirus [71] and adenovirus [72]. Before amplification and expression of the viral genome can take place, it must be made accessible by uncoating. For many viruses, transfer of viral nucleic acid into the nucleus is furthermore required to make use of the genome amplification and transcription machinery of the host cell. Some viruses, including paramyxoviruses, which have their negative strand RNA genome packaged into so-called ribonucleoprotein complexes, harbor the enzymatic machinery necessary for transcription into mRNA [73].

Viral gene expression and genome amplification can be tightly regulated, e.g., by promoter elements and pre-messenger RNA (mRNA) processing. This enables restriction to certain cell types and adaption to changes in the phenotype of the host cell, such as the differentiation of HPV-infected epithelia that is necessary for the virus to complete its life cycle (see reviews [74,75,76]). The presence of immediate early, early and late genes in herpes viruses [77] relates to different phases of infection and ensures efficient replication before cell lysis. The time lapse between viral entry and progeny release, called “eclipse phase”, plays an important role in viral replication kinetics [78].

With regards to the safety aspects of oncolytic virotherapy, such post-entry mechanisms can be exploited to enhance tumor specificity or interfere with viral replication. Introduction of target sites for microRNAs (miRNAs) downregulated explicitly in malignant cells can prevent viral replication in normal tissue [79,80,81]. Exchanging furin cleavage sites necessary for cleavage and activation of viral proteins with target sequences for tumor-specific matrix metalloproteinases represents another strategy to prevent off-tumor toxicity [80,82]. Riboswitches, which undergo self-cleavage upon ligand administration, can be applied as off-switches for viral replication [83]. Tumor specificity can also be conferred at the level of protein translation, by inclusion of polyadenylation signals in the mRNA of key viral proteins [84,85]. Depending on the virus, DNA of viral origin might be integrated into the host cell genome, potentially inducing (epi-)genetic deregulation [86] and mutagenesis (reviewed in [87]). In contrast, RNA viruses are obligatory cytoplasmic and cannot integrate into the host cell genome, contributing to a favorable safety profile regarding their use as oncolytic agents.

Once the viral genome is replicated and structural proteins are expressed, assembly of viral progeny is initiated and finally mature particles leave the host cell via budding, pore formation and/or bursting of the cell. Some viruses, including measles virus, can infect neighboring cells without the necessity of forming extracellular particles by exploiting cell–cell contacts [88] and by display of viral attachment and fusion proteins on the host cell surface [89].

Oncolytic virus replication might be regarded in a very basic sense as releasing infectious particles via killing of infected cells, which makes it an appealing target for mathematical modeling. However, as viral infection is a highly complex, virus-specific, multistep process, an understanding of the underlying biology of virus–host interactions is required for careful consideration of appropriate modeling approaches. As a basis for a mathematical model of oncolytic virotherapy, a model of viral infection, uninfected cells becoming infected cells, is first considered.

4. Mathematical Modeling of Infection: Susceptible and Infected Model

Viruses are infectious agents that rely on a living host cell to replicate. Infectious disease modeling has an extensive history in mathematics to simulate viral spread and cytotoxic effects [90]. Host cells are divided into susceptible (uninfected, S) and infected (I) cells, where C (number of total tumor cells) = S + I. (Note that “C = S + I” does not include resistant cells without receptors for viral entry or stromal cells not targeted by OV.) Such SI models evolved from ecological population dynamics, e.g., food chain [91] (predator, prey, and top-predator) first used for infectious diseases (reviewed in [92]) then eventually for viruses as early as 1996 by Nowak and Bangham [93]. In 1995, Gatenby modeled cancer as a population competing with normal cells [94]. These early models have contributed to the evolution of the mathematical modeling of oncolytic virus therapy [36,95]. Nowak, Perelson, and others [12,93,96,97,98] used models that assume constant target cell production. Wodarz first described OV dynamics by adapting basic virus dynamics models, assuming a density-dependent proliferation of target (tumor) cells [99], as well as tumor-specific immune responses (virus-specific cytotoxic T lymphocytes) [99,100]. Others have expanded on this basis such as Dingli [101] and Biesecker [102].

Experimental and clinical observations often present data for the tumor volume, which can be converted into a total number of tumor cells, but typically are not of sufficient resolution to classify subpopulations of infected and non-infected tumor cells. It may be possible to capture viral infection dynamics with an implicit representation of virus particles [104] (Figure 2a), though more complicated models may require an explicit population of virus to be modeled [105] (Figure 2b). It should be noted that a model with an implicit virus representation may introduce an OV treatment by changing a quantity or proportion of uninfected cells to infected cells at the time of treatment (see f0 in supplemental information of [104]). Using an implicit representation limits the number of variables (and parameters) to the model system. However, this formulation is only applicable if the turnover rate of free virus particles is much faster compared to infected cells. Alternatively, the free virus population must be modeled explicitly as a time-dependent variable, which requires the inclusion of more parameters, some of which may be challenging to identify.

Figure 2

Simple susceptible (S) and infected (I) models with implicit (a) and explicit (b) model representations of a virus population. The β term (a) represents the population of dying infected cells. These infected cells release virus particles—not represented—upon death (also called burst [103]), hence, directly affecting the number of susceptible cells, which may become infected, γ(β)S. In explicit virus systems (b), the β term is commonly described as explicitly affecting the number of virus particles (δ) as opposed to the virus being modeled implicitly from its effects on the S and I populations. Consequently, the viral infection term is γVS. The term U(t) simulates the effect of oncolytic virotherapy, a number of virus particles administered at time t. Novozhilov et al. warns that mass-action, the assumption that the rate is proportional to the product of the variables, may only be acceptable when the S and I populations are close in value, S ≈ I [95].

A basic model for simulating viral infection has been previously discussed as a system of ODE [104]. In this simple SI model with implicit virus particles, the uninfected, susceptible tumor cells (variable S) may grow with rate αS (exponential tumor growth). With viral infection rate γ, uninfected cells become infected (variable I), and infected cells die with a fixed rate, β. A schematic of this simple model and the corresponding equations are shown in Figure 2a. Modeling the actual virus population introduces a new variable, V, with an arrival rate (first exposure, therapy, etc.) and intrinsic death. Infection is then explicitly modeled by the interaction of susceptible cells with virus particles, S × V, which occurs with rate γ.

5. Modes of Action in Oncolytic Virotherapy

In recent years, it has become apparent that direct tumor cell destruction via lytic replication is not the only mode of action contributing to the efficacy of oncolytic treatment (Figure 3), and might even be only a minor determinant of treatment success [106]. More importantly, oncolysis represents a form of immunogenic cell death, resulting in the release of potent immune-stimulatory molecules, including cytokines, pathogen-associated and damage-associated molecular patterns (PAMPs and DAMPs, respectively) and tumor-associated antigens (reviewed by Workenhe et al. [107]). This can induce recruitment of cells of both the innate and the adaptive immune system, potentially leading to systemic antitumor immunity. Tumor debulking by lytic replication can be seen as the initial step to reversing immune evasion mechanisms in the tumor microenvironment and evoking tumor-targeting processes [106]. Modes of action of tumor cell killing during the immune response have been summarized by Cassady and colleagues as cytokine-induced apoptosis, cytotoxicity of innate immune cells, and antigen-specific tumor cell lysis by T cells [18].

Figure 3

Multi-mechanistic modes of action of oncolytic virotherapy. Direct lytic replication (A) can be supported by additional antitumor processes in the context of oncolysis, e.g., immunogenic cell death and subsequent innate (B) and adaptive (C) antitumor immune responses, and endothelial targeting (D). This complexity can be exploited and expanded by genetic modifications and combination therapies. One approach is encoding prodrug convertases in the viral genome, which mediate local activation of systemically applied prodrugs (E). CTLs: Cytotoxic T lymphocytes; NK: Natural killer.

Viral infection is detected intracellularly via pattern recognition receptors (PRRs) such as retinoic acid-inducible gene 1 (RIG-1), which induces expression and secretion of pro-inflammatory cytokines and type I IFNs upon binding of viral RNA (reviewed by Barik et al. [108]). This typically leads to an antiviral state in surrounding cells, but can also result in killing of uninfected tumor cells by cytokines such as tumor necrosis factor α (TNFα) [109]. Innate immune effectors encompass natural killer (NK) cells, which are able to induce lysis of tumor cells or virally infected cells (reviewed in [110]), and phagocytic cells such as neutrophils and macrophages.

In the course of phagocytosis or tumor cell lysis, uptake of tumor-associated antigens (TAAs) by antigen-presenting cells (APCs) such as dendritic cells (DCs) is essential for the induction of an adaptive immune response [111,112]. Cross-presentation of such antigens on major histocompatibility complex I (MHC-I) molecules to CD8+ T cells in lymph nodes, followed by cytokine-mediated stimulation, is necessary for activation and expansion of tumor-specific cytotoxic T lymphocytes (CTLs). This is supported by differentiation of CD4+ T cells towards a type 1 T helper (Th1) phenotype, indicated by release of interleukin (IL)-2 and IFN-γ [113] (reviewed by Farrar et al. [114]). Another crucial aspect in terms of timing, efficacy and durability of antitumor immune responses is the development of memory and effector cell subsets [115,116].

Polarization of macrophages in the tumor microenvironment towards a pro-inflammatory (M1) rather than immunosuppressive phenotype (M2) further contributes to tumor destruction [117]. Inflammatory and other adjuvant stimuli provided by viral infection play a major role in promoting such antitumor effector functions. In their absence, innate sensing of tumor cells and TAA presentation will instead induce tolerance mechanisms, e.g., by regulatory T cells (Tregs) [118]. It has been shown that established TAA-tolerance of the immune system can be overcome by an intratumoral injection of a non-replicating, immunogenic adenovirus [119]. An OV-based combination regimen aiming at enhancing immunogenic cell death was also effective in breaking immune tolerance [120].

Based on this understanding, current development of novel oncolytic vectors focuses not only on enhancing lytic replication, but also on harnessing the immune response, for example by the introduction of transgenes encoding cytokines [121,122], checkpoint inhibitors [123], ligands of T cell co-stimulatory receptors [124], bispecific T cell engagers [125,126,127] or tumor antigens [128,129,130], respectively, into the viral backbone. Combination therapies represent another approach to support immune responses to OV treatment, including additional application of cytokines [131], immune checkpoint inhibitors [132,133] or chemotherapeutics [134].

Immune responses to OV treatment are highly complex and can also prevent successful therapy by limiting viral infection, depending on context and timing. For this reason, even immunosuppression might be beneficial prior to OV treatment to enhance viral replication and spread [135,136].

In addition to making use of the immune system, bystander killing of non-infected tumor cells can also be increased by using virus-encoded prodrug convertases [137]. Upon infection of tumor cells, the transgene is expressed at the tumor site. Non-toxic prodrugs are applied systematically and converted locally into a chemotherapeutically active compound, thereby minimizing toxicity to healthy tissue [138].

Targeting of tumor vasculature by oncolytic viruses has been observed [139,140] and was recently explained by suppression of cell-intrinsic antiviral mechanisms via vascular endothelial growth factor (VEGF) in the tumor microenvironment [141]. This could add to the antitumor potency of oncolytic virotherapy by inducing subsequent nutrient deprivation and hypoxia, but also impair viral delivery to the tumor site, warranting careful assessment prior to manipulating the antiendothelial potential of an oncolytic vector to enhance clinical benefit. Sunitinib-mediated inhibition of the VEGF receptor represents a potential approach for combination therapy [136,142].

There is emerging evidence for modes of action of OVT that go beyond the killing of infected cells by lytic replication. Models of viral spread and treatment outcome should appropriately account for these modes to enable accurate predictions. In a mathematical model, these modes first translate into the speed of growth functions, delays due to the viral life cycle, and mathematical descriptions of physical viral tropisms (manner of infection) as discussed in Section 6.

6. Modeling Specific Mechanisms of Action

Wodarz and Komarova have examined fast and slow classes of viral growth with biological interpretations of non-solid (“liquid”) tumors and solid tumors, respectively [104]. Fast viral growth is indicative of a well-mixed system, where there is little to no restriction on the viral infection rate. This represents some in vitro experiments as well as non-solid tumors. The slow class viral growth is necessary when modeling solid tumors due to spatial penetration dynamics. These fast and slow infection terms may be formulated in a variety of ways, which are discussed thoroughly in [104]. Examples of such viral spread terms include (fast) and (slow), where S and I are the number of susceptible uninfected and infected cells, and ε, ε1 and ε2 are empirical coefficients to describe saturation kinetics [143,144]. Wodarz and Komarova have stated that it can be difficult [145] to fit these models to experimental data [3,104,146], such that arbitrary viral expression terms may have no basis of any biological mechanism. That is, when fitting models to experimental data, in vitro experiments of a few days or in vivo experiments of about 30 days, these functions may all fit similarly well, and it is not until well after the time period of experimental data will models differ due to the different infection terms. An intratumoral injection of viral therapy into a solid tumor, for example, places a high concentration of virus particles in one specific location. Infection rate would be highest at the interface of this high virus concentration and the adjacent tumor cells, while almost no infection would take place in the tumor periphery. A spatially explicit partial differential equation (PDE) system may be employed for a more realistic depiction of spatial propagation of the virus after initial infection. An example of such a system is discussed in detail by Jacobsen et al. [147].

For modeling purposes, the time period from viral entry into tumor cell to tumor cell burst (“eclipse phase”) may be important [148]. In cell adhesion assays, for example, time periods that coincide with the viral replication cycle elapse between the time of oncolytic virus treatment of a cell line and observable effect [149]. Such delayed responses can be simulated mathematically by a delay differential equation (DDE) [150]. Further complexities may be added to a mathematical model to explicitly account for the interplay between multiplicities of viral infection and the antiviral states mediated by interferon [148,151,152] as a cellular response to viral infection. Of course, the model would require additional distinctions of antiviral and non-antiviral states [151] for uninfected cells (for example, see [90]).

The potential immune response to an oncolytic virus adds further complexity. In 2011, Eftimie et al. [4] reported a model based on an oncolytic immunotherapy study presented by Bridle et al. [129]. In this study, an adenovirus (Ad) vaccine was first used for immunization [153,154] before a VSV was used for oncolytic treatment against intracranial and systemic tumors (B16-F10 and CT26 cells, established in C57BL/6 and BALB/c mice, respectively). The Ad and VSV were designed to express the same TAA as the tumor cells. Bridle et al. [129] demonstrated that VSV increased a pre-existing antitumor immune response by shifting the immune response from viral antigens to tumor antigens. Further, a distinction between two compartments, lymphatic and peripheral (tumor) tissues, enabled simulation of a systemic effect of the OV on the immune system. Specifically, the recruitment and infiltration of effectors to the peripheral tissue was modeled due to an antiviral response in the lymphatic tissue. Mathematical analysis of model dynamics can identify a possible “tumor only” state (tumor without virus) when: (1) the inactivation rate of peripheral effectors is high; or (2) when the tumor is aggressive (high net tumor growth rate). While these observations may be biologically obvious, such analysis helps validate mathematical descriptions of the complex system as such conclusions are not built in a priori. More complex analyses help identify conditions under which the three equilibrium states of tumor-free, tumor only, and coexistence of tumor with virus are obtained. In this particular example, VSV must persist for a long time to achieve a tumor-free state. The model suggests, again quite intuitively and thus confirming model applicability, that oncolytic viruses with a higher half-life or better replication rate yield increased efficacy.

Bajzer and Dingli as well as Jacobsen and Pilyugin added a syncytia-forming fusion and budding mechanisms to lysis in a PDE model (Figure 4) that may be tailored to a particular viral mode of action [3,101,102,155]. Jacobsen and Pilyugin built upon previous PDE systems from Wu and Friedman et al. [156,157]. These models only allowed budding as a mechanism for viral particle production from syncytia, assuming that no apoptosis occurs from fused cells [158]. Recent reports suggest that this may not be the case [159,160,161,162,163,164,165], depending on the virus; adjustments during model development and fitting can resolve such issues. In their model, Jacobsen and Pilyugin found that an increase in burst size would allow for tumor control. However, different fusion (formation of syncytia) rates, , predict the outcomes of an inert virus or control of tumor growth . The authors hypothesized that at a low yields a governing lysis rate leading to control of the tumor, and at high , syncytia reportedly leads to fusion of all tumor cells and an exponential decay of tumor cells. Fusion rates between these two categories kill neighboring tumor cells but render remaining virus particles ineffective.

Figure 4

Schematic of general virus–host interactions. S: Susceptible cell; I: Infected cell; F: Fused cells; V: Virus. Each arrow can be described mathematically with a rate constant or more complex functions to describe different mechanisms of action (adapted from [147]).

7. Current Developments in Oncolytic Virotherapy

A total of 81 clinical trials for oncolytic viruses were listed on ClinicalTrials.gov [166] compared to only six studies 10 years ago, demonstrating the rapid development of the field. Compared to the first approaches to using viruses for cancer treatment, in some cases by application of infectious body fluids, research has made huge progress due to a deeper understanding of underlying virological processes and the possibility of genetic engineering (reviewed in [167]).

The most prominent example of the clinical translation of oncolytic viruses is talimogene laherparepvec (T-VEC, trade name ImlygicTM), a genetically modified herpes simplex virus type I. In addition to gene knockouts for enhanced tumor specificity and immunogenicity [168], T-VEC encodes the cytokine granulocyte-macrophage colony-stimulating factor (GM-CSF) for increased infiltration and activation of myeloid cells to support antitumor immunity [169]. As the first—and so far only—oncolytic therapeutic, T-VEC has been granted market approval in the United States in 2015, and shortly afterwards also in Europe, for the treatment of advanced melanoma. Investigations of its use in other malignancies and in combination therapies, especially with immune checkpoint inhibitors, are underway (trial examples for breast cancer: NCT02658812; pancreatic cancer: NCT03086642; combination with radiation: NCT02453191; combination with Nivolumab and Pembrolizumab, respectively: NCT02978625 and NCT02965716). The phase III study that led to approval resulted in a significantly enhanced durable response rate and a higher median overall survival in the patient arm treated with T-VEC compared to patients receiving GM-CSF [170]. Importantly, intratumoral injections of T-VEC have led to remissions of uninjected, distant lesions, indicative of the induction of systemic antitumor immune responses. CD8+ CTLs were shown to play a major role in mediating these effects. However, the exact mechanism of action of T-VEC has not yet been fully elucidated [171].

Another prominent example of OV is oncolytic vaccinia virus JX-594 (pexastimogene devacirepvec, Pexa-Vec), which has advanced to a phase III trial for the treatment of hepatocellular carcinoma (NCT02562755). JX-594 is thymidine kinase-deficient, increasing tumor cell specificity, and, like T-VEC, encodes GM-CSF for enhanced antitumor immune activation. In addition, a transgene encoding β-galactosidase allows for analysis of viral replication [172,173,174]. A phase II study in hepatocellular carcinoma showed an acceptable safety profile and dose-dependent survival benefit after intratumoral injection of JX-594 [175]. Viral replication and transgene expression were verified by measuring genome concentrations in blood and by detecting GM-CSF and antibodies against β-galactosidase, respectively [175]. Tumor responses were observed in both injected and uninjected lesions [175].

MV-NIS, an oncolytic measles virus derived from a live attenuated vaccine strain encoding human thyroidal sodium iodine symporter (NIS), entered clinical trials for application in several malignancies including multiple myeloma [176] and ovarian cancer [177]. The NIS transgene allows for imaging of infected cells via intracellular accumulation of 123I and can also be used for enhanced tumor cell killing by radioactive 131I [178,179]. An important aspect of oncolytic measles virotherapy remains the prevalence of neutralizing antibodies against measles virus surface glycoproteins, which cannot be sufficiently addressed by retargeting of the H gene [180]. Although a much-noticed case of durable complete remission of disseminated multiple myeloma has been observed upon systemic treatment with MV-NIS in a phase I study [181], this patient had unique characteristics favoring therapeutic success: a specific gene signature indicating sensitivity to OVT was detected by sequencing [182] while anti-measles virus antibody titers were absent in blood serum. However, there are means to overcome antibody neutralization of oncolytic agents: polymer coating [183], pseudotyping [184,185], complement inhibition [186] and usage of infected cells as virus carriers [187,188,189,190].

In addition, viruses that are non-pathogenic for humans may be used to reduce the prevalence of pre-existing immunity within the population. A genetically attenuated strain of the Maraba rhabdovirus [191], which had originally been isolated from Brazilian sand flies [192], is currently being tested in a prime-boost setting after immunization with an adenovirus (NCT02285816). Both viruses were genetically modified to encode the TAA MAGE A3 in order to evoke a strong tumor-specific immune response rather than boosting anti-viral responses [130,193]. Rodent parvovirus H-1 (H-1PV/ParvOryx) is another example of a non-human host-specific virus currently under clinical investigation as an oncolytic vector. Safety and improved antitumor effects were shown in a glioblastoma trial [194,195], and H-1PV has now entered phase II against pancreatic cancer (NCT02653313).

The amount and diversity of ongoing clinical trials for OVT is encouraging. Current developments in the filed especially highlight the importance of enhancing immunotherapeutic modes of action by means of vector engineering and combination treatment approaches. One such case follows in Section 8.

8. Mathematical Modeling of Oncolytic Virus Treatment with Immunotherapy

Immunotherapy can have different courses of action, including increasing the efficacy of antitumor immunity [196,197], reducing immune inhibitory mechanisms [198,199], or injection of dendritic cells [200] or engineered chimeric antigen receptor (CAR) T cells [201]. A variety of mathematical models have been developed to simulate different immunotherapies [202]. In 2015, Wares et al. modeled B16-F10 melanoma cells with adenovirus and dendritic cell injection in vivo [105]. A general scheme (Figure 5) is shown where dendritic cell injection as immune therapy (IT) is simulated as a direct increase in effector cells.

Figure 5

Explicit virus (V), and susceptible (S) and infected (I) cell model with addition of different immune system compartments (E, effector; IS, other immune system compartments such as antigen presenting cells) and immunotherapy (IT). Here, the effect of dendritic cell injection IT is modeled as an increase in effector cells. U(t) describes the time-dependent introduction of the virus population. Second-order terms and interactions (death/exhaustion) left out for clarity. Note that the viral infection term in [105]. See Figure 2b for comparison.

Appropriate terms and parameters were used including susceptible cell net growth rate, death/exhaustion terms, infection, injections of therapy at specified time points, infected cells (4-1BBL from Ad-ΔB7/IL-12/4-1BBL) activating effector cells, infected cells (IL-12 from Ad-ΔB7/IL-12/4-1BBL) stimulating recruitment of antigen presenting cells (IS). Some terms or portions of terms are functions; for example, the T cell killing rate is expressed as a base T cell killing rate (k0) that can be enhanced (linearly, c) by the presence of cytokine-producing infected cells (I), k(I) = k0 + c. The recruitments of T cells and APCs have a similar functional form of some constant, c, multiplied by the number of infected cells (I) responsible for recruitment, s(I) = c. Of note is that there is no explicit immune response modeled, just the interactions (recruitment) of certain immune components. Second, a conversion from tumor volume (mm3, experimental data) to number of cells is calculated by a conversion ratio of a volume of 1 mm3 to equal 1 million tumor cells. Lastly, all cell populations but uninfected (susceptible) tumor cells are initialized as being 0. Some model parameters are available in the literature including burst size [203], infected cell lysis rate [3,204,205,206], viral decay [207,208], T cell decay [209], maximum fractional T cell killing rate on tumor cells [210], and average time of an APCs to activate a T cell [211,212].

To calibrate unknown parameters with experimental data, a hierarchical approach was used. The control case was first simulated and fit to an exponential growth curve. The derived parameters were then kept constant while individually fitting for each of the individual therapies. Afterwards, those values were used to find parameters for the combination therapy. A Pearson’s r coefficient was calculated to give a rigorous fit, where values close to 1 indicates agreement. All models had a Pearson’s r greater than 0.98 except for the combination therapy model, which had a Pearson’s r = 0.92.

The Ad-ΔB7/IL-12/4-1BBL+DC model was used to predict alternative strategies: increasing dendritic cell dose, alternating Ad and DC injections, varying dose size, varying time between injections and number of doses. Within specifically posed biological and cytotoxic constraints, model simulations suggested that sequencing three injections of Ad followed by three injections of DC would be most effective to shrinking tumor volume. Further, the best order of Ad and DC injections and their respective optimal doses were highly dependent on temporal spacing. However, as appreciated by the authors, such predictions were likely to be model-dependent, emphasizing the need to simulate a variety of mathematical models with increasing and decreasing complexity to evaluate model predictive power and confidence in study conclusions. Model predictions that schedules that begin with viral injections would generally yield better decrease in tumor volume have been partially validated experimentally [213]. Interestingly, a recent robustness evaluation of the Wares model, using experimental data and a fitted mathematical model as inputs for Virtual Expansion of Populations for Analyzing Robustness of Therapies (VEPART) [214], also agrees with model predictions.

As the exact mechanisms and functional dependence of many of the different mechanisms in tumor–immune as well as virus–immune and virus–cancer interactions are yet to be fully deciphered, mathematical modeling in close dialog with experimental data may help identify likely and unlikely mechanistic properties of complex biological systems, which could and should iteratively inform subsequent experimental validation.

9. Discussion and Future Directions

Despite the many breakthroughs in oncolytic virus biology and mathematical modeling, further improvement in several arenas is essential for successful implementation of OVT in clinical reality. Such research fields include the generation of safer and more potent vectors by genetic engineering, the development of high-capacity production pipelines and good manufacturing practice (GMP)-compliant facilities [215], the selection of appropriate oncolytic agents for certain applications by “build[ing] on the strengths of individual virus platforms” [216], and the design of more potent combination therapies [217].

The ideal route of delivery is among the currently most discussed questions in the field and can be exemplified by T-VEC and MV-NIS, which stand for intratumoral injection and systemic administration, respectively [217]. The approval of T-VEC was a huge achievement for the field that might either represent a stepping stone or a barrier for the approval of further oncolytics, which will have to be compared to this benchmark [217].

In a rapidly advancing field where various options for genetic modifications and treatment combinations are becoming available, the necessity for making rational choices of regimen combinations and scheduling becomes apparent. This can be supported by mathematical modeling as an in silico prediction tool for preclinical, experimental treatments and analyses of patient data in clinical trials. A “discovery workflow” has been discussed describing a working relationship between experimental and computational methods [218]. Between genetic engineering and the resulting effects on the immune system, a myriad of combination therapies involving oncolytic virus therapy is inevitable (and ongoing in clinical trials [219]). Oncolytic virus modeling (OVM) is a tool with proven predictive results when data-driven, which has been used to describe viral interactions with immune cells in different tissues and immune cells combined with immunotherapy. This may be a key fundamental feature for future modeling including personal immune profiles when accompanied with immune surveillance data. An iterative process wherein data continuously refines all models and vice versa could tackle the complexities of OVM (Figure 6). While ODE systems may be robust enough to describe spatial interactions, a PDE system may be used to elaborate on several mechanisms of action, developing a better understanding for an oncolytic virus and a specific cancer.

Figure 6

A workflow between experimental and computational data describing an iterative process among in vitro, in vivo, in silico, and clinical models, adapted from [218].

Collaborations between computational and experimental scientists will be instrumental for accurate and predictive in silico modeling for enhancing pre-clinical and eventually clinical efficacies. Examples of successful model development and data fitting have been highlighted herein; satisfactory model predictions, however, remain very much elusive. With sparse experimental data with low biological resolution of spatial and temporal scales of virus and cell dynamics as well as their interactions, mathematical formulations of different action terms may fit the data equally well but their prediction for system behavior in the future or under different experimental conditions may vary widely [104,220]. The integration of mathematical modeling into experimental cancer virology may help design experiments to collect longitudinal, high-resolution biological data to further our understanding of this complex, non-linear biological system. The expertise of several disciplines—immunology, virology, and cancer biology, to name a few—will guide appropriate experiments and allow for rational but shorter steps towards successful oncolytic virus therapies. Additionally, personalized medicine may be a natural extension of these models. The established concept of distinguishing important immune cells to generate an accurate mathematical model is a hopeful glimpse toward a patient-based description of the immune compartment in predicting cancer therapy.