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.

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 [611]. 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, 1520]. 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, 2330]. 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: 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. 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 cells , infected tumor cells , a free virus population , and a chemotherapeutic drug 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 , where β is the virus replication rate measured per day per cells or viruses. The response of the drug to the uninfected and infected tumor is, respectively, modelled by the functions and where and are induced lysis rates caused by the chemotherapeutic drug measured per day per cell. Virus lifespan is taken to be and its production is considered to be 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 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 , (b) an exponential , and (c) a sinusoidal function , 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]):subject to initial concentrations

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 and , respectively, model (1)

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 , , , , and . Taking , the parameters are renamed to become

For simplicity, we drop the bars and equation (1) becomes

, , and , 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:whose solution isimplying 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 (), the system (5) is transformed towith and . The second equation in equation (8) is a first-order linear ordinary differential equation which can easily be solved to givewhere is a constant of integration. The solution to the first equation in equation (8) depends on the infusion function . For a fixed infusion function ϕ,

From the solution of 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. 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 ,where .

From equation (12),where 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 is substituted in equation (8), the resulting differential equations are solved to givewhere

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

Proof. Using Dulac’s criterion ([40]):
Suppose and is continuously differentiable on a simply connected domain . If there exists a real valued function such that has one sign in , then there are no closed orbits in . Using Dulac’s criterion, it is sufficient to show thatConsider

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 , there are two steady states of equation (8): which is locally asymptotically stable provided that and which is locally asymptotically stable provided that ; otherwise, it is unstable.(2)For the exponential drug infusion , equation (8) has two steady states: which is unstable and which is locally asymptotically stable.(3)For the sinusoidal infusion function, there are four steady states of equation (8): , , and which are unstable and which is locally asymptotically stable if and .

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 iswhose roots λ can only be negative if .The characteristic polynomial evaluated at iswhose roots are negative provided that .(2)By letting , equation (8) is turned into an autonomous system:The eigenvalues of the Jacobian matrix for equation (20) evaluated at which are and 0 and at which are and are all negative.(3)Similarly, by letting , equation (8) becomes the autonomous system:The eigenvalues of the Jacobian matrix for equation (20) evaluated at are and 0 and the eigenvalues evaluated at are and 0. For the third steady state to exist, , and the eigenvalues evaluated at this state are , implying that for it to be locally asymptotically stable, ; yet for this to happen, the steady state will not exist. The eigenvalues evaluated at are implying that this steady state is locally asymptotically stable if and .
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 toThe 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 , the system (12) has two steady states: a tumor free cell state which is unstable and an infection tumor free state which is locally asymptotically stable.(2)If , the system (22) has three steady states: the tumor free state and the infected free state which are unstable and a tumor dormant state:which is locally asymptotically stable if and where are coefficients of the characteristic equation.

Proof. (1)The characteristic equation evaluated at isfrom which , and , thus rendering it unstable.The characteristic equation evaluated at isfrom which and which are all negative since .(2)The characteristic polynomial evaluated at the tumor dormant state is whereand are the coordinates of the tumor dormant state. Using Routh–Hurwitz stability criterion, this state will only be locally asymptotically stable if and .

Since the infected tumor-free state is undesirable, the reverse of the condition is necessary for tumor eradication from body tissue. In other words, , 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

It is therefore evident that high virus replication rate and burst size 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:

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 areimplying that this desirable state is locally asymptotically stable if 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, the eigenvalues of the Jacobian matrix evaluated at the infected tumor-free state are 0.403, and , implying that the infected tumor-free state is unstable.(iii)Tumor dormant state: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 , , , and , implying that this state is stable.

With the consideration of an exponential infusion function , the equation (5) are first turned into an autonomous system of differential equations by letting . This system has three steady states:(i)A tumor-free state where all cell concentrations diminish to zero:This state is unstable because the eigenvalues , , , , and α are not all negative.(ii)A state where the tumor grows to its maximum size:(iii)The characteristic polynomial evaluated at this state is equation (25) and this state is locally asymptotically stable if , otherwise it is unstable.

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 Table 1, the eigenvalues are , , 1.054, and , implying that this state is stable.

Similarly, with , one changes equation (5) into an autonomous system of equations by letting . The autonomous system has six steady states:(i)Tumor-free state where all cell concentrations are wiped out of body tissue:The eigenvalues evaluated at this state are , , , , and α implying that it is unstable.(ii)A state where the tumor grows to its maximum size: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:This state is locally asymptotically stable if and because the eigenvalues evaluated at this state areotherwise, 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

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

