# Mathematical Analysis of a Mathematical Model of Chemovirotherapy: Effect of Drug Infusion Method

Authors: Joseph Malinzi
Journal: Computational and Mathematical Methods in Medicine
 (2019)

---

## Abstract

A mathematical model for the treatment of cancer using chemovirotherapy is developed with the aim of determining the efficacy of three drug infusion methods: constant, single bolus, and periodic treatments. The model is in the form of ODEs and is further extended into DDEs to account for delays as a result of the infection of tumor cells by the virus and chemotherapeutic drug responses. Analysis of the model is carried out for each of the three drug infusion methods. Analytic solutions are determined where possible and stability analysis of both steady state solutions for the ODEs and DDEs is presented. The results indicate that constant and periodic drug infusion methods are more efficient compared to a single bolus injection. Numerical simulations show that with a large virus burst size, irrespective of the drug infusion method, chemovirotherapy is highly effective compared to either treatments. The simulations further show that both delays increase the period within which a tumor can be cleared from body tissue.


---

## Body

## 1. Introduction

Tumors possess mechanisms that suppress antitumor activity such as ligands that block natural killer cells and cytotoxic tumor infiltrating cell functions [1]. Greatly because of this, successful cancer treatment often requires a combination of treatment regimens.Nearly all traditional monotherapies, including chemotherapy, surgery, and radiation therapy are not a definite cure for cancer and are highly toxic [2]. Chemotherapy for example, which is the most commonly used regimen, involves the use of medical drugs to lyse cancer cells. These chemotherapeutic drugs circulate in the body and kill rapidly multiplying cells nonselectively, which ultimately results into the destruction of both healthy and cancerous cells [2, 3]. Chemotherapy can thus be toxic to a patient with adverse side effects and can also damage their immune system [2].Presently, combination cancer treatment is a centerpiece of cancer therapy [4]. The amalgamation of anticancer drugs increases efficacy compared to single-drug treatments. Further, anticancer drug combination provides therapeutic benefits such as reducing tumor growth, arresting mitotically active cells, reducing the population of cancer stem cells, and inducing apoptosis [4]. Despite the fact that combination therapy might as well be toxic if one of the agents used is chemotherapeutic, the toxicity is lesser because different pathways would be targeted [4]. Moreover with the use of combination therapy, the toxicity on normal cells can be prevented while concurrently producing cytotoxic effects on cancer cells [4, 5].In the recent past, virotherapy, a less toxic treatment has been identified as a possible cancer remedy [6–11]. Virotherapy involves the use of oncolytic viruses that infect, multiply, and directly lyse cancer cells with less or no toxicity [9]. Their tumor specific properties allow for viral binding, entry, and replication [12]. Oncolytic viruses can greatly enhance the cytotoxic mechanisms of chemotherapeutic drugs [13]. Further, chemotherapeutic drugs lyse fast multiplying cells and, in general, virus infected tumor cells quickly replicate [14].Chemovirotherapy is a combination treatment strategy that involves the use of oncolytic viruses and chemotherapeutic drugs. Recent experimental and mathematical studies have shown that chemovirotherapy is a plausible cancer treatment and leads to enhanced therapeutic effects not achievable when either therapies are independently used [12, 13, 15–20]. Nguyen et al. [12] gave an account of the mechanisms through which drugs can successfully be used in a combination with oncolytic viruses. They however note that the success of this combination depends on several factors including the type of oncolytic virus- (OV-) drug combination used, the timing, frequency, dosage, and cancer type targeted. To date, the best method of OV drug delivery is debatable [21, 22].The main goal of this study is to, thus, consider and compare the efficacy of three drug infusion methods, use mathematical analysis to predict the outcome of OV-drugs combination treatment and determine the effect of drug response and virus infection delays. To this end, we construct a mathematical model in the form of ODEs which we later extend to DDEs to include the virus infection and drug response delays. The model constructed combines elements from existing mathematical models [10, 11, 19, 20, 23–30]. Tian [10] presented a mathematical model that incorporates burst size for oncolytic virotherapy. His study showed that virotherapy is highly effective provided that viruses with large burst sizes are used. Malinzi et al. [19] constructed a spatiotemporal mathematical model to investigate the outcome of chemovirotherapy. Their study suggested that combining chemotherapeutic drugs with oncolytic viruses is more efficient than using either treatments alone. A similar study by Malinzi et al. [20] indicates that chemotherapy alone is capable of clearing tumor cells from body tissue if the drug efficacy is greater than the tumor growth rate. Nevertheless, the study contends that oncolytic viruses highly enhance chemotherapy in lysing tumor cells. The study further postulates that half the maximum tolerated doses of chemotherapy and virotherapy optimize chemovirotherapy, thus answering a very pertinent question in combination cancer therapy.The article is organised as follows: Section2 presents a comprehensive description of the both the ODE and DDE models and the underlying assumptions made in constructing them. In Section 3, the model without delay is analysed. First, without any form of treatment, then with either treatments (that is, with chemo only and virotherapy alone) and with both treatments. The delay model is then analysed in Section 4 and numerical experiments for both the ODE and DDE models are carried out in Section 5. Finally, before concluding in Section 7, a comparison of this study with related works is done is Section 6.


## 2. Model Description

### 2.1. Model without Delay

Time-dependent cell concentrations of uninfected tumor cellsUt, infected tumor cells It, a free virus population Vt, and a chemotherapeutic drug Ct in an avascular tumor localization are considered. The uninfected tumor grows logistically at an intrinsic rate α per day, and the total tumor carrying capacity is K cells in a tumor nodule. The infected tumor cells die off at a rate δ per day. Virus multiplication in the tumor is represented by the function βUtVt, where β is the virus replication rate measured per day per 106 cells or viruses. The response of the drug to the uninfected and infected tumor is, respectively, modelled by the functions δ0UtCt and δ1ItCt where δ0 and δ1 are induced lysis rates caused by the chemotherapeutic drug measured per day per cell. Virus lifespan is taken to be 1/γ and its production is considered to be bδI where b is the virus burst size, measured in number of viruses per day per cell, and δ is the infected tumor cells’ death rate measured per day. Chemotherapeutic drug infusion into the body is modelled with a function gt and that the drug gets depleted from body tissue at a rate λ per day.Drug infusion into the body is simulated using (a) a constant rategt=q, (b) an exponential gt=qexp−at, and (c) a sinusoidal function gt=qsin2at, where q is the rate of drug infusion. The constant a determines the exponential drug decay and period for the sinusoidal infusion. Constant drug infusion 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 [31, 32]. The exponential drug infusion may relate a situation where a cancer patient is given a single bolus and the drug exponentially decays in the body tissue. This form of infusion is not common although it is now used for some drugs, for example, a single dose of carboplatin can be given to patients with testicular germ cell tumors and breast cancer ([33, 34]). The third scenario is possible when a cancer patient makes several visits to a health facility and is given injections or anticancer drugs periodically [35, 36].The assumptions above lead to the following system of nonlinear first-order differential equations (also similarly derived in [11, 19, 20]):(1)U˙t=αUt1−Ut+ItK−βUtVt−δ0UtCt,I˙t=βUtVt−δIt−δ1ItCt,V˙t=bδIt−βUtVt−γVt,C˙t=gt−λCt,subject to initial concentrations(2)U0=U0,I0=I0,V0=V0,C0=C0.


### 2.2. Delay Model

The model is further extended to account for delays as a result of the infection of tumor cells by the virus and responses of the chemotherapeutic drug. In fact, the viruses need time to develop suitable responses when they meet the uninfected tumor cells (e.g., [37]). The drug does not instantaneously kill the cells (e.g., [38, 39]). By denoting the virus and chemotherapeutic response delays as τ1 and τ2, respectively, model (1)(3)U˙t=αUt1−Ut+tK−βUt−τ1Vt−τ1−δ0Ut−τ2Ct−τ2,I˙t=βUt−τ1Vt−τ1−It−δ1It−τ2Ct−τ2,V˙t=bIt−βUt−τ1Vt−τ1−γVt,C˙t=gt−λCt.


