A level-set approach for a multi-scale cancer invasion model

Central to the quest for a deeper understanding of the cancer growth and spread process, the naturally multiscale character of cancer invasion demands appropriate multiscale modelling and analysis approach. The cross-talk between the tissue scale (macro-scale) cancer cell population dynamics and the cell-scale (micro-scale) proteolytic molecular processes along the tumour boundary plays a particularly important role within the invasion processes, leading to dramatic changes in tumour 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(1): 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.


Introduction
Involving a wide range of cross-related processes occurring on several spatiotemporal scales, cancer cell invasion of the human is one of the hallmarks of cancer [32], playing a crucial role in the overall development a spread of a growing malignant tumour. Taking advantage of the heterotypic character of the tumour microenvironment (which includes immuno-inflamatory cells, stromal cell, fibroblasts, endothelial cells, macrophages), complex molecular processes facilitate intense interactions between the cancer cell population and the extracellular matrix (ECM) [23,32,36,37,52]. 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.
Focusing on local tumour progression, the alteration and remodelling of the ECM by the matrix degraded enzymes (MDEs) such as matrix metalloproteinases MMPs or the urokinase plasminogen activator (uPA) play a key role. 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 partial differential equations [2,3,4,5,9,11,15,16,21,28,30,48,49,57] as well as nonlocal integrodifferential systems [6,14,22,29] were derived and proposed to deepen the understanding, validate and create new experimental hypothesis. Furthermore, to capture various heterotypic aspects and related processes within tumour invasion, several multiphase models based on the theory of mixtures [10,13,25,51,59,60] 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 microenvironment. This gives the cancer invasion a moving boundary character, and to that end several levelset approaches were recently proposed to study the tumour progression both in homogeneous environments [26,38,39,40,61] and in complex heterogeneous tissues [41].
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 [56]. 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 [47].
In this work we present a new formulation of the link between the two scales presented in [56]. 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 [7], facilitates the realization of the code.
The level set method was first introduced in [46] for tracking moving interface with complex deformations. This method was developed starting from the notion of weak solutions for evolving interfaces. The main aspect of this method is the fact 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.
Other works have presented a level set approach for moving the tumour interface. In [38] 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 [61] 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 [42] 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 [56] 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 2 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 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.

