Exploring the Interactions of Oncolytic Viral Therapy and Immunotherapy of Anti-CTLA-4 for Malignant Melanoma Mice Model

Oncolytic ability to direct target and lyse tumor cells makes oncolytic virus therapy (OVT) a promising approach to treating cancer. Despite its therapeutic potential to stimulate anti-tumor immune responses, it also has immunosuppressive effects. The efficacy of OVTs as monotherapies can be enhanced by appropriate adjuvant therapy such as anti-CTLA-4. In this paper, we propose a mathematical model to explore the interactions of combined therapy of oncolytic viruses and a checkpoint inhibitor, anti-CTLA-4. The model incorporates both the susceptible and infected tumor populations, natural killer cell population, virus population, tumor-specific immune populations, virus-specific immune populations, tumor suppressive cytokine IFN-g, and the effect of immune checkpoint inhibitor CTLA-4. In particular, we distinguish the tumor-specific immune abilities of CD8+ T, NK cells, and CD4+ T cells and describe the destructive ability of cytokine on tumor cells as well as the inhibitory capacity of CTLA-4 on various components. Our model is validated through the experimental results. We also investigate various dosing strategies to improve treatment outcomes. Our study reveals that tumor killing rate by cytokines, cytokine decay rate, and tumor growth rate play important roles on both the OVT monotherapy and the combination therapy. Moreover, parameters related to CD8+ T cell killing have a large impact on treatment outcomes with OVT alone, whereas parameters associated with IFN-g strongly influence treatment responses for the combined therapy. We also found that virus killing by NK cells may halt the desired spread of OVs and enhance the probability of tumor escape during the treatment. Our study reveals that it is the activation of host anti-tumor immune system responses rather than its direct destruction of the tumor cells plays a major biological function of the combined therapy.


Introduction
Oncolytic virotherapy has shown promising anti-tumoral effects in numerous basic and clinical research. Besides the oncolytic ability to direct target and lyse tumor cells, oncolytic viruses (OVs) also can be designed to selectively replicate to kill tumor cells while preventing their binding to and replication in healthy normal cells through engineered mutations [1]. Nevertheless, infectious virus particles increase with time within the tumor, which may prevent patients' exposing to substantial excess risk during drug delivery [2] and the tumor-selective replication of viruses makes cancer treatment less toxic than standard chemotherapy drugs [3]. Because OVs can be made in the laboratory by modifying their genome, this allows therapeutic exploitation of viruses encoded with transgene proteins to against T cell inhibitory factors such as anti-CTLA-4. Clinical evidence has showed that combining oncolytic viruses with checkpoint inhibitors CTLA-4 in various cancer types, including prostate, and melanoma xenografts, resulted in upregulation of intratumoral T cells and boosting T cell immunity to solid tumors [4,5]. These molecular biotechnology modifications harness the immune system to eliminate tumor cells and provide a new treatment avenue for durable and effective clinical responses in cancer patients [6].
An activation of CD8 + T cells requires two signals. The first signal occurs when a naive CD8 + T cell encounters and interacts with an antigen presenting cell such as a dendritic cell through the T cell receptor [7]. The second signal, the co-stimulation signal, is provided by the interaction between CD28 on the membrane of T cells and B7 on antigen presenting cells [7,8]. However, CTLA-4 (cytotoxic T lymphocyte antigen-4) molecules are expressed by activated T cells such as CD4 + T cells and regular T cells and can out-compete CD28 for binding to B7 and thus dampen the co-stimulatory signals. In 2011, the U.S. Food and Drug Administration approved ipilimumab, an anti-CTLA-4 antibody, for melanoma to increase the mobility of T cells and allow the cytotoxic T cells to continue to destroy cancer cells [9][10][11][12]. Unfortunately, not all patients responded and the long-term treatment efficacy was not satisfactory in solid tumors [13]. Combining the OVT with CTLA-4 blockade could complement the effects of poor tumor targeting through the selective replication ability of OVs in tumor cells and thus dramatically boost anti-tumor therapeutic efficacy [14].
There are several mathematical systems modeling the interactions between OVs and tumor cells along with various immune responses. For example, Mahasa et al. used a model of delay differential equations to describe interactions of normal and tumor cells by considering adaptive immune responses followed by an early viral propagation period. They assumed that OVs infect both host and cancer cell populations and subsequently induce anti-viral immune responses [15]. Eftimie et al. [16] built a mathematical model to study the anti-tumor effect of adenovirus and oncolytic vesicular stomatitis virus, and their interactions with the CD8 + T related immune cells. Their study indicated that cancer treatment could be improved with different types of oncolytic viruses. Senekal et al. built a mathematical model to investigate the dynamical interactions of OV-induced NK cell recruitment during oncolytic virotherapy. Their study suggested that NK response simulated by OV is more efficient at reducing the infected tumor cell population than stimulated by the presence of tumor only. They also suggested other long-term effector cells such as CD8 + T cells should be included in the OV model building [17]. Storey and Jackson developed a spatially explicit hybrid cellular automaton and partial differential equations with a combination of an oncolytic viral therapy and an anti-PD-1 immunotherapy to investigate the influence of spatial location of viral doses and to determine optimal viral dosing to therapeutic efficacy. Their results indicated that the tumor antigenicity level plays a more important role for treatment efficacy than the T cell killing rate [18]. Most of the oncolytic viral therapy models including checkpoint pathway address the discussion of the PD-1/PD-L1 checkpoint in which the immune responses are conveyed by anti-tumor CD8 + T cells.
Deviated from the above mentioned approaches, our model studies the effectiveness of oncolytic virotherapy together with immune checkpoint modulators CTLA-4. We include the effect of cytokines, which is seldom addressed in other models. Although there are more than ten types of immune cells that are known to play a vital role in oncolytic virotherapy, for simplicity, we only consider essential to therapeutic effect of OV such as CTL and NK cells which are addressed by immune cell depletion studies [19]. Moreover, our model construction is based on the population level interactions and does not account for subcellular events and stimulatory pathways. Since CTLA-4 molecules are expressed mostly by activated CD4 + T cells and little by CD8 + T cells, and CD4 + T cells are stimulated by cytokine IFN-γ, their influences are considered in the immune responses. The CTLA-4 is considered as a negative regulator of CD4 + T cell activation, and, as a consequence, can prevent cytokine production. Unlike CD4 + T cells, the anti-CTLA-4 does not have a significantly influence on the proliferation of CD8 + T cells. Chan et al. shows that checkpoint-blocking activity of anti-CTLA-4 on CD4 + T cells permits greater production in CD4 + T cells than in CD8 + T cells [20]. Therefore, the population of CD8 + T cells is not considered to be influenced by the CTLA-4 in our model. Specifically, we distinguish the tumor-specific immune abilities of CD8 + and CD4 + T cells: effector CD8 + T cells directly kill susceptible and infected tumor cells while CD4 + T cells kill tumor cells indirectly through cytokines. T cells are also recruited by immune cells, which become activated when they encounter viruses. The innate immune NK cell population is assumed to be pre-existing within the tumor vicinity and is recruited to the TME and mediates initial OV clearance as in [17]. They are activated due to immunogenic cell death of infected tumor cells and cytokines. Cytokines activated by both susceptible and infected tumor cells and secreted by activated T cells and NK cells conduct the indirect killing. In addition to the tumor-specific immune response, we also take into the account of virus-specific immune response due to the fact that the presence of virons on infected tumor cells also activates anti-viral immune responses. They serve to constrain viral infections by viral lysis free OVs. Our proposed model is supported by the murine experiment of Engeland et al. [5] in which the treatment of murine model of malignant melanoma B16-CD20 using OVT and the immune checkpoint inhibitor CTLA-4 are conducted [5].
The remainder of this paper is organized as follows: In Section 2, we describe our mathematical model incorporating the combined treatment of OVT and anti-CTLA-4. We describe in detail the assumptions considered for each equation in the full model derived, and we parameterize the model parameters. Numerical explorations are provided in Section 3. In particular, Section 3.1 presents model validation and Section 3.2 investigates various treatment protocols for improved therapeutic outcomes. Section 3.3 performs global sensitivity analysis using either the monotherapy of OVT or the adjuvant therapy with anti-CTLA4. The Section 4 presents a brief review and discussion.

