Nonlocal multiscale modelling of tumour-oncolytic viruses interactions within a heterogeneous fibrous/non-fibrous extracellular matrix

1.   Introduction

Cancer invasion is a complex multiscale phenomenon that relies on the structure and composition of the extracellular matrix (ECM). The ECM is an important biological structure that serves as a platform for cellular communication, as well as providing support to surrounding cells and tissues, transducing mechanical signals, and functioning as adhesive substrate ^[1]. This matrix is a highly dynamic structure that controls most fundamental behaviours of cells: from proliferation, to adhesion, migration, differentiation and apoptosis ^[2]. These processes are controlled through the interactions of cells with the components of the ECM: collagen, proteoglycans, elastin and cell-binding glycoproteins ^[2]. The matrix components are continuously deposited, degraded or modified, and thus the ECM is continuously undergoing remodelling. This remodelling process (which involves matrix degrading enzymes that can be secreted by the cancer cells and normal cells) impacts also the evolution of cancer invasion. A fundamental process in the invasive potential of cancer cells is the cell-ECM adhesion (through adhesion molecules such as integrins) ^[3]. In addition, the collective movement of cancer cells is the result of cell-cell adhesion (through adhesion molecules such as E-cadherins) ^[3]. Understanding these cell-cell and cell-matrix adhesion processes is important not only for our understanding of the evolution of cancer and its invasion of the surrounding tissue, but also for our understanding of the efficacy of various anti-cancer therapies.

An emerging effective anti-cancer treatment is the oncolytic virotherapy (OV-therapy). The effectiveness of this treatment lies on viruses selectively infecting and destroying malignant cancer cells without harming the surrounding healthy cells; see ^[4,5,6,7]. The structure and composition of the ECM influences the effectiveness of the OV-therapies, since the ECM components can form a physical barrier that traps the viral particles ^[8]. In particular, experimental studies have shown that the collagen fibres play a very important role in inhibiting viral spread ^[8].

Due to the complexity of the tumour microenvironment, which makes it difficult to understand the interactions between the different components of this environment, mathematical models have been used over the last few decades to answer various questions about these interactions. The great majority of these models are single-scale models, which focus on spatial tumour invasion ^[9,10,11], on tumour oncolytic therapies ^{[12,13,14,15,16,17,18]}, or both ^{[19,20,21,22,23]}. More recently, various multi-scale mathematical models have been derived to reproduce and investigate biological processes that take place at different spatial scales ^{[24,25,26,27,28,29,30,31]}. For example, ^[30] introduced a multi-scale moving boundary model for cancer invasion, which focused on the local interactions between cancer cells and the ECM, via matrix degrading enzymes (MDEs) that act at the micro-scale level of the invading tumour boundary. This study was further extended in ^[31] by considering also the non-local cell-cell and cell-ECM adhesive interactions. This study also considered a heterogeneous ECM population formed of fibres and non-fibrous sub-populations, and investigated the role of fibres (and their re-arrangement at macro-scale and micro-scale) in the evolution of a tumour population. Other local and nonlocal multiscale models were used in ^{[24,25,26,27]} to investigate the interactions between oncolytic viruses and cancer cells.

In this study, we develop further the non-local multiscale approach for modelling cancer-OV interactions introduced in ^[24] (and which focused on the role of cell-cell and cell-matrix adhesive interactions on the spread of OV throughout solid tumours) by combining it with the heterogeneous ECM approach proposed in ^[31], to investigate these tumour-ECM-OV interactions in fibrous ECM; see Figure 1.

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.

2.   Mathematical model

The multiscale moving boundary model used here is based on the two-scale (tissue scale-macro-scale, and cell scale-micro-scale) moving boundary framework introduced in ^[30] and recently applied to nonlocal cell-cell interactions in the context of oncolytic viral therapies ^[24]. Furthermore, in here we explore the dynamic interaction between an invading solid tumour, OV, and a two-component ECM (which was first introduced by ^[31]). This complex dynamic is captured by two interconnected multiscale systems that share the same macro-scale cancer dynamics at the tissue-scale. However, at cell-scale we use two distinct micro-scale dynamics for fibre rearrangement and for cancer invasion boundary, both being linked to the macro-dynamics through two double feedback loops, as illustrated in Figure 2.

2.1.   Macro-scale dynamics

Let us denote by $Y \subset {\mathbb{R}}^{2}$ the macro-scale domain representing the maximal environmental tissue square that we consider in this study. Also, let us denote by $\Omega(t)\subset Y$ the spatial support of the growing tumour within the macro-scale domain, at a time $t\in[0, T]$. Let $c(x, t)$ and $i(x, t)$, $\forall (x, t)\in\Omega(t)\times[0, T],$ represent the spatio-temporal densities of the uninfected cancer cells and the infected cancer cells, respectively. Finally, let $v(x, t)$ and $e(x, t),$ $\forall (x, t)\in Y \times [0, T]$, represent the spatial-temporal densities of the oncolytic virus and the cumulative extracellular matrix, respectively. The later is defined as

