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 [23,43,28,31]. Their tumour specific properties allow for viral binding, entry, and replication [33]. Various studies have investigated combination strategies with OV and chemotherapeutic drugs to optimize both viral oncolysis and the effect of the added therapy [8,33,42]. Oncolytic viruses can greatly enhance the cytotoxic mechanisms of chemotherapeutics [36]. 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 [33,49,3]. For example, Nguyen et al. [33] 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 [26].

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 [16,29,38] 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. [30] 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 [31,53,54,34,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.

---

2.1. Model description and assumptions

To model the effects of the chemo-virotherapy 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 [15,2,40].

3. The model accounts for both virus and tumour specific immune responses [22].

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 [24].

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 [37]). 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):

$\frac{dU}{dt} = \alpha U\left(1-\frac{U+I}{K}\right) - \frac{\beta UV}{K_u +U} -\nu_U U E_T - \frac{\delta_U UC}{K_c +C}$

(1)

$\frac{dI}{dt} = \frac{\beta UV}{K_u +U}- \delta I - \nu_I E_T I -\tau E_V I -\frac{\delta_I IC}{K_c +C}$

(2)

$\frac{dV}{dt} = b\delta I- \frac{\beta UV}{K_u +U}-\gamma V$

(3)

$\frac{dE_V}{dt} = \phi I - \delta_V E_V$

(4)

$\frac{d E_T}{dt} = \frac{\beta_T (U+I)}{\kappa+ (U+I)}-\delta_T E_T$

(5)

$\frac{dC}{dt} = g(t) - \psi C$

(6)

The model variables and parameters, their meaning and base values are summarised in Table 1 and Table 2respectively.

The initial conditions for the model are:

$\label{ICs} U(0) = U_{0}, I(0) = I_{0} , V(0) = V_{0}, E_V(0) = E_{V_0}, E_T(0) = E_{T_0}, C(0) = C_{0},$

(7)

where the constants $U_{0}$, $I_{0}$, $V_{0}$, $E_{V_0}$, $E_{T_0}$, 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 ${\alpha U\left(1-\frac{U+I}{K}\right)}$ represents tumour growth. The term $\beta U V/(K_u+U)$ describes infection of tumour cells by the virus where $\beta$ is the infection rate. Drug effect to the uninfected and infected tumour cells is respectively described by the terms $\delta_U U C/(K_c+C)$ and $\delta_I I C/(K_c+C)$ where $\delta_U$ and $\delta_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 [15,2,40]. The term $\nu_U U E_T$ represents tumour specific immune response where $\nu_U$ is the uninfected lysis rate by tumour specific immune cells. Infected tumour cells have a lifespan $1/\delta$ and are killed by both tumour and virus specific immune cells. These are described by the terms $\nu_I E_T I$ and $\tau E_V I$ where $\delta$ is the natural death rate, $\nu_I$ and $\tau$ are respectively lysis rates of tumour cells by the tumour specific and virus specific immune cells.

In Eq. (3), the term $b \delta I$ represents virus proliferation where $b$ is the virus burst size and $\delta$ is the infected tumour cell death rate. The term $\beta UV/(K_u+U)$ represents loss of free virus due to infection of the uninfected tumour cells. Virus deactivation in body tissue is represented by the term $\gamma V$ where $\gamma$ is the rate of decay.

In Eqs. (4) & (5), $\phi I$ is the virus specific immune production and $\delta_v E_v$ represents its deactivation where $\delta_v$ is the decay rate. The term $\beta_T (U+I)/(k+(U+I))$ represents tumour-specific immune response where $\beta_T$ is the rate of production and $k$ is a Michaelis constant. Tumour-specific immune decay is described by the term $\delta_T E_T$ where $\delta_T$ is the rate of decay.

In Eq. (6), the time dependent function $g(t)$ represents drug infusion into the tumour and $\psi C$ describes drug decay where $\psi$ 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 tumour-free 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 7.

Theorem 3.1.  (ⅰ) There exists a unique solution to the system of equations (1)-(6) in the region $\mathbb{D} = (U,I,V,E_V,E_T,C) \in \mathbb{R}^{6}_+$.

(ⅱ) If $U(0)\geq 0$, $I(0)\geq 0$, $V(0)\geq 0$, $E_V(0)\geq 0, E_T(0)\geq 0$ and $C(0)\geq 0$, then $U(t)\geq 0$, $I(t)\geq 0$, $V(t)\geq 0$, $E_V(t)\geq 0$, $E_T(t) \geq 0$ and $C(t)\geq 0$ for all $t\geq 0$.

(ⅲ) The trajectories evolve in an attracting region $\Omega = \{(U,I,V,E_V,E_T,C)\in \mathbb{D}\mid U(t)+I(t)\leq K,~ V(t) \leq b/\gamma,~ E_V(t)\leq \phi/\delta_V,E_T\leq \beta_T/\delta_T,C(t)\leq C^*\}$, where $C^*$ is the maximum drug concentration and depends on the infusion method.

(ⅳ) The domain $\Omega$ 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 $\Omega$.

---

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:

$\begin{align} \label{eq:no_treatment} \frac{dU}{dt}& = \alpha U\left(1-\frac{U}{K}\right) - \nu_UUE_T, \nonumber \\ \frac{d E_T}{dt}& = \frac{\beta_T U}{\kappa+ U}-\delta_T E_T. \end{align}$

(8)

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,

$\begin{align} \label{eq:steadystates1:a} X_2 = \bigg[ U^*& = \frac{K}{2\alpha}\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg) + \frac{K}{2\alpha}\sqrt{\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg)^2 +4\frac{\alpha^2\kappa}{K}}, \;\;\;E_T^* = \frac{\beta_T}{\delta_T} \bigg( \frac{U^*}{\kappa + U^*}\bigg)\bigg], \end{align}$

(9)

where $M = \nu_U \beta_T/\delta_T$, which is locally asymptotically stable.

