Modelling the spatiotemporal dynamics of chemovirotherapy cancer treatment

ABSTRACT Chemovirotherapy is a combination therapy with chemotherapy and oncolytic viruses. It is gaining more interest and attracting more attention in the clinical setting due to its effective therapy and potential synergistic interactions against cancer. In this paper, we develop and analyse a mathematical model in the form of parabolic non-linear partial differential equations to investigate the spatiotemporal dynamics of tumour cells under chemovirotherapy treatment. The proposed model consists of uninfected and infected tumour cells, a free virus, and a chemotherapeutic drug. The analysis of the model is carried out for both the temporal and spatiotemporal cases. Travelling wave solutions to the spatiotemporal model are used to determine the minimum wave speed of tumour invasion. A sensitivity analysis is performed on the model parameters to establish the key parameters that promote cancer remission during chemovirotherapy treatment. Model analysis of the temporal model suggests that virus burst size and virus infection rate determine the success of the virotherapy treatment, whereas travelling wave solutions to the spatiotemporal model show that tumour diffusivity and growth rate are critical during chemovirotherapy. Simulation results reveal that chemovirotherapy is more effective and a good alternative to either chemotherapy or virotherapy, which is in agreement with the recent experimental studies.


Introduction
Current cancer treatments involve combination therapies such as radioimmunotherapy [29,58], radiovirotherapy [15,59], and immunotherapy combined with targeted therapies [23,41,63], to name but a few. Presently, despite aggressive treatments including combination therapies, most cancers often recur due to their resistance to conventional therapies and limitations to effective therapies [60]. Recently, chemovirotherapy, a combination therapy with chemotherapy and oncolytic viruses, has gained increasing significance in the clinical setting. The essence of using chemovirotherapy is that oncolytic viruses directly lyse tumour cells or deliver genes that make them more susceptible to chemotherapeutic drugs. In fact, many studies have concluded that the combination of oncolytic viruses with chemotherapy acts synergistically mainly because they are mediated by different pathways [1,3,52].
Numerous experimental studies have shown that, with the right combination of oncolytic viruses and standard cytotoxic chemotherapy agents, synergistic interactions occur resulting in enhanced therapeutic effects not attainable by either therapy alone [6,43,45,49]. Chemovirotherapy has been experimentally and clinically tested for the treatment of some cancers such as malignant gliomas [3,61], lymphoma [62], lung cancers [23], and metastatic breast cancer [12]. The study by Alonso et al. [3] showed that the combination of an oncolytic adenovirus (ICOVIR-5), with either everolimus (RAD001) or temozolomide (TMZ) resulted in an enhanced increase in the anti-glioma effect in vitro and in vivo glioma xenograft model. They concluded that the animals' median survival rate has increased by 20-40% and that they remained disease free beyond 90 days after treatment. Another experimental study on malignant glioma by Ulasov et al. [61] suggested that virotherapy and temozolomide chemotherapy, when combined, are able to eradicate malignant glioma. Their results showed that 90% of the treated mice survived beyond the 100 days' mark after being treated. In [62], Ungerechts et al. examined the synergy between a reprogrammed oncolytic virus and two chemotherapeutics in the mantle cell lymphoma (MCL). They investigated the efficacy of different regimens of a reprogrammed measles virus into an oncolytic virus in combination with two chemotherapeutics, fludarabine and cyclophosphamide (CPA) approved for treatment of lymphoma in an MCL xenograft model. Their study showed that CPA administration before virotherapy enhanced oncolytic efficacy. They concluded that a three 23-day courses of triple sequential treatment with CPA, virus and fludarabine treatment resulted in complete regression of the xenografts.
A recent experimental study by Gomez-Gutierrez et al. [23] showed that the chemotherapeutic drug (TMZ) enhances virotherapy in three lung cancer cell lines, concluding that combined therapy of the oncolytic adenovirus (adeAdhz60) and TMZ has a synergistic killing effect in vitro. In [12], Cody et al. presented a summary of oncolytic virotherapy strategies that have been used in a number of combinatorial therapeutic strategies to increase their effectiveness against breast cancer. They particularly reviewed the recent studies in which virotherapy has been combined with chemotherapeutic agents to target metastatic breast cancer. They communicated that all of these studies concluded that the antitumour efficacy of oncolytic viruses can be enhanced by chemotherapeutics agents. In addition, they discussed the challenges facing oncolytic virotherapy of metastatic breast cancer such as enhancing systemic delivery, promoting efficient intratumoral spread (overcoming matrix barriers, diffusion gradients, and poor viral replication), and limiting the antiviral immune response.
In [6], Binz and Ulrich gave a thorough and up-to-date review on the clinical studies in which the concepts of chemotherapy and virotherapy have been combined. They stated that phase II/III clinical trials on combining the adenovirus H101 with the chemotherapeutic drug cisplatin and 5-fluorouracil showed a 40% improvement compared to chemotherapy alone for the treatment of patients with head and neck cancer. They also communicated that these clinical studies have led to an early classification of chemovirotherapeutic combination regimens. In 2005, China approved adenovirus Ad-H101 for the treatment of head and neck cancers after phase III clinical trials showed a higher response rate for Ad-H101 combined with 5-fluorouracil (79-72%) compared to chemotherapy alone (40%) [21]. Ad-H101 became the first oncolytic virus to be approved for the treatment of human cancers. Ten years later, in 2015, the U.S. Food and Drug Administration approved Imlygic (talimogene laherparepvec), a genetically modified live oncolytic herpes virus therapy, as first oncolytic virotherapy for the treatment of melanoma lesions in the skin and lymph nodes [17,18]. Earlier this year, it was also approved in Europe for some inoperable melanomas [30]. Most of these (and other studies not mentioned in this paper) have validated the general principle that oncolytic virotherapy could be harnessed as a component of a multi-modality cancer therapy.
Despite a large body of experimental and clinical work on oncolytic virus therapy, the mechanism of action of chemovirotherapy is still not well understood. Primarily, this is due to the complexity and the high non-linearity of the dynamics between tumour cells, oncolytic viruses, and their microenvironment [43,45]. There is increasing evidence that mathematical models can be helpful in investigating the interactions between tumour cells and viruses. Indeed, various yet limited, mathematical models of oncolytic virus therapies have been developed to examine the emerging properties of these dynamics, suggest optimal treatment strategies as well as to inform the design of new experiments. For a comprehensive review on mathematical models of virotherapy, the reader is referred to excellent papers, reviews, and the references therein in [14,20,25,26,28,[66][67][68][69]. Here, we restrict our literature review to few key papers in modelling virotherapy treatments that have motivated the formulation of our model.
Tian [57] presented a basic mathematical model in the form of ordinary differential equations (ODEs), for virotherapy treatment which incorporated virus burst size. This mathematical model was based on an earlier basic virotherapy models of Wodarz and Komarova [70]. Tian concluded that besides the virus type, the burst size of an oncolytic virus and the tumour size determine the outcome of oncolytic treatment. Crivelli et al. [13] in an analysis of the cell cycle-specific activity of viruses noted that intracellular viral replication and the virus-tumour cells fusion are among the parameters that substantially impact the outcome of virotherapy. These parameters are equivalent to burst size and viral infection rate in our model.
In [36], Malinzi et al. presented a mathematical model to describe the interactions of the tumour-immune-virus dynamics. Their model included local kinetic interaction terms of the tumour and immune cells and a modified functional immune response to account for the saturation of immune cells in a tumour. The results from the study showed that the main virotherapy treatment properties were tumour cell movement and local kinetic interaction terms such as tumour growth and death rates. Unlike Malinzi et al. [36] who examined a one-dimensional spatial domain, in this paper we consider a three-dimensional spatial tumour domain, under the assumption of radial symmetry, thus providing for more realistic tumour geometry. Also, it is important to note that the immune system plays a very important role in tumour cells eradication from body tissue [16,33,39,40]. We, however, for the sake of model tractability, do not consider immune cell interactions.
In this paper, we extend the growing literature on tumour virotherapy models by presenting a mathematical model of chemovirotherapy which builds upon those presented by Tian [57] and Malinzi et al. [36]. To this end, we develop and analyse a mathematical model which combines chemotherapy and virotherapy treatments. The model includes uninfected and infected tumour cells, free virus particles, and chemotherapeutic drug. The aim is to investigate the spatiotemporal distribution of tumour cells and viruses as well as to determine the outcome of chemovirotherapy treatment. Moreover, through travelling wave solutions, numerical simulations, and sensitivity analysis, we intend to establish the key parameters that promote cancer remission during chemovirotherapy treatment.
This paper is organized as follows. In Section 2, we introduce the mathematical model, define the initial and boundary conditions, and re-scale the model. In Section 3, we set the baseline values for the parameters for simulating the model. In Section 4, we present the temporal model analysis and simulations. In Section 5, we determine travelling wave solutions to the spatiotemporal model and simulate the model in the cases of no treatment, chemotherapy, virotherapy, and chemovirotherapy. In Section 6, we perform a sensitivity analysis on the model parameters to determine which parameters have the greatest impact on model outcomes. The conclusion and discussion are presented in Section 7.

Model construction
In this section, we devise a mathematical model that describes an avascular solid tumour growth under chemotherapy and virotherapy treatments in spherical coordinates under the assumption of radial symmetry. We consider an avascular tumour, with a necrotic core containing dead cells with a radius L 0 , in a radially geometric setting with a fixed radius L. We consider a tumour nodule that has grown to its maximum size of about 1-3 mm diameter [40] and just prior to angiogenesis, an assumption necessary for clinical intervention. We model the movement of cells using the Kolmogorov equation [31], a linear diffusion model used for simulating cell movement. The variables used to describe the tumour progress and its interaction with the oncolytic virus and chemotherapy are: where r ∈ [0, L] is the radius of the tumour at a time t ∈ [0, ∞). Figure 1 depicts the interactions between tumour cells, virus, and a chemotherapeutic drug.

Tumour density, U(r, t), and I(r, t)
The following assumptions are made in order to model the evolution of a tumour over time and in space: • The tumour is considered to grow logistically at an intrinsic rate α per day and its carrying capacity K is taken to be 10 6 cells [40]. • The infected tumour density, I(r, t), increases as the oncolytic viruses multiply in the uninfected tumour cells and burst. • Virus infection of the tumour is considered to be of Michaelis-Menten form [64]. • The tumour-virus interaction is one to one, that is, one virus infects one tumour cell which later multiplies into infected tumour cells [57]. • The chemotherapeutic drug is administered as a single bolus and its concentration in the blood stream decays exponentially [55]. • The drug kills the tumour cells in a concentration-dependent manner [19], that is, the drug cytotoxicity increases with increasing drug concentration, asymptotically approaching its maximum.
The uninfected and infected tumour cell densities are governed by the following reaction-diffusion equations: In Equations (1) and (2), the terms D 1 (1/r 2 )(∂/∂r)[r 2 (∂U/∂r)] and D 2 (1/r 2 )(∂/∂r) [r 2 (∂I/∂r)], respectively, represent diffusion in spherical coordinates, under radial symmetry, for the uninfected and infected tumour cells to model tumour movement. D 1 and D 2 are, respectively, the uninfected and infected tumour diffusion constants. The second term αU(1 − (U + I)/K) represents tumour progression, where α is the intrinsic tumour growth rate and K is the carrying capacity. The third term −βUV/(K u + U) models infection of tumour cells by the virus, where β is the infection rate measured per day per number of cells. It is of Michaelis-Menten form to account for saturation in the virus proliferation [36,64]. The terms −δ 0 UC/(K c + C) and −δ 1 IC/(K c + C), respectively, represent chemotherapeutic drug responses to uninfected and infected tumour cells where δ 0 and δ 1 are lysis-induced rates measured per day. K u and K c are, respectively, Michaelis-Menten constants which relate to lysis rates when the virus and the drug are half-maximal. The term βUV/(K u + U), in Equation (2), describes a gain in infected tumour cells due to infection by the virus. The last term in Equation (2), that is −δI, represents the death of infected tumour cells where δ is the death rate.

Virus density, V(r, t)
Virus production is dependent on its burst size, that is, the larger the burst size, the higher the number of viruses produced [57]. Virus production is taken to be bδI, where b is the virus burst size measured per day and δ is the infected tumour cell death rate measured per day [57]. The virus gets deactivated in body tissue at a rate γ per day. The non-linear differential equation describing the virus density is given by: In Equation (3), the term D 3 (1/r 2 )(∂/∂r)[r 2 (∂V/∂r)] represents virus diffusion, where D 3 is the virus diffusion constant. Viruses reproduce by engulfing tumour cells and bursting thus releasing other copies. The term −βUV/(K u + U) describes loss of free virus due to infection of the uninfected tumour cells. The last term in Equation (3), that is −γ I, represents virus decay from body tissue.

Drug density, C(r, t)
Drug penetration into the tumour is modelled using reaction-diffusion. The drug density is governed by the following reaction diffusion equation: Similarly, in Equation (4) the term D 4 (1/r 2 )(∂/∂r)[r 2 (∂C/∂r)] represents chemotherapeutic drug diffusion, where D 4 is the drug diffusion coefficient, C b (t) = ψ e (−kt) is the prescribed drug plasma level in the blood stream, and −μC is the natural drug concentration decay, where μ is the rate of decay. The constant ψ is the rate of drug infusion and k relates to the chemotherapy drug half life, T 1/2 , which is roughly one day [50], and is given by k = ln 2/T 1/2 . Figure 1 displays the tumour-virus-chemotherapy interactions and Table 1 gives a summary of all the parameter descriptions, their baseline values, units, and sources.

Initial and boundary conditions
This section gives an account for the choice of the initial and boundary conditions for the model. A tumour nodule comprises three layers: a necrotic core containing dead cells, and quiescent and proliferating zones comprising living cells [11]. We therefore consider the uninfected tumour density in the necrotic core to be zero and that it increases outwards, that is, in the quiescent and proliferating zones [36]. Initially, the drug and virus densities are considered to lie on the sheath of the tumour in small concentrations. This is because the boundary in which we investigate the tumour-drug-virus interactions is the tumour sheath. We consider no flux boundary conditions for all cell concentrations at r = 0, the centre of the tumour, because it contains dead cells. At the boundary r = L, we consider no flux conditions for the infected and uninfected tumour cells since the tumour is considered to be at avascular growth stage with no small amounts of cells escaping through the tumour sheath to spread into the surrounding tissue [2,40]. Zero flux boundary conditions are assumed for the virus density, at r = L, because of the cancellation of flux from inside and outside the tumour [46]. The chemotherapeutic drug diffuses into the tumour through the outside tumour boundary r = L whose permeability is denoted as ρ and the virus density at the boundary is determined by the tumour-virus interactions. With these assumptions we close off the system (1)-(4) with the following initial and boundary conditions: where L 0 is the radius of the necrotic core.

Model re-scaling
The system (1)-(4) is re-scaled by settingŪ = U/K,Ī = I/K,V = V/V 0 , andC = C/C 0 witht = t/t 0 andr = r/r 0 where U 0 , I 0 , V 0 , and C 0 are initial cell concentrations and the space variable, r, is re-scaled relative to the length of the region under study, that is, r 0 = 1 [40]. The parameters become Dropping the bars for clarity and taking the necrotic core to be of radius 0.2 mm [22] gives a re-scaled model defined by the following parabolic system of non-linear reaction diffusion equations: with the initial and boundary conditions (11) and (12). and

Parameter estimation and baseline values
To date, there is still no available longitudinal data for chemovirotherapy treatment. To simulate our model, we first determine the baseline parameter values from the literature that most correspond to available experimental chemotherapy, virotherapy data, and biological facts. A summary of the parameter descriptions along with their numerical values and references is given in Table 1.

Diffusion coefficients
The diffusion coefficients for the tumour, D 1 and D 2 , were estimated using the formula, D = 4R 2 q, where R and q are, respectively, the tumour radius and growth rate, discussed in Prigogine and Lefever [48]. These values are approximated to be q = 0.1-1 day −1 and R = 4 µm. Therefore, D 1 and D 2 are approximated to be 10 −6 -10 −3 cm 2 day −1 . The virus diffusion coefficient, D 3 , was estimated using Einstein's formula D = κT/(6π Rτ ), where κ is Boltzmann's constant, T is the temperature, R is the radius of a virus cell, and τ is the viscosity coefficient of the medium [40]. The temperature and viscosity coefficients are known to be T = 310 K, and τ = ρ water . We consider the radius of the virus to be a fraction of that for the tumour. Therefore, R = 1 − 4 µm. Virus diffusivity, D 3 , is therefore approximated to be 10 −5 cm 2 day −1 . Since the oncolytic virus is laboratory made to be injected into body tissue, the drug diffusion coefficient is considered to be in the same range as that of the virus.

Tumour growth and carrying capacity
The carrying capacity of the tumour is taken to be 10 6 cells per unit volume because a tumour nodule contains about 10 5 -10 9 tumour cells [54]. The tumour growth rate was taken to be 0.26 day −1 because several experimental tumour growth models estimate it to be in the range of 0.1-1 [5].

Virus production and infection rate parameters
Viral burst size is the number of viruses released from each infected tumour cell. This alludes to the fact that virus emission leads to either cell lysis or cell death. The production of new viruses from infected tumour cells therefore occurs as a burst of a characteristic numbers of viruses, on time scales ranging from minutes to days [9]. Burst sizes for different viruses have a large range corresponding to their genus. The number of oncolytic viruses produced in a day is considered to be in the range 2-1000 [9]. The infected tumour cell death is estimated by Bajzer et al. [4] to be 0.5115 day −1 .

Chemotherapeutic drug parameters
The rate of drug infusion into body tissue is taken to be 20 mg/day since most cancers require high doses of treatment and the drug decay rate is estimated to be 4.17 mg/day, values which concur with cancer pharmacokinetic studies [1,7,44]. Since infected tumour cells multiplication is enhanced by the oncolytic virus replication, the tumour cells lysis is considered to be greater than that for uninfected tumour cells [47]. The half life of several chemotherapeutic drugs is one hour to a month [65]. Therefore, we estimate the value of k to be k = ln 2/T 1/2 = 0.01 to 0.8. To predict tumour cell densities at different time periods without the consideration of space, we first analyse and simulate the model without spatial dynamics which we present in the next section.

Temporal model analysis
In this section, we investigate the temporal model's phase space properties, present its asymptotics and stability, and run some simulations to support analytical findings. Without spatial effects the model, Equations (7)-(10) become:

Phase space properties
The solutions to model (13)-(16) are cell densities and concentrations which would only make sense when they are positive quantities. It is important to first investigate solutions existence and show that trajectories which start non-negative remain non-negative and that the solutions never blow up.
where C * is the maximum value attained by the chemotherapeutic drug concentration.  (16) is a first-order linear differential equation and can be solved using a suitable integrating factor to obtain For Equations (14)- (16), we use Theorem A.4. in [56] which is stated as follows: Theorem: Let R n + = [0, ∞) be the cone of non-negative vectors in R n . Let F : R n+1 + → R n be locally Lipschitz, and satisfy Then, for every x • ∈ R n + , there exists a unique solution of The functions, F(U, I, V) on the left of Equations (14)- (16) and their partial derivatives are continuous on R 4 and (i) In (13), (16), F(U, I, V) ≥ 0 whenever t ≥ 0, (U, I, V) ∈ R n + , C = 0 (2) From Equations (13) and (14), (15), sup The solution to Equation (16) exponentially decreases with time provided that k ≥ μ and the highest possible concentration, C * , is C(t 0 ), where  (17) shows that with time the drug concentration is depleted from body tissue provided that k, a constant which relates to the drug half life, is greater than the rate of drug decay. We next determine the model equilibria and their stability to predict long-term behaviour of the model solutions.

Asymptotics and stability
In order to determine steady-state solutions of the model (13)-(16), the system equations are first transformed into an autonomous system of ODEs by letting y = C b (t) = ψ e −kt such that dC/dt = y − μC and dy/dt = −ky.

Proof:
The characteristic polynomial of the Jacobian matrix evaluated at the tumour-free steady state has both positive and negative real roots, that is, −γ , −μ, −k, α, −δ. This shows that the tumour-free state is unstable.
The other states are either stable or unstable depending on the parameter values, mostly K u , K c , and β. Equation (21) shows that with time, except for the chemotherapeutic drug, all tumour densities co-exist in body tissue. High values of the virus burst size b and the virus infection rate β lead to lower tumour densities. Moreover, Theorem 4.2 biologically implies that without the consideration of space, the temporal model predicts that the tumour may either co-exist with the virus in small concentrations or it may grow to its maximum size depending on several conditions as defined by the parameters in Equation (21).

Temporal model simulations
In this section, we present numerical simulations of the model (13)- (16) to support analytical findings and to predict changes in tumour size in response to either treatments, and to the combination of the chemotherapy with the virotherapy. The fourth-order Runge-Kutta method is used to discretise the system of equations. The dimensional values in Table 1 were re-scaled by taking t 0 = 1/δ to give the following non-dimensional parameter values: α = 0.4027, β = 0.002, δ 0 = 9.76, δ 1 = 11.73, γ = 0.002, K u = 0.01, Initially, the tumour is assumed to have grown to maximum size to necessitate treatment. The fractional concentrations are therefore taken to be U 0 = 1, I 0 = 0, V 0 = 0.1, and C 0 = 0.1.
First, we simulate the model with either treatments, then chemotherapy and virotherapy treatments combined. In this simulation, one unit on the time scale represents two days, whereas one unit on the cell concentrations scale represents 10 6 number of cells. Figure 2 displays chemotherapy treatment with time. The simulation results show that after 40 days, although the tumour decreased to its steady states when the drug is depleted, it is not completely wiped out from the body tissue. Figure 3 shows the temporal variation of the tumour and virus concentrations using virotherapy. The simulation results show that it takes about two days for all tumour cells to be infected and about ten days for all the tumour cells to reduce to undetectable cell counts. Figure 4 depicts the case of the combination of chemotherapy and virotherapy. The simulation results show that it takes less than four days for all tumour cells to be infected and die out completely.
Biologically, these simulation results imply that the treatment of cancer with chemovirotherapy is more effective than chemotherapy and virotherapy alone, which is in accordance with some experimental studies [60,71]. These results also show that it would take a shorter time period to clear all tumour cells from the body tissue which is in line with the experimental study by Weber et al. [66].

Spatiotemporal model analysis
In this section, we analyse the full model (7)-(10) to investigate the spatial dynamics of the cell densities. First, in one dimension and with constant drug infusion, we employ a technique used by Maidana and Yang [34] and Chahrazed et al. [10] to determine the minimum wave speed of the travelling wave solutions to the model (7)-(10). In doing this, we seek to determine the parameters which define the minimum speed of tumour invasion, thus determining the factors which characterize the potential with which a tumour invades body tissue. We later present results from numerical simulations of the model (7)-(10).
By letting y 1 = dU/dζ , y 2 = dI/dζ , y 3 = dV/dζ , and y 4 = dC/dζ , Equations (23)-(26) are transformed into a system of first-order differential equations: and where Y 1 and Y 2 correspond to the equilibrium solutions of Equation (28). Y 1 is the tumour-free equilibrium and Y 2 is the equilibrium solution where either the tumour grows to maximum size or where all cell concentrations co-exist. A travelling wave solution is the trajectory that joins Y 1 and Y 2 . The tumour-free equilibrium must therefore not oscillate, that is a trajectory leaving Y 1 must move to the stable equilibrium. Consequently, in determining the minimum wave speed the eigenvalues of the Jacobian matrix evaluated about Y 1 must have real values. The eigenvalues, λ i , i = 1, 2, . . . 8, of the Jacobian matrix of Equation (28) are: The eigenvalues λ 1 to λ 6 in Equations (29)-(31) are real. Equation (32) therefore determines the minimum wave speed. This is obtained by setting c 2 − 4αφ 1 = 0 which gives c = 2 √ αφ 1 . This speed highlights the parameters which are critical during chemovirotherapy treatment. These parameters are tumour diffusivity and growth rate.

Spatiotemporal model simulations
We now present numerical simulations of the model (7)-(10) with space to investigate tumour cell distributions in space with time. The multi-domain monomial based collocation method [42] and pdepe [53], a finite element-based method in Matlab, are used to determine the numerical solutions. In re-scaling the spatiotemporal model, we take t 0 = 1/D 1 which gives rise to the following non-dimensional parameter values: We simulate the model in the following cases: no treatment, chemotherapy alone, virotherapy alone, and lastly combination of chemotherapy and virotherapy as described in the following sections.

No treatment
To investigate the efficacy of each treatment and their combination, we first simulate Equation (7) which models tumour growth without any form of treatment, that is C = V = I = 0. Figure 5 illustrates the initial tumour distributions. Simulation results show that the tumour density is zero in the necrotic core (i.e. 0 ≤ r ≤ 0.2) and increases towards the tumour sheath. The infected tumour density is zero throughout the tissue since there was no virotherapy treatment. Figure 6 shows the distribution of tumour cells in the tissue without any form of treatment. The simulation results indicate that tumour density increases with time. This can be seen from Figure 6 (a,b), where after 25 days the fractional tumour density rose to 0.998 × 10 6 cells per unit volume from 0.84 × 10 6 cells per unit volume, which is about 16% increase.
Next, we simulate each of the chemotherapy and virotherapy treatments separately to determine how the drug and virus treatment affect the tumour density and spatiotemporal distribution. Figure 7 depicts the fractional distribution of the chemotherapeutic drug in body tissue. Simulation results show that the drug concentration is highest inside the tissue and that the drug concentration is reduced with time. Initially, the fractional drug concentration is maximum outside the tumour and with time it circulates in the tumour and redistributes to have a higher density inside the tissue, that it, maximum concentration is near the core than outside where r ≥ 0.2. Figure 8 (a-c) shows numerical solutions of Equation (7) which simulates the tumour with chemotherapy treatment. The tumour density in the tissue is reduced with time. When comparing the simulation results in Figures 6 and 8, we note that the tumour density was reduced from 0.849 × 10 6 to 0.790 × 10 6 cells per unit volume for a 10-day period treatment, from 0.998 × 10 6 to 0.41 × 10 6 cells per unit volume for a 25-day period treatment, and from 1 × 10 6 to 0.385 × 10 6 cells per unit volume for a 50-day period treatment. This reduction in tumour density is about 7%, 60% and 62% for the 10, 25, and the 50-day periods of treatment, respectively. Figure 9(a), (b), and (c) shows the numerical simulation of tumour treatment with virotherapy after 10, 25, and 50 days respectively. We notice that with time the virus reduces the uninfected tumour cell density from outside the tumour. Initially, we assume that there are no infected tumour cells in the tissue, that is I = 0, and the tumour has grown to maximum size. After a period of 10 days, Figure 9(a) shows that the infected tumour density begins to increase from outside the tissue, that is, 0.8 ≤ r ≤ 1 and the uninfected tumour density begins to reduce throughout the tissue. After 25 days, the infected tumour density and the region in which it lies, increased to about 0.6 × 10 6 cells per unit volume in 0.6 ≤ r ≤ 1. The uninfected tumour density was reduced to zero in this region. After 50 days, the infected tumour density occupied 0.2 ≤ r ≤ 0.6, whereas the tumour density was reduced to zero in the region 0.4 ≤ r ≤ 1.   Figure 10 (a,b) shows the spatiotemporal dynamics of the tumour with both chemotherapy and virotherapy treatments combined. The simulation results are similar to those for virotherapy treatment except that the cell densities are reduced from both sides of the domain. The simulations indicate that with time, the tumour density reduces from both sides of the tissue domain, that is, from r = 0 and r = 1. After a 50-day period of treatment, the uninfected tumour density occupies a narrower domain, 0.2 ≤ r ≤ 0.8, compared to the region, 0 ≤ r ≤ 1, that it occupied after the 10-day treatment period. Simulation results show that, when both chemotherapy and virotherapy treatments are combined, the uninfected tumour density is reduced to a lower value compared to either chemotherapy or virotherapy treatments. For example, for the case of a 50-day period of treatment, we notice that: (a) with chemotherapy alone, in Figure 8(c), the uninfected density is about 0.41 × 10 6 in the whole tumour domain 0 ≤ r ≤ 1, (b) with virotherapy alone, in Figure 9(c), the uninfected density is decreased from its highest peak of 0.9 × 10 6 to zero and lies in a reduced domain 0 ≤ r ≤ 0.4, and (c) with chemotherapy and virotherapy combined, in Figure 10(c), the highest peak of the tumour is about 0.6 × 10 6 and lies in the region 0.2 ≤ r ≤ 0.7. Going by the peak and the domain occupied by the cell densities, this means that combining chemotherapy and virotherapy results in about 20% improvement compared to virotherapy alone and about 37% compared to chemotherapy alone. We further notice that the number of infected tumour cells is highly reduced with chemovirotherapy compared to treatment with virotherapy. For example, after a 50-day period of treatment with virotherapy alone, when comparing Figures 9(c) and 10(c) we notice that the infected tumour cells are 0.5 × 10 6 and lie in the region 0.2 ≤ r ≤ 0.6, whereas with chemovirotherapy they are about 0.2 × 10 6 and covering 0.6 ≤ r ≤ 0.8. That is about 80% difference between the two treatments. The simulations show that chemovirotherapy has a double-edged sword effect, with the ability to reduce tumour cells spatially  throughout body tissue and in a small time period. The viruses infect and directly lyse tumour cells rendering them weaker and prone to lysis by the chemotherapeutic drug [1,3,52]. Moreover, the chemotherapeutic drug directly kills uninfected tumour cells as well.

Combining chemotherapy and virotherapy
These simulation results assert that combining chemotherapy and virotherapy to treat cancer is more effective and a good alternative to either chemotherapy or virotherapy alone, which is in agreement with the latest experimental studies [3,6,8,12,23,24,43,45,49,61,62]. These studies have shown that the right combination of oncolytic viruses and chemotherapy agents resulted in enhanced therapeutic effects not attainable by either therapy alone.

Uncertainty and sensitivity analysis
In this section, we study the dependence of the model solutions on the parameter values. In doing so, we are able to determine a feasible range of parameter values and ascertain the most critical parameters in the model. We employ two techniques which are discussed in detail by Marino et al. [38], that is, Latin hypercube sampling (LHS) for studying uncertainty analysis and partial rank correlation coefficient (PRCC) for calculating the sensitivity analysis indexes of the parameters. We ran the LHS/PRCC analysis with a sample size of 100. The choice of this sample size is due to the fact that PRCC produces accurate results for a lower sample size compared to other techniques like eFAST [38].
Uncertainty and sensitivity analysis were performed on all non-dimensional model parameters in the homogeneous model (13)- (16) with the aim of determining the most sensitive parameters to the model prediction of the cell densities. The parameter baseline values in Equation (22) were varied in the range of 25%. We considered the same initial cell densities as in the homogeneous model simulations. The indices were calculated at t = 5, a time just before all cell densities have reached equilibrium. Figure 11 displays a bar graph of PRCCs plotted against the homogeneous parameter values with the infected tumour density as the baseline dependent variable. The parameters which are significantly positively correlated with infected tumour cell density, at p < 0.01 level of significance, are tumour growth rate α and drug decay μ, while lysis-induced rate of tumour by the drug, δ 0 , is significantly negatively correlated. An increase in the intrinsic production of tumour cells, α, leads to higher numbers of tumour cells, thus leading to an increase in the number of infected tumour cells. An increase in the drug decay, μ, leads to a less drug concentration in the body tissue, thus reducing the effect of the drug on the tumour cells and thereby leading to an increase in the number of infected tumour cells. The higher the drug lysis rate, δ 0 , the lower the number of uninfected tumour cells, thus leading to a decrease in the number of tumour cells infected. Similarly, Figure 12 shows that the highest significantly correlated parameters with the virus density are virus burst size b which is positively correlated and lysis-induced rate of infected tumour by the drug δ 1 which is negatively correlated. An increase in the virus burst size implies an increase in the production of virus copies, thus increasing the virus density. On the other hand, an increase in the lysis rate of infected tumour leads to an increase in their death, resulting in a decrease in the number of virus copies made. Comparing Figures 11 and 12, we notice that b, δ 1 , and ψ have the highest PRCCs and therefore possess a greater influence, compared to other parameters, in determining the infected and virus densities although their ranks differ. We notice that  Figure 11. PRCCs of homogeneous model parameters with the infected cell density as the baseline variable and at the 5th day. All parameter values were varied in 25% of their baseline values in Table 1. The most sensitive parameters are shown to be α, δ 1 , b, and ψ. The p-values of α, δ 1 , δ 0 , b, and ψ are all less than 0.002. Key: alpha = α, beta = β, delta_0 = δ 0 , delta_1 = δ 1 , k u = K u , k c = K c delta = δ, gamma = γ , psi = ψ, and mu = μ.
much as other parameter values, for example, the virus infection rate, β, have a significant p-value, they have less PRCCs, thus rendering them less sensitive to the model predictions compared to α, b, δ 0 , δ 1 , and ψ.

Discussion
In this paper, we developed a mathematical model of the spatiotemporal dynamics of chemovirotherapy treatment to cancer. Unique to this model is the inclusion of both chemotherapy drugs and virotherapy treatments, which has not been investigated before by means of mathematical modelling. The model describes an avascular solid tumour growth under chemovirotherapy treatment in spherical coordinates under the assumption of radial symmetry. We investigate the spatiotemporal distribution of infected and uninfected tumour cells as well as the outcome of chemovirotherapy treatments. We use the minimum wave speed of tumour invasion, numerical simulations, and sensitivity analysis to establish the key parameters that promote cancer remission during chemovirotherapy treatment.
The model was analysed in two main phases. First without the consideration of space and later with the inclusion of a space variable. The homogeneous model solution properties were investigated and a stability analysis was carried out on the steady-state solutions. A minimum wave speed of tumour invasion, in one dimension, was determined by converting the spatiotemporal model into a system of first-order ODEs and analysing its phase space. Numerical simulations for both temporal and spatiotemporal models were carried out using experimental data from the literature summarized in Table 1 Table 1. The bar graph shows δ 1 , b, and ψ as the most sensitive. alpha = α, beta = β, delta_0 = δ 0 , delta_1 = δ 1 , k u = K u , k c = K c delta = δ, gamma = γ , psi = ψ, and mu = μ.
analysis was performed to identify critical parameters of the model. The main results obtained from this study are: (1) Chemotherapy, alone, may not be able to clear all tumour cells from body tissue. This is depicted in Figure 2 which shows that even after 40 days, the tumour concentration was only reduced to its steady state. (2) The virus burst size and infection rates are key factors in determining the success of the virotherapy treatment as was observed in Equation (21) that depicts the homogeneous model equilibriums. The homogeneous model analysis showed that high values of the virus burst size and the virus infection rate lead to lower tumour densities. This is in agreement with the mathematical studies by [13,57]. (3) As shown in Figures 3 & 4 and 9 & 10, chemovirotherapy, as a form of cancer treatment, is more effective than chemotherapy and virotherapy treatments alone and is capable of eliminating tumour cells from body tissue in less than a week. (4) The most critical parameters during chemovirotherapy according to travelling wave analysis in Section 5.1 are: tumour diffusivity and tumour growth rate. The most critical parameters in reference to sensitivity analysis in Section 6 are: virus burst size, tumour growth rate, drug lysis rates, and drug injection rate. Some of these parameters are as well depicted by the steady-state solutions in Section 4. Some of these parameters have already been identified as highly important for successful cancer treatment for in [13,35,57].
Despite the lack of longitudinal data for chemovirotherapy treatment, our analysis and numerical results are in agreement with latest clinical [6,8,51], experimental [60,71], and mathematical studies [13,25,26,36,57] on virotherapy and chemovirotherapy. Tian [57] affirms that tumour concentration can diminish to low undetectable cell counts provided that the burst size is large. Similarly, we have shown in the homogeneous case and sensitivity analysis that virus burst size and infection rates are critical in determining the outcome of virotherapy and chemovirotherapy treatments. It is worth noting that the larger the virus burst size, the higher the number of infected tumour cells, thus rendering them weak to progress and damage body tissue as shown in [57]. Crivelli et al. [13] assert that virus replication and virus-tumour cells fusion should be optimized in order to maximize virotherapy.
The mathematical model presented in this study considers a couple of simplifying assumptions and omits certain biological aspects during tumour invasion mainly the immune response which has been shown to impede tumour development [33,37,40]. Despite this fact, the results of this study affirm that combining chemotherapeutic drugs and viruses to treat cancer is more effective and a good alternative to either chemotherapy or virotherapy, which is in agreement with [6,60]. The parameter sets, highlighted in this study, could be helpful in pointing out pertinent factors to consider in designing treatment protocols for chemovirotherapy treatment. The present model may be viewed as a stepping stone towards understanding the dynamics of chemovirotherapy treatment. This study can be extended to consider a higher spatial domain such as 3D to better describe the tumour geometry. A consideration of the tumour microenvironment, such as the extracellular matrix, would be a good extension of this study to investigate the effect of chemovirotherapy on normal body tissue cells. Our model is limited to a tumour that has grown to maximum size, prior to angiogenesis. Another plausible extension would be to model chemovirotherapy during the angiogenesis stage.