Model Development
We formulate a simple ODE model to investigate the treatment of murine model of malignant melanoma B16-CD2 using OVT and anti-CTLA-4 [5]. The goals of the proposed model are to predict: combining oncolytic virotherapy with immune checkpoint modulators would reduce tumor burden by direct cell lysis and stimulating anti-tumor immunity as described in [5]. The model is constructed to describe the interactions between tumor populations, virus population, tumor-specific CTL response, NK cells, virus-specific immune populations, tumor suppressive cytokine IFN-γ, and the effect of immune checkpoint inhibitor CTLA-4. We consider separately the average temporal changes in the uninfected tumor size (T u (t)) and infected tumor size (T i (t)). For the immune response, we model the evolution of effector cells CD8 + T (X(t)) and CD4 + T (Y(t)) cells, the population of co-inhibitory receptor CTL-associated antigen 4 (CTLA-4) molecules (W(t)) expressed by the CD4 + T cells and CD8 + T cells. The specific or non-specific attack of NK cells during oncolytic virotherapy is specified by the variable (N(t)). To model the evolution of CD4 + T cells, we focus on the effector role of CD4 + T cells. In fact, CD4 + T cells can kill cancer cells through cytokines and chemokines they produced even in the absence of CD8 + T cells and NK cells [21][22][23][24]. Since the level of tumor suppressive IFN-γ is investigated in mice treated with OVs and aCTLA-4 in [5] and tumor eradication or recrudesced after initial regression was found in mice lacking of IFN-γ recipient [25], the effect of IFN-γ (C(t)) is considered within the model. We model the virus population (V(t)) and the total number of anti-viral immune cells (Z(t)) to describe the immune responses to viral infection. The time unit is a day. The cell population has the unit of number of cells while oncolytic virions have the unit of PFU. The cytokines and immune checkpoint have the units of pg/mL and number of molecules, respectively. Prior to the model introduction, we describe the assumptions considered for each equation in the full model. These assumptions are depicted schematically in Figure 1. The summary of model state variables with their respective definitions are presented in Table 1.

