Re-polarisation of Macrophages Within Collective Tumour Cell Migration: A Multiscale Moving Boundary Approach

Cancer invasion of the surrounding tissue is a multiscale process of collective cell movement that involves not only tumour cells but also other immune cells in the environment, such as the tumour-associated macrophages (TAMs). The heterogeneity of these immune cells, with the two extremes being the pro-inflammatory and anti-tumour M1 cells, and the anti-inflammatory and pro-tumour M2 cells, has a significant impact on cancer invasion as these cells interact in different ways with the tumour cells and with the ExtraCellular Matrix (ECM). Experimental studies have shown that cancer cells co-migrate with TAMs, but the impact of these different TAM sub-populations (which can change their phenotype and re-polarise depending on the microenvironment) on this co-migration is not fully understood. In this study, we extend a previous multi-scale moving boundary mathematical model, by introducing the M1-like macrophages alongside with their exerted multi-scale effects on the tumour invasion process. With the help of this model we investigate numerically the impact of re-polarising the M2 TAMs into the anti-tumoral M1 phenotype and how such a strategy affects the overall tumour progression. In particular, we investigate numerically whether the M2→M1 re-polarisation could depend on time and/or space, and what would be the macroscopic effects of this spatial- and temporal-dependent re-polarisation on tumour invasion.


INTRODUCTION
The last few decades have seen a shift in the focus of cancer research: from a research that was focused on individual tumour cells to a research that is now focused on collective cancer cells movement within the tumour microenvironment (TME) and the complex interactions between tumour cells and other types of cells inside the TME [1]. These processes are key for each of the stages of tumour progression, from the early development of the avascular tumour and its local invasion to angiogenesis and subsequent metastasis stages [2,3].
The TME is formed of tumour's vasculature, connective tissue, infiltrating immune cells and the extracellular matrix (ECM). In recent years the ECM has received considerable attention due to its role in cancer evolution and response to therapies [1]. The ECM is a complex network of macromolecules (such as fibrous proteins, water and minerals), which is an essential part of any healthy tissue [4]. To maintain its functionality, the ECM is subject to continuous remodelling via synthesis and degradation. The degradation of ECM is related to the over-secretion of several enzymes, such as the matrix metalloproteinases (MMPs), which are not explicitly tied to cancer cells, but can also be secreted by other stromal cells [5].
One of these stromal cell populations is the macrophages, which can form up to 50% of the tumour mass [6,7]. They are often classified into pro-inflammatory anti-tumour M1like macrophages and anti-inflammatory pro-tumour M2-like macrophages. A specific class of macrophages that is present in TME is called tumour-associated macrophages (TAMs). Since in established tumours TAMs have a phenotype that resembles M2 cells [8], they are in general correlated with poor prognosis in various cancers [9][10][11][12][13][14]. These TAMs are also involved in the active remodelling of ECM [8]. Experimental and clinical studies over the last two decades have shown that re-educating TAMs to exhibit anti-tumour responses, by switching cell phenotype from M2-like to M1-like cells, can lead to successful anti-cancer treatment protocols [8]. An example of such a treatment protocol involving macrophage re-polarisation, which was shown to be successful in human pancreatic cancer, involved the use of agonist anti-CD40 mAb [15]. Finally, an overview of the complex effects of macrophages within cancer progression and their interaction with the TME is provided in Aras and Zaidi [16].
Nutrients (e.g., oxygen and glucose) are essential for any cells to live and function properly, including cancer cells and macrophages. Nutrients are extravasated from the blood flow and diffuse through the ECM. In avascular solid tumours, nutrients diffuse into the tumour through the tumour boundary, and therefore tumour progression leads to the formation of regions inside the tumours with very low nutrient levels. These regions first become hypoxic and then necrotic, and cells in these regions undergo cell death being deprived of nutrients. Hypoxic conditions can also modify the polarisation of macrophages and influence the malignant behaviour of some cancer cells [17]. Biological evidence also suggests that specifically targeting these hypoxic regions and the metabolic activities of cells within these areas might be a beneficial therapeutic option [18,19].
Therefore, understanding the complex dynamic interactions between macrophages and tumour cells (direct antitumour/pro-tumour interactions, or indirect interactions via the degradation/remodelling of ECM), as well as their nutrient consumption patterns, could help us generate new hypotheses regarding the collective cancer cells invasion of surrounding tissue and the mechanisms involved in tumour progression and control.
Over the last decade, there have been substantial advances in mathematical modelling to understand the dynamics of the cancer cell populations as well as their interactions with their surroundings; see [20][21][22][23][24][25][26][27][28][29][30][31][32][33] and references therein. Although these models mainly focused on the interplay between the cancer cells and the neighbouring ECM, the importance of the macrophages during tumour development cannot be overlooked and as a consequence, some mathematical models started to investigate the role of macrophages on the overall tumour progression [26,28,30,[34][35][36][37][38][39]. Initially, these models [37][38][39] focused only on the anti-tumour role of the macrophages, and only later models started to explore also their pro-tumour role [26,28]. Furthermore, most of these initial mathematical models investigated tumour progression only at one spatiotemporal scale [20][21][22]31]. Since the relevance of various biological processes occurring on different scales could not be ignored, more recent models started to capture these multiscale interactions [23,24,29,30,32,40]. However, due the novelty of these multiscale approaches, they have not been extended to capture also the complex roles of macrophages in the tumour development.
In this work, we further extend an existing multi-scale moving boundary modelling platform [29,30,32]. To this end, we build upon a mathematical model [30] that accounts not only for the proteolytic processes occurring at the leading edge of the tumour but also for the fibre and non-fibre components of ECM, as well as for the presence of M2-like macrophages. In the new mathematical model introduced in this study we include the contribution of M1 TAMs to the overall tumour dynamics, by considering the impact of these cells on both macro-and micro-scales. We pay special attention to the dynamics of M1 and M2 cells near the tumour interface, and their interactions with the tumour and with the ECM. We also focus on the M1→M2 tumour-induced polarisation (in the presence/absence of nutrients) and on the M2→M1 re-polarisation as induced by various treatment protocols (e.g., the agonist anti-CD40 mAb mentioned above [15]). In this new study we also, model the potential role of nutrients on the proliferation/death of M1 and M2 macrophages and cancer cells.
With the help of this new extended multi-scale model, we aim to investigate whether an effective macrophage M2→M1 re-polarisation strategy (within a fibrous tissue environment which impacts macrophage dynamics) is spatial and/or temporal dependent. To our knowledge, this question has never been addressed experimentally, even if the spatial heterogeneity of the TME and the ECM is known to impact the polarisation/repolarisation of macrophages, and therefore any potential antitumour treatment involving re-polarisation protocols.
The structure of the paper is as follows. In Section 2.1 we present the extended macro-scale dynamics. Then, in Section 2.2 we outline the two micro-scale processes. We present our numerical simulations in Section 4, by focusing especially on the spatial and temporal dependency of the macrophage repolarisation. At last, we summarise and discuss the results in Section 5.

MULTI-SCALE MODELLING OF THE TUMOUR DYNAMICS
Building on the multi-scale moving boundary framework initially introduced in Trucu et al. [32] and later expanded in Shuttleworth and Truc [29] and Suveges et al. [30], in this work, we explore not only the M2-like macrophages [30] but also the M1 phenotype, by assessing the impact that these bring within the interlinked tissue-scale (macro-scale) and cell-scale (micro-scale) tumour dynamics. Moreover, the biological context of the cancer macro-dynamics is further broadened by considering the presence of the nutrients, such as glucose and oxygen, which are key constituents of the tumour microenvirnment and play an impotan role within overall tumour progression.