Problem description
We present a two-scale model for cancer invasion that links through a double feedback loop the dynamics occurring at two different spatial scales explored by the following two modelling components, namely: 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 at cell-scale. Both scales are considered at the continuum level. We therefore assume that possible stochastic effects, in regions where the continuum assumption is not more valid, can be neglected. Nevertheless, our aim is to derive a flexible numerical framework that would allow to extend the model with a stochastic part (e.g. at the interface of the domain) leading to a hybrid formulation, if needed.
The cancer cells and extracellular matrix are modelled in an invading domain Ω(t) that changes its size and morphology in time during the invasion process within a reference maximal tissue cube Y . Its boundary ∂Ω(t) will also be referred to as interface because this is regarded here as an interface separating the region with zero cancer cells and a given distribution of ECM from a region with a distribution of cancer cells and ECM that satisfies the macroscopic equations. Finally, in appropriate cell-scale neighbourhoods of the interface ∂Ω(t) points, a microscopic problem describing the cross-interface transport of matrix-degrading enzymes is considered and accounted for in order to determine the law for macro-scale boundary movement. Assumption 2.1 (Scale separation) We assume scale separation in space between the two components of the model, i.e. the macroscopic populations of cancer cells and ECM and the microscopic population of matrix degrading enzyme molecules.
The characteristic length L for the macroscopic part of the model relates to the diameter of the cancer region and is considered here as in [4,29], ranging between 0.1cm and 1.0cm. Further, 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 −3 cm, [45]. The ratio between the scales is denoted ε = /L. Therefore, the enzyme population's dynamics can be described in a bundle of microscopic domains εY that are obtained by scaling a reference tissue cubic domain Y by a factor ε and whose union provide a cell-scale neighbourhood for the interfacial points in ∂Ω(t). Due to spatial scale separation, we consider a microscopic problem at each point x ∈ ∂Ω(t) (i.e., at each point of the macroscopic interface) on the corresponding εY micro-domain centred at x. Further assuming for convenience that the maximal cube Y is centred at origin of the space, the micro-scale coordinates y of the microscale problem on a εY centred at x 0 are obtained by appropriate scaling and translation of the macroscopic coordinates, namely as y = x 0 +ε(x−x 0 ). This approach is similar to the heterogeneous multiscale method for interface dynamics presented in [17], but we do not consider an upscaling process.
The first step of this work is to describe the whole coupled problem at the continuum level and then in a second step to discretize it. In particular, we consider a semidiscretization in time by the implicit Euler method and the finite elements method for the discretization in space. The reasons for this choice will be clarified later in the section dedicated to the discretization.
Schematically the two scales are coupled by the following double feedback loop: • The top-down macroscopic-to-microscopic coupling for the microscopic problem at any boundary point x ∈ Ω(t) is done via the source of matrix-degrading enzymes which is induced by the macro-dynamics and is formed as a collective contribution of the cancer cells that arrive during the macro-dynamics within an appropriate distance from x, see ; • The bottom-up microscopic-to-macroscopic coupling is done by defining the velocity of the macroscopic interface in dependence of the enzymes concentration in the microscopic domains, see (9).
Therefore, the rate of cancer cells invasion into the surrounding tissues is defined by the velocity of the interface that depends on the solution of the microscopic enzyme dynamics. The interface is described by the zero-level of a level-set function that has to be initialized with the initial distribution of the cancer cells, see Figure 1 and formula (20), and that is moved by a transport equation. To account for all interactions between the different parts of the problem we introduce the domain Y ⊂ R 2 that is assumed to be sufficiently large such that the complete dynamics happen inside it. This domain is used at the continuum level to define the region in which the transport equation is solved. Later it is used for the discretized problem to define a region where the finite elements become active if their intersection with the cancer region is not empty.
In the following we proceed with the multiscale model description in three parts: macroscopic, microscopic and transport component. At macro-scale, the model considers cancer cells and extracellular matrix (ECM) interaction.

Macroscopic model component
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 [56], the dynamics at macroscopic scale is given by the following PDE system: where D 1 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. The typical values of these coefficients are reported in Table 2. Except for µ 1 , all other coefficients are considered constant. However, due to the direct dependency of the mitotic process on the presence of ECM [18,50,55], the proliferation coefficient µ 1 is a function of v which decreases monotonically with the amount of ECM, exploring the full range of values between 1 and 0, as v varies from its maximal levels, which are ideal for proliferation, to 0-level ECM regions, where cell-death and necrosis occurs due to lack of anchorage and nutrients. It is assumed that the cancer cells are zero outside Ω(t) and that there is no flux of cells through the boundary ∂Ω(t). Furthermore, under the presence of these boundary and initial conditions, for the case of constant proliferation rate µ 1 , the results in [34,54,58] explore the local and global existence of system (1)- (2). The initial distribution of cancer cells c 0 (x) and extracellular matrix v 0 (x) are given in the larger domain Y .
Finally, we note that while this macro-dynamics occurs at tissue-scale, the macroscale cancer invasion is completely controlled by the movement of the tumour interface which is regulated by the macroscopic/microscopic interaction that defines the interface velocity, whose details will be discussed in the next subsections.
At this stage, it is important to note that the macroscale part (1)-(2) of the model can be regarded as a singularly perturbed diffusion-reaction-transport problem of the type with γ = 0. In particular, we observe that the second equation is reaction dominated and, in this context, it is well known that this kind of problems could present one or more boundary layers [35,53]. The thickness of the boundary layers give these problems a multiscale character in its own right. Classical techniques to solve convection/reaction dominated problems are discontinuous and continuous interior penalty methods [8]. Other techniques are streamline-upwind/Petrov-Galerkin formulations [27]. Furthermore, to overcome numerical instabilities other methods have been developed like the Galerkin enriched finite element method (GEM) [24]. A typical problem in case of internal or boundary layers is the missing a priori knowledge of the position of such layers. Nevertheless, using a range of model parameters limited to the values used in literature for our model, we have observed that the solution component v has a mild boundary layer, due to the degradation of the ECM. Therefore, for the simulations shown here no stabilization technique was necessary and the mesh has been refined for accuracy purposes and not for stability issues.
In the next section we introduce the part of the model used to update the domain in time.