Variable Description
T u (t) total number of susceptible (uninfected) tumor cell population T i (t) total number of infected tumor cell population V(t) total number of oncolytic virions X(t) total number of anti-tumor immune cells (CD8 + T cells) total number of anti-tumor immune cells (CD4 + T cells) total number of CLTA-4 protein population Z(t) total number of anti-viral immune cells N(t) total number of natural killer cells • Equation (1) models the uninfected tumor cell population. We assume susceptible tumor cells grow in a logistic fashion with an intrinsic growth rate r u and carrying capacity K t for all tumor cells. We choose logistic growth model for the susceptible tumor population because evidence shows that the growth of tumor is slower when the tumor becomes larger [15,26]. The term β t T u V m v +T u is the viral infection rate of uninfected tumor cells with the maximum rate β t , which turns a susceptible tumor cell into an infected tumor cell. The saturated form of the tumor-virus interaction is considered according to [26]. The killings of tumor cells are modeled by the terms δ x X m x +X T u and δ c T u m t +T u C due to CD8 + T cells and cytokines, respectively. They represent susceptible tumor population killed directly by the effector CD8 + T cells and indirectly by CD4 + T cells through cytokines at a rate of δ x and δ c , respectively. The Michaelis-Menten kinetics form of the tumor and the tumor-specific immune cell interaction is adopted according to [15] and [27]. The saturated form shows the limited property of effector cells and cytokines abilities to lyse tumor cells. NK cells are able to recognize tumor antigens and indiscriminately kill tumor cells. The term d u NT u is used to model the tumor killing by NK cells at a rate d u .
• Equation (2) models the population of infected tumor cells. Similar to susceptible tumor cells, we assume infected tumor cells grow in a logistic fashion with an intrinsic growth rate r i and carrying capacity K t for all tumor cells. Susceptible tumor cells infected by oncolytic virus at a rate β t result in the increase of the infected tumor cells. This cell population dies at a rate a t . Since T cells kill not only virus-free tumor cells but also virus-infected tumor cells, the infected tumor cells are lysed via anti-tumor adaptive immune cells CD8 + T cells (directly killing) with δ x X m x +X T i and by CD4 + T cells (indirectly killing through cytokines) with δ c T i m t +T i C. The killing due to anti-viral adaptive immune cells is modeled by the term δ z ZT i at a rate δ z [15]. The term d i NT i is used to address the killing by NK cells on infected tumor cells at a constant rate d i . • Equation (3) models the virus population. New virus particles are produced upon clearance of an infected cancer cell. The parameter b t is the burst size of viruses released from an infected tumor cell. δ v is the viral lysis by anti-viral immune cells.
The virus decays at a rate γ v and s represents the dose of OV. NK cells not only directly recognize and kill viral-infected cells through their receptors but also provide an antigen-specific adaptive response to viral infections, which represents the first line of defense and a rapid immune response against viral infections [28]. The term d v NV models virus killing of free OV by NK cells at a rate d v . • Equation (4) models tumor-specific CD8 + T cells related adaptive immune responses. The term a x T i +T u h x +T i +T u X models proliferation of CD8 + T cells. Here, a Michaelis-Menten term is used to denote that the anti-tumor immune response is induced by tumor antigens presented on both uninfected and infected tumor cells [15], where h x is the half-saturation constant and a x is the proliferation rate of anti-tumor adaptive immune cells. Since anti-tumor immune response is activated by oncolytic viruses to fight tumor cells, the term a v NV denotes the rate for which T cells are recruited by immune cells through the interactions with viruses at a rate a v . This addresses the T cell activation following encounter with the viruses, rather than by encountering with infected cells [29]. γ x is the natural death rate of anti-tumor T cells. • Equation (5) models tumor-specific CD4 + T cells related adaptive immune responses (e.g., cytokine related immune responses). The evolution of CD4 + T cells takes similar proliferation formation as that of Equation (4) but with different half-saturation constant h y . Since CD4 + T cells are released mainly by IFN-γ and anti-CTLA-4 amplifies CD4 + T cell activation [27], the proliferation of CD4 + T is modeled in the form of a y where a y is the proliferation rate of Th cells. The inhibition of CTLA-4 on CD4 + T cells is modeled by the term (1 + νW) with a measure of inhibition ν. Parameter γ y is the apoptosis rate of Th cells. • Equation (6) models the time evolution of major cytokine IFN-γ. T cells and NK cells produce several cytokines and chemokines that coordinate various immune responses and are the major source of IFN-γ. The terms α x 1+b y W Y, and α n 1+b n W N represent that IFN-γ is activated by both susceptible and infected tumor cells, and secreted by activated CD8 + T [20,30], CD4 + T cells [20,25], and NK cells [31] at a constant rate α x , α y , and α n , respectively. The terms 1 + b x W, 1 + b y W, and 1 + b n W model CTLA-4 engagement that prevents cytokine production and the fact that anti-CTLA-4 therapy results a significant amount of IFN-γ [25,32]. The parameter γ c is the natural degradation rate of IFN-γ. • Equation (7) models soluble proteins CTLA-4. The CTLA-4 molecules are expressed on activated T cells such as effector T cells and regulatory T cells and can out-compete CD28 for binding to B7 and thus dampen the co-stimulatory signals. The CTLA-4 expression rate on a single CD4 + T cell is assumed a constant and is denoted by the parameter r y . Evidence has shown that CTLA-4 is also expressed on CD8 + T cells [20]. The CTLA-4 expression rate on a single CD8 + T cell is denoted by the parameter r x . The natural lost rate of CTLA-4 is denoted by γ w . The blockade rate of immune checkpoint CTLA-4 is presented by the parameter u. • Equation (8) models the virus-specific immune response. The equation describes immune responses to viral infection. The parameter p v is the virus-specific proliferate rate of anti-viral immune cells which become activated due to debris or viral antigens on infected cells [15]. The natural death rate of anti-viral immune cells is denoted by γ z . • Equation (9) models natural killer (NK) cells. NK cells are known to be non-specific and attack "non-self" cells, but NK cells can be specific or non-specific during oncolytic virotherapy [17]. There is a pre-existing NK cell population (s N ) within the tumor vicinity and is recruited to the TME and mediates initial OV clearance [17]. In this study, we assume activation of NK cells is dependent on the contact with OV infected tumor cells, as experimentally observed in Leung et al. [33]. The second term,