Macro-Scale Dynamics
As this work extends the modelling framework introduced in Shuttleworth and Trucu [29], Suveges et al. [30], and Trucu et al. [32], we start this section by introducing some of its key features. Thus, on the macro-scale we focus on the expanding tumour region (t) that progresses within a maximal tissue cube Y⊂R d , for d = 2, 3 and over the time interval [0, T] (i.e., (t) ⊂ Y, ∀t ∈ [0, T]). In this context, at any macroscale spatio-temporal point (x, t) ∈ (t) × [0, T], we consider a mixed cell population consisting of distributions of: (a) cancer cells c(x, t); (b) M1-like macrophages, M 1 (x, t), briefly addressed here as M1 TAM; and (c) M2-like macrophages, M 2 (x, t), which are briefly referred to as M2 TAM. This mixture of cancer cells and macrophages exercise their naturally multiscale dynamics within an extracellular matrix, which, as in Shuttleworth and Trucu [29,41,42], is regarded as consisting of two major phases, namely a fibrous and a non-fibrous one. Specifically, on the one hand, we have the fibre ECM phase, accounting for all major fibrous proteins (such as collagen and fibronectin), whose micro-scale structure enables a spatial bias for withstanding incoming spatial cell fluxes, inducing this way an intrinsic ECM fibres spatial orientation [29,30,41,42]. Therefore, the spatiotemporal distribution of the oriented ECM fibres at the macroscale point (x, t) is described by a vector field θ f (x, t), where representing the amount of fibres at (x, t). Then, on the other hand, besides these fibrous proteins, the ECM also contains many other components such as non-fibrous proteins (for instance amyloid fibrils), enzymes, polysaccharides and extracellular Ca 2+ ions. Hence, in the second ECM phase, we bundle together these constituents and refer to it as the non-fibre ECM phase, and its distribution at each (x, t) ∈ (t) × [0, T] is denoted by l(x, t). Therefore, for compactness, we denote the global fivedimensional tumour vector by u that is given by as well as denoting the total space occupied at (x, t) by ρ(u) and define it as for all (x, t) ∈ (t) × [0, T]. Furthermore, the last component of the macro-scale dynamics is the nutrients density σ (x, t), whose level within the tumour microenvironment is depleted by the invading cancer. Finally, the macrophages are considered to infiltrate the tumour through the outer boundary [30], which is denoted by ∂ o (t) ⊂ ∂ (t) (see Figure 2), and is defined in Appendix B.