Without virotherapy the model (44) in nondimensional form becomes

The system (45) has two steady states: a tumor-free state and a tumor dormant state , just as previously seen. Letting and where are steady states of equation (45), the linearized model about of equation (45) is

The characteristic equation of equation (46) evaluated at the tumor-free steady state is

For , we obtain the same characteristic polynomial as in the ODE case (Theorem 2). For , 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 state is stable if and , i.e., Equation (47) has a negative root and a zero double root.(2) is stable for a sufficiently small , i.e., Equation (47) has negative real roots for .(3)Equation (47) has a pair of purely imaginary roots if and the other root is . Therefore, is unstable.

Proof. (1)If and , thenForm equation (48),Thus, equation (47) has a double zero root.(2)If we denote as the root of the equation (47), the tumor-free state is stable if . By continuity, if is sufficiently small, we still have and the tumor-free state is stable.(3)If equation (47) has only purely imaginary roots, then the roots should be solutions to the equationAssume that , is the root of equation (50). Substituting and separating real and imaginary parts, one getsSquaring both equations and adding givesfrom which if and the critical values of areWithout chemotherapy equation (44) becomesBy letting , and where , and are steady states of equation (54), the linearized model about of equation (54) isThe characteristic equation of equation (55) at isFor , we retrieve the same results as in Theorem 3. For , 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 cells per unit volume. The number of viruses produced per day b is to be in the range [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 , , , and 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 days and 1 unit number of cells.

Figure 1 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 . Nonetheless, the parameter values used do not conform to this condition.

Figure 2 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 , 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 , it takes about 10 days to reduce the whole tumor concentrations to zero while it takes only about 5 days with , 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 3 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 when and this increased to when .

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 and , the virus, and chemotherapeutic delays. Since secondary transcription and viral protein synthesis can be delayed by about six to eight hours [37], was considered to be between 0.001 and 0.01. The chemotherapeutic response delay is higher than the virus response delay, thus 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 in Figure 4(a), it took about 4 days for the whole tumor to clear whereas it took about 8 days with in Figure 4(b). Comparing Figures 4(b) and 4(c), when 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.

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,1520]. 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.

Appendix

A.1. Existence and Uniqueness

Theorem 4. There exists a unique solution to the system of equation (5) in the region .

Proof. The Picard–Lindelöf theorem [45] is used as follows.
Consider the closed interval and the closed ball where T and d are positive, real numbers. Suppose that the functionis continuous and that the partial derivatives in the Jacobian matrix whereexist and are continuous in . Then, there exists a so that the initial value problemhas a unique solution on the interval It is sufficient to show thatexist and are continuous on . The functions (A.4) and (A.5) are polynomials and are therefore continuous on .

A.2. Boundedness and Positive Invariance

Theorem 5. if , , , and , then , , , and for all .

The same idea as in [10] is used to prove positiveness of the model solutions.

Proof. Assuming that the Theory is not true, then there must be a time such that at least one of the solutions becomes zero first. Each possible case is investigated; if first, then . However, from the first equation in system (5), by the uniqueness of the solution, we know that for all . The second equation then becomes and its solution isThe third equation becomes . If you setSimilarly, the fourth equation becomes and its solution iswhich implies that for all .
If first, then , implying that when , since as assumed.
If first, , implying that when , since as assumed.
If first, , so when , since as assumed.
If two solutions are zero (eg., and ) simultaneously at , then following the same steps above, it is trivial to check that the other solutions will be nonnegative for .
If three solutions are zero (eg., , , and ) simultaneously at , it is trivial to check that the other solution will be nonnegative for .
If the four solutions are zero simultaneously at , from the uniqueness theorem, for .

Theorem 6. The trajectories evolve in an attracting region where depends on the drug infusion function used.

Proof. From equation (1), we know that . This implies that .where is an arbitrary constant of integration. For the constant infusion function ϕ, . For , and for , .

Theorem 7. The domain is positive invariant for the model equation (5) and therefore biologically meaningful for the tumor, virus, and drug cell concentrations.

Proof. The proof directly follows from proofs of Theorems 5 and 6.

Data Availability

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

Conflicts of Interest

The author declares that there are no conflicts of interest.

Acknowledgments

This work was partially funded by the University of Kwazulu-Natal and the University of Eswatini.