where $F(x, t)$ denotes the macroscale spatial distribution of the fibre ECM, which accounts for all significant ECM fibres such as collagen fibres or fibronectin fibrils, and $E(x, t)$ denotes the spatial distribution of the non-fibre ECM, which includes all the other non-fibre components of the ECM, i.e., elastin, laminins, fibroblasts, etc. At any spatio-temporal node $(x, t) \in Y \times [0, T]$, the macroscale dynamics for ECM is described by the following equations:

where $\alpha_{c}, \alpha_{i}, \alpha_{cF}, \alpha_{iF} > 0$ are the ECM degradation rates caused by cancer cells subpopulations $c$ and $i$, respectively. Further, $\mu_{2} > 0$ is a fixed remodelling rate. The fibre micro-dynamics will be discussed in details in Section 2.2. Defining $\textbf{s}$ the tumour-ECM vector,

then the volume fraction of space occupied by the tumour and the ECM is given 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.

In this study tumour dynamics consists of three main components: motility, proliferation, and death. Moreover, we describe the spatial fluxes for both cancer cells $(c, i)$ by a combination of a linear diffusion term caused by cells' random walk and a directed migration term due to cell-cell and cell-ECM adhesion. Here we investigate directed cell migration from two perspectives: local migration and non-local migration as a result to cell adhesion ^{[32,33,34,35,36,37]}. Therefore, the dynamics of the uninfected cancer cells subpopulation is given by:

where $D_{c} > 0$ is a constant diffusion coefficient, $\mu_{1} > 0$ is a constant proliferation coefficient, $\varrho > 0$ is a constant rate at which the uninfected cancer population diminishes due to infection by the oncolytic virus $v$, while ${\mathcal{A}}_c(t, x, \textbf{s}, \theta_f)$ is a non-local spatial flux term that describes the cells adhesion process that causes cancer cells to move in a directed manner. In references [24,25] the authors studied the effects of cell adhesion process on cancer-OV interaction, by focusing on cell-cell adhesion and cell-ECM-non-fibres substrate. Here we adopt the modelling concept proposed by ^[31], to consider the essential role performed by the cell-fibres adhesive interaction. Denoting the cell-cell adhesion function by:

where ${\textbf{S}}_{cc} > 0$ and ${\textbf{S}}_ {ci} > 0$ are the strengths of cell-cell adhesion and cross adhesion bonds, respectively that are formed between cancer cells distributed at $x$ and cells at $x+y$. Since the cell-cell adhesion strength $\textbf{S}_{cc}, \textbf{S}_{ci}$ is dependent on the quantity of intercellular $Ca^{2+}$ ions available within the ECM ^[38,39]. We adopt a similar approach as in ^[31] to compute cell-cell adhesion strength, and extend it to determine cross cell-cell adhesion strength as follows:

where $\textbf{S}_{\cdot \cdot} \in \lbrace \textbf{S}_{cc}, \textbf{S}_{ci}, \textbf{S}_{ic}, \textbf{S}_{ii} \rbrace$ with $S_{\cdot \cdot}^{\, max} \in \lbrace S_{cc}^{\, max}, S_{ci}^{\, max}, S_{ic}^{\, max}, S_{ii}^{\, max} \rbrace$ respectively. $S_{\cdot \cdot}^{\, max}$ is fixed and represent the maximum strength of cell-cell adhesive junctions. Within a sensing radius $R > 0$ at time $t > 0$, the non-local adhesive flux is defined as follows:

where ${\bf{B}}(0, R): = \lbrace z\in {\mathbb{R}}^2:|z|\leq R\rbrace$ is a closed ball centred at origin and of radius R, called here the sensing region. $\textbf{S}_{ce}$ is a constant that describes the cell-ECM adhesion strength. $\textbf{S}_{cF}$ is a constant that describes the cell-fibre-ECM adhesion strength. $\chi_{_{\Omega(t)}}(\cdot)$ is the characteristic function of $\Omega(t)$, and $(1-\rho(\textbf{s}))^{+} : = \rm{max} \lbrace (1-\rho(\textbf{s})), 0\rbrace$ is a threshold term to avoid local overcrowding, and ${\textbf{n}}(y)$ denotes the unit radial vector giving by:

with $\parallel {\cdot} \parallel_{_{2}}$ representing the usual Euclidean norm. $\bar{{\textbf{n}}}(y, \theta_f(x+y, t))$ is the unit vector that is dependent on the fibre orientations, defined as follows:

where $\theta_f(x+y, t)$ is the orientation of the fibres at macro-scale, this was first introduced by ^[31]. This orientation is derived by the micro-scale mass distribution of micro-fibres $f(., t)$, in the sense that affects the cell-ECM adhesion and characterises the ECM fibres distributed at the macro-scale location $x\in Y$ namely $F(x, t)$; for detailed mathematical formulas see Section 2.2. Note that $F(x+y, t)$ describes the influence of fibres distributed at $x + y$ on the adhesion of cells at location $x$ (with adhesion strength $S_{cF}$); see Figure 2 for more details. Furthermore, the radially symmetric kernel $\mathcal{K}(\cdot):[0, R]\to [0, 1]$ explores the dependence 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): = x + {\bf{B}}(0, R)$. Since these adhesion junction strengths are assumed to decrease as the distance $r: = \parallel {x-\zeta} \parallel_{_{2}}$ increases, $\mathcal{K}$ therefore is taken here of the form

By summing up the radially distributed adhesive interactions between the cells at $x$, and the cells and ECM at $x+y$ within $y\in {\bf{B}} (x, R)$, the term $\frac{1}{R}$ that appears in the front of the expression (2.7) is simply an interaction range normalisation factor.

Next, we focus on the infected cancer cell subpopulation $i(x, t)$ that emerges within this dynamics due to infections by the OV. We denote by $\varphi_{i}({{\textbf{s}}})$ the effect of the cell adhesion processes that take place either locally (through adhesive interactions between infected cancer cells and ECM, as tumour cells exercise haptotactic movement 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). Mathematically, this can be formalised as

where $\eta_{i}$ is a constant haptotactic rate associated to $i$, while ${\mathcal{A}}_{i}\big(\cdot, \!\cdot, \! {\textbf{s}}(\cdot, \!\cdot) \!\big)$ is a non-local spatial flux term defined as in Eq (2.7). We define ${\mathcal{T}}_i(x+y, t)$ as follows

where $\textbf{S}_{ic}$ and $\textbf{S}_{ii}$ are are the strengths of infected cell-cell adhesion and cross adhesion bonds, and are dependent on ECM as defined in Eq (2.6). Then, ${\mathcal{A}}_{i}\big(\cdot, \!\cdot, \! {\textbf{s}}(\cdot, \!\cdot)\!\big)$ becomes

where $\textbf{S}_{ie}$ is a constant describing the cell-ECM adhesion strength and $\textbf{S}_{iF}$ is a constant describing the cell-fibre-ECM adhesion strength. Thus, the spatio-temporal dynamics of the infected cancer cell subpopulation is governed by the following equation

Here $D_{i} > 0$ is a constant random motility coefficient, and $\varphi_{i}({\textbf{s}})$ represents the directed migration induced by the cell-adhesion processes that corresponds to $i$ and is described in Eq (2.10) Further, the infected population expands at a rate $\varrho$ due to new infections occurring among the uninfected cells, and dies at rate $\delta_{i} > 0$.

Next, for the oncolytic virus spatio-temporal dynamics, we adopt here a similar reasoning as in ^[24,25], and 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 is governed by

Here $D_{v} > 0$ is a constant random motility coefficient, $\eta_{v} > 0$ is a constant haptotactic coefficient, $b > 0$ is a viral replication rate within infected cancer cells, and $\delta_{v} > 0$ is the viral death rate.

Finally, the coupled interacting tumour-OV macro-dynamics is governed by Eqs (2.4)–(2.15) in the presence of initial conditions

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

2.2.   Fibre micro dynamics on the bulk of the tumour

The cancer cells macro-dynamics cause the fibres to undergo a microscopic rearrangement process in addition to the macroscale fibre degradation mentioned in Eq (2.2). We start by defining the micro-fibre domain as $\sigma Y(x) : = x + \sigma Y, \forall x \in Y$ with scale size $\sigma > 0$ following ^[31], in which the fibre ECM rearrangement occurs reflecting on fibre macroscale orientation as illustrated in Figure 2. At any macroscale point $x \in \Omega(t)$, the ECM-fibre phase is described by a macroscale vector field $\theta_{f}(x, t)$, which gives us the amount of fibres distributed at $(x, t)$ and their macroscopic fibres orientation computed by the revolving barycentral orientation $\theta_{f, \sigma Y(x)}(x, t)$ generated by the microscopic mass distribution of microfibres $f(\cdot, t)$ within the micro-domain $\sigma Y(x)$. Based on the approach introduced by ^[31], the macroscale fibre orientation $\theta_{f}(x, t)$ is defined as follows:

where $\lambda(\cdot)$ is the usual Lebesgue measure, $\|\cdot\|_{2}$ represents the usual Euclidean norm, and