Proof. Equating Eqs. (8) to zero yields:

$\frac{\alpha }{K}U^{*^2}+ \bigg(\frac{\alpha\kappa}{K} + M - \alpha \bigg)U^* -\alpha\kappa = 0$

(10)

$E_T^* = \frac{\beta_T}{\delta_T} \bigg( \frac{U^*}{\kappa + U^*}\bigg)$

(11)

From which if $U^* = 0$ then $E_T^* = 0$. The eigenvalues of the Jacobian matrix of the model (8) evaluated at $(0,0)$ are $\lambda_1 = \alpha$ and $\lambda_2 = -\delta_T$. One positive and the other negative. Implying that the tumour free equilibrium is unstable. Compare Equation (10) to $a U^{*^2} + bU^* + c = 0$. We observe that $c: = -\alpha \kappa <0$ implying that $\Delta = b^2-4ac: = \bigg(\alpha - \left( \frac{\alpha \kappa}{K}+M \right)\bigg)^2 +4\frac{\alpha^2\kappa}{K}>0$ and therefore there is only one positive root of the quadratic equation (10) and its given by

$U^* = \frac{-b+\sqrt{b^2-4ac}}{2a}: = \frac{K}{2\alpha}\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg) + \frac{K}{2\alpha}\sqrt{\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg)^2 +4\frac{\alpha^2\kappa}{K}}.$

The characteristic polynomial of the Jacobian matrix evaluated at $X_2$ is given by

$\begin{align} \label{eq:char1:a} P(x) = x^2 + P_1x + P_0. \end{align}$

(12)

We prove in Appendix 7, 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:

$\frac{dU}{dt} = \alpha U\left(1-\frac{U}{K}\right) -\nu_U U E_T - \frac{\delta_U UC}{K_c +C}, \nonumber \\ \frac{d E_T}{dt} = \frac{\beta_T U}{\kappa+ U} -\delta_T E_T, \nonumber \\ \frac{dC}{dt} = q - \psi C.$

(13)

Proposition 2. The chemo-only model (13) has two biologically meaningful steady states: a tumour-free sate, $C_1 = (U^* = 0,E_T^* = 0,C^* = \frac{q}{\psi})$, which is stable if $\alpha K_u > \frac{q(\delta_U-\alpha)}{\psi}$ and a tumour endemic state,

$\label{eq:steadystates1} C_2 = \bigg[ U^* = \frac{K}{2\alpha}\bigg(\alpha - \left( \frac{\alpha \kappa}{K}+L+M \right)\bigg) + \frac{K}{2\alpha}\sqrt{\bigg(\alpha - \left( \frac{\alpha \kappa}{K} +L+M \right)\bigg)^2 +4\frac{\alpha\kappa}{K}(\alpha-L)} \nonumber\\ \;\;\;\;\;\;\;\;\;E_T^* = \frac{\beta_T}{\delta_T} \bigg( \frac{U^*}{\kappa + U^*}\bigg) \;\;\;and \;\; \;C^* = \frac{q}{\psi}\bigg],$

(14)

where

$L = \frac{\delta_U C^*}{K_c + C^*},\;\;\;\; M = \frac{\nu_U \beta_T}{\delta_T}.$

The endemic state, $C_2$, is stable if the conditions (36) in Appendix 7 are satisfied.

Proof. Setting Equations in (13) to zero gives:

$\alpha U^*\left(1-\frac{U^*}{K} \right) - \nu_U U^*E_T^* - \frac{\delta_U U^* C^*}{K_c + C^*} = 0,$

(15)

$\frac{\beta_T U}{\kappa+U^*} -\delta_T E_T^* U^* = 0, \;\;\; C^* = \frac{q}{\psi},$

(16)

from which if $U^* = 0$ then $C_1 = (U^* = 0, E_T = 0, C^* = \frac{q}{\psi} )$. The eigenvalues of the Jacobian of (13) evaluated at $C_1$ are $\lambda_1 = -\psi$, $\lambda_2 = -\delta_T$ and $\lambda_3 = \alpha- \frac{\delta_U q/\psi}{K_u +q/\psi}$ implying that for $\lambda_3$ to be negative, then $\alpha K_u > \frac{q(\delta_U - \alpha)}{\psi}$. If $U^* \neq 0$ then from Equation (15),

$E_T^* = \frac{\beta_T}{\delta_T} \bigg( \frac{U^*}{\kappa + U^*}\bigg) \;\;\text{and}\;\; \frac{\alpha }{K}U^{*^2}+ \bigg(\frac{\alpha\kappa}{K} + M +L - \alpha \bigg)U^* + \kappa (L-\alpha ) = 0$

which when solved, and using the same analysis as in Proposition 14, gives the expressions in Equation (32). 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 ($\alpha < \delta_U$). Otherwise, it would continue to grow uncontrollably.

Since drug dosage, $q$, must not exceed the Maximum Tolerable Dose (MTD), denoted by $q^{MTD}$, we further derive the following result in Proposition 3:

Proposition 3.  1. If $\delta_U < \alpha$, the tumour-free state will always be unstable for any $q < q^{MTD}$.

2. If $\delta_U > \alpha$:

(ⅰ) If $q^{MTD} < \frac{\alpha \psi K_u }{\delta_U - \alpha}: = q^*$, then the tumour-free state is unstable for all $q<q^{MTD}$.

(ⅱ) If $q^{MTD} \geq q^*$, then the tumour-free state is locally asymptotically stable for all $q\in [q^*,q^{MTD}]$.

Proof. 1.From the condition for stability in Proposition 2, if $\delta_U < \alpha$ then $\alpha K_u > \frac{q\psi}{\delta_U - \alpha}$, which violates the condition for stability.

2. If $\delta_U > \alpha$:

(ⅰ) If $q^{MTD} < \frac{\alpha \psi K_u }{\delta_U - \alpha}: = q^*$, yet the drug dosage, $q$, can not be larger than the MTD, $q^{MTD}$, implying that the condition for stability in Proposition 2 is still violated.