represents the stimulation and recruitment of NK cells. A saturation term
T i T i +m n is used to describe the limited effects of NK response by infected tumor cells.
We also assume that NK cell response is further enhanced by lysis, which induces immunogenic cell death (ICD) of infected tumor cells at the rate r n until it attains the maximum capacity K n . Parameter r n is the recruitment/proliferation of NK cells in response to danger signals (such as DAMPs and PAMPs) released during ICD of OV-infected tumor cells. Parameter m n represents the half-saturation constant of NK cells that supports the maximum killing of tumor cells by NK cells. Proliferation of NK cells is also stimulated by IFN-γ [34]. Parameter h n is used to indicate the saturated effects of immune response, and parameter ζ is the stimulation and recruitment of NK cells by IFN-γ. The interactions of NK cells with tumor cells can lead to inactivation of NK cells at a constant rate, δ n that is proportional to their interaction term. In fact, −δ n (T u + T i )N denotes the inactivation of NK cells upon their interactions with tumor cells at the rate δ n . Finally, the last term, γ n N, denotes the natural death rate of NK cells. The formulation of equation of NK cells is derived from [17,34].
The detailed description of model parameters is summarized in Table 2, and the model is given in Equations (1)-(10). Death rate of effector cells ν The measure of CTLA-4 blocking rate on CD4 + T cells It is clear that solutions of (1)-(9) exist and remain nonnegative on [0, ∞) so that the model is biologically feasible. When there is no OVT, s = 0, (1)-(9) always has a unique tumor-free equilibrium E 0 = (0, 0, 0, 0, 0, 0, 0,N), whereN = s n /γ n . It can be shown that E 0 is locally asymptotically stable if r u < d uN and r i < a t . See Appendix A. This indicates that the tumor can be eradicated if it is small and has small growth rates.

Parameter Estimation
In this section, we estimate parameter values from literature. Susceptible tumor cells. Susceptible tumor cells grow logistically with intrinsic growth rate r u = 0.924 day −1 and carrying capacity K t = 3.3 × 10 9 cells. These values are taken from [35]. An uninfected tumor cell becomes infected after infection from oncolytic viruses. The viral infection rate β t = 0.0038 (cells)(PFU −1 )(day −1 ) and the half-saturation constant m v = 1 cells are adopted from [26]. An uninfected tumor cell can be killed by either effector cells or cytokines. The tumor killing rate δ c = 0.2 (cells)(day −1 )(pg/mL) −1 by cytokines with half-saturation constant m t = 10 5 cells are taken from [36]. The halfsaturation constant for the killing by effector cells is m x = 10 3 cells. The lysis rate of tumor cells (infected and uninfected) by immune cells is δ x = 2 day −1 . Both values of m x and δ x are from [26]. The tumor (infected and uninfected) killing rates due to NK cells d u and d i are the same with 8.68 × 10 −10 (days −1 )(cells −1 ) and taken from [15,17].
Infected tumor cells. Infected tumor cells grow logistically with intrinsic growth rate r i = 0.924 day −1 and carrying capacity K t = 3.3 × 10 9 cells. It has a disease related death rate a t . The range of a t is 0.5 − 2.6667 (cell −1 )(day −1 ) rescaled from [29]. We set a t = 1. This value is also adopted from [15] and [16]. In addition to the killings by CD8 + T cells and cytokines, infected tumor cells can be killed by anti-viral immune cells, and this lysis rate is δ z = 1 (cell −1 )(day −1 ) according to [15].
Tumor-suppressive cytokines (IFN-γ). Both CD8 and CD4 Th1 effector T cells are the primary sources of IFN-γ. CD8 + T cells produce copious amounts of IFN-γ in response to activation. The IFN-γ production in CD8 + T cells is central to the generation of Th1 immune responses through NFAT1 protein [39]. The IFN-γ production rate by CD8 + T cells is α x = 9 (pg/mL)(day −1 )(cell −1 )(cell −1 ). The IFN-γ production rate by CD4 + T cells is α y = 9 (pg/mL)(day −1 )(cell −1 )(cell −1 ), while the loss rate of tumor-suppressing cytokines is γ c = 34 day −1 . These values follow from [27,40]. NK cells are also the major source of IFN-γ. The IFN-γ production rate by NK cells is estimated as α n = 0.4 (pg/mL)(day −1 )(cell −1 )(cell −1 ). For different subjects, the parameter range for the measure of inhibition varies from 10 −3 to 1. We estimate the measure of inhibition of CTLA-4 to CD8 + T (b x ), CD4 + T cells (b y ), and NK cells (b n ) based on the model validation. The measure b x of CTLA-4-mediated inhibition on IFN-γ produced by CD8 + T cells has the value b x = 10 −3 molecule −1 . The measure b y by CD4 + T cells has b y = 10 −3 molecule −1 , and the measure b n by NK cells has b n = 10 −3 molecule −1 .
Immune checkpoint CTLA-4. The CTLA-4 express rate on a single CD4 + T cell is assumed to be a constant with r y = 5000 (molecules)(day −1 )(cell −1 ). The degradation rate of CTLA-4 has γ w = 8.3178 day −1 . They are from [27]. It is estimated using the assumption of exponential decay for the degradation of CTLA-4. The CTLA-4 express rate on a CD4 + T cell is estimated in the range of 2500-5000 (molecules)(day −1 )(cells −1 ) under the linear growth assumption. This estimation is due to the fact that the maximum number of CTLA-4 molecular is 10 4 per cell and the maximal protein CTLA-4 expression is being reported at 48-96 hours post-activation in Jaffe [41]. The CTLA-4 express rate on a single CD8 + T cell is estimated as r x = 800 (molecules)(day −1 )(cell −1 ). The estimation is based on the fact that the relative expression of CTLA-4 is higher in CD4 + T cells compared with CD8 + T cells [20]. The inhibition rate, u (dimensionless), of immune cells by checkpoint CTLA-4 is from 0 to 1.
Virus-specific immune cells. The main function of viral-specific immune cells is to kill viruses from the therapy. The birth rate of anti-viral immune cells in response to the presence of viral particles on the surface of infected cancer cells is estimated as p v = 0.6 day −1 with a death rate γ z = 0.13296 day −1 . They are from [15]. The feasible range of the infected cell-mediated proliferation rate due to anti-viral immune response is 0.6-2.5 day −1 from [16].
Natural killer cells. The constant influx of NK cell is s n = 3.2 × 10 3 (cells)(day −1 ), the half-saturation constant of infected tumor cells is m n = 10 4 cells, the recruitment rate of NK cells via ICD by infected cells is r n = 10 −5 day −1 , the carrying capacity for NK cell population is given by K n = 6.63 × 10 10 cells, the inactivation rate of NK cells by tumor cells is δ n = 10 −7 (cell −1 )(day −1 ), and the natural death rate of NK cells is γ n = 4.12 × 10 −2 day −1 . These values are adopted from [17]. The stimulation and recruitment rate of NK cells by IFN−γ is ζ = 0.5 days −1 which is estimated from [42]. The half-saturation constant for NK cell population activated by cytokines is h n = 3 × 10 2 (pg/mL) −1 followed from [34]. Table 3 summarizes model parameter values and value ranges, and values used in simulations.

