Oncolytic viral therapies and the delicate balance between virus-macrophage-tumour interactions

: The success of oncolytic virotherapies depends on the tumour microenvironment, which contains a large number of inﬁltrating immune cells. In this theoretical study, we derive an ODE model to investigate the interactions between breast cancer tumour cells, an oncolytic virus (Vesicular Stomatitis Virus), and tumour-inﬁltrating macrophages with di ﬀ erent phenotypes which can impact the dynamics of oncolytic viruses. The complexity of the model requires a combined analytical-numerical approach to understand the transient and asymptotic dynamics of this model. We use this model to propose new biological hypotheses regarding the impact on tumour elimination / relapse / persistence of: (i) di ﬀ erent macrophage polarisation / re-polarisation rates; (ii) di ﬀ erent infection rates of macrophages and tumour cells with the oncolytic virus; (iii) di ﬀ erent viral burst sizes for macrophages and tumour cells. We show that increasing the rate at which the oncolytic virus infects the tumour cells can delay tumour relapse and even eliminate tumour. Increasing the rate at which the oncolytic virus particles infect the macrophages can trigger transitions between steady-state dynamics and oscillatory dynamics, but it does not lead to tumour elimination unless the tumour infection rate is also very large. Moreover, we conﬁrm numerically that a large tumour-induced M1 → M2 polarisation leads to fast tumour growth and fast relapse (if the tumour was reduced before by a strong anti-tumour immune and viral response). The increase in viral-induced M2 → M1 re-polarisation reduces temporarily the tumour size, but does not lead to tumour elimination. Finally, we show numerically that the tumour size is more sensitive to the production of viruses by the infected macrophages.


Introduction
Oncolytic viral therapies have become one of the most promising therapies for cancer, due to the ability of some viruses (i.e., oncolytic viruses) to replicate inside tumour cells without damaging normal tissue cells [50]. In addition to their ability to selectively replicate inside cancer cells, which leads to the destruction of these cells, oncolytic viruses can also trigger anti-tumour immune responses [34]. However, these anti-tumour immune responses are counterbalanced by anti-viral immune responses which eliminate the virus particles from the body [14]. Therefore, the success of these oncolytic therapies depends on a better understanding of the interactions between immune cells (and in particular cells of innate immunity) and oncolytic viruses.
Macrophages represent the first line of defence against pathogens, and they have been shown to eliminate viruses in an interferon-dependent manner [40,50,52]. However, macrophages can also be infected by viruses, and thus enhance viral dissemination and persistence [37]. The different roles of macrophages on virus elimination and/or persistence might be explained by the heterogeneity of macrophage population. These immune cells can have a variety of polarisation phenotypes, depending on the microenvironment they are in, and on the activation stimuli. The two extreme macrophage phenotypes are represented by the classically-activated anti-tumour and anti-viral M1 cells, and the alternatively-activated pro-tumour and anti-inflammatory M2 cells [31,49]. The M1-like macrophages have been shown to produce anti-viral interferons [49,55]. While many viral pathogens have been shown to activate a M1 polarisation, some viruses benefit from skewing macrophages towards an M2like phenotype [38]. The M2-like macrophages seem to act as reservoirs of replication for many viruses: from the human immunodeficiency virus (HIV) to the human cytomegalovirus (HCMV) [38]. Recent studies have shown that a promising oncolytic virus, the Vesicular Stomatitis Virus (VSV), can infect and replicate inside M2 cells but not inside M1 cells [60].
In addition to macrophages importance in anti-viral immune responses, these innate immune cells have also been found to infiltrate many types of solid tumours, and can represent between 5%-50% of tumour mass [69,73]. The tumour-associated macrophages usually have a M2-like phenotype, as tumour cells educate macrophages towards a phenotype that supports their growth. Since these tumourassociated M2-like cells can be infected by some oncolytic viruses [60,61], a better understanding of the interactions between M1 cells, M2 cells, oncolytic viruses and tumour cells is necessary to improve current treatment approaches [15].
In this theoretical study we focus on one of the most promising oncolytic viruses currently undergoing research, the oncolytic Vesicular Stomatitis Virus (VSV), which naturally infects and replicates inside cancer cells with defects in their antiviral responses (e.g., defects in the IFNγ pathway) [1]. Very recently, this oncolytic virus has been shown to replicate inside M2 macrophages but not inside M1 macrophages [60]. Moreover, the VSV seems to re-program the pro-tumour M2 macrophages towards the anti-tumour M1 phenotype [60], and this re-programming is triggered by the activation of the type-I IFN anti-viral response [58]. In regard to the cancer cell line, the experiments in [25] on VSV-macrophages-cancer interactions focused on a breast cancer cell line: MDA-MB-231. In [30] it was shown that this MDA-MB-131 cell line activates macrophages towards an M2 phenotype. Moreover, this particular cancer cell line is infiltrated by large numbers of macrophages [44]. For all these reasons, in this theoretical study we decided to focus on the same cancer cell line: the MDA-MB-231 cells.
The main goal of this study is to derive and investigate a mathematical model that could help us better understand the complex interactions between macrophages and oncolytic viruses that lead to the elimination/growth of tumours, by tacking into account also the infection of macrophages with oncolytic virus particles. To this end we generalises the model in [4] (which focused only on macrophage-virus interactions) by incorporating also tumour interactions with viruses and with M1 and M2 macrophages. To emphasise the novelty of our model, we note that the majority of mathematical models in the literature focus only on the interactions between tumour cells, oncolytic viruses and anti-tumour/anti-viral immune cells; see for example, [17,19,21,32,33,36,41,46,48] and references therein. Only a very small number of recent mathematical models investigate the infection of macrophages with oncolytic viruses, which can lead to the delivery of these viruses to specific areas of the tumour (e.g. necrotic region) [6].
Given the complexity of the new mathematical model that we propose in this study, it is impossible to focus exclusively on analytical results. Thus, here we combine some analytical results such as the identification of steady states and their linear stability analysis, with computational results, to gain a better understanding of overall model dynamics (e.g., the correlations between different model steady states that cannot be calculated analytically; the changes in these complex steady states as we vary model parameters). We also consider computational approaches to answer the following biological questions: (i) What could be the effect of macrophage polarisation/re-polarisation on tumour-oncolytic virus dynamics? (ii) What could be the effect of VSV infection of M2 cells vs. infection of tumour cells on overall tumour growth/decay? (iii) What could be the effect of VSV replication inside macrophages vs. replication inside tumour cells on overall tumour growth/decay?
Note that the last two questions refer to slightly different aspects of viral cycle: virus entry into the cells (i.e., infection rate, which can be reduced due to physical barriers inside solid tumours) and virus proliferation inside the cells which culminates with cell burst and the release of new virus particles.
The paper is structured as follows. Section 2 focuses on the description of a mathematical models for the interactions between tumour cells, M1 and M2 macrophages and oncolytic VSV particles. In section 3 we present some numerical results for the baseline dynamics of this model, as well model dynamics when we vary the infection rates of macrophages and tumour cells with the VSV particles, as well as the tumour-induced macrophage polarisation rate. We also investigate analytically the steady states and their stability (to get a better understanding of the long-term behaviour of this model), and use numerical approaches to shed some light on the analytical results that are difficult to obtain. We conclude with a summary and discussion in section 4.