Macroscopic time-dependent domain
As mentioned already 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 microdynamics occurring on a small micro-domain εY centred at x. Hence, the velocity field that arises this way on ∂Ω(t) depends on the micro-dynamics MDE molecular distribution m(y, τ ) over an appropriate micro-spatio-temporal domain εY × (0, ∆T ) (which will be detailed in Section 2.3). Therefore, we will denote this velocity field by V (m) and remark at this stage that this establishes a bottom-up feedback link from micro-to macro-dynamics, being directly responsible for the movement of the tissue scale tumour boundary ∂Ω(t).
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: For later purposes, we introduce the notation 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 [20]. In Section 3.5 more details about this point are given.

Microscopic model component
Due to the scale separation Assumption 2.1 we can describe the dynamics of the MDE on a microscopic domain εY defined and centred at each macroscopic interface point x ∈ ∂Ω(t) as follows. As argued in [56], the tumour cells arriving within a certain radius R m from the interface location x ∈ ∂Ω(t) give rise to a source of MDE, which simply represents the collective secretion of matrix degrading enzymes by the cells from the outer proliferating rim of the tumour that get distributed during their macro-dynamics on Thus, mathematically, this source can be formalised as: Therefore, in the presence of source (7), the cross-interface micro-dynamics of MDE molecular distribution m(y, τ ), which takes place on the micro-scale domain εY over a micro-scale time range (0, ∆T ), is governed by the following reaction diffusion equation with ∆T representing here the micro-scale time perspective and serving also later on as natural time splitting step between micro and macro-stages within the computational approach of the multiscale model. Finally, we would like to remark here the special multiscale importance of the microscale source term given via the non-local operator F x,t (c) that is induced by the macro-dynamics and realises a top-down link from macro-to micro-dynamics, enabling this way the entire interface micro-dynamics.
Time splitting At this stage, the three parts of the model, i.e. the macroscopic, the transport and the microscopic components, are fully coupled and the coupling is given at the continuous level.
To define the velocity V used in equation (5) we consider a splitting in time of the overall coupling between the three model components. Therefore, for a given macroscopic time interval ∆T , 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 mathematically described by where |εY | = εY 1 dt and c vel 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 we ultimately obtain the definition of the velocity given in (9), where V (m) is taken as being proportional to this spatio-temporal mean value, with proportionality constant c vel .
In the implementation we use a linearized interface. For each macroscopic cell that is cut by the interface, we consider quadrature points at the interface and assign a velocity vector to each of them, see Figure 2 for a sketch of the interface.

Definition of the computational microscopic problem
To simplify the implementation of the numerical method we consider the microscopic problem in a bundle of boundary microdomains εY , where Y is the reference maximal tissue cube centred at origin, with y :=(y 1 , y 2 ) being the standard local microscale reference system within a given εY and τ denoting always the time at micro-scale.
As will be explained below in Section 3.3, in our finite element approach the macroscopic dynamics will be considered on a appropriately defined macroscopic domain Ω h (t) with a linearized boundary ∂Ω 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 y 1 ) and two edges orthogonal to it (in direction y 2 ) as shown in Figure 2. This orientation of the computational microscopic domain is well defined because it is centered in quadrature points at the boundary, which are never defined at the corners of the piece-wise linear boundary. This simplifies the setting of the microscopic problem. In fact, since we consider the linearized boundary ∂Ω h (t), the right hand side F x,t in equation (8) does not depend on y 1 . In addition, since on the boundaries of the quadrilateral domain no flux conditions are prescribed, it follows that the solution is constant in y 1 direction. Therefore, we can consider the following simplified one-dimensional microscopic problem for the quantity m, which is the integral of m along y 1 (giving the amount of enzyme molecules per unit of length) where F x,t is F x,t integrated over y 1 . Since F x,t and m do not depend on y 1 , the solution m of (8) is the constant extension of the solution m of (10) in y 1 direction.
We introduce now a scaling of the domain to the interval (0, 1) through the following tranformation y 2 = εz with z ∈ (0, 1), then we get after the rescaling the transformed system with note that the coordinate ξ is a macroscopic quantity. Notice furthermore that a solution of (10) is a solution of (11) by m(z, τ ) = m(εz, τ ).