## 2.1. Model without Delay

Time-dependent cell concentrations of uninfected tumor cellsUt, infected tumor cells It, a free virus population Vt, and a chemotherapeutic drug Ct in an avascular tumor localization are considered. The uninfected tumor grows logistically at an intrinsic rate α per day, and the total tumor carrying capacity is K cells in a tumor nodule. The infected tumor cells die off at a rate δ per day. Virus multiplication in the tumor is represented by the function βUtVt, where β is the virus replication rate measured per day per 106 cells or viruses. The response of the drug to the uninfected and infected tumor is, respectively, modelled by the functions δ0UtCt and δ1ItCt where δ0 and δ1 are induced lysis rates caused by the chemotherapeutic drug measured per day per cell. Virus lifespan is taken to be 1/γ and its production is considered to be bδI where b is the virus burst size, measured in number of viruses per day per cell, and δ is the infected tumor cells’ death rate measured per day. Chemotherapeutic drug infusion into the body is modelled with a function gt and that the drug gets depleted from body tissue at a rate λ per day.Drug infusion into the body is simulated using (a) a constant rategt=q, (b) an exponential gt=qexp−at, and (c) a sinusoidal function gt=qsin2at, where q is the rate of drug infusion. The constant a determines the exponential drug decay and period for the sinusoidal infusion. Constant drug infusion 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 [31, 32]. The exponential drug infusion may relate a situation where a cancer patient is given a single bolus and the drug exponentially decays in the body tissue. This form of infusion is not common although it is now used for some drugs, for example, a single dose of carboplatin can be given to patients with testicular germ cell tumors and breast cancer ([33, 34]). The third scenario is possible when a cancer patient makes several visits to a health facility and is given injections or anticancer drugs periodically [35, 36].The assumptions above lead to the following system of nonlinear first-order differential equations (also similarly derived in [11, 19, 20]):(1)U˙t=αUt1−Ut+ItK−βUtVt−δ0UtCt,I˙t=βUtVt−δIt−δ1ItCt,V˙t=bδIt−βUtVt−γVt,C˙t=gt−λCt,subject to initial concentrations(2)U0=U0,I0=I0,V0=V0,C0=C0.


## 2.2. Delay Model

The model is further extended to account for delays as a result of the infection of tumor cells by the virus and responses of the chemotherapeutic drug. In fact, the viruses need time to develop suitable responses when they meet the uninfected tumor cells (e.g., [37]). The drug does not instantaneously kill the cells (e.g., [38, 39]). By denoting the virus and chemotherapeutic response delays as τ1 and τ2, respectively, model (1)(3)U˙t=αUt1−Ut+tK−βUt−τ1Vt−τ1−δ0Ut−τ2Ct−τ2,I˙t=βUt−τ1Vt−τ1−It−δ1It−τ2Ct−τ2,V˙t=bIt−βUt−τ1Vt−τ1−γVt,C˙t=gt−λCt.


## 3. Mathematical Analysis of Model without Delay

In this section, the model without delay (1) is analysed. The variables in system (1) are first rescaled by setting t¯=δt, U=KU¯, I=KI¯, V=V0V¯, and C=C0C¯. Taking V0=K, the parameters are renamed to become(4)α¯=αδ,β¯=βV0δ,δ0¯=δ0C0δ,δ1¯=δ1C0δ,b¯=bKV0,γ¯=γδ,ϕ=qδC0,ψ=λδ,a¯=aδ.For simplicity, we drop the bars and equation (1) becomes(5)U˙t=αUt1−Ut−It−βUtVt−δ0UtCt,I˙t=βUtVt−It−δ1ItCt,V˙t=bIt−βUtVt−γVt,C˙t=ξt−ψCt.ξ
t
=
ϕ, ξt=ϕexp−at, and ϕsin2at, respectively, are the constant, exponential, and sinusoidal infusion functions. For this model to be biologically meaningful, its solutions should be positive and bounded because they represent concentrations. Well-posedness theorems of model (5) are stated and proved in Appendix A.


### 3.1. Model Solutions