(ⅱ) If $q^{MTD} > q^*$, and $q\in [q^*, q^{MTD}]$ then the stability condition if fulfilled.

Biological interpretation: Proposition 3 suggests the followings:

1. For a chemotherapeutic drug that is not highly efficient ($\delta _U < \alpha$) the tumour cannot be eradicated.

2. If the efficacy, $\delta_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:

(ⅰ) a very toxic drug, then one cannot give enough dose to eradicate the tumour

(ⅱ) 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:

$\begin{align} \label{eq:viro-onlymodel} \frac{dU}{dt}& = \alpha U\left(1-\frac{U+I}{K}\right) - \frac{\beta UV}{K_u +U} -\nu_U U E_T, \nonumber \\ \frac{dI}{dt}& = \frac{\beta UV}{K_u +U}- \delta I - \nu_I E_T I -\tau E_V I, \nonumber \\ \frac{dV}{dt}& = b\delta I- \frac{\beta UV}{K_u +U} -\gamma V, \nonumber \\ \frac{dE_V}{dt}& = \phi I - \delta_V E_V, \nonumber \\ \frac{d E_T}{dt}& = \frac{\beta_T (U+I)}{\kappa+ (U+I)}-\delta_T E_T. \end{align}$

(17)

Proposition 4. The virotherapy-only model (17) has three biologically meaningful steady states: a tumour-free state $V_1 = \left(U^* = 0,I^* = 0, V^* = 0, E_T^* = 0, E_V^* = 0 \right)$ which is unstable, a virus-free state

$\begin{align} \label{eq:steadystates3} V_2 = \bigg[ U^*& = \frac{K}{2\alpha}\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg) + \frac{K}{2\alpha}\sqrt{\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg)^2 +4\frac{\alpha^2\kappa}{K}}, \;\;\;I^* = 0, \nonumber \\ &\;\;\;~V^* = 0, E_T^* = \frac{\beta_T}{\delta_T} \bigg( \frac{U^*}{\kappa + U^*}\bigg), \;\;\;~E_V^* = 0 \bigg] \end{align}$

(18)

and a tumour endemic state $V_3 = (U^*,I^*,V^*,E_T^*,E_V^*)$.

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 $\lambda_1 = \alpha$, $\lambda_2 = -\delta$, $\lambda_3 = -\gamma$, $\lambda_4 = -\delta_T$ and $\lambda_5 = -\delta_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.

---

4. 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 7 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

$\begin{align} \label{eq:Ro-constantinfusion} R_0 ~ = \frac{b\beta \delta U^*(K_c + C^*)}{ \gamma \big[(K_c+ C^*) (\nu_I E_T^* +\delta) + \delta_I C^* \big]\big(K_u + U^* \big)}, \end{align}$

(19)

where $U^*$, $C^*$ and $E_T ^*$ are all given in Equation (14).

The proof of Theorem 4.1 can be found in Appendix 7. Next, we calculate elasticity indices of $R_0$ in response to model parameters to ascertain which them are most critical during chemovirotherapy treatment.

---

4.2. Sensitivity of the basic reproduction ratio

Definition 4.2. The Sensitivity and elasticity indices of the basic reproduction ratio, $R_0$, with respect to model parameter $p$ are respectively given by $S_p = \frac{\partial R_0}{\partial p}$ and $e_p = \frac{\partial R_0}{\partial p}\frac{p}{R}$.

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, $\beta$, leads to a $1\%$ increase in $R_0$, whereas a $1\%$ increase in the virus decay rate, $\gamma$, 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 (a) shows the variation of elasticity of $R_0$ with respect to drug dosage. The figure indicates that increasing the amount of drug infused reduces $R_0$ and thus reduces the number of tumour cells. From the figure we can further infer that values of the drug dosage, $q$, from 40 mg/l to 100 mg/l have a less significant impact on $R_0$.

---

4.3. Sensitivity of endemic equilibria

Sensitivity indices of the endemic equilibria is informative about the establishment of the disease. 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

$\begin{align} \label{eq:sensitivity-of-tumour} \Gamma^{U^*+ I^*}_q& = \frac{\partial (U^*+I^*)}{\partial q} \cdot \frac{q}{U^*+I^*} \nonumber \\& = \frac{U^*}{U^*+I^*}\Gamma^{U^*} + \frac{I^*}{U^*+I^*} \Gamma^{I^*}. \end{align}$

(20)

where $\Gamma ^{U^*}$ and $\Gamma ^{I^*}$ are calculated from the formula in Equation (41) derived in Appendix 7. 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 $\alpha = 0.8$ cells/mm$^3$/day and $b = 50$ virions/mm$^3$/day so that $R_0>1$. We calculated $\Gamma ^{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) \equiv 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.

$\frac{dU}{dt} = \alpha U\left(1-\frac{U+I}{K}\right) - \frac{\beta UV}{K_u +U} - \frac{\delta_U UC}{K_c +C} \nonumber \\ \frac{dI}{dt} = \frac{\beta UV}{K_u +U}- \delta I -\frac{\delta_I IC}{K_c +C} \nonumber \\ \frac{dV}{dt} = b\delta I- \frac{\beta UV}{K_u +U}-\gamma V + {\bf{u_1(t)}} \nonumber \\ \frac{dC}{dt} = {\bf{u_2(t)}} - \psi C$

(21)

---

5.1. The optimal control problem

We wish to determine the optimal combination of controls $\left(u_1(t),u_2(t)\right)$ 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:

$\begin{align} \label{eq:objectivefunctional} J(u_1,u_2) = \int_0^{T_f} \bigg[ U(t)+ I(t) + \left(\frac{A_1}{2}u_1^2(t) +\frac{A_2}{2} u_2^2(t)\right)\bigg]dt \end{align}$

(22)

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:

$\begin{align} \label{eq:minimize} J(u_1^*,u_2^*) = \min \big\{J(u_1,u_2) | (u_1,u_2)\in \Lambda \big\} \end{align}$

(23)

where $\Lambda$ is the control set defined by:

$\begin{align} \Lambda = \big\{(u_1,u_2)|u_i ~\text{is measurable with} ~0 \leq u_i(t)\leq u_i^{\text{MTD}}, ~t\in [0,T_f], ~i = 1,2\big\} \nonumber \end{align}$

where $u_1^{\text{MTD}}$ is the maximum number of viruses body tissue can contain and $u_2^{\text{MTD}}$ 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.

Proposition 5. There exists an optimal control pair $(u_1^*,u_2^*)$ with a corresponding solution ($U^*$, $I^*$, $V^*$, $C^*$) to the model system (21) that minimizes $J(u_1, u_2)$ over $\Lambda$.

Proof. The proof of Proposition 5 is based on Theorem 4.1 in Chapter Ⅲ of Fleming and Rishel [17]. 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 $\Lambda$ 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 $\Lambda$.

The right hand sides of the control system (21) are $\mathbb{C}^1$ and bounded below and above (see Appendix 7), 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 $\mathcal{L}$ in the objective functional (21) is convex because it is quadratic in the controls.

---