is the naturally generated revolving barycentral orientation $\theta_{f, \delta Y(x)}(x, t)$ associated with $\sigma Y(x)$ and given by the Bochner-mean-value of the position vectors function $\sigma Y(x) \ni z \mapsto z-x \in {\mathbb{R}}^{N}$ with respect to the density measure $f(x, t) \lambda(\cdot)$. The macroscopic mean-value fibre representation at any $(x, t)$ is then given by the Euclidean magnitude of $\theta_{f}(x, t)$, namely,

At any time $t$ and at any spatial location $x \in Y$, the cancer cells realign the micro-fibres through a microscopic rearrangement process in each micro-domain $\sigma Y(x)$ that is triggered by the combined macroscale spatial flux of both cancer cell subpopulations as follows:

where

where $\bar{\varphi_{i}}({\textbf{s}})$ denotes here the spatial flux triggered by the interactions between the infected cancer cell population and the surrounding cells and ECM, and is given by

For simplicity, denoting the total cancer cell population by $c_{ \rm{total}}(x, t) = c(x, t)+i(x, t)$. The combined flux $\mathcal{F}(x, t)$ acts upon the micro-scale distribution $f(z, t), \forall z \in \sigma Y(x)$ in accordance to the magnitude that the total mass of cancer cells has relative to the combined mass of cells and fibres at $(x, t)$, which is given by the weight

At the same time, the total spatial flux of cancer cells $\mathcal{F}(x, t)$ is balanced in a weighted manner by the orientation $\theta_{f}(x, t)$ of the existing distribution of fibres at $(x, t)$ that is appropriately magnified by a weight that accounts for the magnitude of fibres versus the combined mass of cells and fibres at $(x, t)$ and is given by $(1-\omega(x, t))$. As a consequence, the micro-scale distribution of micro-fibres $f(z, t), \forall z \in \sigma Y(x)$ is therefore acted upon uniformly by the resultant force given by the following rearrangement vector-valued function

In this context, a mass distribution of the micro-fibres $f(z, t)$ on $\sigma Y(x)$ is exercised under the influence of this rearrangement vector $r(\sigma Y(x), t)$, resulting in spatial relocation of micro-fibres on both $\sigma Y(x)$ and its neighbouring micro-domains. Indeed, under the incidence of $r(\sigma Y(x), t)$, a certain fraction of the micro-fibres positioned at a given $z\in \sigma Y(x)$ get transported at new micro-scale position $z^{*}$, given by

where $\nu_{\sigma Y(x)}(z, t)$ is the emerging relocation vector:

Here, $x_{\mathrm{dir}}(z) = \overrightarrow{xz}$ is the barycentric position vector pointing to $z$ in $\sigma Y(x)$ see Figure 3, $f_{max}$ represents the maximum level of fibres that could reside at the micro-location $z\in \sigma Y$ at any given time, and $f^{*}: = f(z, t)/f_{max}$ is the local micro-fibres saturation fraction. Finally, the relocation magnitude in the direction $x_{\mathrm{dir}}(z)+r(\sigma Y(x), t)$ is simultaneously mediated by the ability of micro-fibres to dislocate (which is exercised when these are not at the maximum level) and is attenuated by the level of the micro-fibres mass fraction $f^{*}$ at $z$ in conjunction with the barycentric position defect quantified here by $\left\|r(\sigma Y(x))-x_{\mathrm{dir}}(z)\right\|_{2}$. Therefore, provided that micro-fibres levels at $z$ are not at their maximum level $f_{\max }$ a micro-fibres mass transport from $z$ to the location $z^*$ is exercised, while lower levels of micro-fibres saturations at $z$ combined with a better flux alignment provided by a smaller position defect result in a relocation of the micro-fibres mass in direction $(x_{\mathrm{dir}}(z) + r(\sigma Y(x), t))$ at a greater distance. Finally, this transport occur provided that a enough room is available at position $z$. This being is naturally mediated through the movement probability

that quantifies the capacity that is still available at $z^{*}$, and enables only the amount $f(z, t)P_{ \rm{move}}$ of micro-fibres to be transported to the new location and the rest $f(z, t)(1-P_{ \rm{move}})$ remain at their location.

2.3.   MDEs microscale dynamics at the leading edge and the induced tumour boundary movement

In addition to the fibres micro-dynamics that occurs in the bulk of the tumour, a second type of tumour micro-dynamics is generated by the cell scale activity of the matrix degrading enzymes (MDEs) along the tumour invasive edge ^[40,41]. Indeed, secreted by the cancer cell population from the outer proliferating rim of the tumour, the MDEs are transported across the tumour interface within a cell-scale (micro-scale) neighbourhood of the tumour boundary in the peritumoural region, where they cause degradation to the ECM distribution that they meet, giving rise to further tumour progression. To capture this process mathematically, we adopt here a similar approach to the one developed in ^[30], whereby the entire cross-interface MDEs transport micro-process along $\partial \Omega(t)$ is decomposed into a union of boundary micro-processes that take place on a appropriately constructed covering bundle $\{\epsilon Y\}_{\epsilon Y\in {\mathcal{P}}_{\epsilon}(t)}$ of $\epsilon-$size overlapping micro-domains $\epsilon Y$, whose union form the cell scale neighbourhood ${\mathcal{N}}_{\epsilon}(\partial \Omega(t))$, i.e.,