Model description
To investigate the innate immune responses generated by macrophages roles on the anti-tumour oncolytic viral therapies (with VSV), we extend the mathematical model derived in [4], by considering the presence of a tumour population. Thus we focus on the following interacting cell populations: the density of uninfected tumour cells (T u ), the density of virus-infected tumour cells (T i ), the density of virions (V), the density of M 1 macrophages (which are resistant to VSV infection [61]), the densities of uninfected M 2 macrophages (M 2u ) and the VSV-infected M2 macrophages (M 2i ) [61]); see also 1. The time evolution of these variables is described by the following equations: (2.1f) Equations (2.1a)-(2.1f) incorporate the following biological mechanisms: • The uninfected tumour cells (see Eq (2.1a)), grow at a rate r up to a carrying capacity K 1 . We choose logistic growth since experimental studies have shown that tumour growth slows down as the tumour becomes very large and it depletes the nutrients [29,39]. We assume that the virus particles can infect the tumour cells at a rate β 1 , and the infection term is bilinear (i.e., the infection rate per virus and per uninfected cell is constant; for similar assumptions see [56,57,59,77]). The uninfected tumour cells can be eliminated by M1 macrophages at a rate d u [79]. At the same time, the anti-tumour immune response generated by these M1 cells can be inhibited by the presence of M2 macrophages [67]. Finally, tumour proliferation can be enhanced, at a rate d m , by the presence of M2 macrophages in the tumour microenvironment [2]. • The infected tumour cells (see Eq (2.1b)) are lysed at rate δ i1 following viral replication and cell burst. As for the uninfected cells, we assume that these infected cells can be eliminated by the M1 Figure 1. Graphical description of the possible non-spatial interactions between the tumour cells, oncolytic viruses and M1/M2 macrophages, as given by Eqs (2.1a)-(2.1f). The model was inspired by the experimental studies in [25,45,62], which showed that the VSV infects macrophages (but only the M2 cells, and not the M1 cells) in addition to the tumour cells, and the mathematical modelling studies in [18,19], which focused only on the anti-viral effect of M1 cells and did not consider the infection of macrophages.
macrophages at a rate d i . Again, we assume that the presence of M2 cells can inhibit the anti-viral effect of M1 cells [66]. • The M1 macrophages (see Eq (2.1c)), are activated by uninfected tumour cells at a rate a u 1 , by viral antigens at a rate a v 1 and by virus-infected tumour cells at a rate a i 1 . We assume that the M1 cells can proliferate logistically at a rate p m1 through a self renewal process [31], up to a maximum carrying capacity K 2 . This type of growth depicts experimentally observed cell kinetics [11]. The M1 cells (which resist to VSV infections [61]) can polarise towards a M2 phenotype at a very small constant rate r 0 m1 (in response to cytokines such as IL-4, IL-10 [3] that can be produced by different healthy and immune cells in the microenvironment). We also assume that the presence of the tumour leads to an enhanced M1→M2 polarisation at a rate r u m1 (since tumour cells produce large amounts of TGF-β, which is known to induce a M2 polarisation [27]). The M2→M1 re-polarisation of macrophages is assumed to occur at a very small constant rate r 0 m2 due to cytokines such as IFN-γ, IL-2 [3] produced by different types of healthy and immune cells in the environment. Moreover, it has been experimentally shown in [25] that matrix (M) protein mutant (rM51R-M) VSV could modulate the switch M2→M1 (probably through the induction of IFN-γ response [58]). Furthermore, engineering oncolytic viruses which carry specific cytokines can trigger a macrophages re-polarisation to a M1-like phenotype [28]. Thus, we assume that this enhanced M2→M1 re-polarisation occurs at a rate r v m2 in the presence of the virus. The re-polarisation of M2 cells upon contact with VSV particles is described by a Michaelis-Menten term with constant h v to account for the saturated re-polarising response induced by viruses. Finally, the M1 macrophages have a natural mortality rate d e1 .
• The uninfected M2 macrophages (see Eq (2.1d)) are activated at a rate a u 2 by the tumour cells, proliferate logistically at a rate p m2 , and have a natural mortality rate d e2 . The M1→M2 and M2→M1 polarisation terms have been described above. Finally, the M2 macrophages are predisposed to infections with VSV particles at a rate β 2 [61].
• The infected M2 macrophages (see Eq (2.1e)) are lysed by the replicating viruses at a rate δ i2 . All other terms have been described above. • The oncolytic virus population (see Eq (2.1f)) proliferates when new virus particles are released by infected M2 cells and tumour cells. Parameters b and c describe the number of virus particles released by one infected tumour cell and one infected M2 macrophage, respectively. Moreover, the reduction in the number of virus particles is the result of their elimination, at a rate δ v , by the M1 macrophages. Note that, as discussed in [15], viral clearance may be prevented by the M2 macrophages. Finally, we assume that the virus particles have a natural death rate ω [16]. This last term includes also the virus elimination rate by other innate immune cells (e.g., NK cells [74]) or adaptive immune cells (e.g., T cells [12]) not considered in this study.