5.3. 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 [41]. This principle converts the model system (21) into a problem of minimizing pointwise a Hamiltonian, $\mathcal{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

$\lambda(t) = (\lambda_1(t),\lambda_2(t), \cdots \lambda_4(t) )$

satisfying the following:

$\lambda_1 ^{\prime}(t) = -1 -\lambda_1 \bigg[ \alpha \left(1-\frac{2U+I}{K}\right) - \frac{\beta K_u V }{(K_u +U)^2} - \frac{\delta_U C}{K_c +C} \bigg],\\ \;\;\;\;\;\; \;\;\; -\lambda_2 \bigg[\frac{\beta K_u V}{(K_u +U)^2}\bigg] + \lambda_3 \bigg[\frac{\beta K_u V}{(K_u +U)^2}\bigg],\\ \lambda_2 ^{\prime}(t) = -1-\lambda_1 \bigg[\frac{\alpha U}{K} \bigg] + \lambda_2 \bigg[\delta +\frac{\delta_I C}{K_c+C} \bigg]-\lambda_3 b\delta, \\ \lambda_3 ^{\prime}(t) = \lambda_1 \bigg[\frac{\beta U}{K_u +U}\bigg] -\lambda_2 \bigg[\frac{\beta U}{K_u +U} \bigg] +\lambda_3 \bigg[\frac{\beta U}{K_u +U}+\gamma \bigg],\\ \lambda_4 ^{\prime}(t) = \lambda_1 \bigg[\frac{K_c \delta_u U}{(K_c +C)^2} \bigg] + \lambda_2 \bigg[\frac{K_c \delta_I I}{(K_c +C)^2} \bigg] +\lambda_4 \psi,$

(24)

with transversality conditions

$\begin{align}\label{eq:transversality-conditions} \lambda_i (T_f) = 0, \;\;\; i = 1,2,\cdots 4. \end{align}$

(25)

and optimal controls:

$\begin{align} \label{eq:conditions:a} u_1^* = \min \bigg[u_1^{\text{MTD}},\max \bigg(0,-\frac{\lambda_3}{A_1}\bigg)\bigg],\;\;\; u_2^* = \min \bigg[u_2^{\text{MTD}},\max \bigg(0,-\frac{\lambda_4}{A_2}\bigg)\bigg]. \end{align}$

(26)

Proof. The Lagrangian and Hamiltonian for the optimal control system (21) are respectively given by:

$\begin{align} \label{eq:Lagrangian} \mathcal{L} (U,I,V,C,u_1,u_2) = U +I + \frac{1}{2}\bigg[A_1u_1^2+A_2 u_2^2\bigg] \end{align}$

(27)

and

$\label{eq:Hamiltonian} \mathcal{H} = U+ I + \frac{1}{2}\bigg[A_1u_1^2+A_2 u_2^2\bigg] + \lambda_1 \bigg[ \alpha U\left(1-\frac{U+I}{K}\right)- \frac{\beta UV}{K_u +U} - \frac{\delta_U UC}{K_c +C} \bigg]\\ + \lambda_2 \bigg[ \frac{\beta UV}{K_u +U}- \delta I -\frac{\delta_I IC}{K_c +C} \bigg] + \lambda_3 \bigg[b\delta I- \frac{\beta UV}{K_u +U}-\gamma V + u_1 \bigg] + \lambda_4\big[u_2 - \psi C \big].$

(28)

We thus obtain Equation (24) using the Pontryagin's maximum principle from

$\begin{align} \label{eq:conditions:b} \begin{split} \lambda ^{\prime}(t) = -\frac{\partial \mathcal{H}}{\partial T}\bigg(t,T^*(t),u^*(t), \lambda(t) \bigg). \end{split} \end{align}$

(29)

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 $\mathcal{H}$ is differentiated with respect to $u_1$ and $u_2$ on $\Lambda$, respectively, to obtain

$\begin{align*} 0 = \frac{\partial \mathcal{H }}{\partial u_1}& = A_1 u_1 + \lambda_3,\\ 0 = \frac{\partial \mathcal{H}}{\partial u_2}& = A_2 u_2 +\lambda_4. \end{align*}$

Thus, solving for $u_1^*$ and $u_2^*$ on the interior sets gives

$u_1^* = -\frac{\lambda_3}{A_1}, \;\;\; u_2^* = -\frac{\lambda_4}{A_2}.$

By standard control arguments involving the bounds on the controls, $0\leq u_1(t)\leq u_1^{\text{MTD}},~~0\leq u_2(t)\leq u_2^{\text{MTD}}$, we conclude that:

$\begin{align} \label{eq:OCconditions} u_1^* = \min \bigg[u_1^{\text{MTD}},\max \bigg(0,-\frac{\lambda_3}{A_1}\bigg)\bigg], \;\;\;u_2^* = \min \bigg[u_2^{\text{MTD}},\max \bigg(0,-\frac{\lambda_4}{A_2}\bigg)\bigg]. \end{align}$

(30)

---

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, $\psi$, to be $4.17$ milligrams per day, values which conform to cancer pharmacokinetic studies [10,35]. Since infected tumour cells multiplication is enhanced by the oncolytic virus replication, the tumour cells lysis, $\delta_I$, is considered to be greater than that for uninfected tumour cells $\delta_U$ (see Ref [39]). We considered a virus-specific immune response rate, $\phi$, 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 [27,32]. The tumour-specific production rate was estimated at $0.5$ cells per mm$^3$ per day [19,27]. 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_{v_0} = 100$ cells per mm$^3$, $E_{T_0} = 100$ cells per mm$^3$ and $C_0 = 100$ g/ml with a high percentage of untreated tumour cell count to require treatment.

---

6.2. Chemovirotherapy model

Figures 3 (a) & (b) are respectively plots of the uninfected tumour density for different virus and drug doses. Both figures depict that an increase in the virus burst size ($b$) and drug dose ($q$) reduces the tumour density. In Figure 3 (a), without virotherapy, that is $b = 0$ virions, the tumour is depleted in about 10 days whereas without chemotherapy, that is $q = 0$ milligrams per day, it takes about 20 days.

Figures 4 and 5 (a) & (b) respectively show the effect of virus burst size on the virus and infected tumour cell densities. Both figures show that the densities with increasing burst size and infection rate. From Figure, for example 4, an increase of $b$ from 250 to 500 virions led to an increase of the virus number to about $700$ $virions$ from $500$ $virions$, in the first 100 days.

Figure 5 (a) & (b) show that increasing virus infection rate, $\beta$, reduces the virus density and increases the infected tumour density. For example from Figure 5 (a), an increase in $\beta$ 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 (b).

---

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.

Figures 7 (a) & (b) represent the optimal controls $u_1$ and $u_2$. The figures suggest that the optimal values of the number of viruses and drug concentration that yields efficient combination therapy are half their maximal dosages $u_1^{\text{MTD}}$ and $u_2^{\text{MTD}}$.

---

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 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 [36,8,21,42] confirm this result. In a review of clinical studies by Binz and Ulrich [8], it is stated that phase Ⅱ/Ⅲ 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 Ⅲ clinical trials showed a 72-79% response rate for the combination of the virus with 5-fluorouracil compared to a 40% with chemotherapy alone [20].

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 [30,47,13]). Nonetheless, the right dose of chemotherapy treatment has always been a question and a concern to both clinicians and mathematicians alike [26,18]. 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 [30,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.

---

Appendix A

A.1. Proof of Theorem 3.1(ⅰ) on the existence and uniqueness of model solutions

Proof. The functions on the right hand side of Equations (1)-(6) are $\mathbb{C}^1$ on $\mathbb{R}^n$. Thus, it follows from the Picard-Lindelöf theorem (see [44]) that (1)-(6) exhibit a unique solution.

A.2. Proof of Theorem 3.1(ⅱ) & (ⅲ) on positivity and boundedness

Proof. The solution to Equation (6) is given by

$C(t) = \exp(-\psi t)\left(\int_{0}^t g(s)\exp(\psi s)ds + C(0)\right) \geq 0.$

Let the model equations (1)-(5) be written in the form $x^{\prime} = F(t,x)$. The functions $F(x,t)$ on the right hand side of Equations (1)-(5) have the property of

$F(U,I,V,E_V,E_T,t)\geq 0~\text{whenever}~ x\in [0,\infty )^n,~ x_j = 0,~t\geq 0,$

where $x = (U,I,V,E_V,E_T)$. Thus, it follows from Proposition A.1 in [46] that

$x(t)\in [0,\infty)^n ~\text{for all} ~t\geq t_0 \geq 0~~\text{whenever}~~ x(t_0) \geq 0$

and thus the model solutions $x = (U,I,V,E_V,E_T)$ are positive and bounded for positive initial values $x^0$.

A.3. Proof of Theorem 3.1(ⅳ) on positive invariance

Proof. The proof directly follows from proofs of Theorems 3.1 (ⅰ)-(ⅲ).

---

Appendix B

B.1. Proof of stability of tumour endemic state, $X_2$, in Proposition 1.

Proof. Denote $U^* : = U_p +U_r$ where

$U_p = -\frac{b}{2a} = \frac{K}{2\alpha}\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg),\\ U_r = \frac{\sqrt{b^2-4ac}}{2a} = \frac{K}{2\alpha}\sqrt{\bigg(\alpha - \frac{\alpha \kappa}{K}-M \bigg)^2 +4\frac{\alpha^2\kappa}{K}}.$

The characteristic polynomial of the Jacobian matrix evaluated at $X_2$ is given by

$\begin{align} \label{eq:char1:b} P(x) = x^2 + P_1x + P_0 \end{align}$

(31)

where $P_0 = P_{01}U^* + P_{00}$ and $P_1 = P_{11}U^* + P_{10}$,

$P_{00} = K\kappa\, \left( \alpha\,\delta_{T}\,K-\beta_{T}\,\nu_{u}\,K+\alpha\, \kappa\,\delta_{T} \right), \\ P_{01} = {\frac {K \big( K(\alpha \delta_T - \beta_T\nu_U)^2+\kappa\,{\delta_{T}} ^{2}{\alpha}^{2}+\alpha\,\beta_{T}\,\kappa\,\nu_{u}\,\delta_{T} \big) }{\alpha\,\delta_{T}}}, \\ P_{10} = K\kappa\,\delta_{T}\, \left( \alpha+\delta_{T} \right), \\ P_{11} = K \left( \alpha\,\delta_{T}-\beta_{T}\,\nu_{u}+{\delta_{T}}^{2} \right).$

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