To investigate the efficacy of each treatment and their combination, we first study the dynamics of the system without treatment. Without any form of treatment, model (5) is reduced to only one equation:(6)U˙t=αUt1−Ut,U0=U0,whose solution is(7)Ut=U01−U0exp−αt+U0,limt⟶∞Ut=1,implying that the tumor logistically grows to its maximum fractional size. Next, the model (5) is analysed, with chemotherapy, with virotherapy, and then with both treatments incorporated. We obtain, where possible analytical and time invariant solutions which predict the long term dynamics of the model equation (1).Without virotherapy (Vt=0), the system (5) is transformed to(8)U˙t=αUt1−Ut−δ0UtCt,C˙t=ξt−ψCt,with U0=U0 and C0=C0. The second equation in equation (8) is a first-order linear ordinary differential equation which can easily be solved to give(9)Ct=exp−ψt∫ξtexpψtdt+R,where R is a constant of integration. The solution to the first equation in equation (8) depends on the infusion function ξt. For a fixed infusion function ϕ,(10)Ut=eαt−δ0ϕt/ψ+δ0e−ψt/ψα∫eαt−δ0ϕt/ψ+δ0e−ψt/ψdt+eδ0ψ/U0,Ct=C0−ϕψe−ψt+ϕψ.From the solution ofCt in equation (10),(11)limt⟶∞Ct=ϕψ.Biologically, it can be inferred that with a constant drug infusion and without virotherapy, the tumor is not completely cleared and a certain proportion of the drug remains in body tissue. The tumor clearance depends on the drug induced lysis of the tumor and the drug infusion rate which should be maximized and the tumor growth and drug decay rate which should be minimized.Forξt=ϕexp−at,(12)Ut=eaαt/a−ψ−αψt/a−ψ+cδ0e−ψt/a−ψ−acδ0e−ψt/a−ψψ+δ0ϕe−at/a−ψaa∫eaαt/a−ψ−αψt/a−ψ+cδ0e−ψt/a−ψ−acδ0e−ψt/a−ψψ+δ0ϕe−at/a−ψadt+ecδ0/a−ψ/U0,Ct=C0ψ−ϕψ−ϕe−at+ψta−ψe−ψt,where c=C0ψ−ϕ/ψ.From equation (12),(13)limt⟶∞Ct=0,limt⟶∞Ut=U∗,where U∗ is a fractional tumor cell concentration between 0 and 1. This suggests that with a single dosage infusion of the chemotherapeutic drug with exponential decay and without virotherapy, the tumor cannot be cleared from body tissue. The drug is also completely depleted from the body.Whenξt=ϕsin2at is substituted in equation (8), the resulting differential equations are solved to give(14)Ct=122R−ψ2cos2ateψt+2aψeψtsin2at−4a2+ψ2eψtϕ4a2ψ+ψ3e−ψt,where(15)R=C0+ψ2−ϕ4a2+ψ22ψ4a2+ψ2,limt⟶∞Ct=C∗=fa,ϕ,ψ.This suggests that with time, some drug concentration remains in the body tissue.Theorem 1.
The system (8), with constant infusion, has no periodic solutions for positive Ut and Ct.Proof.
Using Dulac’s criterion ([40]):
SupposeX˙=fx and fx is continuously differentiable on a simply connected domain D⊂R. If there exists a real valued function gx such that ∇⋅gX˙=∇⋅gf has one sign in D, then there are no closed orbits in D. Using Dulac’s criterion, it is sufficient to show that(16)∂∂UgU˙+∂∂CgC˙≠0,∀U,C∈R+2.
Consider(17)gU,C=1UC,∇⋅gX˙=∂∂UgU˙+∂∂CgC˙,=−αC+ξtUC2<0,∀U,C∈R+2.Theorem 2.
The system (8) has at least two steady states for each of the drug infusion functions: (1)
For the constant drug infusion functionξt=ϕ, there are two steady states of equation (8): U=0,C=ϕ/ψ which is locally asymptotically stable provided that δ0ϕ>αψ and U=1−δ0ϕ/αψ,C=0 which is locally asymptotically stable provided that δ0ϕ+αψ>2δ0ϕψ2; otherwise, it is unstable.(2)
For the exponential drug infusionξt=ϕexp−at, equation (8) has two steady states: U=0,C=0,W=0 which is unstable and U=1,C=0,W=0 which is locally asymptotically stable.(3)
For the sinusoidal infusion function, there are four steady states of equation (8): U=0,C=0,W=0, U=1,C=1,W=0, and U=aαψ−δ0ϕ/aαψ,C=ϕ/aψ,W=ϕ/a which are unstable and U=0,C=ϕ/aψ,W=ϕ/a which is locally asymptotically stable if δ0ϕ>aαψ and ϕ<1.Proof.
(1)
It is easy to show that when equation (8) is equated to zero, one obtains two steady states. The characteristic polynomial of the Jacobian matrix for equation (8) evaluated at 0,ϕ/ψ is(18)λ2+−α+δ0ϕψ+ψλ+δ0ϕ−αψ,whose roots λ can only be negative if δ0ϕ>αψ.
The characteristic polynomial evaluated at1−δ0ϕ/αψ,0 is(19)λ2+−2δ0ϕψ+α+δ0ϕψ+ψλ−2δ0ϕψ2+δ0ϕ+αψ,whose roots are negative provided that δ0ϕ+αψ>2δ0ϕψ2.
(2)
By lettingW=ϕexp−at, equation (8) is turned into an autonomous system:(20)U˙t=αUt1−Ut−δ0UtCt,C˙t=W−ψCt,W˙t=−aW.
The eigenvalues of the Jacobian matrix for equation (20) evaluated at 0,0,0 which are −ψ,α, and 0 and at 1,0,0 which are −α,−ψ, and −a are all negative.
(3)
Similarly, by lettingW=ϕsin2t, equation (8) becomes the autonomous system:(21)U˙t=αUt1−Ut−δ0UtCt,C˙t=W−ψCt,W˙t=4aWϕ−W1/2.
The eigenvalues of the Jacobian matrix for equation (20) evaluated at 0,0,0 are −ψ,α, and 0 and the eigenvalues evaluated at 1,0,0 are ψ,−α, and 0. For the third steady state to exist, aαψ≥δ0ϕ, and the eigenvalues evaluated at this state are −aαψ−δ0ϕ/aψ,−4aϕ/4ϕ1−ϕ,−ψ, implying that for it to be locally asymptotically stable, δ0ϕ>aαψ; yet for this to happen, the steady state will not exist. The eigenvalues evaluated at U=0,C=ϕ/aψ,W=ϕ/a are −aαψ−δ0ϕ/aψ,−4aϕ/4ϕ1−ϕ,−ψ implying that this steady state is locally asymptotically stable if δ0ϕ>aαψ and ϕ<1.
Theorems 1 and 2 show that there are no periodic solutions in the dynamics of equation (8) and with a constant drug infusion, the tumor can be eliminated from body tissue by chemotherapy provided that the combination of the chemotherapeutic drug-induced lysis of the tumor and the drug infusion is greater than the combination of the intrinsic tumor growth rate and the drug deactivation rate. The tumor can also be wiped out with a periodic drug infusion provided that the combination of the tumor-induced lysis by the drug and the dosage is greater than the intrinsic tumor growth rate and drug decay rate. With the exponential infusion method, the tumor is not removed from body tissue and may grow to its maximum size.
Without chemotherapy, equation (5) is reduced to(22)U˙t=αUt1−Ut−It−βUtVt,I˙t=βUtVt−It,V˙t=bIt−βUtVt−γVt.
The analytical solutions to system (22) are not easy to obtain. The derivatives of equation (22) are therefore equated to zero to obtain time invariant solutions and investigate their stability by linearizing equation (22) about the steady states.Theorem 3.
(1)
Ifβ+γ>bβ, the system (12) has two steady states: a tumor free cell state 0,0,0 which is unstable and an infection tumor free state 1,0,0 which is locally asymptotically stable.
(2)
Ifbβ>β+γ, the system (22) has three steady states: the tumor free state 0,0,0 and the infected free state 1,0,0 which are unstable and a tumor dormant state:(23)U=γb−1β,I=αγβb−1−γββb−12+γb−1,V=αβb−1−γb−12β2+αγ,b>1,which is locally asymptotically stable if a0,a1,a2>0 and a1a2>a0 where ai are coefficients of the characteristic equation.Proof.
(1)
The characteristic equation evaluated at0,0,0 is(24)λ3+γ−α+1λ2+γ−αγ−αλ−αγ=0,λ−αλ+1λ+γ=0,from which λ1=α, λ2=−1 and λ3=−γ, thus rendering it unstable.
The characteristic equation evaluated at1,0,0 is(25)λ+αλ2+λ1+β+γ+β+γ−bβ=0,from which λ1=−α and λ2+λ1+β+γ+β+γ−bβ=0 which are all negative since β+γ>bβ.
(2)
The characteristic polynomial evaluated at the tumor dormant state isλ3+a2λ2+a1λ+a0=0, where(26)a2=2Aα+Aβ+Cβ−α+γ+1,a1=2A2αβ−Aαβ+Cαβ−Abβ+2Aαγ+Cβγ+2Aα+Aβ+Cβ−αγ−α+γ,a0=−2A2αbβ+2A2αβ+Aαbβ+Cαβγ−Aαβ+2Aαγ+Cβγ−αγ,and A,B,C are the coordinates of the tumor dormant state. Using Routh–Hurwitz stability criterion, this state will only be locally asymptotically stable if a0,a1,a2>0 and a1a2>a0.Since the infected tumor-free state is undesirable, the reverse of the conditionβ+γ>bβ is necessary for tumor eradication from body tissue. In other words, bβ>β+γ, that is, the product of the virus replication rate and their burst size should be greater than the sum of the burst size and virus replication rate. We also notice from equation (23) that(27)limβ⟶∞U=limb⟶∞U=0.It is therefore evident that high virus replication rateβ and burst size b lead to lower tumor cell concentrations. The steady-state solutions of equation (23) involve many parameters, thereby giving rise to large expressions in the conditions for its stability. It is therefore a difficult undertaking to infer biological implications from these conditions. Nevertheless, it can be observed that virotherapy may only succeed in eliminating cancer from body tissue when the virus deactivation rate is very small or even zero and the virus replication rate very high.Next, the model with both treatments is analysed. For a constant drug infusion rateϕ, the system (5) has three steady states;(i)
Tumor-free steady state:(28)U=0,I=0,V=0,C=ϕψ.Here, the tumor and viruses are cleared from body tissue by the coupled treatment and a fraction of the chemotherapeutic drug remains in body tissue. The eigenvalues of the Jacobian matrix evaluated at this state are(29)λ1=−δ0ϕ−αψψ,λ2=−δ1ϕ+ψψ,λ3=−γ,λ4=−ψ,implying that this desirable state is locally asymptotically stable if δ0ϕ>αψ. From this condition, in order to clear a tumor, the combination of the rate at which the drug kills the uninfected tumor cells and the drug infusion must be higher than the tumor growth rate and deactivation of the drug from body tissue.(ii)
Infected tumor-free state:(30)U=1−δ0ϕαψ,I=0,V=0,C=ϕψ.
In this state, the whole tumor is not cleared as a fraction of uninfected tumor cells remain and all the infected ones are cleared by the treatment combination. Using the parameter values in Table1, the eigenvalues of the Jacobian matrix evaluated at the infected tumor-free state are 0.403, 8.13, and 2.598±2.418i, implying that the infected tumor-free state is unstable.
(iii)
Tumor dormant state:(31)U=δ1ϕ+ψγbψb−1−δ1ϕβ,I=Γγδ1ϕ−b−1ψβ,V=Γ,C=ϕψ,where(32)Γ=βδ0δ1ϕ2−bβδ0ϕψ−αβδ1ϕψ−αδ1γϕψ+αbβψ2+βδ0ϕψ−αβψ2−αγψ2b−1β2ψ2+αβγψ2−β2δ1ϕψ.Table 1
Dimensional parameter values.
K
α
β
δ
γ
b
q
λ
δ
0
δ
1
10
6
0.206
0.001
0.5115
0.001
10
5
4.16
0.005
0.006
[41]
[41]
[41]
[41]
[41]
[42]
[23]
[23]
Estimated
EstimatedIt is a difficult undertaking to investigate the stability of this state without substituting parameter values because of the many terms involved. Using the parameter values in Table1, the eigenvalues are −0.031, −0.025, −1.01, and −8.13, implying that this state is stable.With the consideration of an exponential infusion functionξt=ϕexp−at, the equation (5) are first turned into an autonomous system of differential equations by letting W=ϕexp−at. This system has three steady states:(i)
A tumor-free state where all cell concentrations diminish to zero:(33)U=0,I=0,V=0,C=0,W=0.This state is unstable because the eigenvalues−γ, −ψ, −a, −1, and α are not all negative.(ii)
A state where the tumor grows to its maximum size:(34)U=1,I=0,V=0,C=0,W=0.(iii)
The characteristic polynomial evaluated at this state is equation (25) and this state is locally asymptotically stable if β+γ>βγ, otherwise it is unstable.(35)U=γb−1β,I=αb−αβγ−αγ2b−12β2+αb−αβγ,V=αb−αβ−αγb−1β2+αβγ,C=0,W=0,b>1.The eigenvalues evaluated at this steady state are also big expressions and difficult to analyse analytically and extract conditions for stability. With the set of parameter values in Table1, the eigenvalues are −0.1, −8.13, 1.054, and −0.014±0.085i, implying that this state is stable.Similarly, withξt=ϕsin2at, one changes equation (5) into an autonomous system of equations by letting W=ϕsin2at. The autonomous system has six steady states:(i)
Tumor-free state where all cell concentrations are wiped out of body tissue:(36)U=0,I=0,V=0,C=0,W=0.The eigenvalues evaluated at this state are−γ, −ψ, −a, −1, and α implying that it is unstable.(ii)
A state where the tumor grows to its maximum size:(37)U=1,I=0,V=0,C=1,W=0.The condition for stability of this state is same as with the exponential drug infusion case, that is, the state is locally asymptotically stable ifβ+γ>βγ.(iii)
Tumor-free state where some concentration of the drug remains:(38)U=0,I=0,V=0,C=ϕaψ,W=ϕa.This state is locally asymptotically stable ifδ0ϕ>αψ and 1/2<ϕ<1 because the eigenvalues evaluated at this state are(39)δ0ϕ−αψψ,−δ1ϕ+ψψ,−22aϕϕ−ϕ4ϕ1−ϕ,−γ,−ψ,otherwise, it is unstable.(iv)
Infected tumor-free state where all infected tumor cells are wiped but a certain proportion of the uninfected remains:(40)U=aαψ−δ0ϕaαψ,I=0,V=0,C=ϕaψ,W=ϕa,aαψ≥δ0ϕ.(v)
Drug-free state where the chemotherapeutic drug is wiped out of body tissue and proportions of all the other cell concentrations remain:(41)U=γb−1β,I=αb−αβγ−αγ2b−12β2+αb−αβγ,V=αb−αβ−αγb−1β2+αβγ,C=0,W=0,b>1.(vi)
Tumor dormant state:(42)U=δ1ϕ+aψγab−aβψ−βδ1ϕ,I=Γ2γabψ−δ1ϕ−aψ,V=γ2aψ,C=ϕaψ,W=ϕa,where(43)Γ2=a2αbβψ2−abβδ0ϕψ−aαβδ1ϕψ−aαδ1γϕψ−a2αβψ2−a2αγψ2+βδ0δ1ϕ2+aβδ0ϕψabβψ+aαγψ−βδ1ϕ−aβψaβψ.The conditions for stability for the last three states all depend on huge expressions from which it is hard to extract meaningful biological implications. This analysis, however, suggests that with both treatments and using a sinusoidal type infusion, the tumor can be eliminated from body tissue provided that the combination of the drug infusion rate and the lysis rate of the tumor is greater than the combination of the tumor growth rate and rate of drug loss.


