Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 2, Number 1, March 2021, pp.32-54 https://doi.org/10.5206/mase/11087 A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL THOMAS CARRARO, SVEN E. WETTERAUER, ANA VICTORIA PONCE BOBADILLA, AND DUMITRU TRUCU Abstract. The quest for a deeper understanding of the cancer growth and spread process focuses on the naturally multiscale nature of cancer invasion, which requires an appropriate multiscale modeling and analysis approach. The cross-talk between the dynamics of the cancer cell population on the tissue scale (macroscale) and the proteolytic molecular processes along the tumor border on the cell scale (microscale) plays a particularly important role within the invasion processes, leading to dramatic changes in tumor morphology and influencing the overall pattern of cancer spread. Building on the multiscale moving boundary framework proposed in Trucu et al. (Multiscale Model. Simul. 11(2013), 309-335), in this work we propose a new formulation of this process involving a novel derivation of the macro scale boundary movement law based on micro-dynamics, involving a transport equation combined with the level set method. This is explored numerically in a novel finite element macro-micro framework based on cut-cells. 1. Introduction Involving a wide range of cross-related processes occurring on several spatio-temporal scales, cancer cell invasion in human tissue is one of the hallmarks of cancer [36], playing a crucial role in the overall development of a growing malignant tumour. Taking advantage of the heterotypic character of the tumour microenvironment (which includes immuno-inflamatory cells, stromal cell, fibroblasts, endothe- lial cells, macrophages), complex molecular processes facilitate intense interactions between the cancer cell population and the extracellular matrix (ECM)[29, 36, 39, 40, 58]. These interactions lead to a cascade of specific developmental patterns and behaviours of the growing tumours, most notable stages including the degradation of the ECM, the local progression of the tumour, followed by the tumour angiogenesis process and the subsequent metastatic spread of the cancer cells in the human body. The alteration and remodelling of the ECM by the matrix degraded enzymes (MDEs) such as ma- trix metalloproteinases MMPs or the urokinase plasminogen activator (uPA) play a key role in local tumour progression. Alongside cell-adhesion and multiple taxis processes (including haptotaxis and chemotaxis), the matrix degrading enzymes processes degrade various components of the surrounding ECM that leads to further tumour progression. However, as the full mechanisms involved in these complex processes is yet to be deciphered biologically, over the past two decades or so cancer invasion received extensive mathematical modelling attention, in which systems of reaction-diffusion-taxis par- tial differential equations [3, 7, 8, 9, 13, 15, 19, 20, 25, 33, 35, 49, 50, 54, 55, 74] as well as nonlocal integro-differential systems [10, 18, 26, 34] were derived and proposed to deepen the understanding, validate and create new experimental hypothesis. Furthermore, to capture various heterotypic aspects Received by the editors 15 October 2020; accepted 11 March 2021; published online 15 March 2021. 2000 Mathematics Subject Classification. 65M60, 35K57, 35Q92. Key words and phrases. cancer invasion; multiscale modelling; level set; FEM; cut-cells. AVPB was funded by the Heidelberg Graduate School of Mathematical and Computational Methods for the Sciences (HGS MathComp), founded by DFG grant GSC 220 in the German Universities Excellence Initiative. 32 A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 33 and related processes within tumour invasion, several multiphase models based on the theory of mix- tures [14, 17, 21, 31, 56, 61, 69, 77, 78] were derived (by exploring the mass and momentum balances as well as the inner multiphase constitutive laws). A particularly important role in cancer invasion is played by the MDEs (such as the MMPs) that are secreted from the outer proliferating rim and released within the tumour peritumoural microenvi- ronment. This gives the cancer invasion a moving boundary character, and to that end several level set approaches were recently proposed to study the tumour progression both in homogeneous environments [32, 43, 44, 45, 80] and in complex heterogeneous tissues [46]. Despite recent advances, the multiscale modelling of the processes involved in cancer invasion remains an open problem. Although this is a truly multiscale process, most mathematical models were offering a one-scale perspective, whether that is from a purely macro scale (tissue scale) or an exclusively micro scale (cell scale) stand point. However, recently a novel 2D multiscale moving boundary modelling platform for cancer invasion was proposed in [73]. This explores in an integrated manner the tissue- scale cell population dynamics and relevant cell scale molecular mechanics together with the permanent link between these two biological scales. This addresses directly the dynamics of the MDEs proteolytic processes occurring at the tumour boundary (i.e., at the invasive edge of the tumour) that are sourced from within outer proliferating rim of the tumour and facilitate the complex molecular transport and ECM degradation within the peritumoural region. The tissue-scale progression of tumour morphology is captured here in a multiscale moving boundary approach where the contribution arriving from the cell scale activity to the cancer invasion pattern is realised by the micro scale MDEs dynamics (occurring along the tumour invasive edge), which, for its part, is induced by the cancer macro-dynamics. This was recently applied to the extended context in which, rather than the MMPs dynamics, the uPA is considered as the proteolytic system, and has led to biologically relevant results [53]. In this work we present a new formulation of the link that connects the two scales of the tumour dynamics that was considered in the initial multiscale moving boundary modelling approach presented in [73] as well as in the multiscale cancer invasion modelling developments that followed [4, 5, 6, 71, 52, 62, 63, 64, 65, 67, 72]. In particular, the movement of the tumor boundary is here defined by a velocity field instead of a displacement of the interface. Therefore, the moving interface is defined and implemented in a different way. Furtheremore, the discretization is based on finite elements instead of finite differences as in the previous work. The new model is based on a level set approach in which the moving domain is defined as the zero level of a level set function. The reason for this choice of the problem setting is twofold: on one hand, all components of the model can be described by partial differential equations at the continuum level allowing the complete separation between modelling and discretization; on the other hand, it is better suited for an extension to the three-dimensional case since the formulation of all components of the model is dimension independent and the use of a dimension independent implementation of the discretization, like in our case using the finite element method (FEM) package deal.II [11], facilitates the realization of the code. The level set method was first introduced in [51] for tracking moving interface with complex defor- mations. This method was developed starting from the notion of weak solutions for evolving interfaces. The main aspect of this method is that an interface or a domain boundary is defined through the embedding of the interface as the zero level set of a higher dimensional function. Furthermore, the velocity of the interface is also embedded to the higher dimensional function. We avoid handling a sharp interface, i.e. a lower dimensional manifold in the computational domain, but the velocity needs to be extended from the interface to the rest of the domain. While the original setting, with the sharp interface, poses several numerical difficulties due to its Lagrangian approach, the later setting, using an Eulerian approach, can exploit techniques developed for hyperbolic problems. 34 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU Other works have presented a level set approach for moving the tumour interface. In [43] the authors use a level set function to define the boundary of a tumour mass and extend the velocity orthogonally to the interface using a filter technique to damp numerical noise coming from the extension procedure. The velocity at the interface is defined as a function of the gradient of a computed quantity (the pressure). This work nevertheless does not link different model scales. In [80] an adaptive finite element combined with a level set approach is used to solve a model that considers tumour necrosis, neo-vascularization and tissue invasion. The model is composed of a continuum part and a hybrid continuum-discrete part. The velocity of the interface is the cell velocity. Therefore, the velocity does not need to be extended into the neighbourhood of the interface. In [47] a level set approach with a ghost-cell method is applied to tumour growth of glioblastioma. The velocity of the interface depends on solutions of linear and nonlinear equations with curvature-dependent boundary conditions. Since the velocity is only defined at the interface the authors extend it beyond the interface and use a narrow band/local level technique to update the interface velocity and level set function only in the vicinity of the interface. Our approach uses a level set method with an extension of the velocity. While the coupling between the macroscopic and microscopic scales was originally introduced in [73] considering a Lagrangian approach to move the nodes of the discrete approximation of the interface, we introduce here a continuous link between the two scales defined by the velocity at the interface at the continuum level. This changes the formulation of the multiscale coupling that goes with the definition of the velocity starting from heuristic arguments. The scope of this work is to present as a whole the new formulation of the tumour invasion model explaining the possible advantages that this approach can have for future developments and the numerical aspects that need further attention and further development. The paper is organized as follows. In Section 1 we state the problem setting and describe the different components of the model: the macroscopic and microscopic components and the description of the moving boundary. In Section 3.2 we introduce the weak formulation of the model which is needed for the approximation of the continuum problem with a finite element method. We introduce the discretization of the problem using cut-cells for the approximation of the cancer region domain. In Section 4 we present some numerical results showing the interplay of the different parts of the multiscale model. Finally we present an outlook and some concluding remarks in Section 5. 2. The two-scale tumour dynamics We present a two-scale model for cancer invasion that involves a double feedback loop to link the dynamics occurring at two different spatial scales explored by the following two modelling components: a macroscopic component describing the population of cancer cells and extracellular matrix at tissue-scale and a microscopic component describing the dynamics of a generic matrix-degrading enzyme molecular population whose cell-scale takes place at the leading edge of the tumour. Both scales are considered at the continuum level, and we assume that possible stochastic effects (in regions where the continuum assumption is not valid) can be neglected. Nevertheless, besides the establishment of a new approach to linking the scales, our aim includes also the derivation of a flexible numerical framework that would allow to extend the model with a stochastic part (e.g. at the interface of the domain) that would require to a hybrid formulation. The cancer cells population mixed with the extracellular matrix density exercises its macro-scale interacting dynamics within a tumour invading domain Ω(t) that changes its size and morphology in time during the invasion process within a reference tissue domain Y . Its boundary ∂Ω(t) will also be referred to as interface because this represents the tumour interface separating the healthy region with zero cancer cells from the region with a distribution of cancer cells. A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 35 Assumption 2.1 (Scale separation). We assume here scale separation in space and time between the tissue-scale tumour macro-dynamics (involving cancer cells population mixed with ECM) and the cell-scale molecular boundary micro-dynamics (involving matrix degrading enzymes). The model is of two-scale nature, since, due to the different physical dimensions of tissue and tumour constituents that are involved in the process (such as of cells, ECM, and matrix degrading enzymes), cancer invasion is genuinely a multiscale biological process [75], which in its most basic setting involves a cancer cell population mixed with ECM density (for the macroscopic part) and a population of matrix degrading en- zymes (for the microscopic part) that have their dynamics occurring at separated spatial scales (namely, at macro- and micro- scale, respectively). Hence, the natural assumption of scale separation enables us to derive the multiscale model and to address the dynamics at each of the two scales on correspondingly different scale domains. Finally, the characteristic length L for the macroscopic part of the model re- lates to the diameter of the cancer region and is considered here as in [8, 34], ranging between 0.1cm and 1.0cm. Furthermore, the characteristic length ` of the microscopic part related to the region where the matrix degrading enzymes are spatially transported is considered to be of the order of 10−3cm, [48]. The ratio between the scales is denoted ε = `/L, and so in the non-dimensional setting of the model we propose here ε will represent the cell-scale size (i.e., the micro-scale size). During the macro-dynamics, the cancer cels from the outer proliferating rim secrete matrix degrading enzymes (such as MMPs), providing this way a source of MDEs along the leading edge of the tumour [36, 75], i.e., in the cell-scale proximity of the tumour boundary. Once secreted, these MDEs cell-scale exercise a diffusive transport across the tumour interface in the peritumoural region, causing degradation of the ECM that they meet in the immediate proximity of the tumour and this way enabling further tumour progression [36, 75]. Proceeding in a similar manner as in [73], the boundary micro-dynamics of the matrix degrading enzymes population is explored here on a bundle of micro-domains εY , of cell-scale size � > 0, whose union provides a cell-scale neighbourhood for the interfacial points in ∂Ω(t). Using the scale separation, we explore the micro-dynamics within a cell-scale neighbourhood of each of the boundary points x ∈ ∂Ω(t) (i.e., at each point of the macroscopic interface) that is enabled by acorresponding εY micro-domain centred at x. In brief, adopting here a multiscale modelling perspective similar to the one proposed in [4, 5, 6, 71, 52, 62, 63, 64, 65, 67, 72, 73], the coupling of two-scale dynamics of cancer invasion is captured as follows: • the “top-down” macroscopic-to-microscopic coupling is realised via the source of the matrix- degrading enzymes, which is induced by the tumour macro-dynamics and is formed as a collec- tive contribution of the cancer cells that arrive during the macro-dynamics within an appropriate distance from any boundary point x ∈ ∂Ω(t), see equation (2.8); • the “bottom-up” microscopic-to-macroscopic coupling determines the movement of the macro- scopic tumour interface that is induced by the dynamics of the micro-spatial distribution of MDEs within the cell-scale neighbourhood of the tumour interface (enabled by an appropriate union of boundary micro-domains, see Figure 1) in the form of a velocity field (see equation (2.10)) that drives a transport process for boundary relocation at macro-scale (see equation (2.6)) . The rate of cancer cells invasion into the surrounding tissues is therefore driven by the velocity of the interface that depends on the boundary MDE micro-dynamics microscopic enzyme dynamics. Further- more, the tumour interface is described here by the zero-level of a level set function and its spatial dynamics is governed by a transport equation. To account for all interactions between the different parts of the problem, the tissue macro-domain Y is assumed to be sufficiently large such that the complete dynamics happen inside it. 36 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU In the following we proceed with the multiscale model description in three parts: macroscopic, microscopic and a transport process for the tumour boundary which is induced by the micro-scale. Focussing first on the macro-scale, at this level the model considers cancer cells and extracellular matrix (ECM) dynamic interaction as well as a micro-scale-induced transport process that will ultimately describe the movement of the macro-scale tumour boundary. 2.1. Tumour macro-dynamics. Let c(x,t) and v(x,t) denote the cancer and the extracellular matrix distributions at (x,t) ∈ Ω(t) × (0,T), respectively. Proceeding as in [73], the dynamics at macroscopic scale is given by the following PDE system: ∂c ∂t = Random motility︷ ︸︸ ︷ D1∆c − Haptotaxis︷ ︸︸ ︷ η∇· (c∇v) + Proliferation︷ ︸︸ ︷ µ1 c (1 − c−v), (2.1) ∂v ∂t = − αcv︸︷︷︸ Degradation + µ2(1 − c−v)︸ ︷︷ ︸ ECM Remodelling , (2.2) with boundary condition (D1∇c−ηc∇v) ·n = 0 on ∂Ω × (0,T). (2.3) and initial conditions c(x, 0) = c0(x) on Ω(0) (2.4) v(x, 0) = c0(x) on Ω(0). (2.5) where D1 is the diffusion coefficient for the cancer cells, η is the advection coefficient, µ1 is the proliferation coefficient, α a degradation coefficient and µ2 a coefficient for the remodelling of ECM. All these coefficients are considered constant and their typical values are included in Table 2. It is assumed that the cancer cells are zero outside Ω(t) and that there is no transport of cells through the boundary ∂Ω(t), see boundary condition (2.3). Furthermore, under the presence of these boundary and initial conditions, for the case of constant proliferation rate µ1, the results in [68, 76, 37] explore the local and global existence of system (2.1)-(2.2). The initial distribution of cancer cells c0(x) and extracellular matrix v0(x) are given in the larger domain Y . In the next section we introduce the part of the model used to update the tumour domain in time. 2.2. Two-scale tumour boundary movement: the macro-scale transport process induced by the MDEs boundary micro-dynamics. As mentioned above, the movement of the tumour interface is directly governed by the matrix degrading enzymes (MDE) dynamics occurring in a cell scale neighbourhood of the tumour interface ∂Ω(t). The pattern of degradation of the peritumoural ECM by the advancing front of MDEs drives the invasion of the tumour cells in the surrounding tissues and determines the movement of the tumour boundary ∂Ω(t). Therefore, the movement of the time- dependent macro domain Ω(t) is enabled by a velocity field defined on the points of the interface x ∈ ∂Ω(t), which is determined by the micro-dynamics occurring on a small micro-domain εY centred at x. Hence, the velocity field generated in this way is induced directly by the micro-dynamics of the MDE molecular distribution m(y,τ) over a suitable micro-spatio-temporal domain εY ×(0, ∆T) (which will be detailed in Section 2.3). We denote this velocity field by V (m). Since V (m) is defined only on points at the interface, we consider an extension of the velocity to the whole domain Y . This allows us to describe the cancer region boundary by a level set approach. The interface is defined as the zero-level of the level set function φ which satisfies the following transport equation: ∂φ ∂t + V (m) ·∇φ = 0, in Y × (0,T). (2.6) A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 37 For later purposes, we introduce the notation L0(t) = {x ∈ Y : φ(x,t) = 0} (2.7) for the zero level of the level set function that defines the interface ∂Ω(t). A natural extension of the velocity is the constant continuation of the velocity at the boundary in normal direction [23]. In Section 3.6 more details about this point are given. 2.3. MDEs boundary micro-dynamics. Due to the scale separation Assumption 2.1 we can describe the micro-dynamics of the MDEs within a cell-scale neighbourhood of the tumour interface by exploring this on a bundle of micro-domains εY defined and centred at each macroscopic interface point x ∈ ∂Ω(t). Assuming for convenience that the maximal domain Y is centred at origin of the space, the micro scale coordinates y of the micro-scale problem on a εY centred at x0 are obtained by an appropriate scaling and translation of Y (given by the transformation y = x0 + ε(x−x0)), see Figure 1 for an illustration of these micro-domains on the ∂Ω(t). As argued in [73], collectively, the cancer cells that get to be point on the macroscopic tumour boundary microscopic domains �Y ∂Ω(t) �Y �Y Figure 1. Sketch of micro-dynamics sampling at the macroscopic tumour boundary ∂Ω(t) distributed during their dynamics within a certain radius Rm > 0 from the interface x ∈ ∂Ω(t) secrete matrix degrading enzymes, giving rise this way to a source of MDEs within each boundary micro- domain εY . This micro-scale source MDEs (which is induced by the macro-dynamics) can therefore be formalised mathematically as Fx,t(c)(y) :=   1 |B| ∫ B c(ξ,t) dξ y ∈ εY ∩ Ω(t), 0 otherwise. (2.8) where B := {ξ ∈ Y : ‖ξ−x‖≤ Rm}. Therefore, in the presence of source (2.8), the cross-interface micro- dynamics of MDE molecular distribution m(y,τ), which takes place on the boundary micro-domain εY over a time interval (0, ∆T), is governed by the following reaction diffusion equation ∂m ∂t (y,τ) = D2∆m(y,τ) + Fx,t(c) in εY × (0, ∆T), m(y, 0) = 0 in εY, ∂m ∂n (y,τ) = 0 in ∂εY × (0, ∆T), (2.9) with ∆T> 0 representing here the micro scale time perspective and serving also later on as natural time interval for the coupling between the microscopic and macroscopic stages of the multiscale model. 38 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU Therefore, the pattern of peritumoural ECM degradation caused by the advancing fronts of MDE molecules (which are transported across the tumour interface in the immediate proximity within the appropriate microscale region) gives rise to a boundary velocity that can be described by V (m) = cvel ∆T |εY | ∫ ∆T 0 ∫ εY m∇m dydτ, (2.10) where |εY | = ∫ εY 1 dt and cvel is a tuning scaling factor, see Table 2. Specifically, this form of V is based on the following main considerations: • the term ∇m takes into consideration the assumption that the cancer boundary moves following the gradient with respect to the MDE; • further, by multiplying it by m(y,τ), we are taking into account the influence of the amount of enzymes over their given gradient direction at each spatio-temporal micro-node (y,τ), enabling an appropriate weighting of its “strength” (magnitude); • finally, by considering the average contribution of MDE microdynamics over εY × [0, ∆T ] by simply accounting upon the mean-value in time of the revolving weighted MDE gradient spatial direction [0, ∆T] 3 τ 7→ 1 |εY | ∫ εY m∇m dy, we ultimately obtain the definition of the velocity given in (2.10), where V (m) is taken as being proportional to this spatio-temporal mean value, with proportion constant cvel > 0. 2.4. Schematic summary of our multiscale moving boundary modelling. The new model that we introduced here falls in the class of heterogeneous multiscale models that were developed over the past two decades not only for multiscale biological processes but also for other multiscale processes arising in material science or fluid-structure interactions [1, 4, 5, 6, 22, 71, 27, 28, 30, 52, 62, 63, 64, 65, 67, 72, 73]. Schematically, the two-scale dynamics of our cancer invasion model is coupled across the scales as depicted in Figure 2, and its progression can be summarised in the following three steps: (1) At a time t∗> 0, the macro-scale cancer distribution (i.e., the solution of the macro-dynamics) on the domain Ω(t∗) induces in a non-local manner the enzymatic source for the MDEs boundary micro-dynamics. (2) The boundary micro-dynamics is explored on each boundary micro-domain �Y . The time-space average of the microscale MDEs spatio-temporal distribution over the micro-domain �Y and a fixed but arbitrarily small time range of size ∆T> 0 is used to determine pointwise the tumour interface velocity, which will ultimately result in describing the direction and displacement magnitude of the macroscopic tumour boundary relocation. (3) The interfacial velocity obtained from the boundary MDEs micro-dynamics is then set into a transport equation that finally determines the the position of the tumour boundary at time t∗ + ∆T . Therefore, a new tumour macro-domain Ω(t∗ + ∆T) is defined, and this becomes the new playground for the cancer macro-dynamics, which continues now its evolution on this newly expanded domain, progressing again through the stages described in (1)-(3). 3. Multiscale Computational Approach 3.1. Definition of the computational microscopic problem. We consider the microscopic problem in a bundle of boundary microdomains εY , with y := (y1,y2) being the standard local microscale reference system within a given εY and τ denoting the time at micro scale. As will be explained below in Section 3.5, in our finite element approach the macroscopic dynamics will be considered on an appropriately defined macroscopic domain Ωh(t) with a linearized boundary A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 39 boundary dynamics b o tt o m − u p to p − d o w n transport process for tumor macro−dynamics macro−dynamics MDE micro−dynamics φ(t∗ + ∆T)φ(t∗) Ω(t∗) Ω(t∗ + ∆T) Figure 2. Sketch of the coupling (including the top-down and bottom-up links) be- tween the macro- and micro- dynamics of our two-scale moving boundary model for tumour invasion. ∂Ωh(t). The microscopic dynamics is then explored within a microdomain εY centred at a macro scale boundary point x ∈ ∂Ωh(t) and eventually appropriately rotated so that this is positioned with two edges parallel to the linearized boundary (in direction y1) and two edges orthogonal to it (in direction y2) as shown in Figure 1. This simplifies the setting of the microscopic problem. In fact, since we consider the linearized boundary ∂Ωh(t), the right hand side Fx,t in equation (2.9) does not depend on y1. Furthermore, we would like to note that, due to scale separation, the quantity m does not depend on the macroscopic variable x. Nevertheless, since each microscopic domain is centred at a different point x on the boundary, a potentially different MDEs micro-source is induced by the macro-scale for each micro-dynamics on each microscopic domain �Y . The value of cancer cells concentration c in one single microscopic domain is constant because (due to scale separation) no oscillations in y are considered for c. In addition, since on the boundaries of the quadrilateral domain no flux conditions are prescribed, it follows that the solution is constant in y1 direction. In fact, it is straightforward to show that the 1D solution is also solution of the 2D problem. Due to uniqueness of both problems this constant property is given. Therefore, we can consider the following simplified one-dimensional microscopic problem for the quantity m, which is the integral of m along y1 (giving the amount of enzyme molecules per unit of length) ∂m ∂τ (y2,τ) = D2∆m(y2,τ) + Fx,t(y2) in (0,εL) × (0, ∆T), m(y2, 0) = 0 in (0,εL), ∂m ∂n (y2,τ) = 0 in ∂(0,εL) × (0, ∆T), (3.1) where Fx,t is Fx,t integrated over y1. Since Fx,t and m do not depend on y1, the solution m of (2.9) is the constant extension of the solution m of (3.1) in y1 direction. We introduce now a scaling of the domain to the interval (0, 1) through the following tranformation y2 = εLz with z ∈ (0, 1), 40 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU then we get after the rescaling the transformed system ∂m̂ ∂t (z,τ) = D2ε −2L−2∆m̂(z,τ) + F̂x,t(z) in (0, 1) × (0, ∆T), m̂(z, 0) = 0 in (0, 1), ∂m̂ ∂n (z,τ) = 0 in ∂(0, 1) × (0, ∆T), (3.2) with F̂x,t(z) :=   1 |B| ∫ B c(ξ,t) dξ z ∈ [0, 1/2], 0 otherwise, (3.3) note that the coordinate ξ is a macroscopic quantity. Notice furthermore that a solution of (3.1) is a solution of (3.2) by m̂(z,τ) = m(εz,τ). Remark 3.1 (Limit ε → 0). In case of ε → 0 we have in (3.2) a large diffusion coefficient, therefore a fast redistribution process of the solution occurs, leading to negligible spatial variations of the solution. The only relevant parameter of the problem at the limit becomes the time. The limit problem becomes an ordinary differential equation (ODE). Even if we consider scale separation in this model, we do not consider the limiting case ε → 0. In that case the velocity has to be defined in a different way since the term ∇m becomes the zero vector. The parameter ε in our model has always a finite value bounded below ε ≥ ε cell > 0, where ε cell is assumed here to be a minimal microscale size of the order of a cell length. Thus, using (3.2), we obtain that the velocity can be further expressed as V (m) = cvel ∆Tε2 ∫ [0,ε]2 ∫ ∆T 0 m∇ym dτ dy = cvel ∆Tε2 ε ∫ ε 0 ∫ ∆T 0 m∇ym dy dτ, = cvel ∆T ε ∫ 1 0 ∫ ∆T 0 m̂∇zm̂ dz dτ, (3.4) where m̂ indicates the transformed function on the unit domain. 3.2. Weak formulation of the two-scale tumour invasion model. To describe the model in the setting needed for the FEM we introduce the following weak formulation. We use the notation (·, ·) to define the usual L2 scalar product of Lebesgue square integrable functions. The space H1 is the Hilbert space of square integrable functions with square integrable (weak) first derivative and H∗ is its dual space, i.e. the space of bounded linear functional on H1. Furthermore, we use Bochner spaces of the form U = {u ∈ L2(0,T; H1) : ∂tu ∈ L2(0,T; H∗)} to introduce the weak formulation for the dynamics at each of the two scales. Specifically, we have the following: • at macro-scale: – for the weak formulation of the dynamics of the cancer cells population (2.1) as well as for the dynamics of the ECM (2.2), we consider the space UM = {u ∈ L2(0,T; H1(Ω(t))) : ∂tu ∈ L2(0,T; H∗(Ω(t)))} – for the boundary transport process (2.6), we consider the space UT = {φ ∈ L2(0,T; H1(Y )) : ∂tφ ∈ L2(0,T; H∗(Y ))}; • at micro-scale, for the MDEs micro-dynamics (3.2), we consider: Um = {m ∈ L2(0, ∆T ; H1((0, 1))) : ∂τm ∈ L2(0, ∆T; H∗((0, 1)))} A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 41 Thus, at macro-scale the weak formulation for the tumour macro-dynamics (2.1)-(2.2) is as follows: find (c,v) ∈ UM ×UM such that, for almost all t ∈ (0,T), it satisfies(∂c ∂t ,ϕ ) + ( D1∇c,∇ϕ ) − ( ηc∇v,∇ϕ ) − ( µ1(v) c(1 − c−v),ϕ ) = 0 ∀ϕ ∈ H1(Ω(t)), (3.5a)(∂v ∂t ,ϕ ) + ( αcv,ϕ ) − ( µ2(1 − c−v),ϕ ) = 0 ∀ϕ ∈ H1(Ω(t)), (3.5b) c(x, 0) = c0 in Ω(0), (3.5c) v(x, 0) = v0 in Ω(0). (3.5d) Note that the initial conditions for cancer cells, c0, and ECM distributions, v0, are defined on the larger maximal domain Y and that this formulation implies the natural zero flux condition for the cancer cells and ECM, i.e. ∂nc = ∂nv = 0 on ∂Ω(t). Similarly, the weak formulation of the macro-scale transport equation for interface dynamics (2.6) is: find φ ∈ UT such that, for almost all t ∈ (0,T) it satisfies(∂φ ∂t ,ϕ ) + ( V (m) ·∇φ,ϕ ) = 0 ∀ϕ ∈ H1(Y ), (3.6a) φ(x, 0) = φ0 in Y, (3.6b) with φ0 being the level set function at the initial time. Using the notation ϕ+ := max{ϕ, 0} to define the positive part of a function ϕ, we have that supp(φ+0 ) describes the initial support region of the cancer cells. Finally, at micro-scale, the weak formulation of the MDE micro-dynamics associated with each of the boundary micro-domains �Y is: find m̂ ∈ Um such that for almost all τ ∈ (0, ∆T) it satisfies(∂m̂ ∂τ ,ϕ ) + ( D2ε −2L−2∇m̂,∇ϕ ) = ( F̂x,t,ϕ ) ∀ϕ ∈ H1((0, 1)), (3.7a) m̂(z, 0) = 0 in (0, 1), (3.7b) where F̂x,t is defined as in (3.3) and the natural condition ∂nm̂ = 0 is implicitly defined. 3.3. Discretization. The model is first discretized in time by the implicit Euler method and then discretized in space by the FEM. The discretized system is defined on a regular mesh Mh composed of quadrilateral cells K ∈ Mh of the same dimension. This mesh has the advantage that it can be generated starting from an initial square that is successively refined to achieve a given cell diameter. In this work we use global refinement to generate the mesh. Since the macroscopic domain Ω(t) is time-dependent, the discrete space domain would need to be remeshed at every time step, if a fitted FEM formulation is used. In case of large deformations, the procedure of remeshing has to deal with the possible loss of shape regularity of the mesh. To avoid these complications related to remeshing, we use an unfitted approach by using so called cut-cells. These are a special realization of the FEM as described below. In particular, these are finite elements with shape functions with a support on a subdomain of the cells that is defined by the intersection of the interface with the cells. Let us consider the space of bi-linear polynomials Q1 defined on a unit cell K̂ = [0, 1] 2, i.e. Q1 = span(1,x,y,xy) and the space of linear functions P1 = span(1,x) 42 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU defined on the one-dimensional unit cell K̂ = [0, 1]. The finite element space is defined as UhM (t) = {u ∈ C(Ω(t)) :u|K ◦TK ∈ Q1 if K ∩∂Ω(t) = ∅; u|K∩Ω(t) = ψ|Ω(t), ψ ◦TK ∈ Q1 if K ∩∂Ω(t) 6= ∅}, where TK is a bijective transformation from the unit cell K̂ to the physical cell K. Since the mesh is non-fitted, the shape functions u (in case of cut-cells) are defined only on the portion of cell that is intersected by Ω(t), while outside of Ω(t) they need not to be defined. In Figure 3 the restriction of two shape functions ϕ1 and ϕ2 on the cut-cell is shown. For the transport and microscopic problems 0 1cut−cell ϕ1 ϕ2 Figure 3. Shape functions on a unit cut-cell we use the following spaces UhT = {u ∈ C(Y ) : u|K ◦TK ∈ Q1} and Uhm = {u ∈ C((0, 1)) : u|K ◦TK ∈ P1}. 3.4. Approximation of the solution on the moving domain. In every time step of the time discretization scheme the domains at tn and tn+1, Ω(tn) and Ω(tn+1), are defined by the level set function at the two times tn and tn+1. At the n th time step the solution is known only in Ω(tn). To determine the solution at tn+1 the variation of the domain in time should be taken into account in the formulation of the problem. One possible formulation of the problem would be to consider a reference domain Ω(t0) for a given t0 in each time step and to use a (time dependent) mapping to transform the solution from the domain Ω(tn) to the reference domain Ω(t0) and then to the domain Ω(tn+1). However, if the time step is small enough the combination of these two transformations can be approximated with the identity and, instead of implementing a complicated formulation, an extension of the solution uh(tn) from the old domain Ω(tn) to the new domain Ω(tn+1) could be used. This procedure introduces an approximation error that for small enough time steps can be neglected in comparison to other sources of error (such as space-time discretization of the solution, or splitting error etc.). In the following, we describe the extension used in this work and advice that this should be defined properly depending on the finite element ansatz used. A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 43 0 1st cut 2nd cut extension 1 tn tn+1 Figure 4. One dimensional sketch of the extension of the macroscopic solution in case of a cell cut twice 0 01st cut 1 1 extension 2nd cut tn+1tn Figure 5. One dimensional sketch of the extension of the macroscopic solution in case the interface cuts two neighbour cells at tn and tn+1. Following the above construction, the fully discrete formulation of the macroscopic problem becomes( cn+1h ,ϕ ) Ω(tn+1) + k [( D1∇cn+1h ,∇ϕ ) Ω(tn+1) − ( η cn+1h ∇vh,∇ϕ ) Ω(tn+1) − ( µ1(v n+1 h ) c n+1 h (1 − c n+1 h −v n+1 h ),ϕ ) Ω(tn+1) ] = ( c̃nh,ϕ ) Ω(tn+1) , ∀ϕ ∈ UhM (t n+1), (3.8a) ( vn+1h ,ϕ ) Ω(tn+1) + k [( αcn+1h v n+1 h ,ϕ ) Ω(tn+1) − ( µ2 (1 − cn+1h −v n+1 h ),ϕ ) Ω(tn+1) ] = ( ṽnh,ϕ ) Ω(tn+1) , ∀ϕ ∈ UhM (t n+1), (3.8b) where k = tn+1 −tn is the macroscopic time step, c̃nh and ṽ n h are the extensions from Ω(tn) to Ω(tn+1). In fact, the two components ch(tn) and vh(tn) are defined only in Ω(tn). Therefore an extension in the region Ω(tn+1)\Ω(tn) needs to be defined. We have chosen a continuous extension using the prescribed values c0(x) and v0(x) for x ∈ Y . In particular, we have considered two cases: case (i) the cell where we need to define the extension is cut at time tn and at time tn+1, see Figure 4 and case (ii) the cell is cut at time tn and uncut at time tn+1, see Figure 5. In case (ii) the cut goes to the neighbour cell at time tn+1. In case (i) both components are extended up to the new cut using the values of all degrees of freedom of the considered cell (also those lying outside the domain Ω(t)) with a bilinear nodal interpolation. Note that the bilinear nodal interpolation is justified only within the domain Ω(t) in the cut-cell formulation of the problem, because the integrals in the weak formulation are computed only in the inner part of the cells. In this sense, we are “extrapolating” the values ch and vh outside the region of validity of the finite element interpolation. It is known that the problem formulation using cut-cells can lead to instabilities due to the lost of coercivity outside of the physical domain Ω(tn). To overcome to this 44 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU problem special stabilization techniques can be used, e.g. the ghost penalty method [12]. To implement stabilization techniques extra boundary integrals have to be implemented increasing the complexity of the code. We have therefore decided to use a heuristic solution. We have set a threshold of 1% on the volume to be considered for extrapolation. Cells, whose volume is cut by 99%, are eliminated from the mesh. The value of 1% has been chosen after testing different values and we have seen that also with larger values we got the same boundary, meaning that this threshold is accurate enough. In fact, neglecting the contribution of small cells can be interpreted as a quadrature error. By keeping small the threshold for the suppression, we keep small this quadrature error. Stabilization of the haptotaxis term. Due to the large difference between the diffusion coefficient D1 and the haptotactic coefficient η, see Table 2, the macroscopic problem is transport dominant and therefore must be numerically stabilized. Typically, upwind techniques are used in the finite element framework to stabilize the calculations [41]; see, for example, [66] for an application of these techniques to a chemotaxis problem. In this work, we have chosen a streamline diffusion stabilization [81] that adds an artificial diffusion term only in the direction of the ECM distribution gradient. Hence, the weak formulation (3.5a) of the macroscopic problem becomes(∂c ∂t ,ϕ ) + ( D1∇c,∇ϕ ) − ( ηc∇v,∇ϕ ) + δc ( η∇v ·∇c,∇v ·∇ϕ ) − ( µ1 c(1 − c−v),ϕ ) = 0, (3.9) with a stabilization parameter δc that has been heuristically choosen testing a range of values until the spurious oscillations has been reduced on the choosen refinenemnt level of the computational meshes. The transport equation is defined on Y . Since this is a hyperbolic equation, a suitable discretization is needed. Here we have chosen the streamline diffusion approach for its easy implementation and good performance. We have used an artificial diffusion in the streamline direction scaled with a parameter δ > 0, whose value can be found in Table 2. Therefore, the discretisation of the weak formulation (3.6) for the transport equation is given by( φn+1h ,ϕ ) Y + k ( V nh ·∇φ n+1 h ,ϕ + δ (V n h ·∇ϕ) ) Y = ( φnh,ϕ ) Y ∀ϕ ∈ UhT , φ0h = φ0 ∀x ∈ Y, (3.10) where φ0 is the initial level set function, k the time step and V n h is the discrete velocity defined as V nh := cvel ∆T ε Ix ( Iτ ( m̂h,n∇m̂h,n )) , (3.11) where Ix and Iτ are two quadrature formulas for the approximation of the integral in space and time (see expression (3.4)) and m̂h,n is the discrete solution of the microscopic problem as defined below. Since we assume scale separation, in each point of the macroscopic boundary ∂Ω(t) we need to solve a microscopic problem that defines the local velocity. In the discrete version, we define the microscopic problem in a finite number of points at the interface and we discuss later the issue of how to use these pointwise defined velocities to solve the transport problem. Then the discretisation of weak formulation of the microscopic problem (3.7) reads:( m̂l+1h,n,ϕ ) + kτ ( D2ε −2L−2∇m̂l+1h,n,∇ϕ ) = ( m̂lh,n,ϕ ) + ( F̂x,n,h,ϕ ) ∀ϕ ∈ Uhm, m̂0h,n = 0 in (0, 1), (3.12) where kτ = τn+1 − τn is the time step and m̂lh,n indicates the discrete microscopic solution for the macroscopic step n (note that the right hand side depends on the macroscopic solution at time tn), with l being the time step of the time variable τ, i.e. m̂lh,n = m̂h,n(τl). The term F̂x,n,h is an approximation of F̂x,tn in which c is substituted by its discrete counterpart and the integral over B is A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 45 approximated by a quadrature rule F̂x,n,h(z) :=   IB(c n h) IB(1) z ∈ [0, 1/2], 0 otherwise, (3.13) where IB(·) is a quadrature rule that approximates the integral of the argument over B IB(f) ≈ ∫ B f(ξ) dξ. 3.5. Cut-cells finite element approach in approximating the macroscopic tumour interface. As introduced previously, we discretize in time the tumour macro-dynamics and the transport equation for the tumour boundary relocation. Therefore, in the semidiscrete formulation we have terms that are defined at time t = tn+1 and terms defined at time t = tn. Since the domain is time dependent, the integrals of these terms are defined on different domains. Therefore, in each time step we need to consider two configurations defined by the position of the boundary ∂Ω(t) in two subsequent time steps. In particular, we have to consider the case in which the boundary cuts the same finite element cell in both time steps, see in Figure 6a the cell at the bottom left, and the case in which the boundary cuts one cell at time tn and it goes over to the neighbour cells at time tn+1 leaving the previous cell uncut at time tn+1, see in Figure 6a the cell at bottom right. For the discretized version of the system of equations, we consider the linearized domain Ωh(t), which is defined by the piece-wise linear boundary ∂Ωh(t) := L0,h(t), (3.14) where the linearized zero level L0,h(t) is defined by the polygonal line that connects all intersections of the zero level L0(t) with the mesh cells boundaries shown in Figure 6b. Cell cut twice Ωh(tn+1) Ωh(tn) (a) ∂Ωh(t) ∂Ωh(t) ∂Ω(t) (b) Figure 6. Left: Sketch of the linearized zero level L0,h. Right: Cell cut twice by the interface at two subsequent time steps. Using the linearized boundary ∂Ωh(t) we can apply the quadrature rule described in [16] to integrate the terms of the model on cut-cells with a single cut. Furthermore, if a cell is cut twice, i.e. at time tn and at time tn+1, we apply the previous quadrature rule recursively. 46 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU Approximation of the nonlocal term. The nonlocal term (3.13) is approximated by a quadrature rule. We use a further level set function to define the distance from the macroscopic point x, which determines the domain of integration B. Also in this case we introduce a piece-wise linear approximation of this level set function and integrate the cells using the quadrature rule only for the cell portion contained in B. 3.6. Extension of the micro-scale induced velocity. The velocity induced by the boundary MDE micro-dynamics is determined using formula (3.4). It is computed at the macroscopic point x on the linearized boundary ∂Ωh(t) and then extended to the rest of the domain Ωh(t). In this work, we consider only one point x per finite element cell. This is taken at the midpoint of the segment of the interface L0,h that intersects the cell, as shown in Figure 1. We set the velocity computed at this point x to all cells which center lies at the closest distance from x. Therefore, the velocity is approximated as a piecewise constant function. For finite element cells that lie in the cancer region, i.e. K ∩ Ωh(t) 6= 0, at a distance larger than a prescribed radius of influence ρ> 0, we set the velocity to zero. This enables us to avoid the transport of numerical pollution from the center of the domain due to the singularity of the level set in the point that we take as reference to compute the distance function. This definition of the velocity extension can lead to regularity problems in the transport of the interface if two parts of the boundary approach each other. This happens because the velocity of cells that lie at the same distance from the two approaching boundary parts is not well defined. This is a typical problem in level set approaches that can be overcome using a fast marching method [2]. 3.7. Overall numerical solution process. We sketch the overall solution process underlying the coupling between the different parts of the model. Algorithm 1 Overall solution process 1: Set n = 0 and choose the splitting time step ∆T 2: Set φ0(x), c0(x) and v0(x) in Y 3: Define L0(t n) as in (2.7) and linearize it to get Ωh(t n) 4: Solve macroscopic part (3.8) for (x,t) ∈ Ωh(tn) × (tn, tn + ∆T) 5: Compute F̂x,n,h(z), see (3.13) 6: Solve microscopic part (3.12) for (x,τ) ∈ εY × (0, ∆T) 7: Compute velocity V nh , see (3.11) 8: Extend velocity on all Y 9: Solve transport problem (2.6) for (x,t) ∈ Y × (tn, tn + ∆T) 10: if tn + ∆T = T then 11: stop 12: else 13: Set tn = tn + ∆T 14: goto 3 The macroscopic system is solved with an implicit Euler scheme. At each time step a nonlinear system of the type (3.8a-3.8b) has to be solved. We use an exact Jacobian and no damping for the Newton method, which converges generally in 2 steps to an accuracy lower than 10−6. The linear system arising in each Newton step is solved by the direct solver UMFPACK [24]. The system (3.10) is a linear system and is computationally much cheaper than the macroscopic problem. The direct solver is used in every time step. Finally, the microscopic problem (3.12) is a linear one dimensional parabolic problem solved with an implicit Euler method and the direct solver in each time step. A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 47 4. Numerical results In this section we show some numerical results obtained with the numerical method explained above. We have used the following initial condition for cancer cells c0(x) =   R−‖(x1,x2) − (2, 2)‖2 10 R if ‖(x1,x2) − (2, 2)‖2 < R, 0 otherwise, (4.1) (4.2) where R is the initial radius of the cancer region and the point (2, 2) is the center of the computational domain, see Table 2 for the numerical parameters. For the extracellular matrix distribution we consider two cases: • Case 1: v0(x) = 0.3 sin (3π‖(x1,x2) − (0, 0)‖3) + 0.5. • Case 2: v0(x) = 0.15 (sin (2π‖(x1,x2) − (0, 0)‖3) + sin (2π‖(x1,x2) − (4, 0)‖3)) + 0.75. In Figure 7 are shown the results for Case 1. In figures 7a, 7c and 7e is depicted the distribution of cancer cells at time t = 0, t = 5 and t = 10 respectively, while in figures 7b, 7d and 7f is depicted the ECM distribution at the corresponding time points. The initial tumour mass is located at a valley of the ECM distribution. The white circle in Figure 7b shows the initial cancer cells mass position with respect to the ECM distribution. The coordinates of the center of the initial distribution are in Table 1. We observe that the cancer cells population haptotaxis biases the macro-dynamics of cancer cells towards regions with higher ECM levels, which in turn leads to the formation of stronger MDEs sources for the micro-process taking place on the boundary micro-domains εY situated on the part of the tumour interface facing those elevated ECM regions, leading to a stronger MDEs boundary micro-dynamics that ultimately translates into a larger magnitude velocity that is induced from the micro-scale and is fed into the transport equation which governs the tumour boundary movement, that results into progression of the tumour in those regions. This is a truly multi-scale characteristic of the actual cancer invasion process that our model is able to capture, resulting in this pronounced lobular progression of the tumour. It can be observed at time t = 5 that the cancer cells distribution moves towards the maximum of ECM, and then it follows the path along the two ECM ridges besides the starting position, which confirms in silico the well-known process of durotaxis observed experimentally [38, 42, 57, 59, 60, 70, 79]. In Figure 8 are shown the results of Case 2. In the figures 8a, 8c and 8e is depicted the distribution of cancer cells at time t = 0, t = 5 and t = 10 respectively, while in figures 8b, 8d and 8f is depicted the ECM distribution at the corresponding time points. Also here the initial cacer cells distribution is set at the valley of the ECM distribution, see Table 1. In this case, in the vicinity of the starting position there are four peaks of ECM distribution and the haptotactic term induces a transport of cancer cells distribution towards all four peaks in a nonsymmetric manner, because the initial position is not equidistant to the peaks. Therefore, the results in this figure with changed surrounding tissue structure (induced by the different ECM initial pattern) confirm the same multiscale character of the tumour progression that our multiscale moving boundary model is able to capture. Again here, we observe this pronounced lobular progression underscoring again the durotaxis behaviour observed experimentally for the cancer cell migration. The behavior of the model is consistent with the expected behavior, in which an interplay of the microscopic and macroscopic processes determine the transport of the cancer cell distribution in a nonuniform environment with a nonconstant distribution of the ECM. Specifically, it can be observed that the haptotactic movement of the cancer cell distribution at macro-scale leads higher levels of sources of MDEs that the macro-dynamics induces at micro-scale for the boundary proteolytic processes, which 48 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU (a) (b) (c) (d) (e) (f) Figure 7. Case 1: Distribution of c2ncer cells at time t=0 (7a), t=5 (7c) and t=10 (7e). Distribution of ECM at time t=0 (7b), t=5 (7d) and t=10 (7f). in turn induces a movement of the macroscopic tumourinterface progression in pronounced lobular manner toward the higher values of the ECM distribution. A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 49 (a) (b) (c) (d) (e) (f) Figure 8. Case 2: Distribution of cancer cells at time t=0 (7a), t=5 (7c) and t=10 (7e). Distribution of ECM at time t=0 (7b), t=5 (7d) and t=10 (7f). 50 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU final time T 10 microdynamics time range ∆T 0.1 initial radius of cancer distribution R 0.3 scale factor ε 0.01 diffusion cancer cells D1 0.00043 haptotactic coefficient η 0.2 proliferation coefficient µ1 0.25 ECM remodelling µ2 0.15 degradation α 1.5 diffusion MDE D2 0.001 center of init. cancer distr. Case 1 (2.3,2.2) center of init. cancer distr. Case 2 (2.0,1.9) Table 1. Model parameters in our numerical experiments time step k 0.1 mesh size h 0.0078125 stream-line stabilization transport problem δ 0.5 stream-line stabilization macroscopic problem δc 0.004 total tissue domain Y (0, 4) × (0, 4) Table 2. Details of the numerical setting 5. Conclusions We presented a new formulation of a two-scale model for the simulation of cancer invasion. This included a new derivation of the tumor boundary motion law by considering the contributions of MDE micro-dynamics within a transport equation (2.6) (via the velocity field V (m) that is induced by the micro-scale MDEs processes at tumour interface), the solution of which provides the level set function indicating the cancer boundary progression on the macro scale. For the computational implementation we used an unfitted regular mesh with uniform cell diameters and a cut-cell finite element formulation to avoid the problem of re-meshing in case of large deformations. We have shown that the presented framework is highly flexible to study different aspects of the cancer invasion process. In particular, since it is important to study the interplay between the two scales, the presented implementation allows high flexibility in defining the strength of the coupling via the definition of the velocity field. The effect on the velocity of the tumour boundary from multiple microscopic substrates such as matrix metalloproteinases (MMPs) and urokinasetype plasminogen activators (uPAs) can be studied. Also the level of complexity can be easily increased given the efficiency of the model implementation and the explicit link between the scales. Several numerical aspects are of interest for further work in this framework. The solution of transport equation with finite elements can be substituited by a fast marching algorithm designed for the level set approach as explained in Section 3.6. The stabilization of the macroscopic problem introduced here with a streamline diffusion technique can be improved using (higher order) stabilization techniques based on flux-corrected finite element approaches as indicated in Section 3.4. Furthermore, an adaptive local refinement strategy for moving meshes can be adopted to reduce the computational costs and increase the accuracy in the vicinity of the interface. Finally, we underline the potential of the presented A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 51 method, which allows to go to three-dimensional problems without changing the numerical formulation, thus allowing a significant development of this multiscale modeling framework. The most important extension required for this development is the formulation of cut cells in three dimensions. References [1] A. Abdulle, E. Weinan, B. Engquist and E. Vanden-Eijnden, The heterogeneous multiscale method, Acta Numerica 21 (2012), 1–87. [2] D. Adalsteinsson and J. Sethian, The fast construction of extension velocities in level set methods, Journal of Computational Physics 148 (1999), 2 – 22. [3] J. A. Adam, A simplified mathematical model of tumour growth, Math. Biosci. 81 (1986), 229–244. [4] A. Alsisi, R. Eftimie and D. Trucu, Non-local multiscale approaches for tumour-oncolytic virus interactions, Mathe- matics in Applied Sciences and Engineering 1 (2020), 249–273. [5] T. Alzahrani, R. Eftimie and D. Trucu, Multiscale modelling of cancer response to oncolytic viral therapy, Mathe- matical Biosciences 310 (2019), 76–95. [6] T. Alzahrani, R. Eftimie and D. Trucu, Multiscale moving boundary modelling of cancer interactions with a fusogenic oncolytic virus: The impact of syncytia dynamics, Mathematical Biosciences 323 (2020), 1–22. [7] V. Andasari, A. Gerisch, G. Lolas, A. South and M. A. J. Chaplain, Mathematical modeling of cancer cell invasion of tissue: biological insight from mathematical analysis and computational simulation, J. Math. Biol. 63 (2011), 141–171. [8] A. Anderson, M. Chaplain, E. Newman, R. Steele and A. Thompson, Mathematical modelling of tumour invasion and metastasis, J. Theor. Med. 2 (2000), 129–154. [9] A. R. A. Anderson, A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion, Math. Med. Biol. 22 (2005), 163–186. [10] N. J. Armstrong, K. J. Painter and J. A. Sherratt, A continuum approach to modelling cell-cell adhesion, J. Theor. Biol. 243 (2006), 98–113. [11] D. Arndt, W. Bangerth, D. Davydov, T. Heister, L. Heltai, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin and D. Wells, The deal.II library, version 8.5, Journal of Numerical Mathematics 25 (2017), 137–146. [12] E. Burman, Ghost penalty, Comptes Rendus Mathematique 348 (2010), 1217 – 1220. [13] H. Byrne, M. Chaplain, G. Pettet and D. L. S. Mcelwain, A mathematical model of trophoblast invasion, Appl. Math. Lett. 14 (2001), 1005–1010. [14] H. Byrne and L. Preziosi, Modelling solid tumour growth using the theory of mixtures, Math. Med. Biol. 20 (2003), 341–366. [15] H. M. Byrne and M. A. Chaplain, Modelling the role of cell-cell adhesion in the growth and developement of carci- noma, Math. Comput. Model. 24 (1996), 1–17. [16] T. Carraro and S. Wetterauer, On the implementation of the extended finite element method (xfem) for interface problems, Archive of Numerical Software 4 (2016), 1–23. [17] M. Chaplain, L. Graziano and L. Preziosi, Mathematical modelling of the loss of tissue compression responsiveness and its role in solid tumour development, Math. Med. Biol. 23 (2006), 197–229. [18] M. Chaplain, M. Lachowicz, Z. Szymanska and D. Wrzosek, Mathematical modelling of cancer invasion: The impor- tance of cell-cell adhesion and cell-matrix adhesion, Math. Model. Meth. Appl. Sci. 21 (2011), 719–743. [19] M. Chaplain and G. Lolas, Mathematical modelling of cancer cell invasion of tissue: the role of the urokinase plasminogen activation system, Math. Model. Meth. Appl. Sci. 15 (2005), 1685–1734. [20] M. Chaplain, S. McDougal and A. Anderson, Mathematical modeling of tumor-induced angiogenesis, Annu. Rev. Biomed. Eng. 8 (2006) 233–257. [21] Y. Chen, S. M. Wise, V. B. Shenoy and J. S. Lowengrub, A stable scheme for a nonlinear, multiphase tumor growth model with an elastic membrane, Int. J Num. Meth. Biomed. Eng. 30 (2014), 726–754. [22] L.-T. Cheng and E. Weinan, The heterogeneous multi-scale method for interface dynamics, Contemporary Mathe- matics 330 (2003), 43–54. [23] D. L. Chopp, Another look at velocity extensions in the level set method, SIAM Journal on Scientific Computing 31 (2009), 3255–3273. [24] T. A. Davis, Algorithm 832: Umfpack v4. 3—an unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software (TOMS) 30 (2004), 196–199. [25] N. E. Deakin and M. A. J. Chaplain, Mathematical modelling of cancer cell invasion: the role of membrane-bound matrix metalloproteinases, Front. Oncol. 3 (2013), 1–9. 52 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU [26] P. Domschke, D. Trucu, A. Gerisch and M. Chaplain, Mathematical modelling of cancer invasion: Implications of cell adhesion variability for tumour infiltrative growth patterns, J. Theor. Biol. 361 (2014), 41–60. [27] W. E, Principles of Multiscale Modeling, Cambridge University Press, Cambridge, 2011. [28] W. E and B. Engquist, The heterognous multiscale methods, Commun. Math. Sci. 1 (2003), 87–132. [29] M. Egeblad, E. S. Nakasone and Z. Werb, Tumors as organs: Complex tissues that interface with the entire organism, Develop. Cell 18 (2010) 884–901. [30] S. Frei and T. Richter, Efficient approximation of flow problems with multiple scales in time, Multiscale Modeling & Simulation 18 (2020), 942–969. [31] H. B. Frieboes, F. Jin, Y.-L. Chuang, S. M. Wise, J. S. Lowengrub and V. Cristini, Three-dimensional multispecies nonlinear tumor growth—ii: Tumor invasion and angiogenesis, J. Theor. Biol. 264 (2010), 1254–1278. [32] H. B. Frieboes, X. Zheng, C.-H. Sun, B. Tromberg, R. Gatenby and V. Cristini, An integrated computa- tional/experimental model of tumor invasion, Cancer Res. 66 (2006), 1597–1604. [33] R. A. Gatenby and E. T. Gawlinski, A reaction-diffusion model of cancer invasion, Cancer Res. 56 (1996), 5745–5753. [34] A. Gerisch and M. Chaplain, Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion, J. Theor. Biol. 250 (2008), 684–704. [35] H. P. Greenspan, On the growth and stabiligy of cell cultures and solid tumours, J. Theor. Biol. 56 (1976), 229–242. [36] D. Hanahan and R. A. Weinberg, Hallmarks of cancer: the next generation, Cell 144 (2011), 646–674. [37] T. Hillen, K. J. Painter and M. Winkler, Convergence of a cancer invasion model to a logistic chemotaxis model, Mathematical Models and Methods in Applied Sciences 23 (2013), 165–198. [38] B. C. Isenberg, P. A. DiMilla, M. Walker, S. Kim and J. Y. Wong, Vascular smooth muscle cell durotaxis depends on substrate stiffness gradient strength, Biophysical Journal 97 (2009), 1313–1322. [39] J. A. Joyce and J. Pollard, Microenvironmental regulation of metastasis, Nat. Rev. Cancer. 9 (2009), 239–252. [40] R. Kalluri and M. Zeisberg, Fibroblasts in cancer, Nat. Rev. Cancer. 6 (2006), 392–401. [41] D. Kuzmin, On the design of algebraic flux correction schemes for quadratic finite elements, J. Comput. Appl. Math. 218 (2008), 79–87. [42] C.-M. Lo, H.-B. Wang, M. Dembo and Y.-L. Wang, Cell movement is guided by the rigidity of the substrate, Biophysical Journal 79 (2000), 144–152. [43] P. Macklin and J. Lowengrub, Evolving interfaces via gradients of geometry-dependent interior poisson problems: application to tumor growth, Journal of Computational Physics 203 (2005), 191–220. [44] P. Macklin and J. Lowengrub, An improved geometry-aware curvature discretization for level set methods: application to tumor growth, journal of Computational Physics 215 (2006), 392–401. [45] P. Macklin and J. Lowengrub, Nonlinear simulation of the effect of microenvironment on tumor growth, J. Theor. Biol. 245 (2007), 677–704. [46] P. Macklin and J. Lowengrub, A new ghost cell/level set method for moving boundary problems: Application to tumor growth, J. Sci. Comput. 35 (2008) 266 – 299. [47] P. Macklin, S. McDougall, A. R. Anderson, M. A. Chaplain, V. Cristini and J. Lowengrub, Multiscale modelling and nonlinear simulation of vascular tumour growth, J. Math. Biol. 58 (2009), 765–798. [48] S. Mumenthaler, G. D’Antonio, L. Preziosi and P. Macklin, The need for integrative computational oncology: an illustrated example through mmp-mediated tissue degradation, Frontiers in Oncology 3 (2013): 194. [49] C. F. Ng and H. B. Frieboes, Model of vascular desmoplastic multispecies tumor growth, J. Theor. Biol. 430 (2017), 245–282. [50] C. F. Ng and H. B. Frieboes, Simulation of multispecies desmoplastic cancer growth via a fully adaptive non-linear full multigrid algorithm, Frontiers in Physiology 9 (2018) 821–821. [51] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations, Journal of Computational Physics 79 (1988), 12–49. [52] L. Peng, D. Trucu, P. Lin, A. Thompson and M. A. J. Chaplain, A multiscale mathematical model of tumour invasive growth, Bull. Math. Biol. 79 (2017), 389–429. [53] L. Peng, D. Trucu, P. Lin, A. Thompson and M. A. J. Chaplain, A multiscale mathematical model of tumour invasive growth, Bull. Math. Biol. 79 (2017) 389–429. [54] A. Perumpanani, J. Sherratt, J. Norbury and H. Byrne, Biological inferences from a mathematical model for malig- nant invasion, Invas. Metast. 16 (1996), 209–221. [55] A. Perumpanani, D. Simmons, A. Gearing, K. Miller, G. Ward, J. Norbury, M. Schneemann and J. Sherratt, Extra- cellular matrix-mediated chemotaxis can impede cell migration, Proc. Roy. Soc. Lond. B 265 (1998), 2347–2352. [56] L. Preziosi and A. Tosin, Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications, J. Math. Biol. 58 (2009), 625–656. A LEVEL SET APPROACH FOR A MULTI-SCALE CANCER INVASION MODEL 53 [57] P. P. Provenzano, D. R. Inman, K. W. Eliceiri, J. G. Knittel, L. Yan, C. T. Rueden, J. G. White and P. J. Keely, Collagen density promotes mammary tumour initiation and progression, BMC Med. 6(2008): 11. [58] B.-Z. Qian and J. W. Pollard, Macrophage diversity enhances tumor progression and metastasis, Cell 141 (2010), 39–51. [59] M. Raab, J. Swift, P. D. P. Dingal, P. Shah, J.-W. Shin and D. E. Discher, Crawling from soft to stiff matrix polarizes the cytoskeleton and phosphoregulates myosin-II heavy chain, The Journal of Cell Biology 199 (2012), 669–683. [60] A. Saez, M. Ghibaudo, A. Buguin, P. Silberzan and B. Ladoux, Rigidity-driven growth and migration of epithelial cells on microstructured anisotropic substrates, Proceedings of the National Academy of Sciences 104 (2007), 8281–8286. [61] G. Sciumè, S. Shelton, W. Gray, C. Miller, F. Hussain, M. Ferrari, P. Decuzzi and B. Schrefler, A multiphase model for three-dimensional tumor growth, New Journal of Physics 15 (2013), 1–43. [62] R. Shuttleworth and D. Trucu, Two-scale Moving Boundary Dynamics of Cancer Invasion: Heterotypic Cell Pop- ulations Evolution in Heterogeneous ecm, in Cell Movement Modelling and Applications, eds. M. Stolarska and N. Tarfulea (Birkhäuser, Springer Nature Switzerland AG, 2018), pp. 1-26. [63] R. Shuttleworth and D. Trucu, Multiscale modelling of fibres dynamics and cell adhesion within moving boundary cancer invasion, Bull. Math. Biol. 81 (2019), 2176–2219. [64] R. Shuttleworth and D. Trucu, Cell-scale degradation of peritumoural extracellular matrix fibre network and its role within tissue-scale cancer invasion, Bull. Math. Biol. 82 (2020): 65. [65] R. Shuttleworth and D. Trucu, Multiscale dynamics of a heterotypic cancer cell population within a fibrous extracel- lular matrix, J. Theor. Biol. 486 (2020): 110040. [66] R. Strehl, A. Sokolov, D. Kuzmin, D. Horstmann and S. Turek, A positivity-preserving finite element method for chemotaxis problems in 3d, Journal of Computational and Applied Mathematics 239 (2013), 290–303. [67] S. Suveges, R. Eftimie and D. Trucu, Directionality of macrophages movement in tumour invasion: A multiscale moving-boundary approach, Bull. Math. Biol. 82 (2020), 1–50. [68] Z. Szymanska, C. M. Rodrigo, M. Lachowicz and M. A. J. Chaplain, Mathematical modelling of cancer cell invasion of tissue: the role and effect of nonlocal interractions, Mathematical Models and Methods in Applied Sciences 19 (2009), 257–281. [69] A. Tosin and L. Preziosi, Multiphase modeling of tumor growth with matrix remodeling and fibrosis, Mathematical and Computer Modelling 52 (2010), 969–976. [70] L. Trichet, J. L. Digabel, R. J. Hawkins, S. R. K. Vedula, M. Gupta, C. Ribrault, P. Hersen, R. Voituriez and B. Ladoux, Evidence of a large-scale mechanosensing mechanism for cellular adaptation to substrate stiffness, Proc. Nat. Acad. Sci. 109 (2012), 6933–6938. [71] D. Trucu and M. A. J. Chaplain, Multiscale analysis and modelling for cancer growth and development, in Managing Complexity, Reducing Perplexity, eds. M. Delitala and G. Ajmone Marsan (Springer International Publishing, 2014), pp. 45-53. [72] D. Trucu, P. Domschke, A. Gerisch and M. A. J. Chaplain, Multiscale computational modelling and analysis of cancer invasion, in Mathematical Models and Methods for Living Systems: Levico Terme, Italy 2014, eds. L. Preziosi, M. A. J. Chaplain and A. Pugliese (Springer International Publishing, 2016), pp. 275–321. [73] D. Trucu, P. Lin, M. A. Chaplain and Y. Wang, A multiscale moving boundary model arising in cancer invasion, Multiscale Modeling & Simulation 11 (2013), 309–335. [74] S. Webb, J. Sherratt and R. Fish, Alterations in proteolytic activity at low ph and its association with invasion: a theoretical model, Clin. Experim. Metast. 17 (1999), 397–407. [75] R. A. Weinberg, The Biology of Cancer, Garland Science, New York, 2006. [76] M. Winkler, Boundedness in the higher-dimensional parabolic- parabolic chemotaxis system with logistic source, Communications in Partial Differential Equations 35 (2010), 1516–1537. [77] S. Wise, J. Lowengrub, H. Friebose and V. Cristini, Three-dimensional multispecies nonlinear tumor growth—I. Model and numerical method, J. Theor. Biol. 253 (2008), 524–543. [78] S. M. Wise, J. Lowengrub and V. Christini, An adaptive multigrid algorithm for simulating solid tumor growth using mixture models, Math. Comput. Model. 53 (2011), 1–20. [79] M. H. Zaman, L. M. Trapani, A. L. Sieminski, D. MacKellar, H. Gong, R. D. Kamm, A. Wells, D. A. Lauffenburger and P. Matsudaira, Migration of tumor cells in 3d matrices is governed by matrix stiffness along with cell-matrix adhesion and proteolysis, Proc. Nat. Acad. Sci. 103 (2006), 10889–10894. [80] X. Zheng, S. Wise and V. Cristini, Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method, Bull. Math. Biol. 67 (2005), 211–259. [81] G. Zhou, How accurate is the streamline diffusion finite element method?, Mathematics of Computation 66 (1997), 31–44. 54 T. CARRARO, S.E. WETTERAUER, A.V. PONCE BOBADILLA, AND D. TRUCU Corresponding author. Faculty of Mechanical Engineering, Applied Mathematics, Helmut Schmidt Univer- sity / University of the Federal Armed Forces Hamburg, Holstenhofweg 85, 22043 Hamburg, Germany. E-mail address: carraro@hsu-hh.de Institute for Applied Mathematics, Heidelberg University, Im Neuenheimerfeld 205, 69120 Heidelberg, Germany. Institute for Applied Mathematics, Heidelberg University, Im Neuenheimerfeld 205, 69120 Heidelberg, Germany. Division of Mathematics, University of Dundee, Dundee, DD1 4HN, United Kingdom. E-mail address: trucu@maths.dundee.ac.uk