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

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.


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. ese 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]. e 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][7][8][9][10][11]. Virotherapy involves the use of oncolytic viruses that infect, multiply, and directly lyse cancer cells with less or no toxicity [9]. eir 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][16][17][18][19][20]. Nguyen et al. [12] gave an account of the mechanisms through which drugs can successfully be used in a combination with oncolytic viruses. ey 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]. e 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. e model constructed combines elements from existing mathematical models [10,11,19,20,[23][24][25][26][27][28][29][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. eir 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. e 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. e article is organised as follows: Section 2 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. e 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.

Model without Delay.
Time-dependent cell concentrations of uninfected tumor cells U(t), infected tumor cells I(t), a free virus population V(t), and a chemotherapeutic drug C(t) in an avascular tumor localization are considered. e uninfected tumor grows logistically at an intrinsic rate α per day, and the total tumor carrying capacity is K cells in a tumor nodule. e infected tumor cells die off at a rate δ per day. Virus multiplication in the tumor is represented by the function βU(t)V(t), where β is the virus replication rate measured per day per 10 6 cells or viruses. e response of the drug to the uninfected and infected tumor is, respectively, modelled by the functions δ 0 U(t)C(t) and δ 1 I(t)C(t) 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/c 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 g(t) 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 rate g(t) � q, (b) an exponential g(t) � q exp(−at), and (c) a sinusoidal function g(t) � q sin 2 (at), where q is the rate of drug infusion. e 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]. e 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. is 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]). e 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]. e assumptions above lead to the following system of nonlinear first-order differential equations (also similarly derived in [11,19,20]): subject to initial concentrations

Delay Model.
e 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]). e 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)

Mathematical Analysis of Model without Delay
In this section, the model without delay (1) is analysed. e variables in system (1) are first rescaled by setting t � δt, For simplicity, we drop the bars and equation (1) becomes ξ(t) � ϕ, ξ(t) � ϕ exp(−at), and ϕ sin 2 (at), 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.

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: whose solution is 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).
with U(0) � U 0 and C(0) � C 0 . e second equation in equation (8) is a first-order linear ordinary differential equation which can easily be solved to give where R is a constant of integration. e solution to the first equation in equation (8) depends on the infusion function ξ(t). For a fixed infusion function ϕ, From the solution of C(t) in equation (10), 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. e 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 Computational and Mathematical Methods in Medicine From equation (12), where U * is a fractional tumor cell concentration between 0 and 1. is 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. e drug is also completely depleted from the body. When ξ(t) � ϕ sin 2 (at) is substituted in equation (8), the resulting differential equations are solved to give where is suggests that with time, some drug concentration remains in the body tissue.
e eigenvalues evaluated at implying that this steady state is locally asymptotically stable if δ 0 ϕ > aαψ and ϕ < 1. eorems 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. e 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 e analytical solutions to system (22) are not easy to obtain. e derivatives of equation (22) are therefore equated to zero to obtain time invariant solutions and investigate their stability by linearizing equation (22)   (1) If β + c > 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) If bβ > β + c, 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: which is locally asymptotically stable if a 0 , a 1 , a 2 > 0 and a 1 a 2 > a 0 where a i are coefficients of the characteristic equation. Proof.
(1) e characteristic equation evaluated at (0, 0, 0) is 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 a 0 , a 1 , a 2 > 0 and a 1 a 2 > a 0 .

□
Since the infected tumor-free state is undesirable, the reverse of the condition β + c > bβ is necessary for tumor eradication from body tissue. In other words, bβ > β + c, 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) It is therefore evident that high virus replication rate β and burst size b lead to lower tumor cell concentrations. e Computational and Mathematical Methods in Medicine 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: 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. e eigenvalues of the Jacobian matrix evaluated at this state are 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: 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 Table 1 where It 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 Table 1, 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).
is system has three steady states: (i) A tumor-free state where all cell concentrations diminish to zero: is state is unstable because the eigenvalues −c, −ψ, −a, −1, and α are not all negative.
(ii) A state where the tumor grows to its maximum size: (iii) e characteristic polynomial evaluated at this state is equation (25) and this state is locally asymptotically stable if β + c > βc, otherwise it is unstable.
e 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 Table 1, the eigenvalues are −0.1, −8.13, 1.054, and −0.014 ± 0.085i, implying that this state is stable.
Similarly, with ξ(t) � ϕ sin 2 (at), one changes equation (5) into an autonomous system of equations by letting W � ϕ sin 2 (at). e autonomous system has six steady states: 6 Computational and Mathematical Methods in Medicine (i) Tumor-free state where all cell concentrations are wiped out of body tissue: e eigenvalues evaluated at this state are −c, −ψ, −a, −1, and α implying that it is unstable.
(ii) A state where the tumor grows to its maximum size: e condition for stability of this state is same as with the exponential drug infusion case, that is, the state is locally asymptotically stable if β + c > βc.
(iii) Tumor-free state where some concentration of the drug remains: is state is locally asymptotically stable if δ 0 ϕ > αψ and (1/2) < ϕ < 1 because the eigenvalues evaluated at this state are otherwise, it is unstable. (iv) Infected tumor-free state where all infected tumor cells are wiped but a certain proportion of the uninfected remains: (v) Drug-free state where the chemotherapeutic drug is wiped out of body tissue and proportions of all the other cell concentrations remain: (vi) Tumor dormant state: where e conditions for stability for the last three states all depend on huge expressions from which it is hard to extract meaningful biological implications.
is 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.
For τ ≠ 0, we have a transcendental equation to solve. □

Parameter Values.
e parameter values used were obtained from fitted experimental data for untreated tumors and virotherapy in mice [41]. e tumor carrying capacity, K, is taken to be 10 6 cells per unit volume. e 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].

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 U 0 � 1, I 0 � 0, V 0 � 0.1, and C 0 � 0.1 with a high fractional untreated tumor cell count to necessitate clinical intervention. e 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 � 10 6 number of cells. Figure 1 shows numerical solutions of the chemo-only model (8). e figure shows that despite the drug infusion method, the tumor is not cleared from body tissue. ese 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. 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. ese numerical intimations concur with the analytical results established in Section 3.1. Figure 3 displays the dynamics of model (5)  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. e chemotherapeutic response delay is higher than the virus response delay, thus τ 2 was taken to be between 0.1 and 0.3. e 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. e 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.

Discussion
e 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][16][17][18][19][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). ey investigated the efficacy of different procedures of a measles virus in combination with fludarabine and cyclophosphamide (CPA). eir 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. eir 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]. e duration of cancer treatment depends on several factors including the type of cancer being treated and the patient cells' characteristics. is 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]. e 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.

Conclusion
e 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. e 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 tumorfree situation may be achieved. Numerical simulations for the ODE and DDE models were, respectively, carried out using ode23s and dde23 in MATLAB. e model analysis suggested the following:

A.1. Existence and Uniqueness
Theorem 4. ere exists a unique solution to the system of equation (5) in the region (U, I, V, C) ∈ R 4 + .

Theorem 7.
e domain D is positive invariant for the model equation (5) and therefore biologically meaningful for the tumor, virus, and drug cell concentrations. Proof.
e proof directly follows from proofs of eorems 5 and 6.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e author declares that there are no conflicts of interest.