## 3.1. Model Solutions

To investigate the efficacy of each treatment and their combination, we first study the dynamics of the system without treatment. Without any form of treatment, model (5) is reduced to only one equation:(6)U˙t=αUt1−Ut,U0=U0,whose solution is(7)Ut=U01−U0exp−αt+U0,limt⟶∞Ut=1,implying that the tumor logistically grows to its maximum fractional size. Next, the model (5) is analysed, with chemotherapy, with virotherapy, and then with both treatments incorporated. We obtain, where possible analytical and time invariant solutions which predict the long term dynamics of the model equation (1).Without virotherapy (Vt=0), the system (5) is transformed to(8)U˙t=αUt1−Ut−δ0UtCt,C˙t=ξt−ψCt,with U0=U0 and C0=C0. The second equation in equation (8) is a first-order linear ordinary differential equation which can easily be solved to give(9)Ct=exp−ψt∫ξtexpψtdt+R,where R is a constant of integration. The solution to the first equation in equation (8) depends on the infusion function ξt. For a fixed infusion function ϕ,(10)Ut=eαt−δ0ϕt/ψ+δ0e−ψt/ψα∫eαt−δ0ϕt/ψ+δ0e−ψt/ψdt+eδ0ψ/U0,Ct=C0−ϕψe−ψt+ϕψ.From the solution ofCt in equation (10),(11)limt⟶∞Ct=ϕψ.Biologically, it can be inferred that with a constant drug infusion and without virotherapy, the tumor is not completely cleared and a certain proportion of the drug remains in body tissue. The tumor clearance depends on the drug induced lysis of the tumor and the drug infusion rate which should be maximized and the tumor growth and drug decay rate which should be minimized.Forξt=ϕexp−at,(12)Ut=eaαt/a−ψ−αψt/a−ψ+cδ0e−ψt/a−ψ−acδ0e−ψt/a−ψψ+δ0ϕe−at/a−ψaa∫eaαt/a−ψ−αψt/a−ψ+cδ0e−ψt/a−ψ−acδ0e−ψt/a−ψψ+δ0ϕe−at/a−ψadt+ecδ0/a−ψ/U0,Ct=C0ψ−ϕψ−ϕe−at+ψta−ψe−ψt,where c=C0ψ−ϕ/ψ.From equation (12),(13)limt⟶∞Ct=0,limt⟶∞Ut=U∗,where U∗ is a fractional tumor cell concentration between 0 and 1. This suggests that with a single dosage infusion of the chemotherapeutic drug with exponential decay and without virotherapy, the tumor cannot be cleared from body tissue. The drug is also completely depleted from the body.Whenξt=ϕsin2at is substituted in equation (8), the resulting differential equations are solved to give(14)Ct=122R−ψ2cos2ateψt+2aψeψtsin2at−4a2+ψ2eψtϕ4a2ψ+ψ3e−ψt,where(15)R=C0+ψ2−ϕ4a2+ψ22ψ4a2+ψ2,limt⟶∞Ct=C∗=fa,ϕ,ψ.This suggests that with time, some drug concentration remains in the body tissue.Theorem 1.
The system (8), with constant infusion, has no periodic solutions for positive Ut and Ct.Proof.
Using Dulac’s criterion ([40]):
SupposeX˙=fx and fx is continuously differentiable on a simply connected domain D⊂R. If there exists a real valued function gx such that ∇⋅gX˙=∇⋅gf has one sign in D, then there are no closed orbits in D. Using Dulac’s criterion, it is sufficient to show that(16)∂∂UgU˙+∂∂CgC˙≠0,∀U,C∈R+2.
Consider(17)gU,C=1UC,∇⋅gX˙=∂∂UgU˙+∂∂CgC˙,=−αC+ξtUC2<0,∀U,C∈R+2.Theorem 2.
The system (8) has at least two steady states for each of the drug infusion functions: (1)
For the constant drug infusion functionξt=ϕ, there are two steady states of equation (8): U=0,C=ϕ/ψ which is locally asymptotically stable provided that δ0ϕ>αψ and U=1−δ0ϕ/αψ,C=0 which is locally asymptotically stable provided that δ0ϕ+αψ>2δ0ϕψ2; otherwise, it is unstable.(2)
For the exponential drug infusionξt=ϕexp−at, equation (8) has two steady states: U=0,C=0,W=0 which is unstable and U=1,C=0,W=0 which is locally asymptotically stable.(3)
For the sinusoidal infusion function, there are four steady states of equation (8): U=0,C=0,W=0, U=1,C=1,W=0, and U=aαψ−δ0ϕ/aαψ,C=ϕ/aψ,W=ϕ/a which are unstable and U=0,C=ϕ/aψ,W=ϕ/a which is locally asymptotically stable if δ0ϕ>aαψ and ϕ<1.Proof.
(1)
It is easy to show that when equation (8) is equated to zero, one obtains two steady states. The characteristic polynomial of the Jacobian matrix for equation (8) evaluated at 0,ϕ/ψ is(18)λ2+−α+δ0ϕψ+ψλ+δ0ϕ−αψ,whose roots λ can only be negative if δ0ϕ>αψ.
The characteristic polynomial evaluated at1−δ0ϕ/αψ,0 is(19)λ2+−2δ0ϕψ+α+δ0ϕψ+ψλ−2δ0ϕψ2+δ0ϕ+αψ,whose roots are negative provided that δ0ϕ+αψ>2δ0ϕψ2.
(2)
By lettingW=ϕexp−at, equation (8) is turned into an autonomous system:(20)U˙t=αUt1−Ut−δ0UtCt,C˙t=W−ψCt,W˙t=−aW.
The eigenvalues of the Jacobian matrix for equation (20) evaluated at 0,0,0 which are −ψ,α, and 0 and at 1,0,0 which are −α,−ψ, and −a are all negative.
(3)
Similarly, by lettingW=ϕsin2t, equation (8) becomes the autonomous system:(21)U˙t=αUt1−Ut−δ0UtCt,C˙t=W−ψCt,W˙t=4aWϕ−W1/2.
The eigenvalues of the Jacobian matrix for equation (20) evaluated at 0,0,0 are −ψ,α, and 0 and the eigenvalues evaluated at 1,0,0 are ψ,−α, and 0. For the third steady state to exist, aαψ≥δ0ϕ, and the eigenvalues evaluated at this state are −aαψ−δ0ϕ/aψ,−4aϕ/4ϕ1−ϕ,−ψ, implying that for it to be locally asymptotically stable, δ0ϕ>aαψ; yet for this to happen, the steady state will not exist. The eigenvalues evaluated at U=0,C=ϕ/aψ,W=ϕ/a are −aαψ−δ0ϕ/aψ,−4aϕ/4ϕ1−ϕ,−ψ implying that this steady state is locally asymptotically stable if δ0ϕ>aαψ and ϕ<1.
Theorems 1 and 2 show that there are no periodic solutions in the dynamics of equation (8) and with a constant drug infusion, the tumor can be eliminated from body tissue by chemotherapy provided that the combination of the chemotherapeutic drug-induced lysis of the tumor and the drug infusion is greater than the combination of the intrinsic tumor growth rate and the drug deactivation rate. The tumor can also be wiped out with a periodic drug infusion provided that the combination of the tumor-induced lysis by the drug and the dosage is greater than the intrinsic tumor growth rate and drug decay rate. With the exponential infusion method, the tumor is not removed from body tissue and may grow to its maximum size.
Without chemotherapy, equation (5) is reduced to(22)U˙t=αUt1−Ut−It−βUtVt,I˙t=βUtVt−It,V˙t=bIt−βUtVt−γVt.
The analytical solutions to system (22) are not easy to obtain. The derivatives of equation (22) are therefore equated to zero to obtain time invariant solutions and investigate their stability by linearizing equation (22) about the steady states.Theorem 3.
(1)
Ifβ+γ>bβ, the system (12) has two steady states: a tumor free cell state 0,0,0 which is unstable and an infection tumor free state 1,0,0 which is locally asymptotically stable.
(2)
Ifbβ>β+γ, the system (22) has three steady states: the tumor free state 0,0,0 and the infected free state 1,0,0 which are unstable and a tumor dormant state:(23)U=γb−1β,I=αγβb−1−γββb−12+γb−1,V=αβb−1−γb−12β2+αγ,b>1,which is locally asymptotically stable if a0,a1,a2>0 and a1a2>a0 where ai are coefficients of the characteristic equation.Proof.
(1)
The characteristic equation evaluated at0,0,0 is(24)λ3+γ−α+1λ2+γ−αγ−αλ−αγ=0,λ−αλ+1λ+γ=0,from which λ1=α, λ2=−1 and λ3=−γ, thus rendering it unstable.
The characteristic equation evaluated at1,0,0 is(25)λ+αλ2+λ1+β+γ+β+γ−bβ=0,from which λ1=−α and λ2+λ1+β+γ+β+γ−bβ=0 which are all negative since β+γ>bβ.
(2)
The characteristic polynomial evaluated at the tumor dormant state isλ3+a2λ2+a1λ+a0=0, where(26)a2=2Aα+Aβ+Cβ−α+γ+1,a1=2A2αβ−Aαβ+Cαβ−Abβ+2Aαγ+Cβγ+2Aα+Aβ+Cβ−αγ−α+γ,a0=−2A2αbβ+2A2αβ+Aαbβ+Cαβγ−Aαβ+2Aαγ+Cβγ−αγ,and A,B,C are the coordinates of the tumor dormant state. Using Routh–Hurwitz stability criterion, this state will only be locally asymptotically stable if a0,a1,a2>0 and a1a2>a0.Since the infected tumor-free state is undesirable, the reverse of the conditionβ+γ>bβ is necessary for tumor eradication from body tissue. In other words, bβ>β+γ, that is, the product of the virus replication rate and their burst size should be greater than the sum of the burst size and virus replication rate. We also notice from equation (23) that(27)limβ⟶∞U=limb⟶∞U=0.It is therefore evident that high virus replication rateβ and burst size b lead to lower tumor cell concentrations. The steady-state solutions of equation (23) involve many parameters, thereby giving rise to large expressions in the conditions for its stability. It is therefore a difficult undertaking to infer biological implications from these conditions. Nevertheless, it can be observed that virotherapy may only succeed in eliminating cancer from body tissue when the virus deactivation rate is very small or even zero and the virus replication rate very high.Next, the model with both treatments is analysed. For a constant drug infusion rateϕ, the system (5) has three steady states;(i)
Tumor-free steady state:(28)U=0,I=0,V=0,C=ϕψ.Here, the tumor and viruses are cleared from body tissue by the coupled treatment and a fraction of the chemotherapeutic drug remains in body tissue. The eigenvalues of the Jacobian matrix evaluated at this state are(29)λ1=−δ0ϕ−αψψ,λ2=−δ1ϕ+ψψ,λ3=−γ,λ4=−ψ,implying that this desirable state is locally asymptotically stable if δ0ϕ>αψ. From this condition, in order to clear a tumor, the combination of the rate at which the drug kills the uninfected tumor cells and the drug infusion must be higher than the tumor growth rate and deactivation of the drug from body tissue.(ii)
Infected tumor-free state:(30)U=1−δ0ϕαψ,I=0,V=0,C=ϕψ.
In this state, the whole tumor is not cleared as a fraction of uninfected tumor cells remain and all the infected ones are cleared by the treatment combination. Using the parameter values in Table1, the eigenvalues of the Jacobian matrix evaluated at the infected tumor-free state are 0.403, 8.13, and 2.598±2.418i, implying that the infected tumor-free state is unstable.
(iii)
Tumor dormant state:(31)U=δ1ϕ+ψγbψb−1−δ1ϕβ,I=Γγδ1ϕ−b−1ψβ,V=Γ,C=ϕψ,where(32)Γ=βδ0δ1ϕ2−bβδ0ϕψ−αβδ1ϕψ−αδ1γϕψ+αbβψ2+βδ0ϕψ−αβψ2−αγψ2b−1β2ψ2+αβγψ2−β2δ1ϕψ.Table 1
Dimensional parameter values.
K
α
β
δ
γ
b
q
λ
δ
0
δ
1
10
6
0.206
0.001
0.5115
0.001
10
5
4.16
0.005
0.006
[41]
[41]
[41]
[41]
[41]
[42]
[23]
[23]
Estimated
EstimatedIt is a difficult undertaking to investigate the stability of this state without substituting parameter values because of the many terms involved. Using the parameter values in Table1, the eigenvalues are −0.031, −0.025, −1.01, and −8.13, implying that this state is stable.With the consideration of an exponential infusion functionξt=ϕexp−at, the equation (5) are first turned into an autonomous system of differential equations by letting W=ϕexp−at. This system has three steady states:(i)
A tumor-free state where all cell concentrations diminish to zero:(33)U=0,I=0,V=0,C=0,W=0.This state is unstable because the eigenvalues−γ, −ψ, −a, −1, and α are not all negative.(ii)
A state where the tumor grows to its maximum size:(34)U=1,I=0,V=0,C=0,W=0.(iii)
The characteristic polynomial evaluated at this state is equation (25) and this state is locally asymptotically stable if β+γ>βγ, otherwise it is unstable.(35)U=γb−1β,I=αb−αβγ−αγ2b−12β2+αb−αβγ,V=αb−αβ−αγb−1β2+αβγ,C=0,W=0,b>1.The eigenvalues evaluated at this steady state are also big expressions and difficult to analyse analytically and extract conditions for stability. With the set of parameter values in Table1, the eigenvalues are −0.1, −8.13, 1.054, and −0.014±0.085i, implying that this state is stable.Similarly, withξt=ϕsin2at, one changes equation (5) into an autonomous system of equations by letting W=ϕsin2at. The autonomous system has six steady states:(i)
Tumor-free state where all cell concentrations are wiped out of body tissue:(36)U=0,I=0,V=0,C=0,W=0.The eigenvalues evaluated at this state are−γ, −ψ, −a, −1, and α implying that it is unstable.(ii)
A state where the tumor grows to its maximum size:(37)U=1,I=0,V=0,C=1,W=0.The condition for stability of this state is same as with the exponential drug infusion case, that is, the state is locally asymptotically stable ifβ+γ>βγ.(iii)
Tumor-free state where some concentration of the drug remains:(38)U=0,I=0,V=0,C=ϕaψ,W=ϕa.This state is locally asymptotically stable ifδ0ϕ>αψ and 1/2<ϕ<1 because the eigenvalues evaluated at this state are(39)δ0ϕ−αψψ,−δ1ϕ+ψψ,−22aϕϕ−ϕ4ϕ1−ϕ,−γ,−ψ,otherwise, it is unstable.(iv)
Infected tumor-free state where all infected tumor cells are wiped but a certain proportion of the uninfected remains:(40)U=aαψ−δ0ϕaαψ,I=0,V=0,C=ϕaψ,W=ϕa,aαψ≥δ0ϕ.(v)
Drug-free state where the chemotherapeutic drug is wiped out of body tissue and proportions of all the other cell concentrations remain:(41)U=γb−1β,I=αb−αβγ−αγ2b−12β2+αb−αβγ,V=αb−αβ−αγb−1β2+αβγ,C=0,W=0,b>1.(vi)
Tumor dormant state:(42)U=δ1ϕ+aψγab−aβψ−βδ1ϕ,I=Γ2γabψ−δ1ϕ−aψ,V=γ2aψ,C=ϕaψ,W=ϕa,where(43)Γ2=a2αbβψ2−abβδ0ϕψ−aαβδ1ϕψ−aαδ1γϕψ−a2αβψ2−a2αγψ2+βδ0δ1ϕ2+aβδ0ϕψabβψ+aαγψ−βδ1ϕ−aβψaβψ.The conditions for stability for the last three states all depend on huge expressions from which it is hard to extract meaningful biological implications. This analysis, however, suggests that with both treatments and using a sinusoidal type infusion, the tumor can be eliminated from body tissue provided that the combination of the drug infusion rate and the lysis rate of the tumor is greater than the combination of the tumor growth rate and rate of drug loss.