Numerical Simulations
Since the time period for majority of the experimental studies on tumor-immune system interactions are in the order of few days/weeks after initial injections of oncolytic virus therapy, our study of model (1)-(9) begins by investigating transient behavior of the system as we vary the dosages of immunotherapy and immune responses. We validate the model according to the experimental data of Engeland et al. [5]. In [5], tumor sizes in mice during combined therapy with anti-CTLA4 were recorded. In particular, malignant melanoma B16-CD20 cells were injected subcutaneously into the flank of C57BL/6 mice. When tumors had reached an average volume of 40 mm 3 (4 × 10 7 number of cells), mice were treated with carrier fluid treatment (mock; control) or injected with 2 × 10 6 viruses (OV) in the tumor, or injected with viral particles together with anti-CTLA-4 (OV-aCTLA-4) for 5 consecutive days [5]. They revealed that the combination therapy had a synergistic effect, which enhances host anti-tumor immunity and increases the efficacy of OV treatment alone. Evolutions of tumor volumes for three different treatments on day 18 after implantation are shown in Figure 2c of Engeland et al. [5]. On day 18, the tumor volume reached 1860 mm 3 (1.86 × 10 9 number of cells), 1140 mm 3 (1.14 × 10 9 number of cells), and 400 mm 3 (4 × 10 8 number of cells) approximately for mice treated with mock, OVs, and OV-aCTLA-4, respectively. The tumor volumes for mice treated with OV-aCTLA-4 are significantly lower than the controlled mice of no treatment. In our numerical simulations, tumor volumes are translated to the number of tumor cells through the relation of 1 cm 3 = 10 9 number of cells [46]. The simulations are performed by MATLAB ode15s ODE solver. Treatments with 5 consecutive days are represented by a continuous function. The treatment of oncolytic virus therapy is represented by the parameter s, and the blocking rate of CTLA-4 is represented by the parameter u.

Model Validation
We start our numerical investigation on the dynamics of the model given in Section 2.1.1 by showing that our model output is in reasonable agreement with the data given in Figure 2c of Engeland et al. [5]. For preliminary mathematical analysis of the model, see Appendix A. According to the mice experiment of Engeland et al. [5], the treatment started on day 8 post-implantation and the tumor reached an average number of 4 × 10 7 cells. Therefore, the initial condition for the susceptible tumor is T u (0) = 4 × 10 7 . Following the injection of tumor cells, it is reasonable to assume that the tumor activated primary immune responses. The primary immune responses are represented by the initial number of CD8 + T cells, CD4 + T cells, and NK cells. For simplicity, the initial settings of anti-tumor immune cells of CD8 + T and CD4 + T cells are assumed to be the same with X(0) = Y(0) = 250. All implementations last for 5 consecutive days. Initial settings of other variables are T i (0) = V(0) = C(0) = W(0) = Z(0) = 0 and N(0) = 10 4 from [17].
Validation starts from the case of no treatment. Our simulation results show that the tumor grows up to 1.8630 × 10 9 on day 18 without any treatment. Without any treatment, the tumor size is approximately 1.8 × 10 9 in Engeland et al. on day 18 [5]. For treatment with OV alone, our numerical simulation shows that, if 2 × 10 6 oncolytic viruses are injected into the mice when the tumor reaches an average volume of 4 × 10 7 cells, the tumor reaches 1.1441 × 10 9 cells on day 18 after 5 days of treatment on day 8 as shown in Figure 2. Notice that the average tumor size is approximately 1.14 × 10 9 cells on day 18 in Engeland et al. if monotherapy OV is administered into mice [5].   Table 3.