where the entire enzymatic process and its consequences are explored. Thus, to explore the MDEs enzymatic activity over a small time range $\Delta t > 0$, denoting by $m(y, \tau)$ the micro-scale density of MDEs at the micro-scale spatio-temporal location $(y, \tau)\in \epsilon Y \times [0, \Delta t]$, we have that a micro-scale source of MDEs $\mathcal{G}_{\epsilon Y}:\epsilon Y \times [0, \Delta t] \to [0, \infty)$ is induced naturally as a collective contribution of the cancer cells that are within a small distance $\rho > 0$ from each $y\in \epsilon Y$, and so this can be formalised mathematically as

Here $\lambda(\cdot)$ is the standard Lebesgue measure on $\mathbb{R}^N$, and ${\bf{B}}(z, r): = \lbrace x \in Y: \left\lVert{z-x}\right\rVert_\infty \leq \rho\rbrace$ is the small active tumour region closed to the tumour interface where the cancer cells collectively contribute to the creation of the source of MDEs at $z\in \epsilon Y$ over the time interval $([t_{0}, t_{0}+\Delta t]$. Finally, and $\gamma_{c} > 0$, $\gamma_{i} > 0$ represent constant MDEs secretion contributions of the uninfected and infected cancer cells, respectively. Finally, in the presence of the MDEs source $\mathcal{G}_{\epsilon Y} ({z}, \tau)$, the MDES are assumed here to exercise a diffusive transport within the entire micro-domain $\epsilon Y$, which is mathematically formulated through the following reaction-diffusion equation

Furthermore, as no memory of pre-existing distributions of MDEs are assumed for the enzymatic process, and no molecular transfer is assumed across the boundaries of $\epsilon Y$, the proteolytic boundary micro-dynamics (2.29) that takes place on each $\epsilon Y$ is assumed to take place with null initial conditions and flux-zero boundary conditions, i.e.,

where ${\textbf{n}}_{_{\epsilon Y}}$ is the usual outward unit normal direction on $\partial\epsilon Y$.

Finally, the pattern of ECM degradation within the peritumoural region of the microscale neighbourhood $Ne_{\epsilon}(\partial \Omega(t))$ will correspond to the pattern of significant MDEs transport within each micro-domain $\epsilon Y$. In this context, following the multiscale mathematical modelling approach developed in ^[30], we obtain the law for the macro-scale boundary movement, which is specified in terms of direction and displacement magnitude for the spatial relocation of the tumour interface $\epsilon Y \cap \Omega Y$ for each boundary micro-domain. By determining the law of macro-scale tumour boundary movement, the MDEs micro-dynamics occurring along the tumour interface influences directly the tumour macro-dynamics, completing this way a bottom-up feedback from micro- to macro-scale, contributing significantly to the invasion and spread of tumour within the surrounding tissue ^[30,42].

2.4.   Summary of multiscale model

In summary, the multiscale moving boundary model that we obtained for the tumour-OV-ECM interactions (see Figure 2 for an illustrative sketch) consists of the following parts:

micro-dynamics of fibres triggered by the macro-scale-induced rearrangement flux

the MDEs boundarymicro-dynamics:

The macro-dynamics and the two micro-dynamics are connected through two double feedback loops:

● a tumour bulk top-down link by which the macro-scale cell flux triggers the micro-scale fibres rearrangement

● a bottom-up fibres link by which the micro-scale distribution of fibres naturally induces a spatial orientation that alters the tumour macro-dynamics

● a leading edge top-down link through which the tumour macro-dynamics induces a micro-scale MDEs molecular source in the cell-scale neighbourhood of $\partial \Omega(t)$

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

3.   Brief overview of the computational approach for a couple of tumour-OV-ECM interaction scenarios

The numerical approach and computational implementation of the novel multiscale moving boundary model proposed in this work builds directly on the multiscale moving boundary computational framework initially introduced by ^[30] and further expanded in ^[24,25,31], and so while for full details we defer the reader to these references, in the following we will give a brief summary of this approach.

Macro-scale numerical approach For all the computational simulations we consider the macro-scale tissue domain $Y: = [0, 4]\times [0, 4]$, and we discretise this with a uniform grid $Y^{d}: = \{(x^{1}_{i}, x^{2}_{j})\}_{i, j = 1\dots N}$, with $N = [4/h]+1$, of spatial step size $\Delta x = \Delta y: = h$, with $h > 0$. Correspondingly, the discretised version of the evolving spatial tumour support is denoted by $\Omega^{d}(t)$ (i.e., $\Omega^{d}(t) = Y^{d}\cap \Omega(t)$), with $\partial \Omega^{d}(t)$ standing for the evolving tumour boundary. In brief, for the approximation of the macro-scale dynamics on the evolving tumour domain $\Omega^{d}(t)$, we used a method of lines type approach as proposed initially in ^[30], with the specific off-grid-type approximation for the cell-adhesion fluxes ${\mathcal{A}}_c(\cdot, \cdot, \cdot, \cdot)$ and the novel non-local time-marching scheme that were proposed and developed in ^[31], and adapted here to our particular context of cancer-OV interacting dynamics. Finally, as detailed in ^[30], the discrete domain $\Omega^{d}(t)$ of the progressing tumour is appropriately evolved with additional spatial nodes corresponding to the new locations reached by the invading cancer which are determined by the MDEs micro-dynamics that takes place within a micro-scale neighbourhood of the tumour interface (of micro-scale radius $\epsilon > 0$, see ^[30,42] for complete details).

Approximating the two micro-dynamic processes that occur simultaneously In this model we have two types of micro-scale processes that link to the same macro-dynamics, namely the MDEs boundary micro-dynamics at the leading edge of the tumour and fibre micro-dynamics. Thus, for the MDEs boundary micro-dynamics part, proceeding as detailed in ^[30], the boundary MDEs micro-dynamics is approximated involving central finite differences for the spatial operators on each micro-domain $\epsilon Y$, with a backward Euler time marching scheme. Finally, for the fibres micro-dynamics and their rearrangements, the implementation follows closely the modelling details given in Section 2.2. Thus, for that part, considering the fibres micro-scale taking place on square micro-domains $\sigma Y$ that have their vertices at the dual mesh nodes (i.e., the barycentres of each square pixels formed by the macro-scale grid nodes), we proceed with the evaluation of the spatial flux operators and rearrangement vectors as well as the final numerical inference of the emerging spatial fibres orientation, as detailed in ^[31].

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

Our simulations of the multiscale model explore two distinct scenarios regarding the directed migration of infected cancer cell subpopulation: local and non-local migration at macro-scale (see Eq (2.31b)). Specifically, we consider the following cases:

(1) Local advective flux for the infected cancer cells: in Eq (2.31b) we have $\varphi_{i}({\textbf{s}}) = \eta_{i} \!\nabla\!\! \cdot \!\big(i\nabla e\big)$, and thus the macro-dynamics takes the following form:

Since $\bar{\varphi_{i}}({\textbf{s}}) = i\nabla e,$ the macroscale spatial flux for $\mathcal{F}_i$ in Eq (2.21) becomes

The numerical results for this case are shown in Figures 6–8.

(2) Nonlocal advective flux for the infected cancer cells: in Eq (2.31b) we have $\varphi_{i}({\textbf{u}}) = \nabla \!\! \cdot\! \big(i{\mathcal{A}}_{i}(\cdot, \!\cdot, \! {\textbf{s}}(\cdot, \!\cdot)\!)\!\big)$, and thus the macro-dynamics take the following form:

Since $\bar{\varphi_{i}}({\textbf{s}}) = i{\mathcal{A}}_{i}(x, t, {\textbf{s}}(\cdot, t), \theta_{f}(\cdot, t))$, the macroscale spatial flux for $\mathcal{F}_i$ in Eq (2.21) becomes

The numerical results for this case are shown in Figures 9–11.

Finally, in all our numerical experiments we choose to avoid simple overcrowding by imposing the regime that the cumulated sum of the total volume of cells and ECM to be always kept below $1$. This is obtained by taking unitary volume fractions for each of the involved tissue constituents, i.e., $\nu_{c} = \nu_{e} = 1$, as indicated in Table 1, which obviously leads to a simplified formula for the volume of occupied space, namely: $\rho(\textbf{s}(x, t)) : = e(x, t) + c(x, t)+i(x, t).$

3.2.   Initial conditions

The initial condition for the uninfected cancer cell population, $c(x, 0)$ is chosen to describe a small localised pre-existing tumour aggregation. This is given by

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

with the smooth compact support function $\psi$ given by

We assume no infection at this stage, i.e. zero infected cancer cells ($i(x, 0)$):

Motivated by several biological works underscoring the importance of the highly heterogeneous character of the fibrous microenvironment where tumours develop ^[43,44,45], we adopt the following in silico initial conditions the two-phase ECM considered in this modelling. Specifically, the initial condition for the non-fibre ECM density, $E(x, 0)$, is given by an arbitrarily chosen heterogeneous pattern described by the following equations (as in ^[31])

and is shown in Figure 4(c). Here, we have

Furthermore, the initial condition for one micro-scale fibre domain $\sigma Y(x)$ is shown in Figure 4(d), and it is repeated for all macro-scale locations. To determine the density of the fibres at any macro-point $x$, we integrate the corresponding fibre-micro domain $\sigma Y(x)$. Due to visibility reason, we avoid presenting the pattern of the microscale fibres on the macroscale. For the baseline simulations presented here, we choose the ratio of fibres to non-fibres components of ECM at $20\%:80\%$.

Finally, the initial condition for the OV population ($v(x, 0)$) is chosen to describe one single injection in the middle of the tumour aggregation, as in ^[24,25]:

where

which is smoothed out on the frontier of the viral density support $\Gamma_{v}: = \partial \{x\in Y\, |\, v^{0}(x) > 0\}$.

4.   Results

The numerical results shown in this section are computed using the parameter values listed in Table 1, which we refer to as 'baseline parameters' for ease of reference. Whenever we change these parameters, we clearly specify the new values we use for the simulations.

We start in Section 4.1 by investigating numerically the impact of fibre-ECM local approach for the infected cancer cells $i$ used to describe the cell-cell and cell-matrix adhesion flux on cancer-OV interaction. Then, in Section 4.2, we investigate the impact of varying the amount of fibres in the ECM ($\mathcal{R}_F$). Next, in Section 4.3 we increase the cell-fibre-ECM adhesion strength for some of the cases studied in the previous sections. In Section 4.4 we investigate numerically the impact of fibre-ECM non-local approach for the infected cancer cells $i$ used to describe the cell-cell and cell-matrix adhesion flux on cancer-OV interaction and comparing it to a variety of distinct adhesion strengths. In section 4.5 we investigate the impact of increasing the amount of fibre in ECM for the non-local system. Finally, in Section 4.6 we investigate the impact of different cross cell-cell adhesion strengths.

4.1.   Baseline dynamics for the model with local flux of infected cells

In Figures 5 and 6 we present the numerical simulations obtained for the system (3.1), with the baseline parameters in Table 1, under the assumption that the directed movement of infected cancer cells is described by a local flux term. While Figure 5 shows a time-series progression of the spatial distributions of uninfected cancer cells, infected cancer cells, and oncolytic virus, Figure 6 shows the spatial distribution of all macroscopic variables only at final simulation time $75 \Delta t$. Figure 6 includes also the distributions of non-fibre ECM and full ECM, which, due to lack of visible differences between the two, will be omitted in the next figures, and only results of the oriented fibre-ECM field will be shown in the next sections.

4.2.   The effect of varying the amount of fibres in the ECM

In Figure 7 we show simulations of system (3.1) using the parameters in Table 1, when we vary the ratio $\mathcal{R}_F$ of fibres and non-fibres components of ECM. In Figure 7(a) (first column) we consider $\mathcal{R}_F = 30\%:70\%$, in Figure 7(b) (second column) we consider $\mathcal{R}_F = 35\%:65\%$, and in Figure 7(c) (third column) we consider $\mathcal{R}_F = 40\%:60\%$. We note that as we increase the ratio $\mathcal{R}_F$, the amount of fibres increases forcing the uninfected cancer cells to migrate away from the centre of the domain, where they were initially located. This also leads to higher OV and overall lower tumour spread.

4.3.   The effect of increasing cell-fibre ECM adhesion strengths

In Figure 8 we show simulations of system (3.1) for the parameters in Table 1, but with different $\mathcal{R}_F$ and different cell-fibre ECM adhesion strengths. In Figure 8(a) (first column) we consider $\mathcal{R}_F = 20\%:80\%$, in Figure 8(b) (second column) we consider $\mathcal{R}_F = 30\%:70\%$, and in Figure 8(c) (third column) we consider $\mathcal{R}_F = 40\%:60\%$. Moreover, in all these sub-panels we take $S_{cF} = 0.5$ (compared to $S_{cF} = 0.2$ in Figure 7). We note that increasing cell-fibre ECM adhesion strength leads to spatial pockets of very high cancer density (i.e., $\max \; c(x, 75\Delta t) = 1$ in Figure 8(c), versus $\max \; c(x, 75\Delta t) = 0.55$ in Figure 7), and even a better spatial cancer spread.

4.4.   Dynamics of the model with nonlocal flux of infected cells

In this section we investigate numerically not only the impact of non-local advection fluxes for the infected cancer cells, but also the effect of varying the cell-cell and cell-matrix adhesion strengths. In Figure 9(a) we show the baseline dynamics of the system (3.3) (i.e., dynamics obtained with the baseline parameters in Table 1). In Figure 9(b) we keep most of the parameter the same, with the exception of $S_{cF} = S_{iF} = 0.3$. In Figure 9(c) we keep again most of the parameter the same, with the exception of $S_{ie} = 0.001$. In this case we note that varying the strengths of nonlocal cell-fibre interactions for the uninfected or infected cells does not have a significant impact on tumour structure.

4.5.   Increasing the amount of fibres for the model with nonlocal flux of infected cells

In Figure 10 we investigate numerically the effect of increasing the amount of fibres in the ECM from $\mathcal{R}_F = 20\%:80\%$ (the baseline case shown in Figure 9(a)) to $\mathcal{R}_F = 30\%:70\%$ (here) and varying cell-cell and cell-matrix adhesion strengths. More precisely, in sub-panels (a) we have $S_{cc} = 0.1$, $S_{ce} = 0.5$, while in sub-panels (b) we have $S_{cc} = 0.05$, $S_{ce} = 0.001$. The simulation results do not show significant differences between the two cases. However, if we compare these results in Figure 10(a) (nonlocal advection of infected cells) with those in Figure 7(a) (local advection for infected cells) we observe that accounting explicitly for nonlocal interactions leads to a better spread of cancer cells along the ECM fibres (so the cancer cells "sense" better the direction of the fibres).

4.6.   Cross adhesion strength

Finally, in Figure 11, we investigate numerically the effect of varying the cross adhesion strengths. In sub-panels (a) we assume that the cell-cell adhesion strengths for uninfected cancer cells ($S_{cc}$ and $S_{ci}$) are lower than the cell-cell adhesion strengths for infected cancer cells ($S_{ic}$ and $S_{ii}$). In sub-panels (b) we make the reversed assumption: the cell-cell adhesion strengths for the uninfected cells are higher than for the infected cells. We see that in this second case the tumour spreads faster through the domain (spread which is helped also by a stronger $S_{ce} = 0.5$, compared to case (a) where $S_{ce} = 0.001$).

5.   Conclusions

In this study we extended a nonlocal multi-scale moving boundary model proposed in ^[24] for oncolytic virotherapies, by considering cancer cell interactions with a heterogeneous ECM formed of fibrous and non-fibrous components. With the help of this model, we investigated numerically the impact of assumptions of local vs. non-local interactions between the infected cancer cells and uninfected cells and/or ECM. (It is known that the ECM and its components constitute a physical barrier in the spread of OVs ^[8], but it is not clear what are the interactions between the infected cancer cells and the environment, i.e., other cells and ECM.)

The numerical simulations showed that for the case of local advection of infected cells (see Figures 7 and 8), the ratio $\mathcal{R}_{F}$ of fibre versus non-fibre components of the ECM combined with the strength of cell-fibre ECM adhesion played an important role in the extent of spatial spread of tumour cells. In particular, very large $\mathcal{R}_{F}$ ratios seemed to cause an accumulation of OVs at particular positions in space (see Figures 7(a)–8(c), where tumours exhibited less spatial spread for large $\mathcal{R}_{F}$). Larger $\mathcal{R}_{F}$ ratios also led to a more pronounced tumour cell movement along the ECM fibres. This pronounced directional movement of cancer cells along the fibres could be obtained also for lower $\mathcal{R}_{F}$ ratios, but by considering explicit non-local interactions for the advection of infected cells (see Figure 10(a) versus Figure 7(a)). In this case, the spread of the virus inside tumour seemed to depend on the strength of cell-cell interactions, with larger $S_{cc}$ causing higher viral accumulation at specific positions in space (Figure 11(b)). The nonlocal interactions between the infected cells and the environment (i.e., other cells and ECM components) played an important role in tumour and viral spatial spread only when the magnitude of these interactions was very high (see Figure 11(a)).

All these numerical results suggested that for complex multi-scale biological problems it is difficult to discern between the different biological mechanisms behind the infected cells distribution (i.e., viral infection patterns), and directional versus omni-directional invasion of infected and non-infected tumour cells into the surrounding tissue. This uncertainty is also because many model parameters (and parameter ranges) used for the simulations were unknown and were estimated (i.e., "guessed").

The parametrisation of these types of multi-scale moving-boundary models, to make them more quantitatively relevant, is an open problem which involves different aspects. First, multi-scale data is not always available, and when such data exists it is the result of multiple experiments performed in different conditions for different biological systems and therefore it might not be at all relevant. Second, to apply the results to clinically-relevant situations one needs to perform an uncertainty and sensitivity analysis---which again, is an open problem for such multi-scale moving-boundary models. To make some progress with this issue, one needs to have multi-scale data from one single experimental setting, parametrise the model using this data, and then make further computational predictions that could be tested experimentally.

This work can be further generalised to investigate various aspects of the interactions between oncolytic virus particles and the tumour microenvironment: from the importance of combining OV therapies with other classical therapies such as chemotherapies, to the investigation of go or grow hypothesis within a two-phase heterogeneous ECM.

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.