Remark 2.1 (Limit ε → 0)
In case of ε → 0 we have in (11) 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 (11), we obtain that the velocity defined by problem (8) via equation (9) can be further expressed as where m indicates the transformed function on the reference domain Y .

Weak formulation and discretized model
To describe the model in the setting needed for the FEM we introduce the following weak formulation.

Weak formulation
We use the notation (·, ·) to define the usual L 2 scalar product of Lebesgue square integrable functions. The space H 1 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 H 1 . Furthermore, we use Bochner spaces like U = {u ∈ L 2 (0, T ; H 1 ) : ∂ t u ∈ L 2 (0, T ; H * )} to introduce the weak formulation of each subproblem. In particular, we consider the functional space for the transport component, the space for the macroscopic component and U m = {c ∈ L 2 (0, T ; H 1 ((0, 1))) : ∂ t c ∈ L 2 (0, T ; H * ((0, 1)))} for the scaled microscopic component. The weak formulation problem of the macroscopic part is: find the pair (c, v) ∈ U M × U M such that it satisfies for almost all t ∈ (0, T ) Note that the cancer cells c 0 and ECM distributions v 0 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. ∂ n c = ∂ n v = 0 on ∂Ω(t).
The weak formulation of the transport equation is: find φ ∈ U T such that for a.a. t ∈ (0, T ) it satisfies with φ 0 being the level set function at the initial time. Using the notation ϕ + := max{ϕ, 0} to define the positive part of a function ϕ in its domain of definition, we have that supp(φ + 0 ) describes the initial support region of the cancer cells. The weak formulation of the microscopic problem is: find m ∈ U m such that for a.a. τ ∈ (0, ∆T ) it satisfies where F x,t is defined as in (12) and the natural condition ∂ n m = 0 is implicitly defined.

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 M h composed of quadrilateral cells K ∈ M h of the same dimension. This discretization has in a natural way the property of shape regularity, i.e.
there exists a constant κ such that max where h K and ρ K are respectively the cell diameter and the diameter of the largest ball inscribed into K. Furthermore it has the advantage that the mesh can be generated starting from an initial quadrilateral that is refined successively to reach a given cell diameter.
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 Q 1 defined on a unit cellK = [0, 1] 2 , i.e. Q 1 = span (1, x, y, xy) and the space of linear functions defined on the unit one-dimensional cellK = [0, 1]. The finite element space is defined as where T K is a bijective transformation from the unit cellK to the physical cell K. The functions v are called pull-back functions and are the transformations to the real coordinates of the polynomial functions defined on the unit cells.
In our case, we consider only a translation and a scaling of the unit cells.
Since the mesh is non-fitted, the functions v (in case of cut-cells) are defined only on the portion of cell that is intersected by Ω(t). On the portion of the cell that lies outside of Ω(t) the shape functions need not to be defined. In Figure 3 a one dimensional sketch is shown, where one can see how the restriction of the shape functions ϕ 1 and ϕ 2 on the cut-cell is defined. Note that the shape functions are not modified, we just use a restriction of them on the part of the cell that belongs to the domain. The Lagrangian formulation remains the conventional one for continuous finite elements. The degrees of freedom are defined on the cell nodes as usual. Therefore, although the cells are arbitrarily cut by the boundary ∂Ω(t), the convergence rate of this finite element formulation is not affected by the position of the cut, leading to optimal discretizations. This result can be directly derived by the convergence proof shown in [33].
For the transport and microscopic problems we use the following spaces wherec n h andṽ n h are the extensions from Ω(t n ) to Ω(t n+1 ). In fact, the two components c h (t n ) and v h (t n ) are defined only in Ω(t n ). Therefore an extension in the region Ω(t n+1 ) \ Ω(t n ) needs to be defined. We have chosen a continuous extension using the prescribed values c 0 (x) and v 0 (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 t n and at time t n+1 , see Figure 4 and case (ii) the cell is cut at time t n and uncut at time t n+1 , see Figure 5. In case (ii) the cut goes to the neighbour cell at time t n+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 c h and v h outside the region of validity of the finite element interpolation. We have observed that this "extrapolation" gives good results only for cutcells that are not too small. In case of small cut-cells, the extrapolation leads to wrong values that introduce instabilities, visible as large peaks in the solution, that destroy the convergence of the method. 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 active mesh. We are currently studying the reason for this behavior of small cut-cells and let the rigorous discussion of this issue for a subsequent work since it goes beyond the scope of this paper.
The system of equations is solved with the following initial conditions The transport equation is defined on Y . Since this is a hyperbolic equation, a suitable discretization is needed. We choose the stream-line diffusion approach for its easy implementation and good performance. We have used an artificial diffusion in the stream-line direction scaled with a parameter δ > 0 whose value can be found in Table 2. On general meshes, the stream-line diffusion stabilized formulation of the transport equation converges in the L 2 norm with the rate h 3/2 with bilinear finite elements. Nevertheless, several authors show an optimal convergence rate of h 2 on regular meshes, see for example [62]. The dynamics of the cancer cells and ECM are solely defined by the velocity at the boundary and the initial distributions c(x, 0) and v(x, 0).

The discrete transport equation is
where φ 0 is the initial level set function, k the time step and V n h is the discrete velocity defined as where I x and I τ are two quadrature formulas for the approximation of the integral in space and time (see expression (13)) 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 point wise defined velocities to solve the transport problem. The weak formulation of the microscopic problem reads: where we have used the notation m l h,n to indicate the discrete microscopic solution for the macroscopic step n (note that the right hand side depends on the macroscopic solution at time t n ), with l being the time step of the time variable τ , i.e. m l h,n = m h,n (t 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 approximated by a quadrature rule where I B (·) is a quadrature rule that approximates the integral of the argument over B

Approximation of the interface and cut-cells
As introduced previously we discretize the problems in time by a time step method. Therefore, in the semidiscrete formulation we have terms that are defined at time t = t n+1 and terms defined at time t = t n . Since the domain is time dependent, the integrals of these terms are defined on different domains. Therefore, as explained above, 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  7 the cell at the bottom left and the case in which the boundary cuts one cell at time t n and it goes over to the neighbour cells at time t n+1 leaving the previous cell uncut at time t n+1 , see in Figure 7 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 where the linearized zero level L 0,h (t) is defined by the polygonal line that connects all intersections of the zero level L 0 (t) with the mesh cells boundaries ∂K as shown in Figure 6a. In Figure 6b a part of the actual mesh is shown. Furthermore, a part of the actual level zero isoline and the cut-cells are shown. Note that the irregular cells are not finite element cells, but are shown here only for visualization purpose. By the linear approximation of the boundary we introduce an integration error in the solution of the weak formulation. It can be shown that this error converges to zero with the order of the discretization error. Therefore, the integration error can be interpreted as a perturbation of the underlying Galerkin method that does not change the convergence rate of the method. Using the linearized boundary ∂Ω h (t) we can apply the quadrature rule described in [12] 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 t n and at time t n+1 , we apply the previous quadrature rule recursively.

Approximation of the nonlocal term
The nonlocal term (18) is approximated by a quadrature rule. We use a second level-set function to define the distance from the macroscopic point x that is used to define the domain of integration B. Also in this case we  introduce a piece-wise linear approximation of this second level-set function and we apply recursively the quadrature rule on the neighboring cells.

Extension of the velocity
The velocity is determined using formula (13). It is computed at the macroscopic point x on the linearized boundary ∂Ω h (t) and then extended to the rest of the domain. In this work, we consider only one point x per cell. This is taken at the midpoint of the segment of the interface L 0,h that intersects the cell, as shown in Figure 2.
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 cells that lie in the cancer region, i.e. K ∩ Ω h (t) = 0, at a distance larger than a prescribed radius of influence ρ (see Table 2) we set velocity zero. This procedure avoids 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 [1]. Fast marching algorithms are often used in the context of level set approaches. Important issues in this context are the initialization of the level set function [19], its re-initialization as a signed distance function [31,43,44] and the possibility to update the extension of the velocity using an accurate characteristics reconstruction in case the velocity changes rapidly in time [20].

Solution process
We sketch the overall solution process underlying the coupling between the different parts of the model.
The macroscopic system is solved with an implicit Euler scheme. At each time step a nonlinear system of the type (14a-14b) 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 a direct solver.
The system (15) is a linear system and is computationally much cheaper than the macroscopic problem. A direct solver is used in every time step.
Finally, the microscopic problem (17) is a linear one dimensional parabolic problem solved with an implicit Euler method and a direct solver in each time step.

Numerical results
In this section we show some numerical results obtained with the numerical method explained above. We have used the following initial conditions for cancer cells and ECM where R is the initial radius of the cancer region and the point (4, 4) is the center of the computational domain, see Table 2 for the numerical parameters. We first consider the configuration described in Table 1. In particular, we consider a constant value for the proliferation coefficient where µ * 1 is a constant shown in Table 1. It can be observed in Figure 8 that the cancer cells spread in the surrounding environment with a preferential path along the regions where the  Figure 9, the transport of cancer cells following the gradient of ECM is limited to the vicinity of the boundary ∂Ω(t). The ECM shows some boundary layers at the interface with regions where the value of ECM outside of Ω(t) is larger. The layers are shown in red in the electronic version of the manuscript.
In the second numerical test, we consider the ECM−dependent case for the proliferation parameter µ 1 that was briefly introduced in Section 2.1. In this context, the proliferation parameter µ 1 (v) is given here by which explores the proliferation conditions offered by the ECM, progressing smoothly from the worst conditions due to lack of ECM to the optimal conditions offered by abundant ECM density.  Figure 10: Comparison of cancer distributions: (a) for a constant proliferation parameter µ 1 = µ * 1 ; and (b) for ECM−dependent proliferation parameter µ 1 (v) given in (23). The black line is the contour of the level set c ≤ 0.5.
In Figure 10 we compare on the left side the distribution of cancer cells at time t = 10 computed with constant µ 1 = µ * 1 and on the right side the distribution computed with the function µ 1 as in (23). The black line shown on both pictures is the isoline for the level c(10, x) = 0.5. While a similar morphology of the tumour boundary ∂Ω(10) is observed, the effect of a decreasing µ 1 with a decreasing v results in the spatial spread of cancer cells of high distribution levels (above 0.5) being much more reduced in Figure  10(b) than in Figure 10(a), which corresponds to the case of constant µ 1 .

Conclusions
We have presented a new formulation of a two-scale model to simulate cancer invasion. This included a new derivation of the tumour boundary movement law by accounting on the MDE microdynamics contributions within a transport equation whose solution provides the level-set that indicates the new macro-scale cancer boundary.
At the core of the proposed numerical approach for the proposed model stands the moving boundary method based on a combination of the levelset approach and cut-cells. The later are based on a continuous Lagrangian finite element formulation that on one side has high flexibility and accuracy properties and on the other side has revealed some instabilities issues that have been solved with a modification of the formulation. In particular, we have suppressed the cut-cells that were below a certain volume threshold. This numerical aspect is important and we are investigating the reason for such behaviour. Nevertheless, the neglection of the contribution of small cells can be interpreted as a quadrature error. By keeping the threshold for the suppression small, we keep this quadrature error small.
For the computational implementation, we have used an unfitted regular mesh with uniform cell diameters to avoid the problem of remeshing in case of large deformations.
A further important aspect for the numerical solution of the problem, that might be studied in a future work, is the singularly perturbed character of the macroscopic problem. Even if, in the considered configurations, we have not observed instability problems, it is of importance on its own to study possible stabilization techniques in combination with cut-cells.
We have shown that the presented framework is highly flexible to study possible variations of the model. 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. In conclusion, we underline the potential of the presented method, that allows to go to three dimensional problems without changing the numerical formulation, enabling this way a major development of the multiscale modelling framework introduced in [56]. The major extension needed for this development is the formulation of cut-cells in three dimensions.