## 4. Mathematical Analysis of Delay Model

The nondimensionalised delay model is(44)U˙t=αUt1−Ut+tK−βUt−τ1Vt−τ1−δ0Ut−τ2Ct−τ2,I˙t=βUt−τ1Vt−τ1−It−δ1It−τ2Ct−τ2,V˙t=bIt−βUt−τ1Vt−τ1−γVt,C˙t=ϕt−ψCt.Without virotherapy the model (44) in nondimensional form becomes(45)U˙t=αUt1−Ut−δ0Ut−τ2Ct−τ2,C˙t=ϕt−ψCt.The system (45) has two steady states: a tumor-free state T0=0,ϕ/ψ and a tumor dormant state 1−δ0ϕ/αψ,0, just as previously seen. Letting Z1t=Ut−U∗ and Z4t=Ct−C∗ where U∗,C∗ are steady states of equation (45), the linearized model about U∗,C∗ of equation (45) is(46)Z1˙t=αZ1t−2αU∗Z1t−δ0U∗Z4t−τ2−δ0C∗Z1t−τ2,Z4˙t=ϕt−ψZ4t.The characteristic equation of equation (46) evaluated at the tumor-free steady state is(47)fλ=λ+ψλ−α+e−λτ2δ0ϕψ=0.Forτ2=0, we obtain the same characteristic polynomial as in the ODE case (Theorem 2). For τ2≠0, equation (47) is a transcendental equation and therefore has infinitely many roots and also makes it nontrivial to determine these roots. Nonetheless, the following is noticed:Lemma 1.
(1)
The tumor-free stateT0 is stable if τ2=ψ/δ0ϕ and α=δ0ϕ/ψ, i.e., Equation (47) has a negative root −ψ and a zero double root.
(2)
T
0 is stable for a sufficiently small τ2, i.e., Equation (47) has negative real roots for 0≤τ2≤τ20.
(3)
Equation (47) has a pair of purely imaginary roots ±iω if αψ≤δ0ϕ and the other root is −ψ. Therefore, T0 is unstable.Proof.
(1)
Ifτ2=ψ/δ0ϕ and α=δ0ϕ/ψ, then(48)dfdλ=2λ−α+δ0ϕψe−λτ2−λτ2δ0ϕe−λτ2ψ+ψ−τ2δ0ϕe−λτ2,d2fdλ2=2−2τ2δ0ϕψe−λτ2+τ22δ0ϕλψ+1.
Form equation (48),(49)dfdλλ=0=0,d2fdλ2λ=0=ψδ0ϕ>0.
Thus, equation (47) has a double zero root.
(2)
If we denoteρτ2+iωτ2 as the root of the equation (47), the tumor-free state is stable if ρ0<0. By continuity, if τ2>0 is sufficiently small, we still have ρτ2<0 and the tumor-free state is stable.
(3)
If equation (47) has only purely imaginary roots, then the roots should be solutions to the equation(50)λ−α+e−λτ2δ0ϕψ=0.
Assume thatλ=iω, ω>0 is the root of equation (50). Substituting λ=iω and separating real and imaginary parts, one gets(51)−α+δ0ϕψcosωτ2=0,ω−δ0ϕψsinωτ2=0.
Squaring both equations and adding gives(52)ω2+α2−δ02ϕ2ψ2=0,from which ω=±δ0ϕ/ψ+αδ0ϕ/ψ−α if αψ>δ0ϕ and the critical values of τ2 are(53)τ2k=1ω+arccosαψδ0ϕ,k=0,1,2,….
Without chemotherapy equation (44) becomes(54)U˙t=αUt1−Ut−βUt−τ1Vt−τ1,I˙t=βUt−τ1Vt−τ1−It,V˙t=bIt−βUt−τ1Vt−τ1−γvt.
By lettingZ1t=Ut−U∗, Z2t=It−I∗ and Z3t=Vt−V∗ where U∗,I∗, and V∗ are steady states of equation (54), the linearized model about U∗,I∗,C∗ of equation (54) is(55)Z1˙t=αZ1t−2αU∗Z1t+αI∗Z1t−αU∗Z2t+βV∗Z1t−τ1+βU∗Z3t−τ1,Z2˙t=−βV∗Z1t−τ1+βU∗Z3t−τ1−Z2t,Z3˙t=bZ2t+βV∗Z1t−τ1+βU∗Z3t−τ1−γZ3t.
The characteristic equation of equation (55) at 1,0,0 is(56)λ+αλ+1λ+γ+βe−λτ1−bβe−λτ1=0.
Forτ1=0, we retrieve the same results as in Theorem 3. For τ≠0, we have a transcendental equation to solve.