$\begin{align*} U_r > -\frac{P_{00}}{P_{01}}-U_p \Leftrightarrow U_r ^2 - \left( \frac{P_{00}}{P_{01}} +U_p\right)^2 > 0. \end{align*}$

With the help of Maple, we find that

$U_r^2 - \left( \frac{P_{00}}{P_{01}} +U_p\right)^2 = \alpha\,\beta_{T}\,{\kappa}^{2}\nu_{u}\,\delta_{T}\, \bigg( K^2 (\alpha \delta_T - \beta_T\nu_U)^2 +2\,K\kappa\,{\delta_{T}}^{2}{\alpha}^{2}\\ \;\;\;\;\; \;\;\;\;\; \;\;\;\;\;\;\;\;\;\; \;\;\;\;\; \;\;\;\;\;\;\;\; +2\,K\alpha\,\beta_{T}\,\kappa\,\nu_{u}\,\delta_{T}+{ \alpha}^{2}{\kappa}^{2}{\delta_{T}}^{2} \bigg) > 0.$

Therefore $P_0 >0$.

For $P_1$, we discuss the following cases:

(ⅰ) If $P_{11}> 0$ then $P_1>0$.

(ⅱ) If $P_{11} <0$ then $P_1>0$ if and only if $U^* < -P_{10}/P_{11}$ which is equivalent to

$\begin{align*} U_r^2 - \left(\frac{P_{10}}{P_{11}} + U_p \right)^2 < 0. \end{align*}$

Using Maple, we find that

$U_r^2 - \left(\frac{P_{10}}{P_{11}} +U_p \right)^2 = \beta_{T}\,\kappa\,\nu_{u}\,\delta_{T}\, \bigg( K(\alpha\,\delta_{T}\, - \beta_{T}\,\nu_{u}\,+{\delta_{T}}^{2})-\alpha \kappa(\alpha+\delta_{T}) \bigg) <0.$

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