Nutrient
In this study we focus on an avascular tumour mass, whose growth is supported by an influx of nutrients that diffuse through the outer boundary of the tumour. To incorporate this aspect into our mathematical model, we consider the overall influx of a generic nutrient through the outer tumour boundary ∂ o (t) by using Dirichlet boundary conditions. Although different cell types uptake the supplied nutrients at different rates, for simplicity here we assume that all present cell populations (cancer cells, M1 and M2 TAMs) uptake nutrients at the same constant rate d σ > 0. The spatial transport of the nutrients is modelled by a diffusion term with constant coefficient D σ > 0. Since this diffusion occurs more rapidly than cell diffusion (i.e., cell random walk), we use a quasi steady state reaction-diffusion equation (similar to the one for instance in Macklin et al. [27]) for the generic nutrients σ (x, t): Here, σ nor is the normal level of nutrients along the tumour interface, which is considered here to be a constant.
Since cells require nutrients to function properly, here we introduce four smooth bounded effect-functions that we use to model the effects of the nutrients on the different cell functions that considers the following critical nutrients levels: • the necrotic minimal threshold σ n > 0 (i.e., if σ ≤ σ n this leads to necrotic tumour cell death in that area [43]; • the intrinsic nutrient level sufficient for cells to function properly: σ p > 0; • the normal level of nutrients: σ nor (threshold value that is used also as the Dirichlet boundary condition in Equation 3).
Hence, we have the following relationship between these three values, namely: σ nor > σ p > σ n . Starting with the effect of nutrients on cell proliferation, we first assume that this effect is relatively similar for cancer cells and both macrophage populations, in the sense that very low nutrient level impedes cell proliferation and extremely high nutrient level cannot increase cell proliferation rate above a certain maximum. We consider a maximal proliferation enhancement rate p,max > 0 which corresponds to nutrient levels σ ≥ σ p . Also, we assume no proliferation in the necrotic regions. Thus, we define the proliferation effect-function: where (σ , ·, ·, ·) is a generic cosine function that describes a smooth transition between the two extrema, i.e., for any level σ n < σ < σ p . Thus, (σ , ·, ·, ·) is defined by where min is the minimum, max is the maximum and L controls the phase shift of the cosine function. We illustrate the proliferation effect-function, defined in Equation (4) in Figure 1.
As opposed to cell proliferation, we distinguish between the death of cancer cells and the death of macrophages, since cancer cells resist death [3] while macrophages do not. Moreover, experimental studies have shown that macrophages can die in the absence of nutrients [44]. Also, it is known that following hypoxia-induced chemo-attractant signals, macrophages position themselves into pre-necrotic regions to clear out dead cells in the necrotic regions [45]. Therefore, we consider a maximal enhancement death rate d,max > 0 in necrotic regions for both populations, and while we assume no death for cancer cells, we consider a minimal level of macrophage enhancement death rate dM,min > 0 in regions where the nutrient level is σ ≥ σ p . Hence, cancer cell death and macrophage death effect-functions are defined as (using the transition function Equation 5): These death effect-functions are illustrated in Figure 1.
Finally, it was shown experimentally [46] that the nutrient level affects the macrophages polarisation rate giving rise for instance to hypoxia-inducible factors generated by the link with the TME. Therefore, here we consider also a polarisation effectfunction by which we ensure a smooth transition between levels of maximal enhancement rate M,max > 0 in necrotic regions and levels of minimal enhancement rate M,min > 0 in regions with normal nutrient amounts. Using the smooth transition function (Equation 5), we define the polarisation effect-function This function is illustrated in Figure 1.

Dynamics of both M1 and M2 TAMs
Focusing now on the macrophage population, specifically to the two extreme phenotypes, M1 and M2 TAMs, in this paper we focus on the dynamics of the two TAMs populations exclusively inside the tumour domain (t), the evolution of the macrophages distribution in the surrounding tumour stroma being beyond the scope of this current work. Hence, since macrophages are recruited to tumour sites as an immune response through the peritumoral vasculature, here this influx is represented by a source of M1 TAMs that is localised along the outer tumour boundary ∂ o (t), which is enabled by the immediate activation of macrophages into M1 TAMs as they enter the tumour [30]. For simplicity, we assume that both profile and maximal magnitude M 0 > 0 of this source are identical along the tumour interface, and so this influx term is given by Here, ψ ρ is the standard mollifier with appropriately chosen range ρ > 0, " * " is the convolution operator [47] and χ ∂ o (t) is the characteristic function of the outer boundary ∂ o (t).
Experimental studies have shown that resident macrophages can proliferate in situ, which can lead to the accumulation of tumour-associated macrophages inside various tumours [48]. Moreover, as in Suveges et al. [30], we recognise that the stiffness of the ECM plays a role in the proliferation of the macrophages [49] and that according to several biological studies [50][51][52], cancer cells trigger the proliferation of macrophages by producing survival and proliferation factors. Hence, we denote by µ MF > 0 the macrophages proliferation enhancement rate due to the fibres while the proliferation effect-function p (σ ) defined in Equation (4) accounts for the effect of nutrients on the macrophage proliferation. For simplicity, both the fibre enhancement rate and the effect-function are assumed to remain unchanged for both phenotype. Thus, we formulate the proliferation laws for the M1 and M2 TAMs as respectively. In Equation (9), µ M is the positive baseline proliferation rate while the term (1−ρ(u)) + : = max(0, 1−ρ(u)) ensures that there is no overcrowding.
Experimental studies have shown that macrophages' death can be induced by nutritional starvation [44]. Thus, for both M1 and M2 TAMs, we consider a natural death rate d M > 0 that is regulated by the available nutrients through the death effectfunction dM (σ ) introduced in Equation (6b), and so the death terms for each of the two phenotypes are defined by for the M1 and M2 TAMs populations, respectively.
Due to the versatility of the macrophages, their phenotype can be switched from one to another [53,54]. In the present work, we focus on two factors that drive the polarisation of M1 TAMs into M2 TAMs, which are detailed as follows. On the one hand, cytokines secreted by the cancer cells were shown [55,56] to trigger the polarisation process. On the other hand, the nutrient level was also shown [46] to affect this process. As a consequence, we describe the polarisation of M1 TAMs to M2 TAMs by where p 12 > 0 is a constant proliferation rate, and M is the polarisation effect-function defined in Equation (7). Further, in vitro, it has been demonstrated [57] that the M2-like macrophages can be re-polarised back into the M1 phenotype which may be a viable strategy against tumour development.
To that end, we explore here mathematically the possibilities of the re-polarisation strategy through a re-polarisation term of the form Here, p 21 > 0 is the constant re-polarisation rate, t p > 0 is the activation time and χ p(t,Rp) is the characteristic function of the repolarisation domain p (t, R p ) that is defined in Appendix C and illustrated in Figure 2. This re-polarisation term Equation (12) allows us to examine whether or not we would need to account for spatio-temporal dependencies through the domain p (t, R p ) and activation times t p > 0 in order to obtain an effective re-polarisation strategy. The motility of both macrophages phenotypes is driven both by random and directed movement. Based on recent biological evidence [49], increased stiffness of the substrate leads to an increase in macrophages' speed, aspect explored in our modelling through a diffusion enhancement that corresponds to with the level of ECM fibres. To that end, we consider a stiffnessdependent macrophage diffusion coefficient D M (u) of the form where D M > 0 is the baseline macrophage diffusion rate, and D MF > 0 is the diffusion enhancement rate due to the presence of fibres. On the other hand, besides random movement, macrophages also exercise directed migration due to a both adhesive interactions with the surrounding cells and the ECM as well as an underlying cross-talk between themselves and the cancer cells. A similar "non-local flux term" to the one introduced in Suveges et al. [30] is used here to explore the complex interactions of the cells distributed at x ∈ (t 0 ) with other cells within a sensing region B(0, R), and this accounts for: (1) cell-cell TAMs self-adhesion [58]; (2) nutrients level mediated movement [59]; and (3) the contribution of the cancer cells to the directional movement of the macrophages [60][61][62][63]. Specifically, the contribution of the cancer cells to the directional movement of the macrophages account not only for the biological evidence that cancer cells can bind themselves to TAMs [60] but also for the fact that cancer cells can attract TAMs [61][62][63] by secreting various chemokines. Although we neither model explicitly the involved chemokine activities within this cross-talk nor the chemo-attractant activities involved with the nutrients, here we appropriately account for both of them through the following non-local flux term: where R represents the radius of the sensing region B(0, R).  (14) for the gradual weakening of these different adhesions as we move away from the centre x within B(x, R), we use a radially symmetric kernel K(·) that is given by where ψ(·) is the standard mollifier. Moreover, in Equation (14), [1 − ρ(u) + ensures that overcrowded tumour regions do not contribute to macrophage migration and n(·) is the unit radial vector given by Thus, aggregating now all these cell movement aspects explored in Equations (8)- (14), the dynamics of the two distinct macrophages phenotypes are mathematically formulated as where S M 1 M > 0 and S M 2 M > 0 are the self-adhesion strengths of M1 and M2 TAMs, respectively.

Dynamics of the Cancer Cell Population
The third cell population that we consider at macro-scale is the cancer cell population. Crucially important for cancer development and invasion, the cancer cell proliferation is a complex process that is regulated by several processes involving nutrients and macrophages. From the modelling perspective, while we consider the proliferation process as being of logistic type [64][65][66], we explore the influence of nutrients and macrophages as follows. On the one hand, similar to both TAMs populations, we consider the proliferation effect-function p (σ ) defined in Equation (4) to explore the influence of the available nutrients on the rate of cancer cell proliferation. On the other hand, biological evidence shows that while M2 TAMs promote cancer cell proliferation [67], M1 TAMs inhibits this [68]. Thus, expanding here the proliferation law introduced in Suveges et al. [30] by accounting for the negative effect of M1 TAMs, we obtain leading to the following proliferation law: where µ c > 0 is a baseline proliferation rate that is being regulated by the available nutrients, being enhanced by the M2 TAMs at a rate µ cM 2 > 0 and at the same time weakened by the presence of the M1 TAMs at a rate µ cM 1 > 0. Again, here the term (1 − ρ(u)) + ensures that there is no overcrowding. Besides proliferation, it is well known that cancer cells resist death [3,69]. However, due to the peritumoral vasculature as well as the excessive degradation of the ECM, the efficiency of the nutrients delivery significantly reduces inside the tumour, leading to necrosis [70]. In addition, numerous studies have shown [71][72][73][74][75] that classically activated M1-like macrophages can produce significant amounts of pro-inflammatory cytokines, and thereby have the ability to kill cancer cells. To that end, we assume here a baseline death rate d c > 0 that is regulated not only by the cancer cell death effect-function dc (σ ) introduced in Equation (6a), but also by the M1 TAMs at a rate d cM 1 > 0. This results in the following mathematical representation of the cancer cell death process, namely Similar to the macrophages, for the cancer cell population we also account for the diffusion enhancement that the spatial distribution of ECM fibres enables [76][77][78][79][80][81][82][83][84]. Furthermore, the random movement of the cell population is also affected by the presence of both macrophage populations. While in general, the M2 TAMs were shown to promote cancer cell motility, Afik et al. [85] recent biological evidence [68,86] indicates that the M1 phenotype has a negative effect on the cancer cell motility. Therefore, the diffusion coefficient for the random movement of the cancer cells can be formulated mathematically as where D c > 0 is a baseline diffusion rate, D cF > 0 is the ECM fibres enhancement coefficient, D cM 1 > 0 represents the weakening effect due to the presence of M1 TAM, and D cM 2 > 0 accounts for the positive motility effect due to the presence of M2 TAM. Besides random motility, the directed movement of the cancer cells induced by various adhesion mediated processes [60,61,[87][88][89][90] is a central player in cancer invasion within the oriented fibrous environment. To that end, extending here on the modelling approach proposed in Suveges et al. [30] to include the interactions of cancer cells with both families of macrophages, i.e., M1 and M2 TAM, we have that the non-local spatial flux that drives the directed movement is given in this case as: where R, n(·) and K(·) are the same as in Equation (14). Further, in Equation (19) n(·, ·) is the unit radial vector biased by the orientation of the fibres, i.e., as illustrated in Figure 3. Moreover, in Equation (19) S cM > 0 represents the strength of the adhesion relationship between the cancer cells and M1 and M2 TAMs, S cF > 0 is the strength of the cell-fibre ECM adhesion [91] and S cl > 0 corresponds to strength of adhesion between the cancer cells and the non-fibre ECM phase (that includes for instance amyloid fibrils, which can support cell-adhesion processes [92]). Furthrmore,as high level of extracellular Ca +2 ions (which form one of the constituents of the non-fibre ECM phase) are necessary for cell-cell adhesion [93,94], proceeding as in Shuttleworth and Trucu [29,41,42], and Suveges et al. [30] the cancer cells self-adhesion coefficient S cc is taken here as where S max > 0 and S min > 0 correspond to maximum and minimum levels of Ca +2 ions. Therefore, S cc smoothly increases from a minimal to a maximum value in order to fully explore the varying strengths of cell-cell adhesion.
Thus, using Equations (16)-(19) the spatio-temporal dynamics of the cancer population c(x, t) is expressed as

Two-Phase ECM Macro-Scale Dynamics
Besides the cancer cells, both macrophage phenotypes contribute to the degradation of the ECM by secreting proteolytic enzymes [95][96][97][98][99] (e.g., various classes of matrix metalloproteinases). To that end, we extend the dynamics of the fibre, and non-fibre ECM components used in Suveges et al. [30] by incorporating the effects of the M1 phenotype. Thus, the dynamics of the non-fibre l(x, t) as well as the fibre ECM F(x, t) are formalised as where β lc , β lM 1 , β lM 2 are the positive degradation rates of the non-fibre ECM phase due to the cancer cells, M1 and M2 TAMs, respectively. Similarly, β Fc , β FM 1 , β FM 2 are all positive and describe the degradation rates of the fibre component of the ECM due to the cancer cells, M1 and M2 TAMs, respectively. Finally, in Equation (22)

The Full Macro-Scale Dynamics
In summary, using Equations (3), (15), (21), and (22) the nondimensional macro-scale dynamics is given by the following coupled PDEs in the presence of appropriate initial conditions (such as those specified in Equation (51)) with zero-flux boundary conditions for c, M 1 , M 2 , l and F, as well as Dirichlet boundary condition (Equation 3) for the nutrients σ .

Processes on the Micro-Scales and Links Between the Scales
As the process of cancer invasion is truly a multi-scale phenomena [2], the macro-scale dynamics is tightly linked together with several micro-scale processes. Among the microscale processes of important for cancer invasion, of main interest for us in this work are the micro-scale rearrangement of ECM fibre micro-constituents as well as the cell-scale proteolytic processes that take place at the leading edge of the tumour. In the following, we outline the details of these two micro-dynamics as well as the naturally occurring double feedback loop that links them to the tumour macro-scale dynamics (Equation 23).

Fibres on the Micro-Scale and Their Bottom-Up and Top-Down Links to Macro-Dynamics
Following [29], the macroscopic oriented ECM fibres are represented not only through their amount F(x, t), but also via their spatial bias that characterise their ability to withstand incoming forces. Indeed, both of these characteristics of the ECM fibres are induced by the distribution of the non-dimensional micro-fibres f (z, t) on a cell-scale micro-domain δY(x) : = x+δY of appropriate micro-scale size δ > 0, and are captured via a vector field representation [29] θ f (x, t) of the ECM fibres which is defined by where λ(·) is the Lebesgue measure in R d and θ f ,δY(x) (·, ·) is the revolving barycentral orientation with respect to the measure f (z, t)λ(·) given by Shuttleworth and Trucu [29] Finally, the second characteristic, namely the ECM fibres amount, is given by the Euclidean norm of the vector field θ f (x, t), i.e., and so it describes the mean-value of the micro-fibres distributed on δY(x). Therefore, the micro-scale naturally links with the macro-scale since the representation of the ECM fibres on the macro-scale is induced by the micro-scale fibre distribution, and so we refer to this as the fibres bottom-up link.
On the other hand, the macro-scale spatial fluxes, generated by the tumour macro-dynamics (Equation 23), trigger a rearrangement of the ECM fibres micro constituents on each micro-domain δY(x). Indeed, the collective migration of the cancer cells, M1 TAMS, and M2 TAMs lead naturally to the emergence of the associated spatial fluxes F c , F M 1 , and F M 2 given by respectively. The combined action of these fluxes upon the ECM fibres distributed at (x, t) ∈ t × [0, T] is felt uniformly by its constituent micro-fibres f (z, t) distributed on the associated micro-domain δY(x), consequently inducing a micro-fibres rearrangement vector similar to the ones proposed in Shuttleworth and Trucu [29] and Suveges et al. [30] of the form which triggers a spatial redistribution of the micro-fibres in δY(x). In Equation (25), the non-linear weights ω c , ω M 1 , ω M 2 and ω F are the associated mass factions of cancer cells, M1 TAMs, M2 TAMs, and ECM fibres at (x, t), and so these are , , , respectively. Ultimately, the rearrangement vector (Equation 25) induces a relocation vector ν δY(x) (z, t), and as a result we can appropriately calculate the new positions of any micro-scale node z ∈ δY(x) that are given by In Figure 4, we illustrate the micro-fibre rearrangement process by considering a typical example of these vectors while for further details we refer the reader to Shuttleworth and Trucu [29] and Suveges et al. [30]. Finally, this rearrangement of the ECM fibres at micro-scale triggered by the emergent macro-scale spatial fluxes (F c , F M 1 , and F M 2 ) establishes the fibre top-bottom link. In Figure 5, we illustrate the fibres bottom-up and top-down links, connecting the macro-scale and the ECM fibre micro-scale.

MDE Boundary Micro-Dynamics and Its Connections to Macro-Dynamics
The second micro-process that we consider is the proteolytic molecular processes which are driven by the secretion of matrixdegrading enzymes (MDEs) (such as matrix-metalloproteinases) and take place along the leading edge of the tumour. Indeed, besides cancer cells, both M1 and M2 TAMs within the proliferating outer rim of the tumour secrete MDEs [95][96][97][98][99], and this way the tissue-scale dynamics induces a source for cellscale proteolytic activity. Once secreted, the MDE molecular population exercises a spatial transport across the tumour interface within a cell scale neighbourhood of the tumour boundary, causing degradation of the szrperitumoral ECM, and ultimately leading to changes in the tumour morphology [2]. Thus, following here the approach introduced in Trucu et al. [32], we explore the emerging spatio-temporal molecular MDEs micro-dynamics on an appropriate cell-scale neighbourhood of the tumour interface ∂ (t) enabled by the union of a covering bundle of cubic micro-domains {ǫY} ǫY∈P(t) , with each ǫY being of micro-scale size ǫ > 0. This allows us to decompose the MDE micro-dynamics occurring on ǫY into a union of MDE micro-processes occurring on each ǫY [32].
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org  Therefore, considering the tumour evolution over a time perspective [t 0 , t 0 + t], for an arbitrary instance t 0 ∈ [0, T], and of appropriate micro-scale range t > 0, on any of the microdomains ǫY ∈ P(t 0 ) we denote by m(y, τ ) the spatio-temporal distribution of MDEs at micro-scale point (y, τ ) ∈ ǫY × [0, t]. In this context, at any spatio-temporal (y, τ ) ∈ (ǫY ∩ (t 0 )) × [0, t], a source of MDEs arises as a collective contribution of the cancer cell and both macrophage populations from the outer proliferating rim of the tumour that are situated within a distance γ h > 0 from y ∈ ǫY. Hence, denoting this micro-scale MDE source by h(y, τ ), this can be formalised mathematically via the non-local expression where B(y, γ h ) : = {z ∈ Y| y − z ∞ ≤ γ h } denotes the · ∞ ball with appropriately chosen radius γ h > 0 and 0 < ρ < γ h is a small mollification range which smooths out the source function h(·, ·). Further, in Equation (26) h is given by where α c > 0, α M 1 > 0 and α M 2 > 0 are constant secretion rates of the MDEs by the cancer cells, M1 and M2 TAMs respectively. As the MDE micro-source is naturally induced by the macroscale, this establishes a MDE top-down link between the tumour macro-dynamics and MDE-micro-dynamics occurring at the tumour interface. Finally, under the presence of the MDE source h(·, ·), the MDE micro-dynamics is given by where D m > 0 is the constant diffusion coefficient and n is the outward normal vector. Finally, as the patterns of significant degradation of the peritumoural ECM correspond the regions where significant levels MDEs are transported during the microdynamics, the micro-dynamics (Equation 27) enables us to capture to changes in tumour morphology by determining the direction of movement and magnitude of the displacement of invading cancer in the surrounding tissue [32]. Thus, the MDE micro-process induces changes in the shape of the tumour boundary ∂ (t 0 ), affecting directly the macro-dynamics, which is continued on a modified tumour macro-domain (t 0 + t), for further details we refer the reader to Trucu et al. [32]. Hence, a MDE bottom-up link between the proteolytic boundary micro-dynamics and macro-scale cancer dynamics is this way established. In Figure 5, we illustrate both MDE top-down and bottom-up links that connect the macro-scale and the MDE micro-scale.

NUMERICAL APPROACHES
In this section, we describe the numerical approach developed to solve the tumour macro-scale dynamics (Equation 23). First, to solve the quasi-steady nutrients σ equation (23f), we use the usual successive over-relaxation method with relaxation parameter ω = 0.5 and tolerance of 10 −5 . For the rest of the dynamics (Equations 23a-23e), we use the method of lines approach to discretise the system (Equation 23) first in space and then for the resulting ODEs, we use a non-local predictorcorrector scheme [29]. In this context, we carry out an accurate approximation of the two distinct spatial operators, namely the diffusion and adhesion operators, by using fast convolutiondriven approaches. Specifically, while for the diffusion, we construct a convolution-based second-order central difference scheme [30], for the adhesion operators, we formulate a convolution-driven fifth-order weighted essentially non-oscillatory (WENO5) finite difference scheme. Furthermore, to approximate the adhesion integrals (Equation 14 and 19), we adopt the numerical approach used in Suveges et al. [30]. Here, we omit the technical details of this method since it follows identical steps, the sole difference being that within the formulas, we use the total macrophage population, i.e., M1 + M2 TAMs instead of only the M2 TAMs population. For the same rationale, we also skip the details of the numerical method used in the MDE micro-scale, and for further details, we refer the reader to Suveges et al. [30].
To facilitate the description of the numerical approaches, let us introduce some basic notations. First, we discretise the maximal tissue cube Y ∈ R 2 uniformly in each direction with spatial step size x = y. Therefore, the discretisation of Y can be represented through the macro-spatial locations {(x i , y j )} i,j=1..N , with N = L/h L + 1. Then the discretised global tumour vector at any time-step n is denoted by u n = [c n , M n 1 , M n 2 , l n , F n ] ⊺ , while the discretised nutrient field is denoted by σ n . Furthermore, we denote the diffusion coefficient functions and adhesion integrals for cancer cell and TAMs by (D c ) n , (D M ) n , (A c ) n and (A M ) n , respectively. Finally, with these notations, we are able to rewrite the macro-scale dynamics (Equation 23) in the following convenient form where

Diffusion Operators: ∇ · [D(u)∇u]
Starting with the discretisation of the diffusion operators ∇ · [D(u)∇u] in Equation (28), during the computations, we first detect whether a spatial node (i, j) is inside or outside the expanding tumour domain (t 0 ) via an indicator function I(·, ·) : X × X → {0, 1}, with X = {1, ..., M} that is defined by Similarly, we detect any boundary nodes using a boundary indicator function defined by where ∂ (t 0 ) is boundary of (t 0 ), and for convenience we also define the interior indicator function by These two functions, given in Equations (30) and (31), enable us to split the domain into two parts, namely to boundary and strictly inside parts. The motivation behind this is to use two different computational technique on these parts which eventually reduces the computational cost. Hence, while for any interior node we can use the same discrete universal numerical operator (convolution), for a boundary point we need to apply unique operators due to the zero-flux boundary condition and the continuously changing tumour domain (that may result in an irregular tumour shape). Therefore, to achieve the reduction in computational cost, we accordingly split the spatial operators ∇ · [D(u)∇u] into two components as well, and so at any spatial node (x i , y j ) ∈ (t 0 ) ⊂ Y, the diffusion operators can be represented as otherwise. (32) In this context, the usual two dimensional second-order central difference scheme for a non-constant diffusion operator is given by where u n i,j = [c n i,j , M n 1 i,j , M n 2 i,j , l n i,j , F n i,j ] ⊺ . Further, we observe that Equation (33) can be equivalently expressed by sum of discrete convolutions and so for all interior node (x i , y j ) the scheme is given by where * is the discrete convolution and • is the Hadamard product [101]. Also, in Equation (34) each K k A , with k = 1 . . . 4 describe the average between the point (i, j) and one of its immediate neighbour and so they are defined by Moreover, in Equation (34) K k F and K k B with k = 1, 2 are induced by the forward and backward differences, respectively. Hence, they are defined in both direction i (if k = 1) as well as in direction j (if k = 2) by However, for boundary nodes, we cannot use the form Equation (34) due to the imposed zero-flux boundary condition and the continuously changing tumour domain. This is because the calculation of the diffusion operators at the boundary nodes involves the approximation of the values at any node that does not belong to the tumour domain, i.e., for any node (x i , y j ) / ∈ (t 0 ). In some cases, due to the irregular domain, such values may not be uniquely defined because multiple nodes can require the value of the same ghost/outside node, however, with different values. Consequently, for any boundary node, instead of the convolutional form Equation (34), we instead use the usual form Equation (33) and symmetrically reflect the values of the interior nodes to the ghost/outside nodes in a node by node fashion.

Adhesion Operators: ∇ · F (u)
The other spatial operators that contribute to the motility are the adhesion operators ∇ · F (u) in Equation (28). The procedure to approximate these differential operators is based on the standard WENO5 scheme proposed in Jiang and Shu [102] and Liu et al. [103]. However, it was shown [104] that these standard WENO5 schemes suffer from slight post-shock oscillations. As a consequence, we adopt here the modified WENO5 scheme proposed in Zhang and Shu [104], which uses modified smoothness indicators and is referred to as the ZSWENO scheme.
Similarly to the diffusion operator case, here we split the domain into two parts as well. However, ZSWENO schemes induce larger stencils compared to the second-order central difference scheme and so, we need to split the domain (t 0 ) into two parts differently. We refer these two parts as the inside and layer parts. First, the former one is detected by using an inside indicator function I I (i, j) that we define by where I is defined in Equation (29) and K I is given by which is induced by the fifth-order ZSWENO stencils. Then the layer part of the domain is detected by a layer indicator function formulated as Similarly to the diffusion, using these indicator functions (Equations 37 and 38), we split the adhesion operator into two parts and so at any spatial node (x i , y j ) ∈ (t 0 ) this operator is represented as Frontiers in Applied Mathematics and Statistics | www.frontiersin.org Therefore, for the inside operator ∇ · F (u) I , we use the usual conservative form whereF i+ 1 2 ,j ,F i− 1 2 ,j ,Ĝ i,j+ 1 2 andĜ i,j− 1 2 are the numerical fluxes at (x i+ 1 2 , y j ), (x i− 1 2 , y j ), (x i , y j+ 1   2   ) and (x i , y j− 1 2 ), respectively. Also in Equation (40), for compact notation,F andĜ denotes the x and y components of the vector field F (·), respectively.
For brevity, we will only focus on definingF i+ 1 2 ,j andF i− 1 2 ,j , and note that corresponding calculations forĜ i,j+ 1 2 andĜ i,j− 1 2 follows identical steps. In this context, we split the x component of F (·) into two partsF where dF + (u)/du > 0 and dF − (u)/du ≤ 0. Then, we define these parts by using the popular Rusanov-type flux splitting method [105] that is given bŷ where we approximate α : = max |dF(u)/du|, i.e., the spectral radius of the Jacobian generated byF(·) in a Jacobian-free manner, detailed in Section 3.3. To this end, let us denote bŷ F +  [102][103][104]106]. Hence, they are given bŷ where ω ± k,i± 1 2 ,j are the non-linear weights that we will define later, F ± k,+ are the ENO approximations given bŷ and similarly,F ± k,− are defined bŷ Finally, in Equations (43) and (44)F ± (u) are given by the Rusanov-type flux splitting method defined in Equation (41).
Since our aim is to reduce the computational cost and so for this we seek to use convolution, we observe that indeed these ENO fluxes (Equations 43 and 44) can be equivalently represented by discrete convolutions, i.e., where K ± k,+ and K ± k,− are the induced vectors from Equations (43) and (44), and for completeness they are defined in Appendix D.
Let us now shift our attention to the non-linear weights ω ± k,i± 1 2 ,j that we used to construct the ZSWENO approximation in Equation (42). Following again the usual procedure [103], we define these weights as where we take the usual p = 2, ǫ W = 10 −6 values [102], and define the optimal weights ω o k,+ and ω o k,− as: are the ZSWENO smoothness indicators [104], and the explicit formulae first for IS ± k,i+ 1 2 ,j are given by and then for IS ± k,i− 1 2 ,j , we have Once again, in Equations (47) and (48) the positiveF + (u) and negativeF − (u) parts are defined by using the Rusanov flux splitting Equation (41). Furthermore, we can observe that these smoothness indicators (Equations (47), 48) can be equivalently expressed in terms of discrete convolutions, i.e., where the appropriately induced vectors K ± k,+ and K ± k,− are defined in Appendix D. This completes the description of the convolution-driven ZSWENO scheme for the inside differential operators ∇ · F (u) I , i.e., the approximation of ∇ · F (u) for all inside nodes.
However, since the adhesion operators ∇ · F (u) were split into two parts in Equation (39), it remains to describe the ZSWENO scheme for the layer operator ∇ · F (u) L , accounting for all nodes that are considered to be in the layer part of the domain. As in Section 3.1 for the boundary diffusion operators, here we also need to appropriately approximate the value of any point that is located outside of the tumour domain (t 0 ). To this end, we symmetrically reflect the values of the interior nodes to the ghost nodes in a node by node fashion due to the irregular tumour domain. This ultimately enables us to use the standard nonconvolutional ZSWENO scheme (Equations 43, 44, 47, and 48) instead of the convolutional forms Equations (45) and (49) to approximate the layer operators ∇ · F (u) L .

Approximation of the Propagation Speed α
Due to the complexity of F (·), calculating its Jacobian dF(u)/du is extremely time consuming. Hence, to find the largest eigenvalue α, we rather skip the exact evaluation of the Jacobian and choose to approximate the propagation speed in a Jacobianfree manner. For this, let us first define the eigenvalue problem by where J = dF(u)/du is the Jacobian, v is an eigenvector, and λ is an eigenvalue of J. Then in order to find the largest eigenvalue α, we solve (Equation 50) iteratively with convergence tolerance of 10 −14 and with a random initial guess for v. Hence, similarly to, for instance, a Jacobian-free Newton-Krylov method [107, 108], to find λv in Equation (50) we evaluate the Jacobian-vector product Jv instead of J, which is proved to be a significantly less time-consuming task. To that end, the approximation of Jv is carried out via the first-order Taylor series expansion, and so it is given by where ǫ p is a small perturbation parameter. Since the precision is limited in the evaluation of the flux F (·), a good choice to evaluate this small parameter ǫ p is given by Knoll and Keyes [108] where ǫ mach > 0 is the machine precision.

NUMERICAL RESULTS
In this section, we present the numerical results for our multi- l(x, 0) = min 1 2 which are illustrated in Figure 6A. Here the white curves indicate the boundary of the tumour domain ∂ (0). Besides these macroscale initial conditions, we also illustrate the initial state of one micro-scale fibre domain δY(x) in Figure 6B, which is repeated for all macro-scale points. Also, the ratio between the fibre and non-fibre ECM phases are assumed to be 20% : 80% for all simulations. Finally, all presented simulations corresponds to time 50 t. The baseline parameters values are provided in Appendix A, and any alteration from these values will be stated accordingly.

Spatial Dependency of the Re-polarisation
First, we investigate numerically the effects of changing the repolarisation domain p (t, R p ) used in Equation (12), defined in Appendix C and illustrated in Figure 2. Hence, here we study whether the success of a M2→M1 re-polarisation strategy against the tumour is dependent on the spatial domain p (t, R p ), specifically on R p =the distance from the outer boundary ∂ o (t). For this let us use the radii R p ∈ {0, x, 2 x, 3 x, 4 x} for the re-polarisation domain p (t, R p ), as well as let us start the process at time 0, i.e., we take t p = 0 in Equation (12). To further study these spatial effects, we also consider multiple tissue conditions by changing the tissue environment controller β ∈ {0.75, 0.7875, 0.825} (see Trucu et al. [32]).
To compare the resulting tumours, we measure the overall tumour mass and spread at final time 50 t, using the total density of the cancer cells as well as the area of the tumour. Ultimately, this enables us to quantify the changes, resulted from the modification of the re-polarisation domain p (t, R p ). Specifically, the overall tumour mass is measured by integrating the cancer cell density c(x, 50 t) and the overall tumour spread as given by the area of the tumour domain (50 t).

Baseline Cases
In Figure 7 we present the numerical results for the baseline case characterised by the absence of M2→M1 repolarisation. Specifically, in Figures 7A-C we investigate the no-repolarisation case in the context of different tissue environment controllers: Figure 7A β = 0.75, Figure 7B β = 0.7875, and Figure 7C β = 0.825. As we can see, at time 50 t, the tumour is enlarged and has invaded some of its surroundings, and within the tumour domain (50 t) we observed heterogeneity in all three cell populations (cancer cell, M1 and M2 TAMs). We also notice that near M2 TAMs accumulations sites, the density of cancer cells also seems to be higher. In contrast, higher M1 TAMs density usually corresponds to both low M2 TAMs as well as low cancer cell populations. The movement and spatial distribution of tumour and immune cells are influenced directly and indirectly by the ECM fibres. For illustrative purposes, the rearranged fibre structure is portrayed by a four-fold coarsened oriented ECM fibres field in Figure 7. The peritumoral degradation of the two-phase ECM (caused by the cancer cells and both TAMs) allows the tumour to expand and spread to the neighbouring tissues, resulting in some tumour fingering patterns that vary with the controller β: higher β leads to more tumour fingering. This induces an irregular tumour domain, as well as the formation of "islands" inside the tumour, which corresponds to areas of initially low ECM density, i.e., the ECM level was too low in such areas to support tumour movement. Finally, in Figure 7, we also present the level of available nutrients. Hence, we can see that since the nutrients are only supplied through the outer boundary o (t), the initial normal level of nutrients is significantly depleted (by the three cell populations) inside the tumour. This can lead to hypoxia and then eventually create a necrotic core (not modelled explicitly in this study).

The Impact of M2→M1 re-Polarisation
In Figures 8A-C Figure 8 we see that the re-polarisation of M2 TAMs leads to an increase in the M1 TAMs population inside the tumour, leading to several accumulation sites located further away from the leading edge. In contrast, the presence of M2 TAMs inside of the tumour is decreased compared to Figure 7, and these immune cells now accumulate only in a small neighbourhood of the boundary (because we do not re-polarise M2 TAMs into M1 TAMs in a R p = x neighbourhood of the boundary; see Figure 2). Since both macrophage populations interact with the cancer cells, we also see some differences in the cancer cell population and these differences depend on the degradation level of ECM as controlled by the value of β. Using two measures (spread and mass), we first focus on the β = 0.75 case and compare Figures 7A, 8A: re-polarisation leads to an  ≈ 11% increase in tumour spread as well as an ≈ 20% decrease in tumour mass. Then for β = 0.7875 (Figures 7B, 8B) we observe an ≈ 5% reduction in spread and an ≈ 34% decrease in mass. Finally, for β = 0.825 (Figures 7C, 8C), the tumour spread is increased by ≈ 15% and the tumour mass is reduced by ≈ 31%. Hence, although re-polarising M2 TAMs into M1 TAMs does not show an overall improvement in terms of tumour spread, it significantly impeded the migration of the cancer cell mass from hypoxic regions to the proliferating rim as well as helped reducing the mass by killing the cancer cells. However, since R p = x, the M2 TAMs are able to accumulate in a x neighbourhood of the interface and so we can still observe a moderate density of cancer cells near the boundary.

Tumour Spread/Mass Changes With Respect to R p
In Figure 9 we vary R p and present the changes in the dimensionless tumour area (spread) and mass resulted by changing the radius R p . Specifically, in Figures 9A-C we again consider the three previously used environment controllers β = 0.75, β = 0.7875, β = 0.825. Moreover, since various types of cancers secrete MMPs at different rates [109], here we also consider high, medium and low MDE secretion rates. Finally, we vary the radii for the re-polarisation domain, R p ∈ {0, x, 2 x, 3 x, 4 x}, while we assume that the re-polarisation process starts at time 0, i.e., we take t p = 0. The left panels of Figure 9 show the changes in the dimensionless variable for tumour areas/spreads with respect to R p , while the right panels show the changes in the dimensionless variables for tumour mass with respect to R p . Here, the blue circles correspond to high ECM degradation rates (β lc = 3.0, β Fc = 1.5, α c = 0.625), the orange diamonds correspond to medium degradation rates (β lc = 2.0, β Fc = 1.0, α c = 0.42), and the red crosses correspond to low degradation rates (β lc = 1.0, β Fc = 0.5, α c = 0.21). Comparing first the changes due to varying the tissue controller β, we can see a clear overall decrease in both tumour spread and mass as we increase β. This was seen also in Figures 7, 8, where a more prominent tumour fingering pattern was present as we increased β which resulted in a decrease in tumour spread. Furthermore, in Figures 9A-C the overall behaviour of the tumour spread does not change significantly as we vary the radius of the re-polarisation domain R p . Hence, even though the proteolytic molecular processes at the leading tumour edge change (via the MDE source Equation 26) following the M2→M1 re-polarisation, we cannot see a clear trend in tumour spread. This might be because of the complex interactions between tumour and infiltrating immune cell populations: although we see an increase in the M1 TAMs populations near the boundary [that also secrete more MDEs than the M2 cells [110]], the overall proteolytic molecular processes at the leading edge of the tumour might not change too much, which would mean similar tumour spread. Second, the rearrangement of the micro-fibre distribution also affects tumour spread, since the amount of fibres that are being re-located near the leading edge is dependent on the fluxes generated by the different cell populations. Therefore, our model suggests that merely re-polarising the M2 TAMs into the anti-tumoral M1 phenotype might not be enough to stop tumour spread.
On the other hand, in Figures 9A-C we see that the tumour mass is greatly reduced in all of the presented cases. These results also show that for the largest reduction we need to re-polarise the M2 TAMs inside the whole tumour domain (i.e., R p = 0, which means that the re-polarisation domain is p (t, 0) = (t)). Hence, the presence of M2 TAMs in the proliferating rim may be enough to maintain the tumour mass and to induce tumour spread by leading to a moderate presence of cancer cells in the proliferating rim. In conclusion, the results in Figure 9 suggest that re-polarising M2-like macrophages to the M1 phenotype have the ability to reduce tumour mass, but we may need a supplementary strategy primarily focusing on tumour spread or tumour stroma in order to restrain tumour development.

Temporal Dependency of the Re-polarisation
Let us now concentrate on the temporal aspects of a possible re-polarisation strategy. Since in the previous section, we have shown that the M2→M1 re-polarisation impacts predominantly tumour mass, and since this effect was the strongest when R p = 0, here we investigate the effect of changing t p on tumour mass when we re-polarise the M2 TAMs within the whole tumour domain. Thus, by varying t p used in Equation (12) we aim to investigate the effectiveness of such strategy when we introduce the re-polarisation of M2 TAMs at time t p > 0. This is crucial since tumours are only detectable above a certain size, and so a potential treatment that uses re-polarisation cannot be started at time t p = 0. To this end, in Figure 10 we present the change in the dimensionless tumour mass with respect to the repolarisation start time t p . As before, Figure 10A corresponds to β = 0.75, Figure 10B to β = 0.8785 and Figure 10C to β = 0.825. For each of these sub-cases, the blue circle corresponds to high ECM degradation rates by cancer cells, the orange diamond to medium degradation rates, and the red cross to low degradation rates. The results in Figure 10 suggest that the lowest tumour mass can be achieved by starting the repolarisation as soon as possible, as one would expect. However, even at moderate times, for instance at 25 t, the tumour mass can be controlled at relatively low values. All these tumour values obtained following re-polarisation at different times t p are smaller than the tumour value obtained for no re-polarisation (see the last point on the horizontal axis in Figure 10). This indicates that a re-polarisation based strategy is not only viable when it is introduced at a very early stage, but it could also be effective even against a more advanced and larger tumour.

In the Absence of Nutrients
So far, we have considered the presence of nutrients and the related effect-functions within the macro-scale dynamics (Equation 23). However, here, we also investigate how the absence of these components influence the overall tumour dynamics. Hence, for the rest of this section, we take σ (x, t) ≡ 0, ∀(x, t) ∈ Y × [0, T] and p (σ ) ≡ 1, dc (σ ) ≡ 1, dM (σ ) ≡ 1 and M (σ ) ≡ 1, ∀σ and S Mσ = 0. For reasons of comparability, we use the same scenarios as before, and so   Figure 11, we present the case of no re-polarisation and no nutrients. Although comparing Figures 7 and 11 we can observe some expected tumour morphology changes, the main difference that we would like to emphasise is that in the absence of nutrients, all three cell populations form a more homogeneous distribution. For example, while in Figure 11 there is only one cancer cell accumulation region (about (2, 2) where we centred our initial condition), in Figure 7 this is not the case, and rather cancer cells are more heterogeneously distributed throughout the tumour domain (50 t) creating various accumulation regions. Furthermore, we can also see an M2 TAMs population that penetrate the tumour deeper in the presence of nutrients. This is because in this case in Figure 11 there are no nutrients mediated movement (since S Mσ = 0), and so the combination of macrophage-macrophage and macrophage-cancer adhesions is not adequate to guide the M2-like macrophages away from the boundaries in most cases. Hence, these results suggest that the nutrient mediated movement play an important, nonnegligible role in guiding macrophages towards necrotic regions, which is consistent with biological experiments [59]. Moreover, in Figure 12, we introduce M2→M1 re-polarisation (with p (t, x) and t p = 0), and so we can compare it to the results of Figure 8 where we presented the same scenario but in the presence of nutrients. By comparing these results, we can observe more homogeneous cell populations (cancer, M1 and M2 TAMs) as well as the same phenomena as before but with the M1 TAMs population (due to the re-polarisation), i.e., in the absence of macrophage-nutrient relationship (since S Mσ = 0) M1-like macrophages are accumulated closer to the boundary and are less efficient in penetrating the tumour. This may be of particular interest for future mathematical models due to recent advances regarding macrophage-mediated drug delivery [111], an approach that can be used not only for tumours but for other diseases as well.

first, in
To investigate the effects of nutrients on the tumour spread and mass behaviour that we have seen in Figure 9, here for comparability, we also present the same scenarios but in the absence of nutrients in Figure 13. Hence, we consider the three previously used environment controllers β = 0.75, β = 0.7875, β = 0.825 with the three ECM degradation rates (high, medium and low) and present both the tumour spread and mass for each case. Although with some expected changes in the values, in Figure 13 we detect a similar behaviour as in Figure 9 that the proteolytic molecular processes at the leading edge of the tumour do not change drastically, i.e., varying the radius of the re-polarisation domain R p does not affect tumour spread substantially. Furthermore, we find consistent results for tumour mass as well, since we achieve the greatest reduction in the case of R p = 0 and then as we increase the radius R p , tumour mass also increases. Therefore, our results suggest that although nutrients play an important role in the movement of both macrophage phenotype, it does not affect the features and properties of the spatial dependency of the re-polarisation, i.e., both in the presence and in the absence of nutrients the tumour spread and mass react to changes in the re-polarisation radius R p in a similar fashion.

CONCLUSIONS
In this study, we have further developed and extended a multiscale moving boundary framework for cancer invasion [29,30,32] by also considering the dynamics of the anti-tumoral M1 TAMs as well as the nutrients. On one hand, we took into account the nutrients since every cell requires them to live and function properly, and so their presence is indispensable. On the other hand, we focused on the classically activated M1like macrophages since they are known to be capable of killing cancer cells. Moreover, since macrophages are one of the most abundant immune cells in the tumour micro-environment and their plasticity enables them to switch phenotypes, they are prime candidates to assist in the fight against cancer. To this end, we investigated how the re-polarisation of the M2 TAMs into M1 TAMs can affect cancer development, by focusing especially on the macrophage populations near the leading edge of the tumour. Specifically, we studied the spatial aspect of the M2→M1 re-polarisation through the definition of a repolarisation domain p (t, R p ), and the temporal aspect via the starting re-polarisation time t p used in Equation (12). Finally, while in this work we considered the death of macrophages induced through nutritional starvation, this is only one aspect in the wider picture concerning death and survival of macrophages within the necrotic region, which so far has not been fully explored experimentally. Future mathematical modelling studies will explore other death and survival aspects involved concerning macrophages in necrotic regions.
To propose new hypotheses, we first introduced a macroscale quasi-steady reaction-diffusion equation for the nutrients where we considered the spatial transport to be described by diffusion with constant-coefficient as well as a linear uptake rate. To account for the effect of nutrients on the different cell functions, we defined four effect-functions that we used for the rest of tumour dynamics. Furthermore, we introduced another macro-scale equation, describing the behaviour of the M1 TAMs where the motility of the M1 phenotype is represented both by random and directed movements. The rest of the equation involves an influx term, a nutrient-dependent proliferation and death laws, as well as a nutrient-dependent polarisation and a re-polarisation terms that describe the switch between the two phenotypes. Similarly to the M2 TAMs, the M1 phenotype also secrete MDEs that can degrade the ECM. Therefore, this M1 phenotype directly contributes to the re-arrangement of the micro-fibres constituents, as well as serving as a source of MDEs for the proteolytic processes that occur on the invasive edge of the tumour.
We used this extended model to explore the possibilities of macrophage re-polarisation to depend on the spatial domain as well as on time. First, in Figures 7-9 we presented the result of the spatial dependency by varying the re-polarisation domain p (t, R p ). We concluded that even though the tumour spread does seems to be affected much by the M2→M1 repolarisation (which may be expected as biological studies [110,112,113] have shown that M1 TAMs located in the stoma can promote cancer progression), the tumour mass can be reduced significantly. Therefore, even though we may need additional strategies directly targeting tumour spread, the tumour mass can potentially be reduced by the re-polarisation of M2 TAMs to M1 TAMs.
Finally, since tumours are only detectable above a certain size and therefore the M2→M1 re-polarisation is usually applied at later stages in tumour development, in Figure 10 we investigated the temporal dependency of M2→M1 re-polarisation. There, we showed that while the smallest tumour mass can be obtained when the re-polarisation starts as soon as possible, in some cases, it is possible to keep the tumour under control even when we repolarise at later times. For example, in Figure 10A, for high ECM degradation levels, tumour mass did not change when t p ≤ 25 t. However, tumour mass slowly increased as we delayed the repolarisation time for medium and low ECM degradation. All these theoretical results suggest that macrophages re-polarisation protocols (through various immunotherapies, such as the use of agonist anti-CD40 mAb [15]) might have to be adapted to the tumour environment and the degradation levels of the ECM.
We emphasise that due to the complexity of this modelling framework, which forced us to work with a nondimensionalised model, the results in this study are only qualitative. Nevertheless, they provide new theoretical hypotheses regarding the possible roles of macrophage re-polarisation (within specific regions of the solid tumours, and within specific time scales) in the context of immunotherapies for cancer. Furthermore, while in this study we considered a generic solid tumour rather than a specific one, the findings here are in principle applicable to specific tumour types such as gliomas as we have done it also in our recent publication [114].
To conclude, we suggest that in addition to the re-polarisation of M2 TAMs, we also need additional strategies targeting tumour spread or tumour stroma in order to fully stop the tumour from advancing. In this context, it is worth mentioning that we started modelling the effect of vasculature (existing and newly growing blood vessels) on the transport of nutrients and immune cells into the tumour; but this will form the topic of a future publication. Moreover, beyond the current modelling approach, considering the time evolution of the vasculature would also allow us to gain a more comprehensive insight into the effects of nutrients on other immune cells (such as T cells), the interplay between different cells [such as cancer, TAMs and cancer-associated fibroblasts [115,116]] and additional surplus amount of cytokines secreted by other immune cells in the tumour microenvironment. Furthermore, modelling an evolving network of blood vessels would enable us to use it as a supply of macrophages (and other cells) or create more complex hypoxic and necrotic regions that would induce further heterogeneity in all cell populations, ultimately leading to more realistic scenarios.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.