Non-local multiscale approaches for tumour-oncolytic viruses interactions

. Oncolytic virus (OV) therapy is a promising treatment for cancer due to the OV selective ability to infect and replicate inside cancer cells, thus killing them, without harming healthy cells. In this work, we present a new non-local multiscale moving boundary model for the spatio-temporal cancer-OV interactions. This model explores an important double feedback loop that links the macro-scale dynamics of cancer-virus interactions and the micro-scale dynamics of proteolytic activity taking place at the tumour interface. The cancer cell-cell and cell-matrix interactions are assumed to be nonlocal, while the cell-virus interactions are assumed local. With the help of this model we investigate computationally various cancer treatment scenarios involving oncolytic viruses (i.e., the eﬀect of injecting the OV inside the tumour, or outside it). Moreover, we investigate the eﬀect of diﬀerent cell-cell and cell-matrix interaction strengths on the success of OV spreading throughout the tumour, and the eﬀect of constant or density-dependent virus diﬀusion coeﬃcients on viral spread.


Introduction
Oncolytic virotherapy (OV-therapy) represents an emerging treatment approach for cancers, which involves a naturally or genetically modified viral strain that selectively infects and damages malignant cancer cells without causing harm to the surrounding healthy cells [21,30,40,41]; see Figure 1.The capacity of viruses to kill cancer cells has been acknowledged for almost a century [21], but the detailed mechanisms behind tumour-virus interactions have started to be understood only over the last few decades.These oncolytic viruses can enter cancer cells either through receptor binding, or through fusion with plasma membrane, and once inside they replicate by taking advantage of pathways disregulated inside cancer cells, or by genetic modification of these cells [21].Due to high viral replication, the infected cancer cells burst open releasing new viral particles in the environment, which can infect new cancer cells.Although the oncolytic viruses can enter also healthy cells, they do not usually replicate inside these cells and thus do not kill them.Recent advances in the understanding of tumour-virus interactions at molecular and cellular levels led in 2015 to the approval by the United States of America Food and Drug Administration (FDA) of the first genetically engineered OV (a Herpex Simplex Virus) that can be used for treating melanoma patients [40].
Even if multiple oncolytic viruses (OVs) are currently under clinical development [47,21], the OVtherapies still have some limitations [21].These limitations are the result of many factors: from premature virus clearance due to circulating antibodies and various immune cells, to physical barriers inside tumours (e.g., interstitial fluid pressure, extracellular matrix (ECM) deposits, tight inter-cellular junctions, or stromal cells that hinder virus movement), and various cellular defences against viruses (e.g., anti-viral innate immunity) [2,20,21,51].To overcome these limitations, numerous experimental and clinical approaches are currently being considered: from genetically manipulating natural OVs Figure 1.Caricature description of the interactions between oncolytic viruses (OVs) and healthy and cancerous cells: the OVs replicate inside cancer cells, leading to their lysis; the OVs do not usually replicate inside healthy cells.
to include additional features for improving safety and efficiency [47], to modifications of the immune responses to enhance virus replication and tumour lysis, or modifications of the physical barriers (e.g., via ECM degradation) to improve virus spread [15,51,57].
Understanding the interactions between cancer cells and ECM during cancer invasion could shed a light on the biological mechanisms that might help OV spread by overcoming the physical barriers inside the tumour microenvironment (as OVs are usually administered during the late stages of cancer, when cancer invasion is an important process).It is worth mentioning here that cancer invasion involves several significant steps: reduction or loss of cell-cell adhesion, enhanced cancer cell adhesion to the extracellular matrix, secretion of matrix degrading enzymes (MDEs) leading to extracellular matrix degradation, which facilitate the movement of cancer cells [13].However, the degradation of the ECM can help also OV spread [51], and therefore a better understanding of the cancer-ECM-OV dynamics might help improve current OV therapies.
Mathematical models can describe, in an abstract way, the interactions between the various components of the tumour microenvironment [17,42].They can also be used to perform in silico (i.e., computational) experiments when in vivo or in vitro experiments are too expensive or too difficult to be performed.There are currently many mathematical studies that focus on cancer growth and its interactions with the oncolytic viruses; see, for example [14,18,19,23,26,35,32,52,55] and the references therein.While most of these models focus on the temporal dynamics of OVs-tumour interactions (mainly because of the availability of temporal data for tumour growth in the presence/absence of various OV therapies), recent advances in tumour imaging have led to the recent development of models investigating the spatio-temporal dynamics of tumour-OV interactions [7,29,33,46,56].The large majority of these temporal and spatio-temporal models focus tumour-OV interactions at single spatial and/or temporal scale.However, since cancer invasion is a complex multiscale phenomenon occurring across different spatial and temporal scales, it is important to capture and model these multiscale processes also for tumour-OV interactions.The last two decades have seen the development of various mathematical models for multiscale cancer dynamics [1,34,39,48,49,50].Nevertheless, there are not many multiscale mathematical models for oncolytic viral therapies and tumour-viral interactions; among the very few multiscale models we mention those in [3,4,37,38].These models focus on local cancer cell dynamics and local interactions between cells, ECM and viruses.In particular, the models assume that cells move randomly and in a directed manner in response to local densities of other cells or ECM.
In this study we develop further the multi-scale approach for modelling cancer-OV interactions introduced in [3,4], and consider non-local cell-cell and cell-matrix interactions for cancer dynamics.We are particularly interested in investigating the role of cell-cell and cell-matrix adhesion on cancer cell movement and degradation of ECM -since this further impacts the spread of OVs.These cell-cell and cell-matrix adhesion forces are modelled with the help of non-local terms, which describe long-distance interactions with other cells or with various components of the ECM, which occur when cells extend long pseudopods to migrate [6,10].These non-local adhesive interactions create an adhesive flux that bias significantly the direction of cell movement [5,22,16].With the help of this new nonlocal multiscale model for cancer-OV interactions we investigate the impact of cancer cell-cell and cell-matrix adhesion to the spread of OV through the solid tumour.
In Section 2 we describe the new multiscale mathematical model.In Section 3 we describe the numerical approach used to discretise the macroscale and microscale equations, while in Section 4 we present the results of the numerical simulations.We conclude in Section 5 with a summary and discussion of the results.

Mathematical Model
Building on the multiscale moving boundary model for cancer invasion introduced in [50] and the multiscale model of cancer response to OV proposed in [3], in this work we address the multiscale nature of tumour-OV interactions through a new non-local modelling approach.This modelling captures the complex tumour-virus-ECM dynamics by exploring the genuine multiscale character of this interaction.Indeed, the tissue-scale (macro-scale) dynamics of cancer cells, ECM and the virus particles, is investigated here in conjunction with the crucial cell-scale (micro-scale) proteolytic activity of the MDEs occurring along the invasive edge of the tumour, while exploring the key double-feedback loop that links these processes across the two spatial scales.This double feedback loop is realised via a top-down link and a bottom-up link, as illustrated schematically in Figure 3.The top-down link leads to a MDE source for the cell-scale proteolytic micro-dynamics occurring at the leading edge of the tumour.In turn, this micro-dynamics induces a bottom-up link that describes how the tumour boundary moves.
2.1.Macro-Scale Dynamics.At the macro-scale, on a maximal tissue cube Y ⊂ R N (N = 2, 3), we explore the coupled dynamics of a population density of uninfected cancer cells, c(x, t), and an ECM density, e(x,t), which interact with a population density of oncolytic virus v(x, t), resulting in an emergent density of infected cancer cells i(x, t).Denoting, for convenience, the total population of infected and uninfected cancer cells by as illustrated also in Figure 2, the evolving tumour domain Ω(t) ⊂ Y represents the support of c(•, t), where the combined distributions of uninfected cancer cells c(x, t) and infected cancer cells i(x, t) exercise their dynamics and interact with the oncolytic virus density v(x, t) within the supporting density of ECM e(x, t).Thus, since the available macro-scale space is a key condition for this complex dynamics to take place, denoting by u the total tumour vector, namely u (x, t) = (c(x, t), e(x, t)) T ,  The bottom-up link focuses on the boundary relocation triggered by the microdynamics.All the notation in this figure were chosen to be consistent with the ones used in [50], and in this context, x * Y represents a "midpoint" boundary location on Y ∩ Ω(t 0 ), as determined in [50], whose spatial relocation is representative for the choreographic movement exercised by the localised tumour interface Y ∩ Ω(t 0 ) over a small time interval [t 0 , t 0 + ∆t].
then the volume fraction of space occupied by the tumour is given as ρ(x, t) ≡ ρ(u(x, t)) := ν e e(x, t) + ν c c(x, t) + i(x, t) . (2.1) Here ν e represents the fraction of physical space occupied by the ECM and ν c is the fraction of physical space occupied by c and i.
Throughout this study we assume that cancer cells proliferate logistically.Moreover, we assume that the uninfected cell population exercises a random movement that is biased by emerging cell adhesion fluxes (enabled by the underlying cell-cell and cell-matrix processes [8,11,28,31,54]).This cancer cell population can be infected by the oncolytic virus at a rate > 0. We can recast these assumptions mathematically through the non-local equation where D c > 0 is a constant diffusion coefficient, µ 1 > 0 is a constant proliferation coefficient, while A(t, x, u(t, •)) is a non-local term accounting for cell adhesion processes for which we adopt a formulation similar to the one proposed in [5,16,22,43].Indeed, at each spatio-temporal point (x, t), the cell adhesion flux captures cumulative influence of the cell-cell and cell-matrix adhesion junctions established between the cells distributed at (x, t) and the cells and the ECM distributed at any spatial position y within a maximal sensing region B(x, R), of radius R > 0, and so this can be mathematically formulated as Here χ Ω(t) (•) is the characteristic function of Ω(t), the term (1 − ρ(u)) + := max{(1 − ρ(u)), 0} enables the avoidance of local overcrowding, and n(y) denotes the unit vector originating from x and pointing to x + y that is given by: with • 2 representing the usual Euclidean norm.Furthermore, the radially symmetric kernel K(•) explores the fact that the adhesion junctions established with both cells and ECM distributed at y ∈ B(x, R) decrease as the distance r := x − y 2 increases, and so this is taken here as Finally, S cc > 0 and S cv > 0 are the strengths of cell-cell and cell-ECM adhesion bonds that are established between the cancer cells distributed at x and the cells and ECM phase distributed at x + y, respectively.While the cell-matrix adhesion S cv is considered to be constant, the cell-cell adhesion strength S cc models the ability of the cell distributed at x to establish cell-cell adhesive junctions with the cells distributed at the other locations y ∈ B(x, R), which depends on the amount of intercellular Ca 2+ ions available within the ECM [24,27].Adopting a similar approach to the one in [44], we assume here that S cc is dependent on the ECM density and it takes the form with S max representing the maximum strength of cell-cell adhesive junctions, as shown in Figure 4.
Finally, as the non-local adhesion flux is effectively formed by summing up the radially distributed adhesive interactions between the cells at x and the cells and ECM distributed at x + y with y ∈ B(0, R), and since the influence of the interaction distance is accounted for through the radially symmetric diffusion kernel K, the term 1 R that appears in the front of the expression is simply here an interaction range normalisation factor.The infected cancer cell population that emerges due to the OV presence exercises also a spatiotemporal dynamics in which the spatial redistribution per unit time is driven by local random motility that is biased by the haptotactic movement against ECM gradients.Simultaneously, the infected cell population expands due to the viral infection of uninfected cancer cells at rate > 0, and reduces due to infected cells dying at rate δ i > 0. Therefore, the dynamics of the infected cancer cell population can be formulated mathematically as where D i > 0 is a constant random motility coefficient, and Further, the ECM is degraded by both uninfected and infected cancer cell populations and is remodelled within the limit of available space.Thus, its dynamics can be captured mathematically by the following equation, where µ 2 > 0 is a constant ECM remodelling rate, while α c and α i are the ECM degradation rates caused by the uninfected and infected cancer cells populations, respectively.Finally, the virus is assumed here not only to diffuse but also to exercise a "haptotactic-like" spatial transport towards higher ECM levels.Further, while exercising the spatial movement, the oncolytic virus replicates itself at a given rate and is eliminated from the system not only by being untaken by the uninfected cells during the infection process, but also via natural decay.Hence, its dynamics can be formulated mathematically as where is the viral decay rate, and D v (c) > 0 is a diffusion coefficient that may be dependent on the presence of the cancer population as only in the presence of cancer the virus exercise spatial transport.Thus, in summary, the non-local model that we obtained for the macroscale dynamics is as follows: ) ) This macro-dynamics will be addressed in the results section both for the case of constant and for the full cancer-dependent viral diffusion coefficient, D v (c), which here is considered to be of the Hill-type saturated form 2.2.Micro-Scale Dynamics and Its Top-Down and Bottom Up Connection to Macro-Scale.
While exercising their macro-scale dynamics, the cells from both the uninfected and infected cancer populations that arrive close to the tumour interface (within the outer proliferating rim of the tumour) secrete matrix degrading enzymes (MDE) such as the matrix metallo-proteinases [25,53], enabling this way the formation of an enzymatic source for a cell-scale (micro-scale) proteolytic micro-dynamics that takes place along the invasive edge of the tumour.In the presence of this source, a cross-interface micro-scale MDE spatial transport occurs within a micro-scale neighbourhood of thickness > 0 of the tumour boundary, which we formally denote by N (∂Ω(t)).The pattern of propagation of the advancing front of MDEs within the peritumoural region N (∂Ω(t)) \ Ω(t) corresponds to the areas of significant ECM degradation, which will ultimately enable us to determine (from the micro-scale MDE dynamics) the law of the macro-scale tumour boundary movement.
Adopting here the modelling approach introduced in [50], the micro-scale neighbourhood N (∂Ω(t)) is described by the union of a covering bundle of − size overlapping micro-domains { Y } Y ∈P (t) , as illustrated in Figures 2-3.These micro-domains Y "sit on the interface" and capture relevant parts of both inside and outside regions of the tumour where the proteolytic activity takes place, enabling us to decompose and explore the proteolytic dynamics on Y as a union of micro-dynamic processes occurring on each Y ∈ P (t).
At any instance in time t 0 > 0, on each micro-domain Y ∈ P (t 0 ), a source of MDEs emerges at every micro-scale location z ∈ Y ∩Ω(t 0 ) as a collective contribution of all the cells (both infected and uninfected) from the tumour outer proliferating rim that arrive during their dynamics within a distance r > 0 from z.Thus, this MDEs source induced by the macro-dynamics at micro-scale realises a key "macro-to-micro" top-down link within the overall tumour dynamics.Mathematically, over a small time interval of length ∆t > 0 and at any micro-scale spatial location z ∈ Y , this MDEs source is therefore given as where λ(•) is the standard Lebesgue measure on R N , B(z, r) := {x ∈ Y : z − x ∞ ≤ r}, and γ c > 0 and γ i > 0 represent the contributions (per unit time) of the uninfected and the infected cancer cells to the formation of the micro-source of MDEs, respectively.Once secreted, these matrix degrading enzymes exercise a diffusion transport process within the entire Y , and so, to describe mathematically their spatio-temporal micro-dynamics, let us denote by m(z, τ ) the MDEs distribution at (z, τ With this notation, under the presence of the micro-source (2.13), the rate of change of the MDEs distribution per unit time is therefore given by ∂m ∂τ where z ∈ Y , τ ∈ [0, ∆t].Furthermore, as we assume no memory of previous molecular proteolytic dynamics within Y and that the molecular dynamics takes place entirely on Y , for (2.14) we assume zero initial conditions as well as no-flux boundary conditions, and so these can be summarised mathematically as m(z, 0) = 0, where n is the outward unit normal.
Finally, during the micro-dynamics (2.14)-(2.15), the MDEs transported across the interface in the peritumoural region Y \ Ω(t 0 ) interact with ECM distribution that they encounter, resulting in degradation of its constituents.The pattern of this ECM degradation caused by the advancing front of MDEs corresponds to the choreographic direction of movement of the tumour interface, and following the derivation in [50], a direction of movement η Y and a displacement magnitude ξ Y is this way established for the movement macro-scale cancer boundary.These key characteristics for the macro-scale boundary movement that are determined at micro-scale realise therefore a "micro-to-macro" bottom-up feedback link that is crucial within the overall tumour dynamics, with direct impact on the changes in tumour morphology.This is represented back at macro-scale through the movement of appropriately defined boundary mid-points x * Y uniquely associated to each Y , which exercise their relocation in the direction η Y and with the displacement magnitude ξ Y to the new spatial positions x * Y , as illustrated in Figure 3, leading to the expansion of the tumour boundary Ω(t 0 ) to an evolved domain Ω(t 0 + ∆t) where the multiscale dynamics is continued.

Computational Approach
As discussed previously, the implementation of the novel multiscale moving boundary model require a number of extensive computational steps, which were built within the multiscale moving boundary computational framework initially introduced by [50] and further expanded in [44].
3.1.Macro-scale Computations on the Expanding Tumour Domain.For the computational approach, we consider that the tumour growth and tumour-OV interaction takes place within a maximal tissue domain Y ⊂ R 2 , which is taken here to be Y = [0, 4] × [0, 4].This macro-scale domain Y is discretised with a uniform grid {(x s , x p )} s,p=1...M of spatial step size ∆x = ∆y := h, see Table 1.Further, considering a uniform discretisation {t l } l=0...k of time step δt > 0 for any given small macroscale time interval [t 0 , t 0 + ∆t], with k := ∆t/δt, let us denote by c l s,p , i l s,p e l s,p , A l s,p , and v l s,p the discretised values of c, i, e, A, and v at any spatio-temporal node ((x s , x p ), t l ), respectively.
As the expansion of the macro-scale tumour interface ∂Ω(t) is induced by the boundary MDEs micro-dynamics, to address the macro-dynamics on a continuously evolving domain Ω(t), our macroscale scheme follows the approach introduced in [44] and summarised in Appendix A. Focusing now on the discretisation of the spatio-temporal operators in (2.10a)-(2.10d)over any given small macro-scale time interval [t 0 , t 0 + ∆t], we solve the macro-dynamics (2.10) via a finite differences Method of Lines type scheme, using the non-local predictor-corrector method for time marching introduced in [44] (and detailed for completeness in Appendix C), while the involved spatial operators are discretised via a combination of local "on grid" and non-local "off grid" approaches that we will detail as follows.
For spatial operators of the same types appearing in (2.10a)-(2.10d)we adopt identical steps in their discretisation.Indeed, except for the adhesion flux term A(t, x, u(t, •)) which will be addressed in Section 3.2, the diffusion and hapto-taxis terms ∇•(D c ∇c), D i ∆i + η i ∇•(i∇e), and η v ∇•(v∇e) (with constant coefficients D c , D i , η i , η v ), are all discretised through a combination of mid-point and central differences schemes and are detailed in Appendix B. Finally, for the discretisation on the entire Y of the oncolytic virus diffusion operator involved in (2.10d), where the cell population dependent diffusion D v (c) defined in (2.12) is involved, we proceed as follows.
Using second-order midpoint rule accompanied by central differences applied at "virtual midpoints" (s ± 1 2 , p) and (s, p ± 1 2 ) corresponding to ((x s , x p ) + (x s±1 , x p ))/2 and ((x s , x p ) + (x s±1 , x p ))/2, for every s, p = 1 . . .M and l = 1 . . .k, we denote D l v,s,p := D v (c l s,p ) and we have that the diffusion coefficient values at these midpoints are given by , and while the central differences for the derivatives of the virus distribution at these midpoints are Therefore, the approximation of the diffusion part in equation (2.10d) becomes: Finally, we note that the derivation of the discretisation of this term for the case when D v (c) is constant is only a particular of the approach exposed above.

Adhesive Flux
Computation.An important part of our mathematical model is the non-local adhesive flux term A(t, x, u(t, •)), which addresses the effects of cell-cell and cell-ECM adhesion of the uninfected cancer cell population and that requires a special numerical treatment.To that end, we follow an approach initially introduced in [43] and further enhanced in [44], which involves a series of off-grid computations on a particular decomposition of the sensing region B(x, R).Specifically, at any given spatio-temporal node ((x s , x p ), t l ), we decompose the sensing region B((x s , x p ), R) as a union of 2 m+(j−1) annulus radial sectors S 1 , ..., S q , (3.4) which are obtained by intersecting each annulus j ∈ {1, ..., s} annuli with a corresponding dyadically increasing number of 2 m+(j−1) uniformly distributed radial sectors of B((x s , x p ), R), while the circle at the centre is considered of a numerically negligible radius, as shown in Figure 5(a).Further, for each ν ∈ {1, . . ., q}, denoting by b Sν the centroid of the sector S ν , from Figure 5(b) we observe that these are distributed so that the position vectors indicating their associated normalised radial directions enables are across the sensing region in a balanced manner, leading to a robust approximation of the adhesion flux.On each annulus sector S ν we calculate the mean values of all the macro-scale densities of both uninfected infected as well as ECM namely: ) respectively.Furthermore, for each centroid b Sν we evaluate the associated radial unit vector denoted by n ν := n(b Sν ) (that points from the centre of the sensing region to b Sp ) and given as Thus, finally, the approximation of the adhesion flux A(t, x, u(t, •)) at the spatio-temporal node ((x s , x p ), t l ) is given by where u l b Sν = [ω l Sν ,c , ω l Sν ,i , ω l Sν ,e ] T , we have that ρ(u l b Sν ) is a mean-value approximation of the tumour volume fraction given in (2.1) associated to the centroid b Sν .Figure 6 shows the simplified case of "duro-taxis" that is mimicked by the cell-adhesion flux approximated in (3.9)where only the cell-ECM adhesion is considered while assuming there the absence of cell-cell-adhesion (i.e., S cc = 0) and ignoring the tumour volume fraction term.We observe that the adhesion vector field follows regions of high ECM density, being directed out of low density into high density.

3.3.
Approximating the Micro-Dynamics and Its Impact on the Macro-Scale Tumour Boundary Movement.Each of the boundary micro-domain Y ∈ P (t 0 ) are considered here to be squares of micro-scale size := 2h (where h is the step size used for the discretisation of the macro-domain), which are centred at each of the interface point (x s , x p ) from the discretisation of the macroscopic tumour boundary ∂Ω(t 0 ).Hence, we have that Y := B • ∞ ((x s , x p ), ) for some boundary point (x s , x p ),, where • ∞ is the usual ∞−norm, and so we simply have that Y is given by the union of the four squares on the macro-grid that have (x s , x p ) as a vertex.Thus, adopting a similar approach as this stage as in [50], this enables us to use bilinear shape functions on each of these four squares (that have (x s , x p ) as a vertex) alongside midpoint averaging to evaluate the MDE source on the entire microdomain Y given in (2.13).Next, we solve MDEs system (2.14) by using backward Euler in time combined with central differences for the spatial discretisation.After finding the MDE distribution m(z, τ ) on Y × [0, ∆t], we then follow precisely the modelling and computational approach and steps introduced by [50] to determine the direction and magnitude for the boundary movement as well as to proceed with the actual expansion of the tumour boundary ∂Ω(t 0 ) to its new spatial configuration ∂Ω(t 0 + ∆t) that emerged over during multiscale tumour evolution over the time interval [t 0 , t 0 + ∆t].

Computational Results
As mentioned in the previous section, the simulations are run on a spatial domain Y = (0, 4) × (0, 4).The initial condition for the uninfected cancer cells population (c) is given by (see also Figure 7 (a)): where ψ γ : R N → R + is the standard mollifier that acts within a radius γ << ∆x 3 from ∂B((2, 2), 0.5−γ) to smooth the characteristic function χ B((2,2),0.5−γ) .This way we assure that the uninfected cancer cell population occupies the region Ω t0 (0) = B((2, 2), 0.5) positioned in the centre of Y .This mollifier ψ γ was introduced in [50]: where ψ is the smooth compact support function given by For the infected cancer cells (i), we assume that there are no such cells at the beginning of the simulations: The initial condition for the ECM (e) is given by the following distribution that was previously introduced in [44] (see also Figure 7(b)): where c(x, 0) is the initial condition of uninfected cancer cells given by equation (4.1), and Finally, we consider two types of initial conditions for OV (v), describing current approaches for OV administration: intravenous and/or intra-tumoural administration.The two cases considered here are: (1) one OV dose is injected outside the tumour (via intravenous administration), and ( 2) two OV doses are injected: one inside the tumour (via intra-tumoural administration) and one outside the tumour (via intravenous administration); see Figure 7(c),(d).
(1) The one-dose initial condition for v(x, t), shown in Figure 7(c), is given as follows where ) Let Γ be the boundary of v(x, 0).In order to smooth the characteristic function θ(v), we take the average of the boundary points of v(x, 0) ∈ Γ, as follows v(x, y) = 1 8 i j v (x + ih, y + jh) − v(x, y), ∀x, y ∈ Γ, where i, j = −1, 0, 1. (4.9) (2) The two-dose initial condition, shown in Figure 7(d), is defined as follows where To smoothen the characteristic function θ(v), we are using the same boundary averaging mentioned previously in eq.(4.9).
The numerical results presented in this Chapter are obtained with the parameter values described in Table 1.Whenever we vary these parameters, we state clearly the new values we use for the simulations.

4.1.
Results of Baseline Parameters.In this section we investigate the effect of haptotactic cell and virus movement (as given by parameters η i and η v ) when one OV dose is administered (Figure 8) or two OV doses are administered (Figure 9): • One OV dose.In Figure 8(a) we show the dynamics of our multi-scale model in the absence of haptotactic terms for the infected cells and virus particles (η i = η v = 0), while in Figure 8(b) we show the dynamics of this model in the presence of such haptotactic terms (where η i = η v = 0.0285, as given in Table 1).We can see that, for the parameter values used in these simulations, The fraction of physical space occupied by the ECM [44] ν c 1 The fraction of physical space occupied by cancer cells [44] γ c 1 MDEs secretion rate by uninfected cancer cell [45] γ i 1.5 MDEs secretion rate by infected cancer cell [45] D m 0.0025 MDE diffusion coefficient [39] h 0.03125 spatial discretisation step size in both x and y directions [44] there is no difference in the spatial distribution of uninfected cancer cells between panels (a) and (b).However, the addition of haptotactic movement impacts the spatial distribution of infected cancer cells and OVs, leading to more localised viral infections.Since the change in ECM density follows the uninfected cancer cells distribution, we choose to show the ECM separately in Figure 8(c), where we illustrate its distribution at two macro-micro stages: I) stage 50; II) stage 75.• Two OV doses.In Figure 9(a) we show the dynamics of our multi-scale model in the absence of any haptotactic movement for the infected cells and virus particles, while in Figure 9(b) we show what happens when we consider haptotactic movement.As for the one-dose case, the haptotactic terms do not affect the spatial distribution of the uninfected cancer cells.The inclusion of the haptotaxis leads to slightly higher densities of infected cells and virus particles,  1.
compared to the case of no haptotaxis.Also the presence of haptotactic terms seems to lead to more localised viral infections.
If we compare the one-dose OV case (Figure 8) with the two-dose OV case (Figure 9) we observe, as expected, that the injection of two OV doses leads to better elimination of uninfected cancer cells through a much better viral replication.We also observe that the injection of two OV doses leads to a faster movement of outside virus particles inside the tumour (e.g., in Figure 8 the outside viruses is still entering the tumour at stage 50, while in Figure 9 the outside viruses have already entered the tumour by stage 50).
4.2.The Impact of Adhesion Strength.In this Section we investigate the impact of cell-cell and cell-matrix adhesion strengths on the spread of OV, and the elimination of tumour cells.
• Constant vs. density-dependent S cc .In equation (2.6), we defined the cell-cell adhesion strength (S cc ) as a function of the ECM density, with a maximum of S max (see also Figure 4).In Figure 10 we investigate the assumption that this cell-cell strength is a constant: S cc = S max .As with the previous figures, in panels (a) we show model dynamics without haptotactic cell movement, while in panels (b) we show model dynamics with haptotactic cell movement.Comparing Figure 8 (where S cc = S max ) with Figure 10 (where S cc = S cc (e)), we notice a difference in cancer invasion, with a slightly larger invasion range when S cc = S max .However, S cc seems to have no impact on the densities of cells and virus particles.Regarding the role of haptotactic terms (i.e., panels (a) vs. (b)), we see again that these terms have no effect on the uninfected tumour cells, but they impact the spatial distributions of infected tumour cells and virus particles, which become more localised.• Different adhesion strengths: S cc > S ce vs. S cc < S ce .In Figure 11, we vary the cell-cell and cell-matrix adhesion strengths from the baseline values listed in Table 1 (i.e., max(S cc ) = S max = 0.1, S ce = 0.5) to max(S cc ) = S max = 0.5, S ce = 0.01, to investigate the case where S cc > S ce .In this Figure 11 we also investigate the effect of increasing the ECM degradation on overall cancer invasion, by changing the baseline rate α c = 0.15 to α c = 2. Comparing Figure 8(a) (where S cc < S ce ) with Figure 11 we observe that high cell-cell adhesion (S cc > S ce ) leads to a very-high density tumour mass which does not spread too much.Moreover, as expected, the increase in the ECM degradation rate α c leads to an almost complete degradation of ECM in the regions with high tumour densities (see the ECM distribution in Figures 11  and 8(c)).

4.3.
Non-Constant OV Diffusion Coefficient.In all previous simulations we assumed that the diffusion coefficients were constant.However, in equation (2.12) we modelled the assumption that the virus diffusion coefficient might depend on the uninfected cancer cell density (D v (c)).In Figure 12 we show numerically the dynamics of our multi-scale model (with haptotactic cell and virus movement) when we assume that D v (c) is given by equation (2.12) with α v = 0.1.In panels (a) we show cancer-OV dynamics when we inject one-dose of OV, while in panels (b) we show model dynamics when we inject two doses.Comparing Figure 12(a) with Figure 8(a), we see that the assumption of a density-dependent diffusion allows the virus to spread faster through the whole tumour.However, in the absence of a higher virus proliferation, the density of infected cancer cells is lower than for the baseline case.A similar conclusion can be obtained when comparing Figure 12(b) with Figure 9(a).

4.4.
Multiple Virus Shots at Different Times.In this final section, we explore the effect of multiple virus injections when these are applied at the same location and with the same dose levels (as the initial viral condition) but at different times during tumour evolution.To capture this multiple temporal viral injections at the same initial location, we amend the virus equation (2.10d) by adding the following viral source term  1.
where, where p represents the number of shots distributed at different times, v(x, 0) are OV initial viral conditions, and Figure 11.Simulations of system (3) investigating the case S max > S ce .We show the densities of cancer cells, OVs (top two rows) and ECM (bottom row) at two macromicro stages (50 and 75), when we ignore the haptotactic movement of infected cells and virus particles (i.e., η i = η v = 0).We consider only one dose of OV as initial condition (see eq. (4.6)).Most of the parameters values are as in Table 1, with the exception of max(S cc ) = S max = 0.5, S ce = 0.01, α c = 2.
represents the time distribution of OV shots which are localised in time in the neighbourhood of the time instances {t 1 , . . ., t p }, each shot having a total time span of 2θ, i.e., occurring on the time interval [t i − θ, t i + θ], for a very small θ > 0, so that consecutive injection intervals would not overlap.Thus, with this amendment, equation (2.10d) becomes We explored this numerically by considering virus injection of the same spatial distribution but administered at two different times t 1 and t 2 during the tumour evolution.These injection times were chosen during the 25th and 50th stages of the tumour evolution as follows:  1.
and so the viral injections scheduling could therefore be captured by g v (x, t) = rt1 (t) + rt2 (t) while taking θ := 3δt in (4.14).Figure 13 shows that when we give multiple virus injections at the same position but at different times, we obtain better tumour control than the case where there is one viral injection at one position in space (shown in Figure 8b).Moreover, as shown in Figure 13, increasing the viral infection rate gives better tumour control results (i.e., the middle of the tumour is eliminated by the viral therapy).as well as for the adhesion fluxes while the central differences at the virtual nodes (s, p ± 1 2 ) and (s ± 1 2 , p) are given by: for c: for i:  Hence, the discretisation of the spatial operator ∇ ∆y .

(B.7)
Denoting now by F l s,p the discretised value of the flux F := D c ∇c − cA(t, x, u(t, •)) at the spatiotemporal node ((x s , x p ), t l ), we observe that the discretisation of ∇ given in equation (2.10a) can therefore be equivalently expressed in a compact form as (∇ • F) where the spatial discretization ∇ • F s,p is given by equation (B.8) but applied to the spatial flux F, and ρ(ū s,p ) is the volume fraction defined in equation (2.1) evaluated for the discrete vector ūs,p = [c s,p , īs,p , ēs,p ], ∀s, p = 1...M .On the time interval [t l , t l+1 ], we first predict c at t l+ 1  2 using an explicit method as follows s,p we calculate the corresponding predicted flux Fl+ 1 2 at t l+ 1 2 .Then, we construct a non-local corrector that involves the average of the flux at the active neighbouring spatial locations {(x s , x p±1 )}, {(x s±1 , x p )}, {(x s±1 , x p−1 )}, {(x s±1 , x p+1 )} ∩ Ω(t 0 ). (C.3) Now, let us denote the set of indices corresponding to these active locations by N , we have that the corrector flux is calculated as and then we compute the flux F l+

Figure 2 .
Figure 2. Schematic diagram showing the cubic tissue region Y .In yellow we illustrate the growing cancer cluster Ω(t 0 ), and the small red cube Y belongs of the family of micro-domains Y covering the tumour boundary ∂Ω(t 0 ).

Figure 3 .
Figure 3. Schematics showing the interactions between the macroscale (tissue-level) dynamics of solid tumours and the proteolytic MDEs microscale dynamics which controls tumour boundary relocation.The top-down link describes the interactions between various cells and their production of MDEs for the proteolytic micro-dynamics.The bottom-up link focuses on the boundary relocation triggered by the microdynamics.All the notation in this figure were chosen to be consistent with the ones used in[50], and in this context, x * Y represents a "midpoint" boundary location on Y ∩ Ω(t 0 ), as determined in[50], whose spatial relocation is representative for the choreographic movement exercised by the localised tumour interface Y ∩ Ω(t 0 ) over a small time interval [t 0 , t 0 + ∆t].

Figure 4 .
Figure 4. Graphical description of the cell-cell adhesion strength as given by equation (2.6).

Figure 5 .
Figure 5. (a) Sensing region B(0, R) approximated by the annulus radial sectors with the barycentre b Sp associated with each sector S p highlighted in red and its centroid points highlighted with a blue dot.(b)The normalised direction of each centroid point corresponding to sub-panel (a); note that these directions are not overlapping, enabling a balanced discretised covering of the radial vector field across the sensing region.

Figure 6 .
Figure 6.Plot of the non-local cell-ECM adhesion flux given by (3.9), while ignoring there the cell-cell adhesion and the effect of overcrowding.

Figure 7 .
Figure 7. Initial conditions used for the numerical simulations: (a) uninfected cancer cells density, as described by equation (4.1);(b) ECM density, as described by equation (4.5); (c) OV density (one initial dose), as described by equation (4.6);(d) OV density (two initial doses), as described by equation (4.10).The white curve indicates the tumour boundary (uninfected and infected cells).

Figure 8 .
Figure 8. Simulations of system (3) using the parameters inTable1 with one-dose OV initial condition (see eq.(4.6)).Here we show the cell and virus distribution at two different micro-macro stages: 50 and 75.(a) No haptotactic behaviours (η i = η v = 0); (b) With haptotactic behaviours (η i = η v = 0.0285).In sub-panels (c) we show ECM density when we consider haptotactic cell movement for I) stage 50, II) stage 75.(We don't show the ECM for the case without haptotactic cell movement, since the ECM follows the uninfected cancer cells distribution, which is the same for both cases.)

Figure 9 .
Figure 9. Simulations of system (3) investigating the injection of two doses of OV (see the initial condition given by eq.(4.10)).We show the densities of uninfected and infected cancer cells, and the OV densities at two macro-micro stages (50 and 75).(a) Without haptotactic movement (η i = η v = 0); (b) With haptotactic movement (η i = η v = 0.0285).The rest of parameter values are as in Table1.

Figure 10 .
Figure 10.Simulations of system (3) investigating the case of constant S cc : S cc = S max .We plot the macroscopic densities of uninfected and infected cancer cells, as well as the densities of OVs at two macro-micro stages (50,75), when we inject two-dose OVs (inside and outside the tumour; see eq. (4.6)).(a) Without haptotactic behaviours (η i = η v = 0); (b) With haptotactic behaviours (η i = η v = 0.0285).The rest of the parameters values are as in Table1.

Figure 12 .
Figure 12.Simulations of system (3) investigating the role of density-dependent viral diffusion D v (c) coefficient (given by eq.(2.12) with α v = 0.1).As for other cases before, we ignore the haptotactic movement of infected cancer cells and virus particles (η i = η v = 0).We show the densities of uninfected and infected cancer cells, and densities of OVs at three macro-micro stages (25, 50, and 75), when we inject: (a) one-dose OV (see initial condition given by eq.(4.6));(b) two-dose OV (see initial condition given by eq.(4.10)).The rest of parameter values are as in Table1.

Figure 13 .
Figure 13.Simulations of system (3) using the parameters in Table 1 with multipledose OV initial conditions (see eq. (4.13)).Here we show the cell and virus distribution at two different micro-macro stages: 50 and 75.(a) multiple OVs injections (b) multiple OVs injections with increased rate of viral infection = 0.158.

Table 1 .
Baseline parameters values used for our multiscale computations.

Table
1  2at t l+12 , in which we use the values of cl+ 1 2 s,p .Next, we use it to initiate the predictor-corrector steps described above on this new time interval t l+12 , t l+1 .Therefore, the predictor step obtained as follows ], ∀s, p = 1...M .Finally, we correct these values at t l+1 with the same non-local trapezoidal-type corrector as described in equation (C.5), here involving the corrector flux calculated as average of the predicted flux values Fl+1 (corresponding to the predicted values cl+1 ) at the active neighbouring locations given in equation (C.3)Thus, this last corrector step gives us ultimately the values that we accept at t l+1 , namely p ], ∀s, p = 1...M .