$P(x) = \newcommand{\Bold}[1]{\bf{#1}}x^{3} + \left(\frac{2 \, U \alpha}{K} - \alpha + \delta_{T} + \frac{C \delta_{U}}{C + K_{c}} + \frac{U \beta_{T} \nu_{U}}{{\left(U + \kappa\right)} \delta_{T}} + \psi\right) x^{2}$

(32)

$\;\;\;\;\; +\bigg(\frac{2 \, U \alpha \delta_{T}}{K} - \alpha \delta_{T} + \frac{C \delta_{T} \delta_{U}}{C + K_{c}} + \frac{2 \, U \beta_{T} \nu_{U}}{U + \kappa} - \frac{U^{2} \beta_{T} \nu_{U}}{{\left(U + \kappa\right)}^{2}}$

(33)

$\;\;\;\;\; + \frac{2 \, U \alpha \psi}{K} - \alpha \psi + \delta_{T} \psi + \frac{C \delta_{U} \psi}{C + K_{c}} + \frac{U \beta_{T} \nu_{U} \psi}{{\left(U + \kappa\right)} \delta_{T}}\bigg) x$

(34)

$\;\;\;\;\;+ \frac{2 \, U \alpha \delta_{T} \psi}{K} - \alpha \delta_{T} \psi + \frac{C \delta_{T} \delta_{U} \psi}{C + K_{c}} + \frac{2 \, U \beta_{T} \nu_{U} \psi}{U + \kappa} - \frac{U^{2} \beta_{T} \nu_{U} \psi}{{\left(U + \kappa\right)}^{2}} = 0.$

(35)

The conditions for stability of the endemic state, $C_2$, therefore follow from the Routh-Hurwitz criterion and are stated as:

$a_2: = \frac{2 \, U^* \alpha}{K} - \alpha + \delta_{T} + \frac{C^* \delta_{U}}{C^* + K_{c}} + \frac{U^* \beta_{T} \nu_{U}}{{\left(U^* + \kappa\right)} \delta_{T}} + \psi > 0, \nonumber \\ a_1 : = \frac{2 \, U^* \alpha \delta_{T}}{K} - \alpha \delta_{T} + \frac{C^* \delta_{T} \delta_{U}}{C^* + K_{c}} + \frac{2 \, U^* \beta_{T} \nu_{U}}{U^* + \kappa} - \frac{U^{*^2} \beta_{T} \nu_{U}}{{\left(U^* + \kappa\right)}^{2}} \nonumber \\ \;\;\;\;\;\;\;\;\;\;\;\; + \frac{2 \, U \alpha \psi}{K} - \alpha \psi + \delta_{T} \psi + \frac{C^* \delta_{U} \psi}{C^* + K_{c}} + \frac{U^* \beta_{T} \nu_{U} \psi}{{\left(U^* + \kappa\right)} \delta_{T}} > 0, \nonumber \\ a_0: = \frac{2 \, U^* \alpha \delta_{T} \psi}{K} - \alpha \delta_{T} \psi + \frac{C^* \delta_{T} \delta_{U} \psi}{C^* + K_{c}} + \frac{2 \, U^* \beta_{T} \nu_{U} \psi}{U^* + \kappa} - \frac{U^{*^2} \beta_{T} \nu_{U} \psi}{{\left(U^* + \kappa\right)}^{2}} > 0 \nonumber \\ \;\;\;\;\;\;\;\; \text{and} ~ a_2a_1>a_0.$

(36)

---

Appendix C

C.1. Proof of conditions $A_1 -A_5$ for the derivation of $R_0$

Proof. Denote $x = (x_1,x_2,\cdots x_5)^T = (U,I,V,E_V,E_T)^T$ with $x_i \geq 0$ be the set of cell compartments divided into infected cells $i = 1,\cdots,m$ and uninfected cells $i = m+1,\cdots, n$. Define $X^0$ to be the set of virus or infection free states:

$X^0: = \{x\geq 0: x_i 0,~ \text{for all} ~ i = 1,\cdots, m \}.$

Let $\mathcal{F}_i (t,x)$ be the input rate of newly infected cells in the $i^{\text{th}}$ compartment. Let $\mathcal{V}_i^+(t,x)$ and $\mathcal{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

$\frac{dx_i}{dt} = \mathcal{F}_i(t,x)-\mathcal{V}_i (t,x).$

We now proceed to prove the assumptions $A_1 -A_5$:

$A_1$: For each $1\leq i \leq n$, the functions $\mathcal{F}_i(t,x)$, $\mathcal{V}_i^+$ and $\mathcal{V}_i^-$ are non-negative and continuous on $\mathbb{R} \times \mathbb{R}^n_+$ and continuously differentiable with respect to $x$. From Equations (1)-(6)

$\label{eq:fv} \mathcal{F}_i = \left( \begin{array}{c} \frac{\beta U V}{K_u + U } \\ 0 \end{array} \right), \;\;\mathcal{V}_i^+ = \left( \begin{array}{c} \delta I + \nu_I E_T I + \frac{\delta_I I C}{K_c +C} \\ \gamma v \end{array} \right), \;\;\mathcal{V}_i^- = \left( \begin{array}{c} \tau E_V E_V I \\ b \delta \end{array} \right).$

(37)

Clearly, $\mathcal{F}_i$, $V_{\left. {} \right\rangle i}^ +$ and $\mathcal{V}_i^-$ satisfy $A_1$.

$A_2$: If $x_i = 0$, then $\mathcal{V}_i = 0$. Clearly, from (37), if the state variables are all equal to zero, then $\mathcal{V}_i = 0$.

$A_3$: $\mathcal{F}_i = 0$ for $i>m$. This can be seen from the compartments $U$, $E_V$ and $E_T$.

$A_4$: If $x\in X^0$, then $\mathcal{F}_i(x) = \mathcal{V}^+_i (x) = 0$ for $i = 1,\cdots , m$. At $X^0$, $V$ and $I$ are zero, implying that $\mathcal{F}_i = \mathcal{V}^+_i = 0$.

$A_5$: If $\mathcal{F}_i = 0$, then all the eigenvalues of $D \mathcal{V} (X^0)$ have positive real parts. This can be ascertained from Equation (37).

C.2. Proof of Theorem 4.1

Proof. From Equation (37),

$F = \left( \begin{array}{cc} 0& \frac{\beta U}{K_u + U } \\ 0 &0 \end{array} \right), \;\; V = \left( \begin{array}{cc} \delta +\nu_IE_T + \frac{\delta_I C^*}{K_c + C}+\tau E_v&0 \\ b\delta &\gamma \end{array} \right),$

$~ FV^{-1} = \frac{1}{\gamma(\delta + \nu_I E_T^* +\frac{\delta_IC^*}{K_c + C^*})}\left( \begin{array}{cc} \frac{b \delta \beta U}{K_U + U } & \frac{b\beta U^*}{(K_u + U^*)} \\ 0&0\end{array} \right).$

Thus

$\begin{align} \label{eq:Ro-proof} R_0 = \rho(FV^{-1}) = \frac{b\beta \delta U^*(K_c + C^*)}{\gamma \big[(K_c+ C^*) (\nu_I E_T^* +\delta) + \delta_I C^* \big]\big(K_u + U^* \big)}. \end{align}$

(38)

---

Appendix D

D.1. Elasticity indices of ${R_0}$

$\begin{align} \begin{split} e_b & = 1\\ e_{\beta} & = 1\\ e_{\gamma} & = -1\\ e_{\delta} & = \newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\Gamma_1\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma}{{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K b \beta}}\\ e_{\alpha} & = \newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\Gamma_2\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \alpha \gamma}{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K b \beta \delta}\\ e_{\delta_U} & = \newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\Gamma_3\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \delta_{U} \gamma}{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K b \beta \delta}\\ e_{\delta_I} & = \newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\delta_{I} q}{K_{c} \delta \psi + \delta q + \delta_{I} q}\\ e_{K_c} & = \newcommand{\Bold}[1]{\mathbf{#1}}\frac{\Gamma_4 K_{c} \psi q}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} {\left(K_{c} \psi + q\right)}} \\ e_{K_u} & = \newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(K_{c} \psi + q\right)} K_{u} \alpha}{K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q}\\ e_{\psi} & = \newcommand{\Bold}[1]{\mathbf{#1}}\frac{\Gamma_4 K_{c} \psi q}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} {\left(K_{c} \psi + q\right)}} \\ e_q & = \newcommand{\Bold}[1]{\mathbf{#1}}-\frac{\Gamma_4 K_{c} \psi q}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} {\left(K_{c} \psi + q\right)}}, \end{split} \end{align}$

(39)

where

$\Gamma_1 = \frac{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)}^{2} K b \beta \delta}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)}^{2} \gamma}\\ \;\;\;\;\;\;\;\;\;\;\;\;- \frac{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K b \beta}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma}\\ \Gamma_2 = \frac{{\left(K K_{c} \psi + K_{c} K_{u} \psi + K q + K_{u} q\right)} {\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K b \beta \delta}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)}^{2} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma}\\ \;\;\;\;\;\;\;\;\;\;\;\;- \frac{{\left(K_{c} \psi + q\right)}^{2} K b \beta \delta}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma} \\ \Gamma_3 = \frac{{\left(K_{c} \alpha \psi + \alpha q - \delta_{U} q\right)} {\left(K_{c} \psi + q\right)} K^{2} b \beta \delta q}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)}^{2} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma}\\ \;\;\;\;\;\; \;\;\;\;\;\;- \frac{{\left(K_{c} \psi + q\right)} K b \beta \delta q}{{\left(K K_{c} \alpha \psi + K_{c} K_{u} \alpha \psi + K \alpha q + K_{u} \alpha q - K \delta_{U} q\right)} {\left(K_{c} \delta \psi + \delta q + \delta_{I} q\right)} \gamma}\\ \Gamma_4 = K K_{c}^{2} \alpha^{2} \delta_{I} \psi^{2} + K_{c}^{2} K_{u} \alpha^{2} \delta_{I} \psi^{2} + K_{c}^{2} K_{u} \alpha \delta \delta_{U} \psi^{2} + 2 \, K K_{c} \alpha^{2} \delta_{I} \psi q \\ \;\;\;\;\;\;\;\;\;\;\;\;+ 2 \, K_{c} K_{u} \alpha^{2} \delta_{I} \psi q + 2 \, K_{c} K_{u} \alpha \delta \delta_{U} \psi q - 2 \, K K_{c} \alpha \delta_{I} \delta_{U} \psi q + K \alpha^{2} \delta_{I} q^{2} \\ \;\;\;\;\;\; \;\;\;\;\;\;+ K_{u} \alpha^{2} \delta_{I} q^{2} + K_{u} \alpha \delta \delta_{U} q^{2} - 2 \, K \alpha \delta_{I} \delta_{U} q^{2} + K \delta_{I} \delta_{U}^{2} q^{2}.$

D.2. Elasticity indices of endemic equilibrium.

Consider a system of differential equations dependent on a parameter $p$

$\begin{align} x^{\prime} = f(x,p) \end{align}$

(40)

where

$f(x,p) = \left( \begin{array}{c} f_1(x,p) \\ f_2(x,p) \\ .\\ .\\ .\\ f_n(x,p) \end{array}\right), \;\;x = \left( \begin{array}{c} x_1 \\ x_2 \\ .\\ .\\ .\\ x_n \end{array}\right).$

At equilibrium point $x^* = x^*(p)$ satisfies $f(x^*(p),p) = 0$. Differentiating with respect to $p$, one obtains

$\;\;\;\frac{\partial f}{\partial x_1} \left(x^*(p),p \right)x^{\prime}_1(p) + \frac{\partial f}{\partial x_2} \left(x^*(p),p \right)x^{\prime}_2(p)\\ +\cdots + \frac{\partial f}{\partial x_n} \left(x^*(p),p \right)x^{\prime}_n(p) + \frac{\partial }{\partial p} \left(x^*(p),p \right) = 0.$

This implies that

$\begin{align*} \left( \begin{array}{cccccc} \frac{\partial f_1}{\partial x_1}&\frac{\partial f_1}{\partial x_1} & . &. &.&\frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_2} &\frac{\partial f_2}{\partial x_2} & .&.&.&\frac{\partial f_2}{\partial x_n} \\ .&.&.&. &.&.\\ .&.&.&. &.&.\\ .&.&.&. &.&. \\ \frac{\partial f_1}{\partial x_n}&\frac{\partial f_n}{\partial x_2} & . &. &.&\frac{\partial f_n}{\partial x_n} \end{array}\right) \left( \begin{array}{c} x^{\prime}_1 (p) \\ x^{\prime}_2 (p)\\ .\\ .\\ .\\ x^{\prime}_n (p) \end{array}\right) = -\left( \begin{array}{c} \frac{\partial f_1}{\partial p} \\ \frac{\partial f_2}{\partial p}\\ .\\ .\\ .\\ \frac{\partial f_n}{\partial p} \end{array}\right) \end{align*}$

$\begin{align*} \left( \begin{array}{c} x^{\prime}_1 (p) \\ x^{\prime}_2 (p)\\ .\\ .\\ .\\ x^{\prime}_n (p) \end{array}\right) = -J^{-1}_{x^*_p} \left( \begin{array}{c} \frac{\partial f_1}{\partial p} \\ \frac{\partial f_2}{\partial p}\\ .\\ .\\ .\\ \frac{\partial f_n}{\partial p} \end{array}\right). \end{align*}$

Multiplying both sides with the diagonal matrix, $diag(p/x^*_i)$, matrix

$\begin{align*} \begin{pmatrix} \frac{p}{x^*_1}&& \\ & \ddots&\\ &&\frac{p}{x^*_n} \end{pmatrix} \left( \begin{array}{c} x^{\prime}_1 (p) \\ x^{\prime}_2 (p)\\ .\\ .\\ .\\ x^{\prime}_n (p) \end{array}\right) = -\begin{pmatrix} \frac{p}{x^*_1}&& \\ & \ddots&\\ &&\frac{p}{x^*_n} \end{pmatrix} J^{-1} \left( \begin{array}{c} \frac{\partial f_1}{\partial p} \\ \frac{\partial f_2}{\partial p}\\ .\\ .\\ .\\ \frac{\partial f_n}{\partial p} \end{array}\right). \end{align*}$

The sensitivity indices of the steady state variables at equilibrium $x^*_i$ with respect to a parameter $p$ are therefore given by

$\begin{align} \label{eq:derived-sensitivity} \left( \begin{array}{c} \Gamma^{x^*_1}_p \\ \Gamma^{x^*_2}_p \\ .\\ .\\ .\\ \Gamma^{x^*_n}_p \end{array}\right) = -K J^{-1} \frac{\partial f}{\partial p} \left(x^*(p),p \right). \end{align}$

(41)

---

Acknowledgments

Joseph Malinzi was jointly supported by the University of Pretoria and DST/NRF SARChI Chair in Mathematical Models and Methods in Bioengineering and Biosciences.

Amina Eladdadi and K.A. Jane White would like to acknowledge and thank the UK-QSP Network (Grant-EP/N005481/1) for their financial support to attend the QSP-1st Problem Workshop and collaborate on this research.

Torres was supported by FCT through CIDMA, project UID/MAT/04106/2013, and TOCCATA, project PTDC/EEI-AUT/2933/2014, funded by FEDER and\break COMPETE 2020.

---