Remark 1.
Because many experimental studies on the proliferation rates of macrophages do not distinguish between the M1 and M2 cells [11], throughout most of this study we will assume that p m1 = p m2 =: p m .

Remark 2.
Many modelling studies have included infected tumour cells in the carrying capacity terms (see [35,48]). Since in our study we did not see any significant changes in the model dynamics when we considered the infected cells (see Figure 14 in Appendix B we have decided to ignore the infected cells in the carrying capacity terms for both tumour and macrophages (see Eqs (2.1a), (2.1c) and (2.1d)).
In the next sub-section we discuss in more detail the parameter values that we use to parametrise model (2.1).

Parameter values
In Table 1, we summarise the parameter values used throughout this theoretical study. Column 3 in Table 1 shows the dimensional values of the parameters (with units in column 4) -as obtained from the published literature (see discussion below) or estimated values/ranges. However, to avoid numerical problems caused by very large parameter values (see K 1 , K 2 in column 3 in Table 1) and very small parameter values (see β 1 , β 2 in column 3 in Table 1), we decided to rescale the cell populations by their carrying capacities (see also Appendix A): In addition, we rescaled the total virus population by V max = 10 11 , which is assumed to be the maximum viral load that does not cause neurotoxicity (see also Appendix A): V * = V/V max . This rescaling of variables leads to a rescaling of 11 parameters: = 1000 mm 3 [71]. We can assume that at this size the tumour has likely depleted the organism of nutrients (so the tumour has reached its carrying capacity). Moreover, based on the study in [24] we assume that 1cm 3 tumour tissue contains maximum 10 9 tumour cells. Thus, for this study we choose a tumour carrying capacity of K 1 = 10 9 cells/cm 3 . In this study we consider a tumour proliferation rate r ∈ (0.47, 0.93), with an average of r = 0.62. • In [22] the authors suggested that the doubling time of macrophages is around 27 hrs. In [81] the authors estimated the doubling time of untreated murine macrophage-like RAW264.7 cells to be between 18-22 hrs, while cells stimulated with bacterial lipopolysaccharide (LPS) had an estimated doubling time of 35 hrs. In [72], the authors estimated that M1 macrophages have a doubling time between 23.86 hrs and 28.97 hrs. In [11] it has been indicated that the average doubling time of macrophages is between 20-30 hrs. Therefore, from all these experimental studies we deduce that the doubling time of macrophages is likely between 18-35 hrs, which corresponds to proliferation rates between 0.4-0.9/day. For simplicity, through this study we choose the baseline proliferation rates p m1 = p m2 = 0.57/day, although we varied these rates within the interval (0.4, 0.9). • In regard to the M1/M2 macrophages natural death rates, various modelling studies used different values. For example, in [20] it was assumed that d e1 = d e2 = 0.02/day, the same as in [47]. On the other hand, in [18] the authors considered d e1 = d e2 = 0.2/day, as approximated from the experimental study in [78]. A recent experimental paper [31] suggested that the half-life of M1 pro-inflammatory murine M1 macrophages is between 18-20 hr (corresponding to a death rate between ln(2.0)/20hr and ln(2.0)/18hr), while the half-life of anti-inflammatory murine M2 cells is between 5-7 days (corresponding to a death rate between ln(2.0)/7days and ln(2.0)/5days). Therefore, for this study we choose d e1 ∈ (0.83, 0.93) and d e2 ∈ (0.099, 0.138). • It is known that in breast cancer, up to 40% tumour size can be represented by macrophages [69].
Therefore, throughout this study we assume that the carrying capacity for macrophages is K 2 = 40%K 1 = 4 × 10 8 . • In regard to the lysis rate of infected macrophages, Rager et. al [63] showed that two macrophages cell lines (clones J774.16 and C3C, derived from the murine reticulum cell sarcoma J774) were completely lysed by the virus within 1-2 days after infection. Thus, in this study we assume a death rate of δ i2 ∈ ln(2.0) 2 , ln(2.0) 1 = (0.35, 0.69)/day. For the simulations we use an average value of δ i2 = 0.52.
• In regard to the VSV burst size, [63] showed that each productively infected macrophage was able to produce viral progeny of at least 1000PFU. Moreover, Zhu et. al [80] showed that each virusinfected tumour cell produced between 50 to 8000 progeny virus particles. For our numerical simulations we assume that the average burst size of infected tumour cells is b = 2500, which is the same value as in [18]. For the burst size of infected macrophages we assume c = 2000. • The VSV death rate varies between different mathematical studies: e.g., ω = 2/day in [18,46]), or ω ∈ (1 − 2.56)/day in [21]. This is because experimental studies have shown that the intracellular half-lives of non-replicating wild type and mutant strains of VSV can vary between 5.3 hrs and Mathematical Biosciences and Engineering Volume 18, Issue 1, 764-799. 18 hrs, which translates into a death rate between ln(2)/5.3hr = 3.13 days and ln(2)/18hr= 0.92 days [16]. In [26] the authors have calculated the half-life of a replicating VSV strain between 2-15.9 hrs, which translates into a death rate between 1.04-8.31/day. In this study, we choose an average death rate of ω = 2/day. • In regard to the M1→M2 and M2→M1 baseline polarisation rates, the theoretical studies in [20,75] used r 0 m1 ∈ (0.05, 0.09) and r 0 m2 ∈ (0.05, 0.08). In [18] the authors used the baseline values of r 0 m1 = 0.001 and r 0 m2 = 0.01. Here, we assume that these two parameters vary between (10 −5 , 10 −1 ). • For the virus-induced re-polarisation rate r 2 m2 we used the same values as in [18], and thus we consider the baseline value r v m2 = 0 (meaning that by default there is no virus-induced re-polarisation).
• The rest of the parameter values (i.e., h m , h v , δ i1 , β 1 , β 2 , δ v , r u m1 , g 2 ) used in this study for the numerical simulations are listed in Table 1. The ranges of these parameters, as well as their baseline values are "guessed", since we could not find any references for these parameters.
To investigate the impact of these guessed parameters on the dynamics of uninfected tumour cells (T u ), we perform a global sensitivity analysis using the Latin Hypercube Sampling/Partial Rank Correlation Coefficient (PRCC) approach [5,51]. In Figure 2 we plot the PRCC index for each model parameter, and we see that among these guessed parameter values only β 1 and r u m1 have an impact on uninfected tumour cells T u (since β 1 ≈ −0.8, and r u m1 ≈ −0.4). Here we simulated tumour dynamics for 14 days; simulating tumour dynamics for longer time (e.g., 60-80 days) reduces the impact of r u m1 ; but β 1 still has an important role on tumour dynamics. In a similar manner we can investigate the sensitivity of M1 and M2 macrophages to model parameters; see Figures 15 and 16 in Appendix E. Both uninfected macrophage populations (M 1 and M 2 ) are most sensitive to two guessed parameters, h u and r v m2 (whose upper range is also guessed in Table 1).
Remark 3. The study in [24] suggested that the detection level for human tumours is between 10 7 −10 9 cells. We can assume that for in vivo murine experiments, the detection level is between 10 6 − 10 7 cells (since we are expecting to see tumour growth at the injection site). After the rescaling cell variables -see Appendix A -we assume that the detection level for the new rescaled tumours will be between 10 −3 − 10 −2 . We will use this information in our discussion of numerical results.