Treatment Protocol for Improved Therapeutic Outcomes
In the above experiment, none of the treatment protocol can suppress the tumor completely on day 18. We continue to investigate whether tumors can be completely killed (the number of tumor cells is less than one) on day 18 post implantation with the same treatment schedule but various dosages. We either vary the amount of oncolytic virus or the blockade rate of CTLA-4 or both. With the same initial settings and obstructing rate of CTLA-4, u = 0.76, if the dosage of oncolytic virus is s = 9 × 10 6 on day 8 as the tumor size reaches to 4 × 10 7 and the treatment lasts for 5 days, the susceptible tumor size is 2.2660 × 10 7 on day 18 post-implantation. If the dosage of oncolytic virus increases to s = 2 × 10 7 , the susceptible tumor size is reduced to 1.2457 × 10 6 . As the dosage of oncolytic virus further increases to 9 × 10 7 , 6 × 10 8 , and 2 × 10 9 , the tumor size is 5.8229 × 10 3 , 160.9773, and 42.7871, respectively. The tumor can be eradicated (T u = 0.9755) on day 18 if s = 7.2 × 10 9 . See Figure 3. Note that the amount of oncolytic virus ranging from 10 6 to 10 9 is adopted in our numerical investigations. On the other hand, 10 9 is used in [15] and 10 7 to 10 8 is implemented in [16], whereas 2 × 10 6 is infused in the mouse experiment of Engeland et al. [5]. The virus inoculum is often manipulated in clinical trials in the orders of magnitude (10 3 -10 10 ) [15]. Thus, the level of oncolytic virus is safe in our simulation.  Table 3.
Next, we examine the influences of the blocking rate of CTLA-4 to the tumor size on day 18 by fixing the amount of oncolytic virus at s = 2 × 10 6 and only varying the CTLA-4 blocking rate u. When the CTLA-4 blocking rate is varied between u = 0.93 and u = 0.9315, the corresponding tumor size is 4.9361 × 10 5 and 8.1300 × 10 3 , respectively. However, if u = 0.9321, the tumor size dropped to T u = 0.0080 on day 18. The tumor can be eradicated on day 18. See Figure 4. It shows that the CTLA-4 blocking rate is very sensitive when it is above a threshold and there exists a critical CTLA-4 blocking rate for tumor eradication.
Following from the above results, the level of oncolytic virus and the rate of blockade of CTLA-4 have deterministic influences on the treatment outcomes. With the above suggested treatment protocols, either the therapeutic outcomes can be improved or the tumor can be eradicated on day 18 after implantation. Observe that the level of oncolytic virus or the CTLA-4 blocking rate has to be high to eradicate the tumor on day 18. The high level of dosages may result in immense immunopathology and cause immune-related adverse events. Therefore, we continue to explore the possible treatment protocol of effective therapy with reduced dosages. If u = 0.87, and s = 5 × 10 9 , the tumor size is 0.2893 on day 18. If u = 0.818, and s = 6 × 10 9 , the tumor size is 0.9669 on day 18. If u = 0.77, and s = 7 × 10 9 , the tumor size is 0.8133 on day 18. The tumor is possible to be suppressed on day 18 if the above treatment protocols are administered following from our numerical results.  The therapy is administered on day 8 and holds for 5 days. The same dosages of oncolytic virus as in Figure 2, s = 2 × 10 6 , is applied to OV-aCTLA-4 treatment. Other parameter values are given in Table 3.

Parameter Sensitivity Analysis
We start by performing global sensitivity analysis of mono-therapy of OVT and also of combined treatment of OVT and anti-CTLA-4 to replicate in silico of a virtual experimental trial with 500 different mice. Sensitivity analysis is used to identify primary parameters that influence treatment efficacy. The majority of range of parameter values are estimated based on one over ten to twice the baseline values as in [27], and otherwise using estimates in the literature when available. The range of values of each model parameter is shown in Table 3. We perform sensitivity analysis using partial rank correlation coefficient (PRCC) analysis and Latin hypercube sampling (LHS) [47,48]. The sensitivity indices of the PRCC vary between -1 and +1, which measure nonlinearly but the strength of monotonic relation between output variables and parameters of interest. The global sensitivity analysis is obtained at the endpoint t = 10 days. This, in particular, is correspondent to 18 days postimplantation, in agreement with the model validation. Note that the dummy parameter does not belong to the model parameters. The model parameters with sensitivity indices less than or equal to that of the dummy parameter are considered not significantly different from zero [47].

Monotherapy of Oncolytic Virus
In this subsection, we investigate the model results when the immune checkpoint inhibitor is not applied. Figure 5 shows the PRCC analysis result in the scenario where each parameter in the model is considered. For this case, the parameters with the strong correlation to the susceptible tumor size are tumor growth rate r u , with P(r u ) = 0.6663, tumor killing rate by immune cells δ x , with P(δ x ) = −0.5705, half-saturation constant of the tumor killing rate by effector cells m x , with P(m x ) = 0.5195, decay rate of tumorsuppressing cytokines γ c , with P(γ c ) = 0.5057, and tumor killing rate by cytokines δ c , with P(δ c ) = −0.4167. See Figure 5.

