Enhancement of chemotherapy using oncolytic virotherapy: Mathematical and optimal control analysis

Oncolytic virotherapy (OV) has been emerging as a promising novel cancer treatment that may be further combined with the existing therapeutic modalities to enhance their effects. To investigate how OV could enhance chemotherapy, we propose an ODE based model describing the interactions between tumour cells, the immune response, and a treatment combination with chemotherapy and oncolytic viruses. Stability analysis of the model with constant chemotherapy treatment rates shows that without any form of treatment, a tumour would grow to its maximum size. It also demonstrates that chemotherapy alone is capable of clearing tumour cells provided that the drug efficacy is greater than the intrinsic tumour growth rate. Furthermore, OV alone may not be able to clear tumour cells from body tissue but would rather enhance chemotherapy if viruses with high viral potency are used. To assess the combined effect of OV and chemotherapy we use the forward sensitivity index to perform a sensitivity analysis, with respect to chemotherapy key parameters, of the virus basic reproductive number and the tumour endemic equilibrium. The results from this sensitivity analysis indicate the existence of a critical dose of chemotherapy above which no further significant reduction in the tumour population can be observed. Numerical simulations show that a successful combinational therapy of the chemotherapeutic drugs and viruses depends mostly on the virus burst size, infection rate, and the amount of drugs supplied. Optimal control analysis was performed, by means of Pontryagin's principle, to further refine predictions of the model with constant treatment rates by accounting for the treatment costs and sides effects.

1. Introduction. Combination therapy strategies have shown significant promise for treating cancers that are resistant to conventional modalities. Oncolytic virotherapy (OV) is an emerging treatment modality that uses replication competent viruses to destroy cancers [22,43,27,30]. Their tumour specific properties allow for viral binding, entry, and replication [32]. Various studies have investigated combination strategies with OV and chemotherapeutic drugs to optimize both viral oncolysis and the effect of the added therapy [8,32,41]. Oncolytic viruses can greatly enhance the cytotoxic mechanisms of chemotherapeutics [35]. Furthermore, chemotherapeutic drugs lyse fast multiplying cells and, in general, virus infected tumour cells quickly replicate [5]. In most combination studies with chemotherapy drugs, apoptosis was increased but viral replication was not enhanced [32,49,3]. For example, Nguyen et al. [32] gave an account of the mechanisms through which drugs can successfully be used in a combination with oncolytic viruses. They, however, noted that the success of this combination depended on several factors including the type of oncolytic virus (OV)-drug combination used, the timing, frequency, dosage, and cancer type targeted. A review of recent clinical studies by Binz and Urlrich [8] shows the classification of possible combination of drugs and oncolytic viruses. Oncolytic virotherapy is still ongoing clinical trials [49,3], and some virus-drug combos, for example talimogene laherparepvec, have been approved for the treatment of melanoma [51].
While, there is a growing body of scientific research showing that combination therapies are the cutting edge for cancer treatment, the design of an optimal protocol remains an open question. The chemo-viro therapy combination is no exception. Despite a noticeable success in the clinical and experimental studies to investigate and characterise the treatment of cancer with chemotherapeutics-virus combinations, much is still not understood about chemovirotherapy. Some of the main open questions in designing an optimal chemovirotherapy treatment are figuring out the right dose combination, the most effective method of drug infusion, and the important treatment characteristics [25].
Mathematical modelling and optimal control theory have over the years played an important role in answering such important research questions which would cost much to set up experimentally. Several mathematical studies including [15,28,37] have addressed the dynamics of oncolytic virotherapy treatment. For example, Ursher [50] gave a summary of some mathematical models for chemotherapy. Tian [47] presented a mathematical model that incorporates burst size for oncolytic virotherapy. His analysis showed that there are two threshold burst size values and below one of them the tumour always grows to its maximum size. His study affirmed that a tumour can be greatly reduced to low undetectable cell counts when the burst size is large enough. A recent study by Malinzi et al. [29] showed that that chemotherapy alone was not able to deplete tumour cells from body tissue but rather in unison with oncolytic viruses could possibly reduce the tumour concentration to a very low undetectable state. Other similar mathematical models on virotherapy include [30,53,54,33,1].
In this paper, we address the question on: "what is the optimal chemotherapeutic drug and virus dosage combination for the elimination of tumour cells in body tissue?" To this end, we develop a mathematical model and optimal control problem that account for the combination of cancer treatment using chemotherapy and virotherapy. The paper is organised as follows. Section 2 is devoted to the model description and the underlying assumptions. In Section 3, we analyse three sub-models: without-treatment model, chemotherapy-only model, and the oncolytic virotherapy-only model. The analysis of the whole model, that is, the chemotherapeutic drug combined with the virus therapy is detailed in Section 4. In Section 5, we investigate an optimal control problem with a quadratic objective function to determine the optimal dosage combination of the chemotherapy drug and the virotherapy. The numerical simulations illustrating our theoretical results are presented in Section 6. We conclude with a discussion in Section 7.
2. Model Formulation. In this section, we propose and formulate our new mathematical model describing the growth of an avascular solid tumour under the effect of both chemotherapy and oncolytic virotherapy treatments. The model considers six state variables summarized in Table 1. We first start by stating the underlying model assumptions in Subsection 2.1, followed by the model equations and their detailed description in Subsection 2.2.