## 5. Numerical Simulations

### 5.1. Parameter Values

The parameter values used were obtained from fitted experimental data for untreated tumors and virotherapy in mice [41]. The tumor carrying capacity, K, is taken to be 106 cells per unit volume. The number of viruses produced per day b is to be in the range 10−1000 [42]. Drug infusion and decay rates q and λ agree with cancer pharmokinetic studies [33, 34].


### 5.2. Simulation Results

Numerical simulations of the models (5) and (45) are presented, first with monotherapies followed with combination treatment. In all simulations, unless stated otherwise, initial concentrations are considered to be U0=1, I0=0, V0=0.1, and C0=0.1 with a high fractional untreated tumor cell count to necessitate clinical intervention. The equations were integrated using a Runge–Kutta fourth-order scheme and implemented in MATLAB. It is worth noting that the scale for the time and concentrations is, respectively, 1 unit ≈2 days and 1 unit =106 number of cells.Figure1 shows numerical solutions of the chemo-only model (8). The figure shows that despite the drug infusion method, the tumor is not cleared from body tissue. These numerical solutions agree with the analytical results obtained in the previous sections that chemotherapy on its own may not clear all tumor cells in body tissue, and the tumor grows to its maximum size and the drug concentration decays to zero. Section 3 revealed that total tumor clearance from body tissue can possibly be achieved if δ0ϕ>αψ. Nonetheless, the parameter values used do not conform to this condition.Figure 1
Solutions of the model without virotherapy: equation (5) showing a variation of fractional concentrations with time, using (a) a constant, (b) an exponential, and (c) a sinusoidal drug infusion. The initial cell concentrations are U0=0.8 and C0=0.2.
(a)
(b)
(c)Figure2 shows the dynamics of the viro-only model. It is clear from this figure that virotherapy alone could possibly clear all tumor cells from body tissue provided that the virus infection rate is high and with a large virus burst size. Figures 2(a) and 2(b) show a variation of the fractional tumor and virus concentrations against time for different values of the virus replication rate. It is noticed that with a small virus infection rate for example β=10−6, it takes a longer time to clear the tumor cells. Figures 2(b) and 2(c) show a variation of the fractional concentrations with two different burst sizes. From these figures, we notice that when b=10, it takes about 10 days to reduce the whole tumor concentrations to zero while it takes only about 5 days with b=100, implying that a high virus burst size yields a quick recovery with virotherapy treatment. These numerical intimations concur with the analytical results established in Section 3.1.Figure 2
Solutions of model (12) without chemotherapy showing a variation of fractional concentrations with nondimensional time using low and high values of the virus replication rate β and burst size b, that is, (a) β=10−6, b=10; (b) β=10−3, b=10; (c) b=10, β=10−6; and (d) b=100, β=10−6.
(a)
(b)
(c)
(d)Figure3 displays the dynamics of model (5) with both treatments. The numerical results are similar to those of Figure 2, only that with both treatments, it takes a shorter time to bring the tumor cell concentrations to zero. High values of the virus infection rate and burst size lead to tumor clearance in a shorter time period. Both Figures 2 and 3 show that an increase in the virus multiplication rate and burst size increase the infected tumor cells concentration. For example, comparing Figures 3(a) and 3(b), the number of infected tumor cells was about 0.15×106 when b=15 and this increased to 0.35×106 when b=25.Figure 3
Solutions of model (3) with both treatments showing a variation of fractional concentrations with time using high and low virus burst sizes, that is, b=2 and b=5 and with different drug infusion functions, that is, (a) constant, (b) exponential, and (c) sinusoidal.
(a)
(b)
(c)Model (44) was simulated using dde23 in MATLAB. Figure 4 displays the numerical simulations for the delay model (44) with low and high values of τ1 and τ2, the virus, and chemotherapeutic delays. Since secondary transcription and viral protein synthesis can be delayed by about six to eight hours [37], τ1 was considered to be between 0.001 and 0.01. The chemotherapeutic response delay is higher than the virus response delay, thus τ2 was taken to be between 0.1 and 0.3. The figure shows that when both delays are increased, the time it takes for cell concentration solutions to converge is slightly increased although they converge at the same steady states as without the delay. For τ1=0.001 in Figure 4(a), it took about 4 days for the whole tumor to clear whereas it took about 8 days with τ1=0.01 in Figure 4(b). Comparing Figures 4(b) and 4(c), when τ2 was increased from 0.001 to 0.3, the time it took for the whole tumor to be cleared was increased from six to eight days. Initially there are oscillations for high values of the chemotherapeutic delay although the cell concentrations converge at the same steady states, just as the case with no delays. Figure 5 is a close up form of Figure 4 to display oscillations caused by the virus and drug delays. The oscillations only occur in the initial stages of treatment but later fade away. Nonetheless, the results in Figures 4 and 5 suggest that it is imperative to design viruses and drugs which are highly responsive in order to minimize these delays.Figure 4
Solutions of delay model (15) showing a variation of fractional concentrations with time using (a) τ1=0.001, τ2=0.1; (b) τ1=0.01, τ2=0.1; (c) τ1=0.001, τ2=0.001; and (d) τ1=0.001, τ2=0.3.
(a)
(b)
(c)
(d)Figure 5
Close up images of Figures4(c) and 4(d) showing oscillations in cell concentrations.
(a)
(b)