Results: steady states and their stability
We start investigating the dynamics of model (2.1) by focusing first on the long-term (asymptotic) behaviour of this model, i.e., the steady states and their stability. Due to the complexity of model (2.1), we combine analytical results with numerical simulations to gain a better understanding of this long-term dynamics. These steady-state results will be further used in section 3.2, where we will investigate numerically the transient and asymptotic behaviour of model (2.1), to explain the oscillatory and chaotic dynamics observed numerically in this model.

Steady states and their stability
Steady states. In this section we investigate the long-term behaviour of the model (2.1) by discussing all possible steady states. This steady-state investigation will allow us to understand better activation rate of M1 macrophages in response to tumour antigen [18] proliferation rate of M1 cells [11] 0.4 − 0.9 (0.57) day −1 0.4 − 0.9 (0.57) r 0 m1 baseline M1→M2 re-polarisation rate in response to antiinflammatory cytokines in the microenvironment [18]  Each parameter is sampled randomly 1000 times from the parameter ranges shown in Table 1, using a uniform distribution. We simulate tumour dynamics for 14 days (while the tumour exhibits a transient behaviour), with a time step of 1 day. The PRCC index varies between −1 and +1, and the largest PRCC index (in absolute value) corresponds to the parameter with respect to which the model outcome is most sensitive. We see that T u is most sensitive to r, followed by β 1 and r u m1 .
the numerical simulations in sections 3.2.1-3.2.2. The equilibria of (2.1a)-(2.1f) satisfy the following equations: As mentioned in the previous section, since we do not have data which differentiates between the proliferation rates for M1 and M2 cells, we assume that pm1 = pm2 := pm (see Table 1). Under these assumptions, model (2.1) exhibits two general types of equilibria: tumour-free steady states (summarised in Proposition 1) and tumour-present steady states (summarised in Proposition 2). ( given implicitly by the following equation (which is plotted in Figure 3(a) for the baseline parameter values given in Table 1): , and V * given implicitly by the following equations (the last two being plotted in Figure 3(b) for the baseline parameter values given in Table 1): From Eq (3.3b) it is clear that this T F steady state exists only when M * 2u > ω c β 2 . Proposition 2. Model (2.1) can exhibit the following two types of tumour-present steady states: where T * u , M * 1 and M * 2u are given implicitly by the following equations (also plotted in Figure 3(c), after Eq (3.4a) was substituted into (3.4b)): , where the states are given implicitly by the following equations (also plotted in Figure 3(d) after substituting T * u , T * i and V * in (3.5d)): In regard to Figure 3, where we plot the implicit expressions describing the various steady states, we note that for T V F, V F and T IV steady states, an increase in M1 cells is associated with an increase in the uninfected M2 cells. This is true also for the steady state T F (see Eq (3.3b)). Since the steady states T V F, T F, V F and T IV are too complex to calculate their closed-form expressions, in Figure 4 we show how they change as we vary two model parameters: (a) the tumour infection rate β 1 , (b) the M2 macrophages infection rate β 2 . We see that a decrease in β 1 below β 1 = 729.9 or an increase in β 2 above β 2 = 5.9 leads to the bifurcation of a TIV state (which contains tumour) from the TF state (with no tumour).
Stability of steady states. In the following we investigate analytically the stability of two of the steady states summarised in Proposition 1 -the simplest states, with no tumour and no virus. For the more complex steady states we will have to use numerical approaches to investigate their stability. We emphasise that this stability analysis will help us understand the numerical results presented in sections 3.2.1-3.2.2.
Proposition 3. The tumour-free/virus-free/macrophages-free state (T IV F) and the tumour-free/virus-free/macrophage-present state (T V F) have the following stabilities: (i) The TIVF state is always unstable. (ii) The TVF state is asymptotically stable provided that the following inequalities hold true:    Table 1. Here "BP 1 " denotes the bifurcation point where T IV steady state bifurcates out of the T F steady state as we decrease β 1 below β 1 = 729.9, and "BP 2 " denotes the bifurcation point where T IV bifurcates out of T F as we increase β 2 above β 2 = 5.9.
• Condition λ 1 < 0 leads to the following inequality between the steady state values for M1 cells and uninfected M2 cells: • Condition λ 5,6 < 0 holds true if c 11 + c 22 < 0 and c 11 c 22 − c 12 c 21 > 0. When p m1 = p m2 := p m these last two inequalities are equivalent to This completes the proof.
Remark 4. Note that the stability of TF, VF, and TIV steady states is very difficult to investigate analytically, due to the complexity of the model and the fact that all these states are defined implicitly; see Eqs (3.3)-(3.5d). Even the stability conditions (3.6) for the TVF state are not straightforward to understand, since they depend on the implicit expression (3.2) connecting M * 1 and M * 2 . To address the analytical issues highlighted in the above Remark, in Figure 5 we investigate numerically the stability of steady states exhibited by model (2.1), as we vary the two infection parameters: β 1 , β 2 ∈ (1, 1000). Since the TIVF state is always unstable, we do not show it here. We note that TVF and VF steady states are always unstable when varying β 1 and β 2 . An interesting dynamics is observed for the TF state when we vary β 2 (see panels (c)(ii), (d)(ii) and (e)). As β 2 increases towards β * 2 = 10.5 the larger amplitude limit cycle around the unstable TF point becomes unstable, and a smaller amplitude limit cycle emerges. This stable limit cycle evolves around the unstable TIV point.  , to emphasise the two Hopf bifurcation points "HB1" and "HB2", and the two limit cycles with different amplitudes that exchange stability at β 2 = 10.5. The rest of parameter values are presented in Table 1.

Numerical results
Numerical approach. To approximate numerically the ODE system (2.1), which might become stiff due to small and large parameter values (see Table 1), here we use an implicit 3rd order Rosenbrock method [65].
Initial conditions. The initial conditions for the time-evolution of model (2.1), summarised also in Table 2, were chosen to replicate some of the experimental conditions for tumour-macrophage-VSV dynamics from [25]. In [25], VSV was added to a plate with 5×10 5 MDA-MB231 human breast cancer cells, at a multiplicity of infection (MOI) of 1 or 10 PFU/cells. Here, we assume an average MOI of 2, and thus consider V(0) = 10 6 PFU. The experiment in [25] considered the same number of macrophages and tumour cells (i.e., 5 × 10 5 ). Since in this study we investigate the activation/proliferation of M1 and M2 cells, as well as the M1→M2 polarisation (as tumour progresses) and M2→M1 re-polarisation (triggered by the VSV), here we start with much lower levels of M1 and M2 macrophages, and with more M1 than M2 cells; see also Table 2. Table 2. Summary of initial conditions used for the numerical simulations of system (2.1). Remark 5. We emphasise that there are very few experimental studies that investigate the interactions between the tumour cells, macrophages and oncolytic viruses -VSV in particular [25,42,43]. Most of the experimental studies existent in the literature record the dynamics of macrophages in general, without distinguishing the anti-viral/pro-viral roles of this heterogeneous population of cells [42,43]. Among the very few experimental studies we found to distinguish the anti-tumour/pro-tumour and anti-viral/pro-viral roles of macrophages, we used [25] to generate the initial conditions for our mathematical model. Note that some of the experimental results in [25] (which did not include tumour explicitly, but focused on macrophages infection with VSV) were published very recently in [60]. This scarcity of experimental results did not allow us to parametrise model (2.1) using one single experimental study, and thus many parameters were "guessed"; see section 2.1.

Baseline model dynamics
We start the investigation into the dynamics of model (2.1) by showing in Figure 6 the baseline dynamics obtained for the parameter values listed in Table 1. In sub-panel (a) we show the time evolution of model variables, while in sub-panel (b) we graph the percentage of macrophages in the system (since various studies consider the percentage of macrophages as a predictor for tumour relapse [9]). We see that the increase in viral load leads to tumour elimination by day t = 16. The increase in M2 cells for t ∈ (20, 30) leads to an un-detected growth in tumour population. This growth is faster and becomes detectable for t > 50, and is the result in a decrease in the M1 cell population around t ≈ 50. The role of M1 and M2 macrophages in tumour relapse can be seen more clearly in sub-panel (b).   Table 1, with the initial conditions listed in Table 2

Numerical results: transient and long-term behaviour as we vary model parameters
In the following we investigate numerically the dynamics of model (2.1) as we vary the rate of tumour cells infection with VSV (β 1 ) and the rate of M2 cells infection with VSV (β 2 ). We also investigate numerically the effect of varying the macrophages polarisation (r u m1 ) and re-polarisation (r v m2 ) rates (to address question (i) in the Introduction), and the effect of varying the viral bursts sizes c and b (to address question (ii) in the Introduction).
• In Figure 7(a),(b) we see that decreasing tumour infection rate β 1 = 10 → β 1 = 1 (while fixing β 2 = 500) leads to a faster tumour relapse and a large tumour size. In Figure 7(c),(d) we see that decreasing the macrophage infection rate β 2 = 10 → β 2 = 1 (while fixing β 1 = 500) leads to tumour elimination in the presence or absence of virus-immune oscillations. • In Figure 8 we investigate the oscillatory dynamics observed above for small β 2 . To this end, we fix β 1 = β 2 = 1.4 and vary the tumour-induced macrophages polarisation rate r u m1 . We see that a decrease in this rate leads, in the long term, to a transition from regular oscillations (panels (a)) to chaotic oscillations (panels (c)). This chaotic dynamics occurs during tumour relapse; tumour seems to be eliminated before day t = 35, but relapse after day t = 50 in (b) and after day t = 150 in (c).    Table 1.  Table 1.
• In Figure 9 we investigate again the oscillatory dynamics observed above for small β 2 , and how it changes as we vary the virus-induced macrophages re-polarisation rate r v m2 . As we increase this rate, from r v m2 = 0.01 in sub-panel (a) to r 2 m2 = 10 in sub-panel (c) we see a loss in oscillations and a reduction/elimination of virus population. The virus particles are eliminated by the M1 cells, which cannot eliminate also the tumour cells.
• In Figure 10 we investigate the effect of varying the number of VSV particles released from infected M2 cells. We see that an increase in macrophages burst size from c = 1 (sub-panels (a)) to c = 80 (sub-panels (c)) leads to an increase in a transient immune-virus oscillatory dynamics,   Table 1. which occurs when tumour is not detectable. In the long term, the tumour always approaches a steady state T * u , which increases with the ratio b/c (see panel (d)). For b/c < 1 this asymptotic tumour size decreases fast towards zero. For b/c > 1 this asymptotic tumour size increases at a much slower rate. These results suggest that the non-linearity in model dynamics plays a role in tumour evolution when the same number of virus particles are released by the infected tumour cells and infected macrophages. In particular, for the parameter values investigated here, larger macrophages burst sizes (c) have a bigger impact on tumour reduction compared to the tumour burst sizes (b).
The above results, on the effect of β 1 , β 2 , r u m1 , r v m2 on tumour dynamics, are summarised in Table 3.  Table 1. . Periodic and chaotic oscillations. We have seen in Figure 8(c) that, for some specific parameter values, model (2.1) can exhibit chaotic dynamics. To have a better understanding of the short-term and long-term tumour behaviour corresponding to this case (where tumour seems to be eliminated periodically), in Figure 11 we graph the chaotic attractor corresponding to Figure 8(c), together with all the fixed points of model (2.1) (which will be discussed in more detail in section 3.1. We see that the tumour (T * u + T * i ) is decreased towards zero due to the unstable tumour-present fixed points. The unstable tumour-free (TF) steady state resides in the chaotic attractor, while all other fixed points are outside the attractor. For comparison with the above chaotic attractor, in Figure 12 we show a periodic attractor. where the tumour is decreased to lower values but detectable values. The unstable tumour-free (TF) steady state resides in the middle of the periodic orbit, with the other fixed points (all unstable) being outside this orbit. Finally, we summarise these different asymptotic dynamics in Figure 13, where we show the regions in the (β 1 , β 2 ) parameter space where we can find stable steady states (TF or TIV), stable limit cycles, and chaotic attractors.

Summary and discussions
In this paper, we considered a mathematical modelling and computational approach to investigate the complex interactions between tumour cells, an oncolytic virus -the Vesicular Stomatitis Virus (VSV) -and the innate immune response generated by the M1 and M2 macrophages, which can be found in large numbers inside different solid tumours. The novel aspect of this model is the possibility of infection of M2 cells with this oncolytic virus (as shown experimentally in [60]).
To understand the dynamics of this new (complex) mathematical model, we first focused on the steady states and their stability. Analytical results combined with numerical simulations were used to show how the steady states varied when we changed parameters related to cells' infection rates (i.e., β 1 , β 2 ; see Figures 4 and 5). We mention here that a global sensitivity and uncertainty analysis using the classical LHS-PRCC approach identified β 1 as an parameter to which the uninfected tumour cells (T u ) are most sensitive to (see Figure 2); however, throughout this study we also wanted to compare the roles of β 1 and β 2 on the overall short-term and long-term tumour dynamics. Then, we used numerical simulations to investigate the transient and long-term dynamics of this model, as we varied the infection rates (β 1 , β 2 ; see Figure 7), the tumour-induced M1→M2 polarisation rate (r u m1 ; see Figure 8), the virus-induced M2→M1 re-polarisation rate (r v m2 ; see Figure 9), and the viral burst size for tumour cells and macrophages (b, c; see Figure 10). Overall, numerical simulations showed a variety of asymptotic dynamics: from stable steady states, to stable limit cycles and chaotic oscillations (when all steady states were unstable). Given that for the parameter values used in this study (many of which were guessed, due to a lack of experimental studies on both tumour and macrophage infections with oncolytic VSV) the regular and irregular oscillations were observed also above the tumour detection level (see Remark 3). To our knowledge, this is not a biologically realistic behaviour for aggressive breast cancers, as the majority of experimental in vitro and in vivo studies on the highly aggressive MDA-MB-231 cell line show exponential cancer growth [10,71,76], in the absence of any external treatment that could be periodically administered (and thus force a periodic decay/re-growth of tumour). Nevertheless, the results are interesting from a dynamical systems point of view, since they allow us to understand the whole range of possible dynamics exhibited by model (2.1).
Returning to the biological questions we raised at the end of the Introduction section, our computational results have suggested the following biological hypotheses: (i) The increase in the macrophage tumour-induced polarisation rate r u m1 leads to a fast tumour relapse (see Figure 8): from relapse on day t ≈ 170 for r u m1 = 0.01 to a relapse on day t ≈ 70 for r u m1 = 0.1

Mathematical Biosciences and Engineering
Volume 18, Issue 1, 764-799.  Table 1): (a) Densities of cell/virus populations vs. time (for t ∈ [700, 1000] days); (b) 3D phase plane plot showing macrophages vs. virus vs. tumour populations (for t ∈ [700, 1000] days); (c) (i) Phase plane for macrophages-virus plane, (ii) Phase plane for macrophages-tumour plane, (for t ∈ [700, 1000]). In addition, we show the fixed points TIVF, TVF, VF, TF and TIV (corresponding to these parameter values.) . and continuous tumour presence when r u m1 = 1.0. The increase in the macrophage virus-induced re-polarisation rate r v m2 can lead to a tumour reduction (but not elimination) for as long as the virus is present in the system (see Figure 9). (ii) Changes in the infection rate of macrophages, β 2 (for large β 1 , which lead to tumour elimination) seems to lead to a transition between a state characterised by oscillations in immune response to a state characterised by an immune-only steady state (see Figure 7 (c),(d)). We do not know yet if such oscillations in immune response (which take place over multiple days) are realistic from a biological point of view. The increase in the infection rate of tumour cells, β 1 (for large β 2 ) leads to a delay in tumour relapse (see Figure 7  particles from the lysed cells, we note that b > c is associated with a fast tumour relapse after a temporary reduction in tumour size. In contrast, c > b is associated with an increase in tumour relapse. Moreover, the tumour seems to be more sensitive changes in macrophages burst size c compared to tumour burst size b (see Figure 10).
Since the VSV-infected macrophages do not release immediately the virus particles, being able to act as reservoirs of VSV gRNA [54,68], our current work in progress focuses on the impact of delayed release of VSV particles from the infected macrophages. Also, since infected M2 macrophages can deliver oncolytic viruses to hypoxic areas of the tumour [53] another current research direction focuses on the migration of macrophages that can be infected with VSV, and their interactions with the tumour cells.
To conclude this discussion, we emphasise that the mathematical models derived to investigate specific aspects of tumour-immune interactions and different immunotherapies for cancer are becoming more and more complex, thus acknowledging the complex biological interactions that made it impossible until now to find a cure for cancer. But this increased model complexity also forces us to rely more and more heavily on computational approaches to understand model dynamics, and to make biologically-testable predictions.  Table (1).
Appendix C: Linear stability of steady states for model (2.1) For the proof of Proposition 3 we use the Jacobian matrix associated with system (2.1a)-(2.1f). At a generic equilibrium point, this matrix is: