From 3D models to FDM 3D prints: experimental study of chemical treatment to reduce stairs-stepping of semi-sphere profile

1.   Introduction

The go or grow (GOG) hypothesis or the migration-proliferation dichotomy, proposes that cell proliferation and cell migration are two temporally exclusive events: cells either migrate or proliferate and they can periodically switch between proliferative and migratory states ^[1]. Studies on the GOG hypothesis are conflicting, with some studies supporting and confirming this hypothesis in vitro for different cancer cell types (e.g., glioma, melanoma, or breast cancer) ^{[2,3,4,5,6,7,8,9]} and other studies challenging this hypothesis ^[1,10,11]. For example, Tektonidis et al. ^[9] presented a computational data-driven study of in vitro glioma invasion based on three experimental papers, and concluded that the GOG mechanism combined with self-repulsion and a density-dependent phenotypic switch is mandatory to duplicate the experimental results ^[7,12]. On the other hand, Corcoran et al. ^[1] used time-lapse video-microscopy to monitor directional migration, invasion and mitosis of cancer cells, and concluded that in medulloblastoma cell lines there is no evidence to support the GOG hypothesis. More precisely, their results suggested that migrating and non-migrating cell lines have similar mitotic activities. Similarly, Garay et al. ^[10] tested this hypothesis on 12 mesothelioma, 13 melanoma and 10 lung cancer cell lines using time-lapse video-microscopy, and their results also contradicted the concept of GOG hypothesis.

Over the last decades, mathematical models have been used in addition to experimental studies to shed some light on this go or grow hypothesis. Many of these models are discrete models of cellular automata (CA) type or lattice gas cellular automata (LGCA) type ^{[7,13,14,15,16,17]}. There are also various continuum models described by partial differential equations ^{[11,18,19,20,21,22,23,24,25]}. These continuum models can either model separately the migrating and proliferative cancer cell sub-populations ^{[20,21,22,24]}, or can model them via a single equation for one cancer population that can move for some time instances and proliferate for other time instances ^[11,18]. Those models that consider separate sub-populations of migrating and proliferative cancer cells incorporate also transition (i.e., switching) rates between these two sub-populations. These transition rates can be either constant ^[20,26] or density dependent: they can depend for example on the total density of cells ^[20,22], on the concentration of integrins bound to extracellular matrix (ECM) fibres ^[23], on the density of ECM fibres ^[25], on the level of oxygen ^[19,24].

It should be mentioned here that some GOG models have been reduced to simpler forms; for example, in ^[20] the authors assumed that proliferating cell sub-population could stop proliferating and thus their go-or-growth model was reduced to a go-or-rest model.

The majority of continuum mathematical models in the literature for the GOG hypothesis focus mainly on local interactions between cells. However, it is known that cells can mechanically sense and react to the presence of other cells up to $100\mu m$ away ^[27], and thus more and more mathematical models have been recently developed to consider such non-local cell-cell and cell-ECM interactions ^{[28,29,30,31]}. Nevertheless, these nonlocal models do not usually incorporate the GOG hypothesis.

In this study we plan to investigate (for the first time - to our knowledge) the impact of the GOG hypothesis on oncolytic virotherapies. These oncolytic virotherapies are cancer therapies that use oncolytic viruses (OVs), i.e., viruses that replicate inside and destroy cancer cells. Despite some clinical successes with these oncolytic virotherapies (currently a few such viruses are in the late stages of various clinical trials ^[32]) there are still many open questions related to the interactions between oncolytic viruses and tumour cells ^[33]. And, to our knowledge, it is not clear at this moment how the spread of oncolytic viruses through the solid tumours is affected by the GOG hypothesis. In this study, we investigate this particular aspect using a modelling and computational approach, which allows us to test numerically various hypotheses related to cancer-OV interactions.

To this end, we extend our previous non-local multi-scale mathematical model for cancer-oncolytic viruses (OV) interactions ^[34], by considering also the GOG hypothesis. We consider distinct proliferative and migrating cancer cell subpopulations, and assume that they can become infected with an oncolytic virus (see Figure 1). With the help of this new model, we explore numerically different aspects of the GOG hypothesis, as well as the possibility of having local vs. non-local cell interactions and their impact on cancer invasion in the context of oncolytic virotherapies. We need to emphasise here that with the help of this new model we investigate (for the first time - to our knowledge) the impact of the GOG hypothesis on oncolytic virotherapies.

We describe the new multiscale model in Section 2. The computational approach used to simulate numerically this model is described briefly in Section 3. Then, numerical simulations are presented in Section 4. We conclude in Section 5 with a brief summary and discussion of the results.

2.   Mathematical model

Adopting the multiscale moving boundary modelling approach introduced initially in ^[35], in the following we explore the dynamic interaction between an invading heterotypic tumour and an oncolytic virus. Indeed, considering here the go or grow hypothesis ^[1], the invading tumour is assumed to consist of two subpopulations of cancer cells, namely migrating and proliferative, which exercise their dynamics within the surrounding ECM and that can become infected by an oncolytic virus over a time interval $[0, T]$. For $t\in [0, T]$, denoting by $\Omega(t)$ the spatial support of the progressing tumour that evolves inside a maximal tissue cube $Y$ (i.e., $\Omega(t)\subset Y$), for any $x\in \Omega(t)$, let $c_{m}(x, t)$, $c_{p}(x, t)$, and $e(x, t)$ represent the spatial densities of the migrating cancer cells subpopulation, the proliferating cancer cells subpopulation, and the ECM, respectively. Furthermore, denoting here the oncolytic virus density by $v(x, t)$, $\forall x\in \Omega(t)$, as both the migrating and the proliferating cancer cells can become infected by the oncolytic virus $v(x, t)$, let $i_{m}(x, t)$ and $i_{p}(x, t)$ represent the densities of infected migrating cancer cells and proliferating cancer cells, respectively (see Figure 1).

By incorporating here the GOG hypothesis, we expand and generalise the modelling framework for tumour-OV interaction proposed in ^[29,34]. Indeed, building on the multiscale framework introduced in ^[35], we explore these complex cancer-OV interactions by accounting for the interlinked two-scale dynamics that connects the tissue-scale (macro-scale) tumour macro-dynamics with the cell-scale (micro-scale) proteolytic activity of MDEs that occurs along the invasive edge of the tumour. In the following we describe in detail the macro-dynamics, the micro-dynamics, as well as the double feed-back loop that connects the two scales of activity in our new model.

2.1.   Tumour$-$OV interacting macro-dynamics

For each $t\in [0, t]$, and each $x\in \Omega(t)$, denoting the total cancer cell population by $c_{ \text{total}}(x, t)$

and defining the total tumour vector to be ${\textbf{u}} \left(x, t\right) = (c_p(x, t), \, i_p(x, t), \, c_m(x, t), \, i_m(x, t), \, e(x, t))^T$, the volume fraction of space occupied by the tumour can therefore be expressed mathematically as

where $\nu_e$ represents the fraction of physical space occupied by the ECM and $\nu_c$ is the fraction of physical space occupied collectively by all cancer subpopulations.

For the tumour dynamics, we assume that the motility of each of the cancer cells subpopulations is due to a combination of random movement (described by linear diffusion term) and a directed migration due to cell-cell and cell-ECM adhesion. The spatial fluxes triggered by cell adhesion that cause the directed cells migration are considered here both from a local and non-local perspective ^{[36,37,38,39,40,41]}, and in following we will detail their mathematical formulation. For convenience, for each cancer cell subpopulation $c \in \lbrace c_p, \: c_m, \: i_p, \:i_m \rbrace$, we consider a global notation $\varphi_{c}({\textbf{u}})$ describing the effect of the cell adhesion processes either locally, through adhesive interactions between cancer cells and ECM (whereby the tumour cells exercise haptotactic movement ^[42] towards higher levels of ECM), or non-locally, where both cell-cell and cell-ECM adhesive interactions are accounted for within an appropriate cell sensing region. Thus, $\varphi_{c}({\textbf{u}})$ is mathematically formalised as

where, for any given subpopulation $c \in \lbrace c_p, \: c_m, \: i_p, \:i_m \rbrace$, we have that $\eta_{c} > 0$ is a constant haptotactic rate associated to $c$, while ${\mathcal{A}}_{c}(x, t, {\textbf{u}}(\cdot, t))$ is a non-local spatial flux term that is detailed as follows. Indeed, following a similar approach as in ^[28,43,44], $\forall c \in \lbrace c_p, \: c_m, \: i_p, \:i_m \rbrace$, at each spatio-temporal point $(x, t)$ the cell adhesion flux ${\mathcal{A}}_{c}(x, t, {\textbf{u}}(\cdot, t))$ cumulates the strengths of the cell-cell and cell-matrix adhesion junctions that cells from cancer subpopulation $c$ that distributed at $(x, t)$ establish with the other cell subpopulations and the ECM distributed within an appropriate maximal sensing region ${\bf{B}}(x, R)$ of radius $R > 0$. This is formulated mathematically as

here, $\chi_{_{\Omega(t)}}(\cdot)$ is the characteristic function of $\Omega(t)$, while the term $(1-\rho({\textbf{u}}))^{+} : = \text{max} \lbrace (1-\rho({\textbf{u}})), 0\rbrace$ enables the avoidance of local overcrowding. Further, for any $y\in {\bf{B}}(0, R)$, ${\textbf{n}}(y)$ denotes the unit radial vector originating from $x$ and pointing to $x + y\in {\bf{B}}(x, R)$, which is given by

with $\parallel {\cdot} \parallel_{_{2}}$ being the usual Euclidean norm. Moreover, $\textbf{K}(\cdot):[0, R]\to [0, 1]$ is a radially symmetric kernel that explores the dependance of the strengths of the established cell adhesion junctions on the radial distance from the centre of the sensing region $x$ to $\zeta\in {\bf{B}}(x, R)$. Since these adhesion junction strengths are assumed to decrease as the distance $r: = \parallel {x-\zeta} \parallel_{_{2}}$ increases, $\textbf{K}$ therefore is taken here of the form

Furthermore, since in this study we focus only on self-adhesion and we do not consider cross-adhesion bonds between four cancer cells subpopulation, in Eq (2.4) we have that $S_{cc}$ and $S_{ce}$ for any given cancer cell subpopulation $c \in \{ c_p, \: c_m, \: i_p, \:i_m \}$ represent the adhesive interaction strengths for the self$-$cell$-$cell adhesion and cell-ECM, respectively. While the cell-ECM adhesion strength $S_{ce}$ is considered to be a positive constant, the cell-cell adhesion strength $S_{cc}$ explores here the fact that the ability of the cell distributed at $x$ to establish cell-cell adhesive junctions with the cells distributed at the other locations $y\in {\bf{B}}(x, R)$ depends on the amount of intercellular $Ca^{2+}$ ions available within the ECM ^[45,46]. As a consequence, adopting a similar approach to the one in ^[47], we assume here that $S_{cc}$ is dependent on the ECM density and it takes the form

with $S^{\max}_{cc}$ representing the maximum strength of cell-cell adhesive junctions established by the cancer cells subpopulation $c \in \{c_p, \: c_m, \: i_p, \:i_m \}$ . Therefore, the adhesive strengths for cell-cell and cell-ECM adhesion for all four cancer subpopulation can be compactly expressed via the diagonal matrices

respectively.

Finally, in the context of the GOG hypothesis, another important aspect that occurs during the tumour dynamics is the transitions between from the proliferative cancer subpopulation and the migrating one. Adopting a similar form to the one proposed in ^[15,16,26], the transition from proliferative to migrating cancer cells is captured here through the switching term $\lambda_{c_p, c_m}$ that is given by

where $\omega_1$ is the rate of switching from proliferative state $c_p$ to migration state $c_m$, and $\omega_2$ is the rate of switching from migration to proliferative state. On the other hand, the transition from migrating to proliferative cancer cells is expressed through the switching term $\lambda_{c_m, c_p}$ defined as

Thus, for each of the cancer cells subpopulation $c\in\{ c_p, \: c_m, \: i_p, \:i_m \}$, the spatial transport is a combination of random movement (expressed through diffusion) and directed movement due adhesion (explored either locally or non-locally, and represented compactly through $\varphi_{c}({\textbf{u}})$). Furthermore, for the particular case of the migrating and proliferative cells subpopulation $c_{p}$, and $c_{m}$, besides the spatial transport and in addition to their own proliferation (considered here of logistic type ^[48,49]), their dynamics is also affected by the cell population "exchanges" due to proliferative-migrating transitions (i.e., transitions between proliferative and migrating subpopulations) as well as by the presence of the oncolytic virus that is able to infect cells from both populations. Finally, for their part, the infected cancer cells subpopulations, while exercising a spatial transport of the type described above, they contribute to virus replication and die. Therefore, the dynamics for each cancer cell population can be expressed mathematically as follows.

First, the governing equation for the uninfected proliferative cancer cell subpopulation is given by

where $D_{c_p} > 0$ is a constant diffusion coefficient, the term $\varphi_{c_p}({\textbf{u}})$ represents the directed movement triggered by cell-adhesion processes that corresponds to $c_{p}$ and is described in Eq (2.3) for $c = c_p$. Further, $\mu_{p} > 0$ is an intrinsic constant proliferation rate, $\varrho_p > 0$ is the rate at which the oncolytic virus infects the proliferative cancer cell population, and $\lambda_{c_p, c_m}$ is the switching term given in (2.9), representing the process through which migrating cancer cells transition towards proliferative state during the tumour dynamics.

The infected proliferative cancer cell population, $i_p(t, x)$, which emerges within this dynamics due to the OV infection of $c_p$, also exercises a spatio-temporal dynamics that is governed by the following equation

where $D_{i_p} > 0$ is a constant random motility coefficient, and $\varphi_{i_p}({\textbf{u}})$ is the spatial influence of the cell-adhesion processes that is described in Eq (2.3) and corresponds to $i_{p}$. The cancer cell population increases with at rate $\varrho_p$ due to the new infections of the proliferative cancer cells, and decreases at rate $\delta_{i_p} > 0$ due to infected cell death.

Further, since for the migrating cancer cell population we always take into account not only cell-ECM adhesion but also cell-cell self-adhesion, the directed cell migration term $\varphi_{c_m}({\textbf{u}})$ that is defined in Eq (2.3) and corresponds to $c_{m}$ is in this case constantly of the non-local form

where the spatial flux ${\mathcal{A}}_{c_m}(\cdot, \!\cdot, \! {\textbf{u}}(\cdot, \!\cdot)\!)$ is the one defined in Eq (2.4) for $c = c_{m}$. As a consequence, the governing equation for the uninfected migrating cell population is

where $D_{c_m} > 0$ is a constant diffusion coefficient, $\mu_{m} > 0$ is a constant proliferation coefficient, $\varrho_m > 0$ is a constant rate at which the uninfected migrating population diminishes due infection by the oncolytic virus $v$. Further, $\lambda_{c_m, c_p}$ is the switching term given in Eq (2.10) that represents the net transition from the proliferative into the uninfected migrating state that occurs per unit time during the tumour dynamics.

The fourth tumour cell population is the infected migrating cancer cell subpopulation $i_m(t, x)$ that emerges within this dynamics due to infections by the OV, and its spatio-temporal dynamics is governed by the following equation

where $D_{i_m} > 0$ is a constant random motility coefficient, and ${\varphi_{i_m}({\textbf{u}})}$ represents the directed migration induced by the cell-adhesion processes that corresponds to $i_{m}$ and is described in Eq (2.3) for $c = i_m$ Further, the infected migrating population expand at a rate $\varrho_m$ due to new infections occurring among the uninfected migrating cells, and they also die at rate $\delta_{i_m} > 0$.

At the same time, the ECM is degraded by both uninfected and infected cancer cell populations and is remodelled within the limit of available space. Thus, its governing dynamics is given mathematically by

where $\alpha_{c_p} > 0$, $\alpha_{i_p} > 0$, $\alpha_{c_m} > 0$, and $\alpha_{i_m} > 0$ are the ECM degradation rates caused by cancer cells subpopulation $c_p$, $i_p$, $c_m$, and $i_m$, respectively. Further, $\mu_{2} > 0$ is a constant ECM remodelling rate.

Concerning the oncolytic virus spatio-temporal dynamics, we adopt here a similar reasoning as in ^[34], and we assume that the OV motion is described by a random movement that is biased by a "haptotactic-like" spatial transport towards higher ECM levels. Thus, the dynamics of the oncolytic virus that we consider here is governed by

where $D_{v} > 0$ is a constant random motility coefficient, $\eta_{v} > 0$ is a constant haptotactic coefficient, $b_m, b_p > 0$ are a viral replication rates within infected proliferating and infected migrating cancer cells, respectively, and $\delta_{v} > 0$ is the viral decay rate.

Finally, the coupled interacting tumour $-$ OV macro-dynamics is governed by Eqs (2.11)-(2.16) in the presence of initial conditions

while assuming zero-flux boundary conditions at the moving tumour interface $\partial \Omega(t)$.

2.2.   Micro-scale dynamics

During their macro-scale dynamics, the four cancer cells subpopulations that get near the tumour interface (i.e., within the outer proliferating rim of the tumour) are able to secrete matrix degrading enzymes (such as the matrix metalloproteinases ^[50,51]), providing this way a source of MDEs for a cell-scale (micro-scale) proteolytic micro-dynamics that takes place along the invasive edge of the tumour. Indeed, in the presence of this source of MDEs (induced by the tumour macro-dynamics), a cross-interface micro-scale MDEs spatial transport occurs within a micro-scale neighbourhood of the tumour boundary of an appropriate cell-scale thickness $\epsilon > 0$, denoted here simply by ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))$. The areas of significant ECM degradation caused by the pattern of propagation of the advancing front of MDEs within the peritumoural region ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))\setminus \Omega(t)$ will ultimately be explored by the cancer cells that will progress in those regions ^[51], and so precisely these boundary regions (affected by significant ECM degradation) will shape the pattern of tumour progression. Thus, following the modelling approach introduced in ^[35] we depict these regions of significant ECM degradations by exploring the MDEs micro-dynamic processes within ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))$, which enables us ultimately to determine the law of the macro-scale tumour boundary movement.

To formalise these laws of macro-scale boundary movement induced by the boundary MDEs micro-dynamics, we adopt here the approach introduced in ^[35]. Therefore, the micro-scale neighbourhood ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))$ is given here as a union of a covering bundle of $\epsilon-size$ overlapping micro-domains $\{\epsilon Y\}_{\epsilon Y\in {\mathcal{P}}_{\epsilon}(t)}$, namely,

This enables us to decompose the MDEs micro-dynamics on $N_{\epsilon}(\partial \Omega(t))$ by exploring this as a union of micro-dynamic processes occurring on each $\epsilon Y\in {\mathcal{P}}_{\epsilon}(t)$. At any instance in time $t_{0} > 0$, on each micro-domain $\epsilon Y\in {\mathcal{P}}_{\epsilon}(t_{0})$, a source of MDEs appears at every micro-scale location $z\in \epsilon Y \cap \Omega (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 $\rho > 0$ from $z$. Thus, over any small time interval of length $\Delta t > 0$, $[t_{0}, t_{0}+\Delta t]$ and at any micro-scale spatial location $z\in \epsilon Y$, this MDEs source is therefore given as

where $\lambda(\cdot)$ is the standard Lebesgue measure on $\mathbb{R}^N$, the ball ${\bf{B}}(z, r): = \lbrace x \in Y: \left\lVert{z-x}\right\rVert_\infty \leq \rho\rbrace$ is the maximal outer proliferating rim region from where cells that get to contribute to the formation of MDEs source at $(z, \tau)\in \epsilon Y \times [t_{0}, t_{0}+\Delta t]$, and $\gamma_{c_p}$, $\gamma_{i_p}$, $\gamma_{c_m}$, $\gamma_{i_m}$ are all positive constants representing the contributions of the cancer subpopulations of uninfected proliferative cells, infected cancer cells, uninfected migrating cells, and infected migrating cells, respectively.

In the presence of the micro-scale source of MDEs induced from the macro-dynamics on each micro-domain $\epsilon Y$, these matrix degrading enzymes exhibits a diffusion transport process within the entire $\epsilon Y$. Thus, denoting the MDEs distribution at $(z, \tau)\in \epsilon Y \times [0, \Delta t]$ by $m(z, \tau)$, the MDEs micro-dynamics on each $\epsilon Y$ is given by

where $z\in \epsilon Y$, $\tau \in [0, \Delta t]$. Furthermore, since we assume no pre-existing MDEs within $\epsilon Y$ prior to the initiation of the proteolytic micro-dynamics, the MDEs micro-dynamics Eq (2.25g) takes place in the presence of zero initial conditions. Furthermore, we assume the presence of zero-flux boundary condition, namely

where ${\textbf{n}}_{\epsilon Y}$ is the outward unit normal on the frontier of the micro-domain $\partial\epsilon Y$.

2.3.   Bottom-up feedback: the macro-scale tumour boundary movement induced by the boundary MDEs micro-dynamics

During the micro-dynamics Eqs (2.19) and (2.20), the MDEs transported across the interface in the peritumoural region $\epsilon Y\setminus \Omega(t_{0})$ interact with ECM distribution that they encounter, resulting in degradation of ECM constituents. The advancement of MDEs within the peritumoural region $\epsilon Y\setminus \Omega(t)$ lead to a degradation of the ECM in that cell-scale region, and determines the way the macroscopic tumour boundary evolves, leading to the establishment of a boundary movement law. Indeed, following the derivation in ^[35], the MDEs micro-dynamics on each micro-domain $\epsilon Y$ enables us to derive the movement characteristics for the relocation of the macro-scale tumour boundary $\partial \Omega(t)\cap \epsilon Y$, expressing these through the derivation of a direction of movement $\eta_{\epsilon_Y}$ and a displacement magnitude $\xi_{\epsilon Y}$ in that direction for the advancement of $\partial \Omega(t)\cap \epsilon Y$ within the peritumoural region ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))\setminus \Omega(t)$. To simplify the representation, the choreographic movement exercised by the $\partial \Omega(t)\cap \epsilon Y$ over a given time span $[t_{0}, t_{0}+\Delta t]$ is represented back at macro-scale through the movement of the associated boundary midpoint $x^{*}_{_{\epsilon Y}}$ of $\epsilon Y$ (defined topologically with full details in ^[35], and which can be regarded as "the center of $\partial \Omega(t)\cap \epsilon Y$"), as illustrated in Figure 2. For completeness, we briefly outline below the main steps involved in deriving the boundary relocation characteristics that were introduced in ^[35].

On the cell-scale neighbouring bundle ${\mathcal{N}}_{\epsilon}(\partial \Omega(t_{0}))$ of the tumour interface, for each of the boundary micro-domains $\epsilon Y\in {\mathcal{P}}(t_{0})$ at a given a time instance $t_{0} > 0$, we use the regularity property of Lebesgue measure ^[52] to depict the first dyadic decomposition $\{{\mathcal{D}}_{k}\}_{k\in {\mathcal{I}}_{_{\epsilon Y}}}$ of $\epsilon Y$ that has the property that the union of those dyadic cubes ${\mathcal{D}}_{k}$ included in the complement of $\Omega(t_{0})$ approximate to a given global micro-scale accuracy $\delta_{_{\Omega(\cdot)}} > 0$. This is schematically illustrated by the small green squares in Figure 2 that are situated outside the black tumour boundary $\partial \Omega(t_{0})\cap \epsilon Y$. Further, denoting by $y_{k}$ the barycenter of ${\mathcal{D}}_{k}$, we sub-select a sub-family of dyadic cubes $\{{\mathcal{D}}_{k}\}_{k\in {\mathcal{I}}^{*}_{_{\epsilon Y}}}\subset \{{\mathcal{D}}_{k}\}_{k\in {\mathcal{I}}_{_{\epsilon Y}}}$ that consists only of those dyadic cubes that are situated furthest away from the boundary midpoint $x^{*}_{_{\epsilon Y}}$ (corresponding to $\epsilon Y$) with the property that they carry an amount of MDEs above the mean of MDEs transported within the entire preritumoural region $\epsilon Y\setminus \Omega(t_{0})$, hence covering precisely the region of significant ECM degradation caused by MDEs within $\epsilon Y\setminus \Omega(t_{0})$, as illustrated in Figure 2. Thus, by cumulating the contribution to the significant ECM degradation within $\epsilon Y\setminus \Omega(t_{0})$ of all the dyadic cubes $\{{\mathcal{D}}_{k}\}_{k\in {\mathcal{I}}^{*}_{_{\epsilon Y}}}$ while accounting on both their relative spatial location with respect to $x^{*}_{_{\epsilon Y}}$ and the amount of MDEs that they get to carry at time $\tau_{f}: = t_{0}+\Delta t$, we obtain the direction of choreographic boundary relocation (exercised by the $\partial \Omega(t)\cap \epsilon Y$) due to micro-scale MDEs degradation. Therefore, this boundary movement direction is given by the positive direction of the emerging line defined by the position vectors involved $\big\{\overrightarrow{y_{k}x^{*}_{_{\epsilon Y}}}\big\}_{k\in {\mathcal{I}}^{*}_{_{\epsilon Y}}}$ magnified accordingly by the MDEs mass that each dyadic cube in $\{D_{k}\}_{k\in {\mathcal{I}}^{*}_{_{\epsilon Y}}}$, carries, namely

Furthermore, the magnitude of the actual boundary movement in direction $\eta_{\epsilon_Y}$ is appropriately given as a weighted sum of the Euclidean magnitudes of the position vectors $\{\parallel {y_{k}x^{*}_{_{\epsilon Y}}} \parallel_{2}\}_{k\in {\mathcal{I}}^{*}_{_{\epsilon Y}}}$, with the weights accounting on the relative contribution brought to the ECM degradation of each of the corresponding dyadic cubes. Thus the movement magnitude in direction $\eta_{\epsilon_Y}$ is given by

Therefore, as the tumour boundary relocation induced by the micro-dynamics on each $\epsilon Y$ is represented at macro-scale through the movement of the boundary midpoint $x^{*}_{_{\epsilon Y}}$, in the context that enough but not complete ECM degradation occurs within $\epsilon Y\setminus \Omega(t_{0})$ (tissue condition that is detailed and explored in full in ^[35]), we have that $x^{*}_{_{\epsilon Y}}$ exercises a relocation to a new position $\widetilde{x^{*}_{_{\epsilon Y}}}$ that is given by

where

and is illustrated through the dark blue arrow in Figure 2. Thus, a law for macro-scale tumour boundary movement is this way induced by the MDEs micro-dynamics, enabling us to capture the evolution of the tumour boundary over the time interval $[t_{0}, t_{0}+\Delta t]$ from its state $\Omega(t_{0})$ at $t_{0}$ to a new spatial configuration at $\Omega(t_{0}+\Delta t)$ at $t_{0}+\Delta t$. This relocated domain $\Omega(t_{0}+\Delta t)$ allows the initiation of the dynamics on the next time interval $[t_{0}, t_{0}+\Delta t]$ where the tumour-oncolytic virus interaction continues its proceedings.

2.4.   Brief summary of the multiscale model

In summary, the multiscale moving boundary model that we obtained for the tumour$-$OV interaction (schematically illustrated in Figure 2) is structured as follows,

The macro-dynamics and micro-dynamics are connected through a double feedback loop enabled by:

● a top-down link by which the macro-dynamics induces the source for the micro-dynamics given in Eq (2.18).

● a bottom-up link by which the MDEs micro-dynamics induces and determines the law for the macro-scale tumour boundary movement.

3.   Computational approach

The numerical approach and computational implementation of the novel multiscale moving boundary model require a number of steps that build on the multiscale moving boundary computational framework initially introduced by ^[35] and further expanded in ^[34,47].

3.1.   Brief overview of the multiscale computational approach

Macro-scale computations. The maximal macro-scale tissue domain $Y \subset \mathbb{R}^2$, where the tumour-OV interacting macro-dynamics Eq 2.25(a)-(f) takes place, is considered here to be $Y : = [0, 4]\times [0, 4]$ and is discretised uniformly using a spatial step size $\Delta x = \Delta y: = h$, with $h > 0$. Let's denote by $Y^{d}$ the discretised $Y$, i.e., $Y^{d}: = \{(x^{1}_{i}, x^{2}_{j})\}_{i, j = 1\dots N}$, with $N = [4/h]+1$. Further, as the macro-dynamics Eq 2.25(a)-(f) is addressed only on the expanding tumour domain $\Omega(t)\subset Y$, for convenience, for any $t > 0$, we denote the discretised tumour domain by $\Omega^{d}(t)$ (i.e., $\Omega^{d}(t) = Y^{d}\cap \Omega(t)$) and the discretised tumour boundary by $\partial \Omega^{d}(t)$ (i.e., the frontier of $\Omega^{d}(t)$ is $\partial \Omega^{d}(t)$). To carry out the computations exclusively on the expanding tumour, the numerical scheme that we developed here involves a method of lines-type approach combined with a non-local predictor corrector time-marching method introduced in ^[47] (and, for completeness, summarised also in Appendices A and B). Finally, as the tumour progresses, $\Omega^{d}(t)$ is appropriately expanded by activating and including within tumour domain the new points invaded by cancer within $Y^{d}$.

Approximating the micro-dynamics and its top-down and bottom-up links with the tumour$-$OV macro-dynamics. At any instance of time $t_{0}$, we consider that the cell-scale covering bundle $\{\epsilon Y\}_{\epsilon Y\in {\mathcal{P}}_{\epsilon}(t_{0})}$ of the discretised tumour interface $\partial \Omega^{d}(t_{0})$ consists of overlapping squares $\epsilon Y$ of micro-scale size $\epsilon: = 2h$, which are centred at each of the tumour interface spatial node $(x^{1}_{s}, x^{2}_{p})\in \partial \Omega^{d}(t_{0})$, i.e.,

where $\parallel {\cdot} \parallel_{\infty}$ is the usual $\infty-$norm, and ${\bf{B}}_{_{\parallel {\cdot} \parallel_{_{\infty}}}}\!\!(\!(x^{1}_{s}, x^{2}_{p}), \epsilon/2): = \{(z^{1}, z^{2})\in {\mathbb{R}}^{2}\, |\, \parallel {(x^{1}_{s}, x^{2}_{p})-(z^{1}, z^{2})} \parallel_{\infty}\leq \epsilon/2\}$ is the closed ball of radius $\epsilon/2$.} By adopting a similar approach to the one introduced in ^[35], we use using bilinear shape functions to calculate the MDE source given by Eq (2.18) on each micro-domain $\epsilon Y$. To solve MDEs micro-dynamics Eq (2.25g), we use backward Euler in time combined with central differences for the spatial discretisation. After finding the MDE distribution $m(z, \tau)$, with $(z, \tau)\in\epsilon Y\times [0, \Delta t]$, we follow the modelling and computational approach introduced in ^[35] to determine the direction $\eta_{\epsilon Y}$ and displacement magnitude $\xi_{\epsilon Y}$ for the movement of the tumour boundary $\partial \Omega(t_{0})\cap \epsilon Y$ that is captured by each micro-domain $\epsilon Y: = {\bf{B}}_{_{\parallel {\cdot} \parallel_{_{\infty}}}}~~~\!\!(\!(x^{1}_{s}, x^{2}_{p}), \epsilon)$ and is represented trough the movement of its midpoint $(x^{1}_{s}, x^{2}_{p})\in \partial \Omega^{d}(t_{0})$. Finally, we use these movement characteristics induced from the micro-dynamics (i.e., $\eta_{\epsilon Y}$ and $\xi_{\epsilon Y}$, $\forall\, \, \epsilon Y\in {\mathcal{P}}_{\epsilon}(t_{0})$) to proceed with the corresponding global relocation of the macro-scale tumour boundary $\partial \Omega^{d}(t_{0})$ to its new spatial configuration $\partial \Omega^{d}(t_{0}+\Delta t)$, which emerges due to the multiscale tumour evolution over the time interval $[t_{0}, t_{0}+\Delta t]$.

3.2.   Local vs non-local directed migration due to cell adhesion

In our numerical experiments, we explore the multiscale model dynamics on three distinct local and non-local scenarios that we consider within the macro-dynamics Eq 2.25(a)-(f) for the directed migration due to cell adhesion for cancer cell subpopulations $c_p$, $i_p$, and $i_m$. Specifically, we consider the following cases:

1. The cell-adhesion interactions for both the uninfected proliferative subpopulation $c_p$ and for the infected subpopulations $i_p$, and $i_m$ are considered to be local of haptotactic type, i.e., in the coupled macro-dynamics in Eq 2.25(a)-(f) we have $\varphi_{c_p}({\textbf{u}}) = \eta_{c_{p}} \!\nabla\!\! \cdot \!\big(c_{p}\nabla e\big)$, $\varphi_{i_p}({\textbf{u}}) = \eta_{i_{p}} \!\nabla\!\! \cdot \!\big(i_{p}\nabla e\big)$, and $\varphi_{i_m}({\textbf{u}}) = \eta_{i_{m}} \!\nabla\!\! \cdot \!\big(i_{m}\nabla e\big)$. Thus, the macro-dynamics Eq 2.25(a)-(f) is in this case of the form:

with results for this case shown in Figure 4.

2. The cell-adhesion interactions for the uninfected proliferative subpopulation $c_p$ are considered now to be non-local, while the infected subpopulations $i_p$, and $i_m$ are still considered to be local of haptotactic type. Hence, in the coupled macro-dynamics in Eq 2.25(a)-(f) we have $\varphi_{c_p}({\textbf{u}}) = \nabla \!\! \cdot\! \big(c_{p}{\mathcal{A}}_{c_{p}}(\cdot, \!\cdot, \! {\textbf{u}}(\cdot, \!\cdot)\!)\!\big)$, while $\varphi_{i_p}({\textbf{u}}) = \eta_{i_{p}} \!\nabla\!\! \cdot \!\big(i_{p}\nabla e\big)$, and $\varphi_{i_m}({\textbf{u}}) = \eta_{i_{m}} \!\nabla\!\! \cdot \!\big(i_{m}\nabla e\big)$. Thus, the macro-dynamics Eq 2.25(a)-(f) is in this case of the form:

with results for this case shown in Figure 5(a).

3. Finally, all the cell-adhesion interactions for both the uninfected proliferative subpopulation $c_p$ and for the infected subpopulations $i_p$, and $i_m$ are considered to be non-local, i.e., in the coupled macro-dynamics in Eq 2.25(a)-(f) we have $\varphi_{c_p}({\textbf{u}}) = \nabla \!\! \cdot\! \big(c_{p}{\mathcal{A}}_{c_{p}}(\cdot, \!\cdot, \! {\textbf{u}}(\cdot, \!\cdot)\!)\!\big)$, $\varphi_{i_p}({\textbf{u}}) = \nabla \!\! \cdot\! \big(i_{p}{\mathcal{A}}_{i_{p}}(\cdot, \!\cdot, \! {\textbf{u}}(\cdot, \!\cdot)\!)\!\big)$, and $\varphi_{i_m}({\textbf{u}}) = \nabla \!\! \cdot\! \big(i_{m}{\mathcal{A}}_{i_{m}}(\cdot, \!\cdot, \! {\textbf{u}}(\cdot, \!\cdot)\!)\!\big)$. Thus, the macro-dynamics Eq 2.25(a)-(f) is in this case of the form:

with results for this case shown in Figure 5(b).

3.3.   Initial conditions

The initial conditions for the uninfected proliferative cancer cell population, $c_{p}(x, 0)$ is chosen to describe a small localised pre-existing tumour aggregation. This is given by the following equations:

whose plot is shown in Figure 3(a). Here $\psi_{\gamma}:\mathbb{R}^N \rightarrow \mathbb{R}_{+}$ is the usual standard mollifier of radius $\gamma < < \frac{\Delta x}{3}$ given by

where $\psi$ is the smooth compact support function given by

Moreover, we assume that the tumour is detected early enough so that migration is not initiated at the start of these simulations, and thus

Also, since there is no infection at this stage, we assume that both infected proliferative ($i_{p}(x, 0)$) and migrating ($i_{m}(x, 0)$) cancer cells are zero:

Furthermore, the initial condition for the ECM density, $e(x, 0)$, is represented by an arbitrarily chosen heterogeneous pattern described by the following equations (as in ^[47])

and is shown in Figure 3(b). Here, we have

While other types of heterogeneous ECM patterns could be considered (see ^[53]), here we focus our attention to explore cancer-viral dynamics on this particular ECM pattern.

Finally, the initial conditions for the OV population, $v(x, 0)$ is chosen to describe one single injection in the middle of the tumour aggregation, as in ^[34]. This is described by the equation

where

In computations, the initial condition is smoothed out on the frontier of the viral density support $\Gamma_{v}: = \partial \{x\in Y\, |\, v^{0}(x) > 0\}$ via the averaging

4.   Results

The numerical results presented in this Section are obtained with the parameter values described in Table 1 which, for convenience we call them 'baseline parameters'. Whenever we vary these parameters, we state clearly the new values we use for those simulations. Note that these baseline parameters are based on other papers or on our own estimates. For instance, using the GOG hypothesis, we estimated that $D_{c_p}$ is likely much smaller than $D_{c_m}$. Since we could not find an exact value for $D_{c_p}$ ($D_{c_m}$ was assumed to be $0.00035$, as in ^[43]), we arbitrarily estimated $D_{c_p} = 10^{-5}$.

We start in Section 4.1 by investigating numerically the impact of local vs. nonlocal approaches used to describe the cell-cell and cell-matrix adhesion flux. Then, in Section 4.2, we investigate the impact of varying the adhesion strength in the non-local cell flux. Following that we focus on the system of Eq (3.1) without haptotaxis for $c_p$ $(i.e., \eta_{c_p} = 0)$, to investigate the impact of varying different parameters: in Section 4.3 we vary the impact of transition rate between migrating and proliferative cells, in Section 4.4 we vary the impact of OVs infection rate, and in Section 4.5 we vary the impact of OVs replication rate.

4.1.   The importance of cell-cell and cell-matrix adhesion: Local vs non-local approaches

First, we focus on model (3.1), described in detail in Section 3.2. In Figure 4(a) we show the dynamics of our multiscale model in the absence of haptotactic terms for the proliferative uninfected cells $(\eta_{c_p} = 0)$, while in Figure 4(b) we show the dynamics of this model in the presence of such haptotactic terms (where $\eta_{c_p} = 0.00285$, as given in Table 1). We can see that, for the parameter values used in these simulations, there is no difference in the spatial distribution of migrating uninfected or infected cancer cells between panels (a) and (b). However, the addition of haptotactic movement impacts the spatial distribution of proliferative uninfected cancer cells, leading to a more localised cancer cells distribution.

Next, we investigate numerically the dynamics of models (3.2) and (3.3). In Figure 5(a) we show the dynamics of our multi-scale model (3.2) that has two non-local fluxes $c_m$ and $c_p$, while in Figure 5(b) we show what happens when we consider non-local fluxes for all cancer cells subpopulations (i.e., model (3.3)). Clearly there is a difference between these results and the previous one. The cancer cells for the case with two non-local fluxes (i.e., for proliferative $c_p$ and migrating $c_m$ cells) are more invasive, but have lower densities compared to the case where we assume four non-local adhesion fluxes. We can also see that the OV density for the case described by model (3.2) (with two nonlocal fluxes) is almost double compared to the OV density for the case described by model (3.3) (with four nonlocal fluxes).

4.2.   The impact of adhesion strength

In this subsection we investigate the effect of cell-cell and cell-matrix adhesion strengths for the two non-local subpopulations (i.e., model (3.2)) versus the four non-local subpopulations (i.e., model (3.3)). In ^[34] the authors studied the impact of different adhesion strengths in a model with one homogeneous cancer population and showed that when cell-cell adhesion strength was lower than cell-matrix adhesion strength it led to larger tumour spread. Since here we focus on two different cancer cell sub-populations (i.e., migrating and proliferative), we assume that cell-cell adhesion strength is lower than cell-matrix adhesion strength for the migrating cancer cells $c_m, i_m$ (to allow for cell migration), and the other way around for the proliferative cancer cells $c_p, i_p$ (to reduce cell migration). The results of numerical simulations with these different adhesion strengths are shown in Figure 5. In this case we see relatively similar tumour spread patterns for (a) model (3.2) with two non-local sub-populations and for (b) model (3.3) with four non-local sub-populations. The only difference is a slight increase in the density for uninfected cancer cells and a decrease in the density of infected cancer cells for the model in sub-panels (b).