Combined OVT with Anti-CTLA4
We now discuss the model results when combined OVT and the immune checkpoint inhibitor, anti-CTLA-4, are administered. Figure 6 shows the results of PRCC for each parameter in the model. We denote this PRCC byP for the combined treatment scenarios. For this case, the parameters have a strong relationship with the susceptible tumor size are decay rate of tumor-suppressing cytokines γ c , withP(γ c ) = 0.6121, measure of CTLA-4mediated inhibition on IFN-γ of CD4 + T cells b y , withP(b y ) = 0.5760, tumor killing rate by cytokines δ c , withP(δ c ) = −0.5613, IFN-γ production rate of CD4 + T cells α y , witĥ P(α y ) = −0.5270, and tumor growth rate r u , withP(r u ) = 0.3875. Following from the above analysis, parameters γ c , r u , and δ c play important roles on both the sole OVT and the combined therapy with anti-CTLA-4. We are interested in how important parameters γ c and δ c influence the susceptible tumor size by the day 18 post-implantation, if the same treatment protocol as in Figure 2 is applied. Figure 7a shows susceptible tumor size against γ c , the decay rate of tumor-suppressing cytokines. If we increase γ c to twice the baseline value (γ c = 68), the tumor size increases to 1.0220 × 10 9 by day 18. The tumor size is approximately 2.5 times larger than the tumor size when γ c is at the baseline value. If we decrease γ c to one third value in Table 3, γ c = 11.33, the susceptible tumor size is close to zero on day 18. The tumor can be eradicated on day 18. Figure 7b provides susceptible tumor size against δ c , the tumor killing rate by cytokines. If we decrease δ c to one tenth of the value in Table 3, δ c = 0.02, the susceptible tumor size increases to 1.5989 × 10 9 on day 18 post-implantation. The tumor size is approximately 4 times larger than the tumor size when δ c is at the baseline value. If we we increase δ c to 0.67, the tumor size is 5.0596 × 10 −6 . The tumor can be eradicated on day 18 post-implantation. Collectively, parameters related to CD8 + T cell killing play significant roles for OVT only therapy, whereas parameters related to the level of IFN-γ secreted by CD4 + T cells are important for combined OVT with anti-CTLA4. In addition, the tumor growth rate, decay rate of tumor-suppressing cytokines, and tumor killing rate by cytokines are important factors influencing the behavior and fate of the tumor. To pinpoint the parameters that have a major effect on the combined therapy in relation to the sensitivity analysis of monotherapy of OVT, we investigate the rate of PRCCs between the analysis with sole OVT and the analysis of combined therapy. Among significant parameters, the most marked distinction between this part of analysis with the previous analysis of monotherapy of OVT is the parameter b y , which is the measure of CTLA-4-mediated inhibition on IFN-γ produced by CD4 + T cells. With anti-CTLA-4, the PRCC between b y and tumor size isP(b y ) = 0.5760, while it is P(b y ) = 0.3941 when only OVT is applied. Thus, parameter b y exhibits a much stronger correlation with the susceptible tumor size after treatment when the tumor is administered with additional anti-CTLA-4. This indicates that the measure of CTLA-4-mediated inhibition on IFN-γ produced by CD4 + T cells gives more significance to the efficacy of the combined treatment than to the effectiveness of the sole OV therapy. The second and third marked distinction between current analysis and the result with sole OVT are related to the parameters δ c , the tumor killing rate by cytokines, and α y , the IFN-γ proliferation rate by CD4 + T cells, respectively. With anti-CTLA-4, the PRCC values of δ c and α y , with respect to the tumor size areP(δ c ) = −0.5613, andP(α y ) = −0.5270, respectively, while with single dose of OVT, the corresponding PRCC values are P(δ c ) = −0.4167, and P(α y ) = −0.4054, respectively. It is suggesting that the amount of IFN-γ secreted by CD4 + T cells contributes more significantly to the cogency of the combined treatment than to the effectiveness of sole OVT. We also note that the tumor growth rate, r u , the lysis rate of tumor cells by immune cells, δ x , and the half-saturation constant of cytotoxic killing rate by immune cells, m x , are much less significant with combined therapy than with the monotherapy OVT. It is re-emphasized that IFN-γ plays a more important role on eradicating tumor cells with combined immunotherapies.
Particularly, we are interested in the role of oncolytic viruses on the treatment success. To study this effect, we implement different global sensitivity analysis by varying only those parameters associated with the OVs while keeping all other parameters fixed. This simulates an experiment of mice with comparable similar tumors and immune response but treated with various characters of viruses. We realized that the most two significant oncolytic virus-related parameters are the killing rate d v of virions by innate immune cells, and the proliferation rate a v of immune cells due to OVs. Figure 8 shows that the PRCCs for d v and a v wereP(d v ) = 0.9074 andP(a v ) = −0.8801, respectively. It indicates the strong correlation between these two parameter values and the post-treatment susceptible tumor cell population. The left plot in Figure 9a illustrates the tumor size as a function of killing rate of virions by innate immune cells, d v . With the same treatment protocols used in Figure 2, but if d v is increased to 4.8, 19.2, or 33.6, the susceptible tumor size gradually increases to 9.7220 × 10 8 , 9.8844 × 10 8 , or 9.9085 × 10 9 , respectively on day 18. If d v reaches the upper feasible value, 48, in Table 3, the susceptible tumor size grows up to 9.9161 × 10 8 . This strong correlation reflects the fact that host immune resistance is a major obstacle to intravenous delivery of OVs. Particularly, NK cells are known to indiscriminately attack both uninfected and OV-infected tumor cells rapidly [33,49]. This rapid clearance, however, may halt the desired spread of OVs and hence diminish overall OVT efficacy. Consequently, NK cell response should be minimized in order to allow viruses to replicate sufficiently [50,51]. This result is consistent with preclinical and clinical studies [17]. Figure 8. GSA for the oncolytic virus related parameters under combination therapy.
The right panel in Figure 9b plots tumor size against immune cell proliferation rate a v activated by oncolytic viruses. The simulation endpoint is on t = 18 day. As a v decreases respectively to 4 × 10 −6 , 10 −6 , 2 × 10 −7 , or 4 × 10 −9 , the tumor size increases to 1.6255 × 10 8 , 6.4039 × 10 8 , 9.1278 × 10 8 , or 9.9287 × 10 8 , respectively. On the other hand, if a v increases respectively to 6 × 10 −5 , 8 × 10 −4 , 0.048, or 1.14, the tumor size will reduce to 2.0135 × 10 4 , 118.5440, 6.2881, or 0.9339, respectively. The importance of parameter a v reveals that the major role of the combined therapy is its activation of host anti-tumor immune system responses to post-treatment susceptible tumor population rather than to its direct destruction of the tumor cells.

Discussion
In this work, we developed a mathematical model to investigate the treatment of murine model of malignant melanoma B16-CD20 using an immune checkpoint inhibitor anti-CTLA-4 and OVT. The model is constructed to describe the interactions between the tumor populations (susceptible and infected), virus population, tumor-specific CTL response, natural killer cell populations, virus-specific immune populations, tumor suppression cytokines, and the effect of immune checkpoint inhibitor anti-CTLA-4. In particular, we considered the effects of both types of effector cells, namely CD8 + T and CD4 + T cells. In our model, CD8 + T cells play a direct role of tumor killing, whereas CD4 + T cells kill cancer cells through cytokines they produced. T cells are also recruited by immune cells which become activated when they encounter the oncolytic viruses. Another group of immune cells we considered in the model is NK cells. They can be specific or non-specific during oncolytic virotherapy. NK cells become activated due to immunogenic cell death of infected tumor cells and interferons. They kill tumor cells and clear free viruses. Cytokines activated by both susceptible and infected tumor cells, and secreted by activated T cells and NK cells, conduct the indirect killing. CTLA-4 molecules are expressed on CD8 + T and CD4 + T cells and can prevent the cytokine production. Fragments from infected cancer cells activate the anti-viral immunity which subsequently kills infected cells and clear free virus. Our study focused on the transient behavior of the system because the tumor-immune dynamics are displayed only within few days/weeks after the initial injection of OVs in most of the pre-clinic studies as well as in the mouse experiment. Our model was supported by the experimental data of Engeland et al. [5], which showed that mice treated with combined therapy would greatly reduce tumor burden comparing with the ones treated with either mock or only OVT. We took a step further to investigate treatment protocols to improve therapeutic outcomes by studying various amounts of oncolytic virus and rates of blockade of CTLA-4. It was found that the tumor burden can either be reduced or completely eradicated on day 18 post-implementation for the combined therapy if either the level of oncolytic virus or the rate of blockade of CTLA-4 is increased. The blockade rate of CTLA-4 is very sensitive when it was above a threshold and there exists a critical blockade rate of CTLA-4 for tumor eradication. Moreover, to avoid immense immunopathology due to high level of dosages, we also provide suggestions on effective treatment protocols with reduced dosages.
We performed sensitivity analyses with OVT alone and also with the combined therapy of OVT and anti-CTLA-4 to determine parameters that are most significantly impacting the tumor response to treatments. Our study revealed that parameters related to CD8 + T cells killings have a large impact on the treatment outcome with OVT alone, whereas parameters related to IFN-γ secreted by CD4 + T cells strongly influence treatment responses for combined treatment with anti-CTLA-4. In addition, tumor growth rate, decay rate of tumor-suppressing cytokines, and tumor killing rate by cytokines are important factors in influencing the behavior and fate of the tumor treatment. We also observed that the most three substantial differences between analysis with combined OVT and anti-CTLA-4 and the analysis with monotherapy of OVT are the measure b y of CTLA-4-mediated inhibition on IFN-γ of CD4 + T cells, the tumor killing rate δ c by cytokines, and the IFN-γ production rate α y by CD4 + T cells. It suggests that the IFN-γ produced by CD4 + T cells associated with immunity takes part in a much more important role in responding to the auxiliary therapy than to the efficacy of sole OVT. In addition, we performed sensitivity analyses only for parameters directly related to the oncolytic virus. One of the significant OV related parameters is the killing rate d v of virions by innate immune cells, reflecting the fact that host immune resistance is a major obstacle to intravenous delivery of OVs. It indicated that the virus must first survive interactions with antibody neutralization elimination effect in the blood circulatory system [52]. Immunosuppression inhibition such as anti-CTLA-4 can enhance the oncolytic virus therapy. The other significant oncolytic virus related parameter is the proliferation rate a v of immune cells evoked by the interactions with oncolytic virus, suggesting that the principal role of the combined therapy is its activation of host anti-tumor immune system responses to post-treatment susceptible tumor population rather than to its direct destruction of the tumor cells. Similar conclusions can be found in Storey et al. [29]. With combination therapy, the OVs infect susceptible tumor cells, whereas the lifespan of infected tumor cells is much shorter comparing to susceptible tumor cells. The anti-CTLA-4, on the other hand, stimulates IFN-γ production, which in turn kills tumor cells more effectively. Therefore, it is the combined effect of OVs and anti-tumor immunity that facilitates tumor destruction.
Despite the model's ability on explaining how anti-CTLA-4 can improve virotherapy, we acknowledge that the model has some limitations. Our model is based on the population level interactions and does not account for subcellular events and stimulatory pathways. For example, the effect of IFN-γ on tumor gene expression [53], TLR-ligands derived from the OLVs [54], and T cell exhaustion [55] are not incorporated. Despite these limitations, the model outputs using plausible parameter values agree well with the experimental data of England et al.. Our mathematical model provides useful insights on the dynamics of OV-anti-CTLA4. More complex dynamic interactions between OV, tumor cells, CTLA-4 and other types of cytokines and immune responses will be considered in our future research. In addition, the doses and schedule of the treatment in the experiment may not be optimal. The optimal therapeutic protocol will be investigated for the combination therapy to determine how the schedule and dose can produce optimal treatment outcomes in future research. Models (1)-(9) always have a tumor-and virus-free equilibrium E 0 = (0, 0, 0, 0, 0, 0, 0, 0,N), whereN = s n /γ n . The Jacobian matrix evaluated at E 0 is given by where *s are unimportant entries. It follows that E 0 is asymptotically stable if r u < d uN and r i < a t . If one of these or both inequalities is reversed, then E 0 is a saddle point.
Proof. Notice T i (t) ≤ β t V − a t T i and V (t) ≤ b t a t T i − γ v V for t ≥ 0. Consider Observe that (A1) is a linear, cooperative two-dimensional system, and solutions of (A1) converge to (0, 0) by the assumption γ v > b t β t . Since T i (t) ≤ x(t) and V(t) ≤ y(t) for t ≥ 0, the claim is proven.