Model description and assumptions.
To model the effects of the chemovirotherapy on the tumour, we consider the dynamics of the following interacting cell populations: tumour cells (both uninfected and virus-infected cells), immune cells (both virus and tumour specific immune responses), the free virus and the drug concentrations. Based on the discussion above and the scientific literature on the interactions between the uninfected and the virus-infected tumour cells in the presence of the oncolytic virotherapy and chemotherapy treatments, the following assumptions are made in setting up the model: 1. Without treatment, the tumour grows logistically with a carrying capacity K [47]. 2. Virus infection, chemotherapeutic drug response to the tumour, chemokine production and immune cells proliferation are considered to be of Michaelis-Menten form to account for saturation [42,2,39]. 3. The model accounts for both virus and tumour specific immune responses [21]. 4. Virus production is a function of virus burst size and the death of infected immune cells. The number of viruses therefore increases as infected tumour cells density multiplies [47]. 5. We consider both virus and tumour specific immune responses. Virus-specific immune response is proportional to the infected tumour cell numbers whilst the tumour specific immune response is taken to be of Michaelis-Menten form to account for saturation in the immune proliferation [23]. 6. We consider a case where drug infusion per day is constant. The constant infusion rate may relate to a situation where a patient is put on an intravenous injection or a protracted venous infusion and the drug is constantly pumped into the body. This form of drug dissemination is used on cancer patients who stay in the hospital for a long period of time. Higher doses of certain anti-cancer drugs may however lead to hepatic veno-occlusive disease, a condition where the liver is obstructed as a result of using high-dose chemotherapy [12,55]. Another more realistic consideration would be to use periodic infusion [7,48]. This would lead to a periodic system of equations for which standard theory on periodic systems applies (see for example [36]). This case is, however, not dealt with in this work, and will be studied in a future work.

2.2.
Model equations. The proposed mathematical model consists of the following system of six differential equations (1)- (6): The model variables and parameters, their meaning and base values are summarised in Table 1 and Table 2 respectively.
Virus specific immune response cells per mm 3 E T (t) Tumour specific immune response cells per mm 3 C(t) Drug concentration grams per millilitre (g/ml) The initial conditions for the model are: where the constants U 0 , I 0 , V 0 , E V0 , E T0 , and C 0 denote the initial concentrations of the uninfected tumour, virus infected tumour, free virus particles, virus specific immune, tumour specific immune, and the chemotherapeutic drug respectively. They are assumed to be nonnegative to make sense biologically. A detailed description of the model follows.
In Eqs. (1) & (2), the terms αU 1 − U+I K represents tumour growth. The term βU V /(K u + U ) describes infection of tumour cells by the virus where β is the infection rate. Drug effect to the uninfected and infected tumour cells is respectively described by the terms δ U U C/(K c + C) and δ I IC/(K c + C) where δ U and δ I are lysis induced rates. K u and K c are Michaelis-Menten constants and relate to killing rates when the virus and the drug are half-maximal [42,2,39]. The term ν U U E T represents tumour specific immune response where ν U is the uninfected lysis rate by tumour specific immune cells. Infected tumour cells have a lifespan 1/δ and are killed by both tumour and virus specific immune cells. These are described by the terms ν I E T I and τ E V I where δ is the natural death rate, ν I and τ are respectively lysis rates of tumour cells by the tumour specific and virus specific immune cells. In Eq. (3), the term bδI represents virus proliferation where b is the virus burst size and δ is the infected tumour cell death rate. The term βU V /(K u +U ) represents  (5), φI is the virus specific immune production and δ v E v represents its deactivation where δ v is the decay rate. The term β T (U + I)/(k + (U + I)) represents tumour-specific immune response where β T is the rate of production and k is a Michaelis constant. Tumour-specific immune decay is described by the term δ T E T where δ T is the rate of decay. In Eq. (6), the time dependent function g(t) represents drug infusion into the tumour and ψC describes drug decay where ψ is the rate of decay.
3. Mathematical analysis. To better understand the dynamics of the proposed model, we begin by examining the model's behaviour about the steady states. This analysis is crucial for identifying the conditions necessary to achieve a tumourfree state. We first show that the model is well posed in a biologically feasible domain, and then proceed with a stability analysis of the model with constant drug infusion. We then analyse the model in the case of no treatment, with chemotherapy alone, and with oncolytic virotherapy alone. The case of the chemotherapeutic drug combined with the virus therapy is detailed as well.
3.1. Well-posedness. Before we proceed with the mathematical analysis, we need to show that the model is well posed in a biologically feasible domain. The model describes the temporal evolution of cell populations, therefore, the cell densities should remain non negative and bounded. The well-posedness theorem is stated below and the proofs are given in Appendix A.

(iii) The trajectories evolve in an attracting region
where C * is the maximum drug concentration and depends on the infusion method. (iv) The domain Ω is positive invariant for the model Equations (1)- (6) and therefore biologically meaningful for the cell concentrations.
We will analyse the model quantitative behaviour in the domain Ω.
3.2. Without-treatment model. To investigate the efficacy of each treatment, their combination as well as the immune response to tumour cells, we first study the model without any form of treatments. The model system of eqs. (1)- ( 6) without treatment reduces to the following system describing the interactions of the uninfected tumour with the tumour specific immune cells: Proposition 1. The model (8) has two biologically meaningful steady states: a tumour-free state, X 1 = (U * = 0, E * T = 0), which is unstable and a tumour endemic state, where M = ν U β T /δ T , which is locally asymptotically stable.
Proof. Equating Eqs. (8) to zero yields: From which if U * = 0 then E * T = 0. The eigenvalues of the Jacobian matrix of the model (8) evaluated at (0, 0) are λ 1 = α and λ 2 = −δ T . One positive and the other negative. Implying that the tumour free equilibrium is unstable. Compare Equation (10) to aU * 2 + bU * + c = 0. We observe that c := −ακ < 0 implying that and therefore there is only one positive root of the quadratic equation (10) and its given by The characteristic polynomial of the Jacobian matrix evaluated at X 2 is given by We prove in Appendix B, using Routh-Hurwitz criterion, that P 1 > 0 and P 0 > 0 and thus the endemic steady state, X 2 , is locally asymptotically stable.
Biological interpretation: Proposition 1 suggests that, without any form of treatment, the tumour-free state is always unstable and the endemic state is stable implying that the tumour would not be eliminated from the body. Next, the model (1)- (6) is studied with only chemotherapy to investigate the effect of the chemotherapeutic drug on tumour cells.
3.3. Chemotherapy-only model. We consider a case where drug infusion per day is constant, that is, g(t) = q. The model with chemo-only, that is I = V = E V = 0, reduces to: Proposition 2. The chemo-only model (13) has two biologically meaningful steady states: a tumour-free sate, ψ and a tumour endemic state, where The endemic state, C 2 , is stable if the conditions (36) in Appendix B are satisfied.
Proof. Setting Equations in (13) to zero gives: . The eigenvalues of the Jacobian of (13) evaluated at C 1 are λ 1 = −ψ, λ 2 = −δ T and λ 3 = α − δU q/ψ Ku+q/ψ implying that for λ 3 to be negative, then αK u > q(δU −α) ψ . If U * = 0 then from Equation (15), which when solved, and using the same analysis as in Proposition 1, gives the expressions in Equation (14). The characteristic polynomial of the Jacobian matrix evaluated at C 2 is given in Equation (32). It directly follows from Routh-Hurwitz criterion that the tumour endemic state, C 2 , is locally asymptotically stable if the conditions (36) are fulfilled.
Biological interpretation: Proposition 2 suggests that a tumour can be eradicated by chemotherapy from the body tissue if the tumour growth rate is less than the drug efficacy (α < δ U ). Otherwise, it would continue to grow uncontrollably. Since drug dosage, q, must not exceed the Maximum Tolerable Dose (MTD), denoted by q MT D , we further derive the following result in Proposition 3: 1. If δ U < α, the tumour-free state will always be unstable for Proof.
1. From the condition for stability in Proposition 2, if δ U < α then αK u > qψ δU −α , which violates the condition for stability. 2. If δ U > α: (i) If q MT D < αψKu δU −α := q * , yet the drug dosage, q, can not be larger than the MTD, q MT D , implying that the condition for stability in Proposition 2 is still violated. (ii) If q MT D > q * , and q ∈ [q * , q MT D ] then the stability condition if fulfilled.
Biological interpretation: Proposition 3 suggests the followings: 1. For a chemotherapeutic drug that is not highly efficient (δ U < α) the tumour cannot be eradicated. 2. If the efficacy, δ U , which is also the lysis rate measured in number of uninfected tumour cells lysed per mm 3 per day, is high enough but with: (i) a very toxic drug, then one cannot give enough dose to eradicate the tumour (ii) a drug which is less toxic, then one can afford to give a dose which is not larger than the MTD and yet which allows for the condition for stability.
3.4. Virotherapy-only model. To determine the effect of virotherapy on tumour cells, we now study the model with only virotherapy treatment. The model (1)- (6) with only virotherapy treatment, that is, C = 0, reduces to: Proposition 4. The virotherapy-only model (17) has three biologically meaningful steady states: a tumour-free state V 1 = (U * = 0, which is unstable, a virus-free state and a tumour endemic state Proof. The above expressions are obtained when model equations in (17) are equated to zero. The eigenvalues of the Jacobian matrix evaluated at V 1 are λ 1 = α, λ 2 = −δ, λ 3 = −γ, λ 4 = −δ T and λ 5 = −δ V with one positive and the others negative implying that the tumour-free state, V 1 , is unstable.
Biological interpretation: Proposition 4 suggests that virotherapy on its own is not capable of eliminating tumour cells.

Chemovirotherapy.
To assess the effects of the combination therapy of virotherapy and chemotherapy, we use the forward sensitivity index to perform a sensitivity analysis of the virus basic reproductive number and the tumour endemic equilibrium with respect to the chemotherapy key parameters. We now proceed to study the whole model system (1)-(6) with constant drug infusion, that is, g(t) = q. We first derive the model's basic reproduction number, R 0 , and calculate its elasticity indices with respect to the model parameters to identify which of them are most sensitive during tumour infection. Finally, we calculate sensitivity indices of the endemic total tumour density with respect to drug infusion. 4.1. Basic reproduction ratio. A basic reproductive number, in our case, can be defined as the average number of new tumour infections generated by one infected cell, via cell lysis, during virotherapy in a completely susceptible cell population [52]. In general, if R 0 > 1, then, on average, the number of new infections resulting from one infected cell is greater than one. Thus, viral infections will persist in tumour cell populations. If R 0 < 1, then, on average, the number of new infections generated by one infected cell in virotherapy is less than one. This implies that the viral infections will eventually disappear from the cell populations. This threshold can as well be used to depict parameters which are most important during tumour infection. We calculate R 0 using the next-generation matrix approach described in [52]. Firstly, we prove that the conditions A 1 − A 5 in Appendix C are satisfied by the model Equations (1)-(6).
Theorem 4.1. The basic reproduction ratio, R 0 , for constant drug infusion g(t) = q, is given by where U * , C * and E * T are all given in Equation (14). The proof of Theorem 4.1 can be found in Appendix C. Next, we calculate elasticity indices of R 0 in response to model parameters to ascertain which them are most critical during chemovirotherapy treatment.  Table 3 shows the sensitivity and elasticity indices obtained with the use of Sage. From the Table one can note that a 1% increase in virus burst size, b, and infection rate, β, leads to a 1% increase in R 0 , whereas a 1% increase in the virus decay rate, γ, leads to a 1% reduction in R 0 . The indices indicate that R 0 is most sensitive to virus burst size, infection and decay rates. Figure 1 is a graphical representation of the indices. Figure 2  Here, we calculate sensitivity indices of the endemic total tumour density with respect to drug infusion. The gist is to determine the relative change in the tumour equilibria when the amount of drug infused is changed and thus infer the feasible amount of the drug that should be infused. The sensitivity index of the total tumour endemic equilibria U * + I * with response to drug infusion q is given by where Γ U * and Γ I * are calculated from the formula in Equation (41) derived in Appendix D. It is worth noting that the endemic equilibrium of the full chemotherapy model (1)- (6) can not explicitly be determined in terms of the parameters and therefore we numerically computed the equilibrium values. We used high values of the tumour reproduction rate and virus burst size, that is α = 0.8 cells/mm 3 /day and b = 50 virions/mm 3 /day so that R 0 > 1. We calculated Γ U * +I * for several values of q. Table 4 shows selected indices for different amounts of the drug infused with their corresponding basic reproduction ratios. Figure 2 (b) is a plot of total tumour sensitivity index with respect to drug infusion. The figure suggests that the feasible amount of drug that should be infused in a patient lies between 50 mg/l and 100 mg/l. This range of values is similar to those inferred from the elasticity indices of R 0 . It is worth noting that when q becomes larger than 50 mg/l, the impact of chemotherapy on virotherapy (R 0 and endemic equilibria) becomes less significant. This analysis, nonetheless, does not incorporate the cost for example monetary or treatments side effects. These results will further be confirmed using optimal control theory. In the next section, we set up an optimal control problem to explicitly determine the optimal virus and chemotherapeutic drug dosage combination necessary for tumour eradication in body tissue.

5.
Optimal control analysis. In this section, we propose and analyze an optimal control problem applied to the chemovirotherapy model to determine the optimal dosage combination of chemotherapy and virotherapy for controlling the tumour. We set the control variables u 1 (t) and u 2 (t) ≡ g(t) to respectively be the supply of viruses from an external source of the drug dosage, which are then incorporated into the model system's equations (3) and (6) to obtain the following control system. For model tractability, we ignore the immune responses.
5.1. The optimal control problem. We wish to determine the optimal combination of controls (u 1 (t), u 2 (t)) that will be adequate to minimize the total tumour density (U (t) + I(t)) together with the cost of the treatment and negative side effects over a fixed time period. The optimization problem under consideration is to minimise the objective functional: where T f is the termination time of the treatment, subject to the control system (21). The two control functions u 1 (t) and u 2 (t) are assumed to be bounded and Lebesgue integrable. We thus seek an optimal control pair (u * 1 , u * 2 ) such that: where Λ is the control set defined by: is the maximum number of viruses body tissue can contain and u MTD 2 is the maximum tolerated drug dose. These may as well be viewed as the maximum amounts a patient can financially afford. The lower bounds for u 1 and u 2 correspond to no treatment. Here, it is important to note that the balancing factors A 1 and A 2 in the objective functional (22) are the relative measures of both the cost required to implement each of the associated controls as well as the negative sides effects due to the treatment.

5.2.
Existence of an optimal control pair. We examine sufficient conditions for the existence of a solution to the quadratic optimal control problem.
Proof. The proof of Proposition 5 is based on Theorem 4.1 in Chapter III of Fleming and Rishel [16]. The necessary conditions for existence are stated and verified as follows.
(1) The set of all solutions to the control system (21) and its associated initial conditions and the corresponding control functions in Λ is non-empty. (2) The control system can be written as a linear function of the control variables with coefficients dependent on time and state variables. (3) The integrand in the objective functional in Equation (21) is convex on Λ.
The right hand sides of the control system (21) are C 1 and bounded below and above (see Appendix A), thus the solutions to the state equations are bounded. It therefore follows from the Picard-Lindelöf theorem [44] that the system is Lipschitz with respect to the state variables. Thus, condition (1) holds. It can be seen from the control system (21) that the right hand sides are linearly dependent of u 1 and u 2 . Thus, condition (2) also holds. To establish condition (3), we notice that the integrand L in the objective functional (21) is convex because it is quadratic in the controls.

Optimal control characterization.
In this section, we characterize the optimal controls u * = (u * 1 , u * 2 ) which gives the optimal levels for the various control measures and the corresponding states (U * , I * , V * , C * ). The necessary conditions for the optimal controls are obtained using the Pontryagin's Maximum Principle [40]. This principle converts the model system (21) into a problem of minimizing pointwise a Hamiltonian, H, with respect to u 1 and u 2 as detailed below.
Proposition 6. Let T = (U, I, V, C) and u = (u 1 , u 2 ). If (T * (t), u * (t)) is an optimal control pair, then there exists a non-trivial vector λ(t) = (λ 1 (t), λ 2 (t), · · · λ 4 (t)) satisfying the following: with transversality conditions and optimal controls: Proof. The Lagrangian and Hamiltonian for the optimal control system (21) are respectively given by: and We thus obtain Equation (24) using the Pontryagin's maximum principle from The transversality conditions are as given in Equation (25) since all states are free at T f . The Hamiltonian is maximized with respect to the controls at the optimal control u * = (u * 1 , u * 2 ). Therefore H is differentiated with respect to u 1 and u 2 on Λ, respectively, to obtain Thus, solving for u * 1 and u * 2 on the interior sets gives By standard control arguments involving the bounds on the controls, 0 ≤ u , we conclude that: 5.4. The optimality system. In summary, the optimality system consists of the control system (21) and the adjoint system (24) with its transversality conditions 25, coupled with the control characterizations (30). Next, we proceed to solve numerically the proposed model and the optimal control problem.
6. Numerical simulations. In this section, we discuss the numerical solutions of both the chemovirotherapy model equations (1)-(6) and the optimal control problem defined in Section 5.4. We also outline the parameter choices and the initial conditions. We use parameter values in Table 2 to solve model equations and the optimality system. The numerical solutions of the model equations are illustrated using MATLAB, while the optimality system was solved using a fourth order Runge-Kutta iterative method.
6.1. Parameter values and initial conditions. Some of the parameter values were obtained from fitted experimental data for untreated tumours and virotherapy in mice [4] and others from biological facts in the literature. A tumour nodule can contain about 10 5 − 10 9 tumour cells [45]. The carrying capacity is therefore considered to be 10 6 cells per mm 3 . In vivo experiments estimate the intrinsic rate of growth to be 0.1 − 0.8 day −1 [6]. We consider the number of viruses produced per day, b, to be in the range 0 − 1000 virions [11]. The amount of drug infused in body tissue, q, is considered to be 5 milligrams per day and the decay rate, ψ, to be 4.17 milligrams per day, values which conform to cancer pharmacokinetic studies [10,34]. Since infected tumour cells multiplication is enhanced by the oncolytic virus replication, the tumour cells lysis, δ I , is considered to be greater than that for uninfected tumour cells δ U (see Ref [38]). We considered a virus-specific immune response rate, φ, of 0.7 day −1 [9]. Both virus and tumour specific immune decay rates are assumed to be 0.01 day −1 , given the fact that their lifespan is less than 100 days [26,31]. The tumour-specific production rate was estimated at 0.5 cells per mm 3 per day [18,26]. In all our simulations, unless stated otherwise, we considered the initial concentrations U 0 = 10000 cells per mm 3 , I 0 = 100 cells per mm 3 , V 0 = 500 virions per mm 3 , E v0 = 100 cells per mm 3 , E T0 = 100 cells per mm 3 and C 0 = 100 g/ml with a high percentage of untreated tumour cell count to require treatment.    Figure 5 (a), an increase in β from 0 to 1 cells per mm 3 per day led to a reduction by about 100 viruses from 120 whereas it led to an increase by about 10 viruses in the first day, from 0 as seen in Figure 5    6.3. Optimal control problem. Figure 6 is a plot of total tumour density of the optimal control solution to the problem formulated in Section 5. It shows the tumour density being reduced by the combinational therapy treatment to a very low state in less than a week.  Figure 6. Total tumour density with optimal control. The tumour density is reduced in a very short time period.

7.
Conclusion. Successful cancer treatment often requires a combination of treatment regimens. Recently, combined oncolytic virotherapy and chemotherapy has been emerging as a promising, effective and synergetic cancer treatment. In this paper we consider a combination therapy with chemotherapy and virotherapy to investigate how virotherapy could enhance chemotherapy. To this end, we developed a mathematical model and an optimal control problem to determine the optimal chemo-virus combination to eradicate a tumour. We firstly validated the model's plausibility by proving existence, positivity and boundedness of the solutions. We analysed the model in four scenarios: without treatment, with chemotherapy, with oncolytic virotherapy alone, and with combined chemotherapeutic drug and virus therapy. A basic reproduction number for the infection of the tumour cells was calculated to analyse the model's tumour endemic equilibrium. Furthermore, sensitivity and elasticity indices of the basic reproduction number and tumour endemic equilibria with respect to drug infusion were calculated. The optimal control problem was solved using the Pontryagin's Maximum Principle. Numerical solutions to the proposed model were carried out to illustrate the analysis results. Similarly, numerical solutions to the optimal control problem were carried out to identify the optimal dosages for the minimization of the chemo-virus combination. Model analysis and simulations suggest the following results.
The stability analysis showed that a tumour can grow to its maximum size in a case where there is no treatment. It also demonstrated that chemotherapy alone is capable of clearing tumour cells provided that the drug efficacy is greater than the intrinsic tumour growth rate. This can be evidenced in a study by Dasari and Tchounwou [14] and references therein. Furthermore, the result emphasised that if the chemotherapeutic drug's efficacy is high enough but with a toxic drug then the tumour cannot be eliminated. If, however, the drug is not toxic, then a dosage which (i)  is not larger than the maximum tolerated dose can be given while still allowing for a stability condition of the tumour free equilibrium, thus eliminating the tumour.
In the case of virotherapy alone, stability analysis revealed that virotherapy on its own may not clear tumour cells but rather highly enhances chemotherapy in destroying tumour cells. Several clinical and experimental studies for example [35,8,20,41] confirm this result. In a review of clinical studies by Binz and Ulrich [8], it is stated that phase II/III clinical trials on combining adenovirus H101 with a chemotherapeutic drug, cisplatin and 5-fluorouracil, leads to a 40% improvement compared to chemotherapy alone for the treatment of patients with head and neck cancer. Adenovirus Ad-H101 was consequently approved for the treatment of head and neck cancers in China and recently in the United states after phase III clinical trials showed a 72-79% response rate for the combination of the virus with 5fluorouracil compared to a 40% with chemotherapy alone [19].
Sensitivity analysis and elasticity indices of the basic reproduction ratio indicated that successful chemovirotherapy is highly dependent on the virus burst size and its infection rate as well as on the drug infusion and its decay rates. Larger virus burst size and higher replication rates lead to a lower tumour concentration. Whereas, low decay and high doses of the drug that the body can tolerate lead to better treatment results. Sensitivity analysis of both the basic reproduction ratio and model equilibria suggested the feasible drug infusion to be between 50 and 100 milligrams per litre. The effect of virus burst size and infection rates in determining the success of virotherapy have extensively been confirmed in recent mathematical studies (see for example [29,47,13]). Nonetheless, the right dose of chemotherapy treatment has always been a question and a concern to both clinicians and mathematicians alike [25,17]. In this study, optimal control results further indicated that 50% of the maximum drug and virus tolerated dosages (MTDs) optimize chemovirotherapy.
In addition, numerical simulations suggested that with the use of both chemotherapy and virotherapy, a tumour may be eradicated in less than a month. We know that this is not realistic for human cancers and is strictly based on the experimental data used in the numerics. Simulations further showed that a successful combinational therapy of cancer drugs and viruses is mostly dependent on virus burst size, infection rate and the amount of drugs supplied into a patient's body which is in agreement with recent studies [29,47,13].
The mathematical model we developed in this study is a simple one, for example, it considers the cell densities to only be time dependent. It is therefore imperative for further studies to incorporate more biological aspects like spatial variation of the cell concentrations. Another facet would be to investigate the effect of toxicity of both viruses and the drug on normal body tissue. The key to improving combined virotherapy and chemotherapy lies in quantifying the dependence of treatment outcome on immune stimulation. Another extension to this model, would be to include other subpopulations of the immune system. Nevertheless, the results here emphasise the treatment characteristics that are vital in combining drugs and viruses to treat cancer and an optimal drug and virus dosage is suggested.

Appendices
A. Appendix A.

A.1. Proof of Theorem 3.1(i) on the existence and uniqueness of model solutions
Proof. The functions on the right hand side of Equations (1)-(6) are C 1 on R n . Thus, it follows from the Picard-Lindelöf theorem (see [44]) that (1)-(6) exhibit a unique solution. Proof. Denote U * := U p + U r where The characteristic polynomial of the Jacobian matrix evaluated at X 2 is given by where P 0 = P 01 U * + P 00 and P 1 = P 11 U * + P 10 , Using Routh-Hurwitz criterion, the endemic steady state, X 2 , is locally asymptotically stable if P 1 > 0 and P 0 > 0. We notice that P 01 > 0. This means that P 0 > 0 if and only if U * > −P 00 /P 01 which is equivalent to With the help of Maple R , we find that Therefore P 0 > 0.
For P 1 , we discuss the following cases: (i) If P 11 > 0 then P 1 > 0.
(ii) If P 11 < 0 then P 1 > 0 if and only if U * < −P 10 /P 11 which is equivalent to Using Maple R , we find that This implies that P 1 > 0. Therefore the endemic steady state, X 2 , is locally asymptotically stable.
B.2. Proof of stability of endemic state, C 2 , in Proposition 2 Proof. The characteristic polynomial, P (x), of the Jacobian evaluated at C 2 is given by The conditions for stability of the endemic state, C 2 , therefore follow from the Routh-Hurwitz criterion and are stated as: C. Appendix C.
C.1. Proof of conditions A 1 − A 5 for the derivation of R 0 Proof. Denote x = (x 1 , x 2 , · · · x 5 ) T = (U, I, V, E V , E T ) T with x i ≥ 0 be the set of cell compartments divided into infected cells i = 1, · · · , m and uninfected cells i = m + 1, · · · , n. Define X 0 to be the set of virus or infection free states: X 0 := {x ≥ 0 : x i 0, for all i = 1, · · · , m}.
Let F i (t, x) be the input rate of newly infected cells in the i th compartment. Let V + i (t, x) and V − i (t, x) respectively be the input rate of cells by other means and rate of transfer of cells out of a compartment i. Thus the system of Equations (1)-(6) take the form We now proceed to prove the assumptions A 1 − A 5 : A 1 : For each 1 ≤ i ≤ n, the functions F i (t, x), V + i and V − i are non-negative and continuous on R × R n + and continuously differentiable with respect to x. From Equations (1)-(6) (37) Clearly, F i , V + i and V − i satisfy A 1 . A 2 : If x i = 0, then V i = 0. Clearly, from (37), if the state variables are all equal to zero, then V i = 0.