Since it is difficult to measure the adhesion strengths for cells with different phenotypes, in Figure 6 we also investigate numerically what happens with tumour and virus spread patterns when we assume that all cancer subpopulations have similar cell-cell adhesion strengths that are lower than their cell-matrix adhesion strengths. In this case we see that the cancer cells show less spatial spread compare to the case in Figure 5. Moreover, the OV (which has the highest density in the middle of the tumour mass) cannot destroy the tumour in that region; this is more evident in sub-panels (b) (for model (3.3) with four non-local sub-populations), where the level of the virus is also very reduced. In sub-panels (a) (for model (3.2) with two non-local sub-populations) we still see a bit of reduction in tumour size in the middle of the tumour region where the level of OV is similar as in Figure 5.

We conclude from these two numerical studies that the magnitudes of cell-cell and cell-matrix adhesion strengths for different cancer cell phenotypes (here migrating and proliferative cells), combined with their local/non-local character, influence significantly the spread of the virus through the tumour.

In the next three sub-sections we return to model (3.1) without haptotaxis for $c_{p}$ cells (i.e., $\eta_{c_p} = 0$), and investigate the impact of transition rates between migrating and proliferative cell sub-populations, as well as the impact of virus infection and replication rates.

4.3.   The impact of transition rates

We return now to model (3.1) without haptotaxis for $c_{p}$ (i.e., $\eta_{c_p} = 0$), and investigate numerically the impact of varying the transition rates between migrating and proliferating cancer cells, and the differences from the baseline results shown before in Figure 4(a). In Figure 7 we investigate the spatial spread of tumour and virus populations when (a) $\omega_2 = \omega_1 = 0.1$, and (b) $\omega_2 = \dfrac{\omega_1}{4} = 0.025$. We see that decreasing $\omega_2$ leads to an increase in the density of migrating cells ($c_m, i_m$, in sub-panel (b)) compared to the case in Figure 4(a). Moreover, we see an increase in the spatial spread of the tumour between sub-panel (a) and sub-panel (b) where there are more migrating cancer cells.

4.4.   The impact of infection rate

In Figure 8 we investigate the impact of the infection rates of cancer cells by the virus particles for the model (3.1) when we ignore the haptotaxis for $c_p$ (i.e., $\eta_{c_p} = 0$); we compare the results with those in Figure 7(b). In Figure 8(a) we assume that the proliferating cells have a faster infection rates compared to the migrating cells: $\varrho_p = 3\varrho_m = 0.316$. In Figure 8(b) we assume that the migrating cells have a faster infection rates compared to the proliferating cells: $\varrho_m = 3\varrho_p = 0.316$. We see that by increasing the OV infection rate for any cancer subpopulation it leads to an increase viral density, better viral spread and better killing of cancer cells. This cancer-killing effect is slightly more pronounced in sub-panels (a) where $\varrho_p > \varrho_m$.

4.5.   The impact of replicating rate of OVs

In Figure 9 we investigate the impact of varying the OV replication rate for model (3.1) without $c_{p}$ haptotaxis (i.e., $\eta_{c_p} = 0$), and the results are compared to Figure 7(b). When (a) $b_p = 50 \geq b_{m} = 40$, the impact of increasing $b_p$ by a small amount is more evident on the proliferative cancer population $c_p$ and on the virus density. When (b) $b_p = \dfrac{b_m}{8} = 5$ and $b_{m} = 40$, it is clear that decreasing $b_p$ while keeping $b_m$ fixed leads to a very low OV density and a higher density of proliferative cells $c_{p}$ in the middle of the tumour mass compared to the case (a) when $b_{p} > b_{m}$. Compared to the results in Figure 7(b) we deduce that increasing $b_{p}$ leads to a reduction in $c_{p}$ but not $c_{m}$ and an increase in virus $v$ levels, while decreasing $b_{p}$ leads to an increase in both $c_{p}$ and $c_{m}$ and a drastic reduction in virus $v$ levels.

5.   Conclusions

In this study we proposed a new multiscale moving boundary model that considers the local/non-local interactions between cancer cells and ECM, as well as the infections of cancer cells with oncolytic viruses (OV), all in the context of the go or grow hypothesis. This model generalises the previous studies in ^[34] (that focused on nonlocal multi-scale moving boundary models for oncolytic virotherapies in the context of a homogeneous cancer population) and ^[29,58] (that focused on local multiscale moving boundary models for oncolytic virotherapies in the context of a homogeneous cancer population). Here, we consider a heterogeneous cancer cell population formed of two sub-populations: mainly-migrating and mainly-proliferative cells.

Using this new model, we investigated not only the impact of different cell-cell and cell-ECM interaction strengths on the overall spread of cancer cells and OVs (see Figures 4-6), but also the effect of changes in the transition rates between the migrating and proliferative cells (see Figure 7), as well as the effects of varying the infection rates of different cancer cells (Figure 8), and the proliferation rates of viruses inside different cancer cells (Figure 9). First, we have seen that the magnitudes of cell-cell and cell-matrix adhesion strengths for different cancer cell phenotypes (i.e., migrating and proliferative cells), combined with their local/non-local character, influence significantly the spread of the virus through the tumour. Second, we have seen that the killing of cancer cells by the OVs is slightly more pronounced when the proliferating cells have a faster infection rate compared to the migrating cells (i.e., $\varrho_p > \varrho_m$). This suggests that giving the virus during a certain time interval, when the majority of cells in the solid tumour are in a proliferative phase, might eventually lead to better cancer killing. Finally, we have seen that viral replication inside proliferating cells (and viral burst size from these cells) might affect in some cases not only the density of proliferative cells but also the density of migrating cells. For example, when $b_{p} > b_{m}$, the migrating cells do not seem to be greatly impacted; however, for $b_{p} < b_{m}$, the migrating cells are impacted.

To conclude, we emphasise that the heterogeneity of solid cancers (formed of sub-populations of cells with different phenotypes; e.g., mainly-migrating and mainly-proliferative cells) might impact the success of oncolytic therapies, where the virus needs to spread throughout the tumour to be able to eliminate it. We would also like to suggest that one possible explanation for the contradictory experimental results in regard to validity of the go or grow hypothesis for various cancer cell lines (see our discussion in the Introduction) might be related to the differences in the heterogeneity of tumours for these different cell lines. However, this hypothesis will have to be tested experimentally in the future.

Acknowledgments

The first author would like to acknowledge the financial support received from the Saudi Arabian Cultural Bureau in the UK on behalf of Taibah University, Medina, Saudi Arabia.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix

A.   Indicator functions enabling computations exclusively on the tumour moving domain

In order to carry out the macro-scale computation exclusively on the developing tumour $\Omega(t_0)$, for any $t_{0} > 0$, we define the spatial cancer indicator function $\mathcal{I}(\cdot, \cdot):\lbrace1, ..., M\rbrace\times\lbrace1, ..., M\rbrace \rightarrow \lbrace0, 1\rbrace$ given by

Further, in order to identify the immediate outside nodes that are neighbours to the tumour boundary (which are needed during the computation), we define the spatially closest outside neighbour indicator functions $\mathcal{H}_{x, +1}(\cdot, \cdot)$, $\mathcal{H}_{x, -1}(\cdot, \cdot)$, $\mathcal{H}_{y, +1}(\cdot, \cdot)$, $\mathcal{H}_{y, -1}(\cdot, \cdot):\lbrace2, ..., M-1\rbrace \times \lbrace2, ..., M-1\rbrace \rightarrow \lbrace0, 1\rbrace$, defined by

These enable us to exercise the computations for the cancer cells dynamics only on the expanding cancer region $\Omega (t_0)$. Indeed, these indicator functions allow us to take advantage of the zero-flux conditions at the moving boundary $\partial \Omega(t_{0})$ and appropriately assign density values that are required in computations (and that are dictated by the boundary conditions) at each of these closest outside points to the tumour boundary (along each spatial direction). Finally, these closest outside points to the tumour boundary are given by the union of pre-images $\mathcal{H}_{x, -1}^{-1} (\lbrace 1 \rbrace)\cup \mathcal{H}_{x, +1}^{-1} (\lbrace 1 \rbrace)\cup \mathcal{H}_{y, -1}^{-1} (\lbrace 1 \rbrace)\cup \mathcal{H}_{y, +1}^{-1} (\lbrace 1 \rbrace)$.

B.   Finite difference-midpoint involving indicator function

This section explain finite difference-midpoint method used through the paper, into two subsections local flux and non-local flux to accounts for all fluxes considered in this paper.

B.1.   Non-local flux

At any spatial discretised spatial location $(x^{1}_s, x^{2}_p) \in \Omega(t_0)$, and any discretised time $t_{t}: = t_{0}+l\delta t$, $\forall l\in \{0, \dots, [\Delta t /\delta_{t}]\}$, we denote by $c^{l}_{s, p}$ the discretised densities at the $((x^{1}_s, x^{2}_p), t_{l})$ for each of the subpopulations $c\in\{c_{p}, i_{p}, c_{m}, i_{m}\}$. Similar notations we adopt also for the discretised adhesion fluxes at $((x^{1}_s, x^{2}_p), t_{l})$, namely ${\mathcal{A}}_{c, s, p}^l$ stands for the discretised adhesion flux for each of the subpopulations $c\in\{c_{p}, i_{p}, c_{m}, i_{m}\}$. In this context, using the zero-flux conditions across the tumour's moving boundary (accounted here for via the indicators in equations Eqs (A.1) and (A.2)), we are able to cary out computations on the expanding spatial mesh by involving midpoint approximations for the cancer subpopulation densities

as well as for the adhesion fluxes

while the central differences at the virtual midpoint nodes $(s, p \pm \frac{1}{2})$ and $(s\pm \frac{1}{2}, p)$ are given by:

● for $c$:

Hence, the discretisation of the spatial operator $\nabla \cdot \left[D_c \nabla c -c{\mathcal{A}}_{c}(t, x, \textbf{u}(t, \cdot))\right]$ in Eq (2.25c) is obtained by

Denoting now by $\mathcal{F}_{c, s, p}^l$ the discretised value of the flux $\mathcal{F}_{c}: = D_c \nabla c -c{\mathcal{A}}_{c}(t, x, \textbf{u}(t, \cdot))$ at the spatio-temporal node $((x_s, x_p), t_l)$, we observe that the discretisation of $\nabla \cdot \mathcal{F}_{c} = \nabla \cdot \left[D_c \nabla c -c{\mathcal{A}}_{c}(t, x, \textbf{u}(t, \cdot))\right]$ given in Eq (2.25c) can therefore be equivalently expressed in a compact form as

where

B.2.   Local flux

For the case when the spatial flux is given through the local operator $\mathcal{F}_{c}: = D_c \Delta c - \eta_c \nabla \! \cdot \! c\nabla e)$ where $c \in \lbrace c_p, i_p, i_m \rbrace$, we have the discretisation

where

and

and for $e$:

Finally, for the haptotactic-like virus spatial operator $\eta_v \nabla \!\cdot \!(v\nabla e)$ is discretised as

C.   The predictor-corrector

We denote by $H_{{c}}(\cdot, \cdot, \cdot)$ the right-hand side spatial operator for each of the tumour subpopulations Eq 2.25(a)-(d), with the non-transport part of this operator further denoted by $f_{c}({\textbf{u}}, v)$. Then, maintaining the same notation style, at any discretised spatio-temporal node $((x_s, x_p), t_l)$, we have that

On the time interval $[t_l, t_{l+1}]$, we first predict $c$ at $t_{l+\frac{1}{2}}$ using an explicit method as follows

Further, using $\tilde{c}_{s, p}^{l+\frac{1}{2}}$ we calculate the corresponding predicted flux $\tilde{\mathcal{F}}_{c, s, p}^{l+\frac{1}{2}}$ at $t_{l+\frac{1}{2}}$. Then, we construct a non-local corrector that involves the average of the flux at the active neighbouring spatial locations

Denoting the set of indices corresponding to these active locations by $\mathcal{N}$, we we establish a corrector for the spatial flux, given by

which ultimately enables us to use a trapezoidal-type approximation and obtain a corrected value for $c$ at $t_{l+\frac{1}{2}}$ as

Here, $\tilde{\bf{u}}_{s, p}^{l+\frac{1}{2}}$ is the "half-time predicted" tumour vector obtained for the current subpopulation $c$ given by the half-time predicted value $\tilde{c}_{s, p}^{l+\frac{1}{2}}$, while all the other tumour subpopulations and the ECM (that enter as components of ${\textbf{u}}$) being given by their accepted values at $((x_s, x_p), t_l)$.

Then, using

to calculate the flux $\mathcal{F}_{c}^{l+\frac{1}{2}}$ at $t_{l+\frac{1}{2}}$, we initiate the predictor-corrector steps described above on the second half of this time interval, which we now repeat on $[t_{l+\frac{1}{2}}, t_{l+1}]$, and we finally obtain $c_{s, p}^{l+1}$ given through the analogue of Eq (C.5) that results (through the same procedure, described in Eqs (C.2)-(C.5)) on this second half-interval (i.e., on $[t_{l+\frac{1}{2}}, t_{l+1}]$. Finally, the same time-marching steps are adopted also for virus equations Eq (2.25f).