## 5.1. Parameter Values

The parameter values used were obtained from fitted experimental data for untreated tumors and virotherapy in mice [41]. The tumor carrying capacity, K, is taken to be 106 cells per unit volume. The number of viruses produced per day b is to be in the range 10−1000 [42]. Drug infusion and decay rates q and λ agree with cancer pharmokinetic studies [33, 34].


## 5.2. Simulation Results

Numerical simulations of the models (5) and (45) are presented, first with monotherapies followed with combination treatment. In all simulations, unless stated otherwise, initial concentrations are considered to be U0=1, I0=0, V0=0.1, and C0=0.1 with a high fractional untreated tumor cell count to necessitate clinical intervention. The equations were integrated using a Runge–Kutta fourth-order scheme and implemented in MATLAB. It is worth noting that the scale for the time and concentrations is, respectively, 1 unit ≈2 days and 1 unit =106 number of cells.Figure1 shows numerical solutions of the chemo-only model (8). The figure shows that despite the drug infusion method, the tumor is not cleared from body tissue. These numerical solutions agree with the analytical results obtained in the previous sections that chemotherapy on its own may not clear all tumor cells in body tissue, and the tumor grows to its maximum size and the drug concentration decays to zero. Section 3 revealed that total tumor clearance from body tissue can possibly be achieved if δ0ϕ>αψ. Nonetheless, the parameter values used do not conform to this condition.Figure 1
Solutions of the model without virotherapy: equation (5) showing a variation of fractional concentrations with time, using (a) a constant, (b) an exponential, and (c) a sinusoidal drug infusion. The initial cell concentrations are U0=0.8 and C0=0.2.
(a)
(b)
(c)Figure2 shows the dynamics of the viro-only model. It is clear from this figure that virotherapy alone could possibly clear all tumor cells from body tissue provided that the virus infection rate is high and with a large virus burst size. Figures 2(a) and 2(b) show a variation of the fractional tumor and virus concentrations against time for different values of the virus replication rate. It is noticed that with a small virus infection rate for example β=10−6, it takes a longer time to clear the tumor cells. Figures 2(b) and 2(c) show a variation of the fractional concentrations with two different burst sizes. From these figures, we notice that when b=10, it takes about 10 days to reduce the whole tumor concentrations to zero while it takes only about 5 days with b=100, implying that a high virus burst size yields a quick recovery with virotherapy treatment. These numerical intimations concur with the analytical results established in Section 3.1.Figure 2
Solutions of model (12) without chemotherapy showing a variation of fractional concentrations with nondimensional time using low and high values of the virus replication rate β and burst size b, that is, (a) β=10−6, b=10; (b) β=10−3, b=10; (c) b=10, β=10−6; and (d) b=100, β=10−6.
(a)
(b)
(c)
(d)Figure3 displays the dynamics of model (5) with both treatments. The numerical results are similar to those of Figure 2, only that with both treatments, it takes a shorter time to bring the tumor cell concentrations to zero. High values of the virus infection rate and burst size lead to tumor clearance in a shorter time period. Both Figures 2 and 3 show that an increase in the virus multiplication rate and burst size increase the infected tumor cells concentration. For example, comparing Figures 3(a) and 3(b), the number of infected tumor cells was about 0.15×106 when b=15 and this increased to 0.35×106 when b=25.Figure 3
Solutions of model (3) with both treatments showing a variation of fractional concentrations with time using high and low virus burst sizes, that is, b=2 and b=5 and with different drug infusion functions, that is, (a) constant, (b) exponential, and (c) sinusoidal.
(a)
(b)
(c)Model (44) was simulated using dde23 in MATLAB. Figure 4 displays the numerical simulations for the delay model (44) with low and high values of τ1 and τ2, the virus, and chemotherapeutic delays. Since secondary transcription and viral protein synthesis can be delayed by about six to eight hours [37], τ1 was considered to be between 0.001 and 0.01. The chemotherapeutic response delay is higher than the virus response delay, thus τ2 was taken to be between 0.1 and 0.3. The figure shows that when both delays are increased, the time it takes for cell concentration solutions to converge is slightly increased although they converge at the same steady states as without the delay. For τ1=0.001 in Figure 4(a), it took about 4 days for the whole tumor to clear whereas it took about 8 days with τ1=0.01 in Figure 4(b). Comparing Figures 4(b) and 4(c), when τ2 was increased from 0.001 to 0.3, the time it took for the whole tumor to be cleared was increased from six to eight days. Initially there are oscillations for high values of the chemotherapeutic delay although the cell concentrations converge at the same steady states, just as the case with no delays. Figure 5 is a close up form of Figure 4 to display oscillations caused by the virus and drug delays. The oscillations only occur in the initial stages of treatment but later fade away. Nonetheless, the results in Figures 4 and 5 suggest that it is imperative to design viruses and drugs which are highly responsive in order to minimize these delays.Figure 4
Solutions of delay model (15) showing a variation of fractional concentrations with time using (a) τ1=0.001, τ2=0.1; (b) τ1=0.01, τ2=0.1; (c) τ1=0.001, τ2=0.001; and (d) τ1=0.001, τ2=0.3.
(a)
(b)
(c)
(d)Figure 5
Close up images of Figures4(c) and 4(d) showing oscillations in cell concentrations.
(a)
(b)


## 6. Discussion

The results in this study contend with previous experimental and mathematical studies that oncolytic viruses enhance chemotherapeutic drugs in the treatment of cancer [12,13,15–20]. A study by Ungerechts et al. [16] examined the synergy between a reprogrammed oncolytic virus and two chemotherapeutic drugs in the mantle cell lymphoma (MCL). They investigated the efficacy of different procedures of a measles virus in combination with fludarabine and cyclophosphamide (CPA). Their study suggested that that CPA administration before virotherapy enhanced oncolytic efficacy. An experimental study by Ulasov et al. [17] indicated that the combination of virotherapy and temozolomide is capable of eliminating malignant glioma. Their results showed that 90% of treated mice survived beyond the 100 days’ mark after being treated. Another study by Alonso et al. [18] showed that the amalgamation of oncolytic adenovirus (ICOVIR-5) with either everolimus (RAD001) or temozolomide (TMZ) resulted in an enhanced antiglioma effect. Recent mathematical studies by Malinzi et al. [19, 20] assert that combining chemotherapeutic drugs with oncolytic viruses is more efficient than using either treatments alone. In [20], it is indicated that although chemotherapy alone may clear tumor cells from body tissue if drug efficacy is bigger than the tumor growth rate, the use of both OV and drugs leads to enhanced treatment effects.Biologically, the reduction of a tumor to undetectable levels in less than a week is unrealistic in comparison to existing clinical and research studies [43]. The duration of cancer treatment depends on several factors including the type of cancer being treated and the patient cells’ characteristics. This makes it hard to predict the time period to clear a tumor in body tissue. Moreover, a tumor can be reduced to insignificant levels but may later regrow [44]. Nevertheless, this study agrees with the fact that chemovirotherapy is highly likely to bring the tumor to undetectable levels in a short time period just as previously established in [19, 20].The mathematical model considered in this study is built on a couple of simplifying assumptions and thus omits pertinent biological aspects. For example, the drug infusion functions considered are not pragmatic. It is thus imperative to extend this study to consider more realistic drug infusion functions that describe all the important pharmacodynamics properties. Nonetheless, this study indicates that a cancer patient should not be given a single bolus injection as it is less effective compared to periodic or constant drug infusions.


## 7. Conclusion

The aim of this study was to investigate the outcome of the amalgamation of chemotherapy and virotherapy in treating cancer using three different drug infusion methods and to compare the efficacy of using chemotherapy and virotherapy individually.A mathematical model in the form of nonlinear and nonautonomous first-order ordinary differential equations was developed. It was extended into delay differential equations to account for delays as a result of the infection of tumor cells by the virus and chemotherapeutic drug responses. The model’s well-posedness was shown by proving existence, positivity, and boundedness of the model solutions. Analysis of the model was done with each of the treatments and for each of the infusion functions. Exact solutions were determined where possible. Stability of the time invariant solutions was carried out to determine the conditions under which a tumor-free situation may be achieved. Numerical simulations for the ODE and DDE models were, respectively, carried out using ode23s and dde23 in MATLAB. The model analysis suggested the following:(i)
A tumor can grow to its maximum size in case where there is no treatment.(ii)
Chemotherapy alone is capable of clearing tumor cells in body tissue provided that the drug-induced lysis of the tumor and the drug infusion rate are maximized and the drug decay and tumor growth are minimized.(iii)
Constant and periodic drug infusions are more potent than a single bolus injection.(iv)
Successful virotherapy is highly dependent on virus burst size and infection rate.(v)
With the use of both chemotherapy and virotherapy, a tumor may be cleared from body tissue in less than a month.(vi)
Successful chemovirotherapy depends on the virus burst size and replication rate, chemotherapeutic drug lysis, infusion and decay rates, and the method of drug infusion.(vii)
Both the virus and chemotherapeutic response delays increase the period within which a tumor can be cleared from body tissue and thus treatment options should strive to minimize them by designing viruses and drugs which are highly responsive.



---

*Source: 7576591-2019-03-11.xml*