A Hybrid Multiscale Model for Cancer Invasion of the Extracellular Matrix
Nikolaos Sfakianakis, Anotida Madzvamuse, Mark A.J. Chaplain
AA Hybrid Multiscale Model for Cancer Invasion ofthe Extracellular Matrix
Nikolaos Sfakianakis ∗ Anotida Madzvamuse † Mark A.J. Chaplain ‡ Tuesday 29 th May, 2018
Abstract
The ability to locally degrade the extracellular matrix (ECM) andinteract with the tumour microenvironment is a key process distin-guishing cancer from normal cells, and is a critical step in the metastaticspread of the tumour. The invasion of the surrounding tissue involvesthe coordinated action between cancer cells, the ECM, the matrix de-grading enzymes , and the epithelial-to-mesenchymal transition (EMT).This is a regulatory process through which epithelial cells (ECs) ac-quire mesenchymal characteristics and transform to mesenchymal-likecells (MCs). In this paper, we present a new mathematical modelwhich describes the transition from a collective invasion strategy forthe ECs to an individual invasion strategy for the MCs. We achievethis by formulating a coupled hybrid system consisting of partial and stochastic differential equations that describe the evolution of the ECsand the MCs, respectively. This approach allows one to reproducein a very natural way fundamental qualitative features of the currentbiomedical understanding of cancer invasion that are not easily cap-tured by classical modelling approaches, for example, the invasion ofthe ECM by self-generated gradients and the appearance of EC inva-sion islands outside of the main body of the tumour.
Keywords: cancer invasion, multiscale modelling, hybrid continuum-discrete, ∗ Institute of Applied Mathematics, Faculty of Mathematics and Informatics, Heidel-berg University, Heidelberg 69120, Germany, [email protected] † Department of Mathematics, School of Mathematical and Physical Sciences, Univer-sity of Sussex, Brighton BN1 9QH, UK,
[email protected] ‡ Mathematical Institute, School of Mathematics and Statistics, University of St An-drews, St Andrews KY16 9SS, Scotland, [email protected] a r X i v : . [ q - b i o . CB ] M a y oupled partial and stochastic partial differential equations AMS subject classification:
Identified as one of the hallmarks of cancer , [18, 19], cancer invasion (andthe subsequent metastasis) is a complex process involving interactions be-tween cancer cells and the extracellular matrix (the tumour microenviron-ment ) facilitated by matrix degrading enzymes. By its nature, the invasioninvolves the development of and changes to cell-cell and cell-matrix adhe-sion processes. Broadly speaking, during the progression to full malignancy,cancer cells reduce their cell-cell adhesions and gain cell-matrix adhesions(controlled by adhesion molecules such as cadherin ). Coupled with cell mi-gration and proliferation, this enables the local spread of cancer cells intothe surrounding tissue. Any encounter with blood or lymphatic vessels inthe microenvironment initiates the spread of the cancer to secondary loca-tions in the host i.e. metastasis, which accounts for over 90% of deaths dueto cancer, [34, 53].Having been studied in some detail for the past 10-15 years, it has becomeclear that cancer invasion has a certain degree of diversity in its migratorymechanisms and a degree of plasticity in cellular behaviour and properties.The diversity of cancer invasion mechanisms is illustrated schematically inFig. 1 and the plasticity in 2. Both invasions are extensively discussed in[14]. Still, cancer invasion can be broadly classified into two main groups,each differing in the behaviour of how the cells migrate —individual versuscollective migration— and how these are controlled by different intra-cellularmolecular programmes. Accordingly, cancer invasion can be characterised as epithelial or collective invasion whereby clusters or sheets of connected cellsmove en masse, or as mesenchymal or individual invasion whereby singlecancer cells or small numbers of cancer cells actively invade the microenvi-ronment.The plasticity of cancer invasion mechanisms is illustrated schematically inFig. 2, see also [14]. Cancer cells may transition back and forth betweenthe two different invasion mechanisms during the invasion process as theypenetrate into the surrounding tissue. The transition process between thecollective and individual invasion is largely controlled by varying the ex-pression levels of molecules such as integrins, proteases, and cadherins and2igure 1: “Diversity of cancer invasion”. Classification of the various mi-gration and invasion strategies and corresponding types of tumour. As thecomplexity of the tumour increases, so do the expressions of cell-matrix andcell-adhesion molecules ( integrins and cadherins ) and the characterizationof the invasion as individual or collective. Figure adopted from [14]. (PER-MISSION REQUESTED)varying cell-cell communication via gap junctions. This process is known asthe epithelial-mesenchymal transition (EMT), whereas the opposite processis known as mesenchymal-epithelial transition (MET).An alternative invasion mechanism also exists — amoeboid invasion — wherebyindividual cells exhibit morphological plasticity and develop the ability tosqueeze through gaps in the extracellular matrix (ECM), rather than mod-ify/degrade the ECM via matrix degrading enzymes, e.g. urokinase-typeplasminogen activator (uPA), matrix metalloproteases (MMPs), [30, 43].Cancer invasion has also been the focus of mathematical modelling overthe past twenty years or so, beginning with the work of [15]. Since then,many different models and approaches have been formulated, some takingan individual-based (or agent-based) approach e.g. [40, 41, 21, 52, 44],others adopting a continuum approach using systems of partial differentialequations e.g. [10, 39, 8, 9, 1, 12, 11, 37, 26, 45, 13, 38], while others haveadopted a hybrid continuum-discrete approach, e.g. [4, 2].3igure 2: “Plasticity of cancer invasion”. The character of cancer cellmigration changes from collective to individual, following the loss of thecadherin or β multiscale framework approachto the problem. In particular, we describe (cancer) epithelial cells (ECs)by a density distribution and their spatiotemporal evolution by a macro-scopic deterministic model, whereas (cancer) mesenchymal cells (MCs) aremodelled by an atomistic approach, and their spatiotemporal evolution byan individual stochastic model . In this hybrid approach we also include thecoupling dynamics between the two cell types and their corresponding de-scriptions.The macroscopic sub-model for the evolution of the ECs includes primaryprocesses such as diffusion and proliferation of the ECs and the de-differentiation and differentiation of the MCs to ECs. It also includes some basic processesof the MMPs and ECM. The MMPs are assumed to be produced by thecancer cells, to diffuse in the extracellular environment, and to degrade theECM. The ECM, in its turn, is assumed to be non-uniform and is not remod-elled. On the other hand, the MCs are described as a finite set of particles that represent isolated cells or small cell aggregates. This representation is4otivated by the small number in which the MCs appear in the tumour. Forthe spatiotemporal evolution of the particles we assume that they obey an individual stochastic model that includes haptotaxis and random motion foreach particle. The total number of particles in the system varies accordingto the EMT and MET processes between the ECs and MCs.The triggering mechanisms of the EMT and MET are not specifically ad-dressed in this work. Instead, we follow a simplified approach and assumethat the EMT occurs randomly over the ECs with a given constant prob-ability over space and time. In a similar way, we assume that every MCparticle undergoes MET randomly with a given probability rate. Similarly,the description and modelling of specific mechanisms of the ECM dynamicsalso fall outside the scope of this paper. Nevertheless, we assume that theECM is represented by the density of the collagen macromolecules , and istherefore modelled as a spatially non-uniform, immovable component of thesystem. Finally, we assume that the ECM is degraded by the combined func-tion of the cancer cells MMPs complex and, for simplicity, that no matrixreconstruction takes place. In this section we present the main components of the model, the two types(or phases ) of the cancer cells, their properties and interactions. The modelwe propose is a hybrid amalgam of the two phases of the cancer cells thatare described by a continuum density for the DCs and a collection of discreteparticles for the ECs. For the sake of clarity of presentation, the developmentof the model and the corresponding techniques are constrained to the two(spatial) dimensional case only.Therefore, our paper is structured as follows. In Section 2.1 we describe the continuum density submodel of the problem. This is a macroscopic deter-ministic model that addresses the spatiotemporal evolution of the densitiesof the ECs, ECM, and the MMPs. The participation of the DCs in thissubmodel is merely implicit, i.e. the macroscopic model does not dictatetheir time evolution. We then introduce the discrete particle submodel ofthe problem in Section 2.2. This submodel refers to the MCs as particlesand is responsible for the time evolution. In Section 2.3 we describe thetransitions between the two phases of the cancer cells. We address the5ay the MC particles are substantiated from their density formulation viaa density-to-particles process, and how they transition back to density viaan opposite particles-to-density procedure. In Section 2.4 we present thecombined spatiotemporal evolution of the two phases under the prism ofthe EMT and MET processes. We moreover address separately the EMTand MET processes as well as their influence on the phases of the two typescancer cells.
From a macroscopic deterministic scale, we follow the seminal works of[29, 15, 4, 7] and describe the ECs, MMPs, and ECM by their densities. Asthe MCs are primarily described as particles, the model for their evolution isderived in Section 2.2. Still, the MCs appear also in the density formulationas they directly affect the ECs, MMPs, and the ECM. Since the focus ofthis paper is on the combination of the two types/phases of the cancer cellsrather than on the biological applications of the model, we consider onlythe very basic biological dynamics. More detailed and cancer-type specificmodels will be considered in follow up works. Indeed, we mainly assumethat the ECs are transformed to MCs and vice-versa via the MET andEMT processes, and that they proliferate by following a logistic volume-filling constraint as they compete for free space and resources with eachother and with the MCs and the ECM. Furthermore, we assume that theECs diffuse in the environment; although this process is expected to be veryslow.To proceed we denote by Ω ⊂ R the Lipschitz domain of study, and by c α ( x , t ), c β ( x , t ), m ( x , t ), and v ( x , t ), x ∈ Ω and t ≥ α and β the two types (or phases) of cancer cells,the ECs and MCs, respectively.It follows therefore that the equation that controls the evolution of the ECs6eads ∂∂t c α ( x , t ) = D α ∆ c α ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) diffusion − µ EMT α ( x , t ) c α ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) EMT + µ MET β ( x , t ) c β ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) MET + ρ αc c α ( x , t ) (cid:16) − c α ( x , t ) − c β ( x , t ) − v ( x , t ) (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) proliferation , (1a)where µ EMT α ( x , t ) = µ α X A ( t ) ( x ), µ MET β ( x , t ) = µ β X B ( t ) ( x ), with A ( t ) , B ( t ) ⊂ Ω, and D α , µ α , µ β , ρ αc ≥ c β after havingundergone a specific particle-to-density transformation; this is discussed inSection 2.3.The triggering mechanisms of EMT and MET are not the focus of this work.We instead assume a simplified approach where EMT occurs in a randomlychosen domain, denoted by A ( t ) ⊂ Ω in (1a), at a constant rate µ α . In asimilar way, we assume that the MET occurs randomly at every particle;this gives rise in a natural way to the domain B and the MET rate µ B , seealso Sections 2.2 and 2.3 for further details.Both types of cancer cells, ECs and MCs, produce MMPs, which in turndiffuse in the environment (molecular diffusion) and decay with a constantrate satisfy: ∂∂t m ( x , t ) = D m ∆ m ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) diffusion + ρ αm c α ( x , t ) + ρ βm c β ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) production − λ m m ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) decay , (1b)with D m , ρ αm , ρ βm , λ m ≥ ∂∂t v ( x , t ) = − (cid:16) λ αv c α ( x , t ) + λ βv c β ( x , t ) (cid:17) m ( x , t ) v ( x , t ) (cid:124) (cid:123)(cid:122) (cid:125) degradation , (1c)7ith λ αv , λ βv ≥ (advection-)reaction-diffusion (A-)RD system (1a)–(1c) can also bewritten in a more convenient matrix-vector compact form for the numeri-cal treatment formulation, see also Appendix A. In particular, using thenotation w ( x , t ) = (cid:0) c α ( x , t ) , m ( x , t ) , v ( x , t ) (cid:1) T , (1a)–(1c) read w t ( x , t ) = D ( w ( x , t )) + R ( w ( x , t )) , (2)where D ( w ) = D α ∆ c α D m ∆ m and R ( w ) = − µ EMT α c α + µ MET β c β + ρ αc c α (cid:0) − c α − c β − v (cid:1) ρ αm c α + ρ βm c β − λ m m − (cid:16) λ αv c α + λ βv c β (cid:17) mv denote the diffusion and reaction operators, respectively. In the more gen-eral case where chemotaxis or haptotaxis are considered, the correspondingformulation should also include an advection operator.Clearly, cancer invasion models of the form (2) are mere simplifications ofthe biological reality; they are also quite simple in mathematical structure.Nevertheless, their analytical and numerical investigations are challenging,see for example [1, 16, 26, 32, 54]. One of the reasons for this, is theirmixed nature, i.e. the ECs and MMPs obey partial differential equations (PDEs) with respect to time and space, whereas the ECM obeys an ordinarydifferential equation (ODE) with respect to time for every point in space. We are motivated in this description by methods and techniques that havebeen used previously in other scientific fields. One such example is theclassical particle-in-cell (PIC) method which was first proposed in [20] andused among others in plasma physics. A second example is be the smoothed-particle hydrodynamics (SPH) method used in astrophysics and ballistics(see [17] for example). The stochastic nature of the ODEs that the particlesobey is motivated by the seminal work of [47]. For the combination of thetwo cancer cell formulations we are inspired by [6, 24, 31, 48]. In the general case, (1a)–(1c) could include advection as well
8n view of the above, we describe the MCs as a system of N particles that areindexed by p ∈ P = { , . . . , N } , and account for their positions x p ( t ) ∈ R and masses m p ( t ) ≥
0. We allow for their number to vary in time and so weset N = N ( t ) ∈ N .The mass distribution of such system of particles, { ( x p , m p ) , p ∈ P } , isgiven by ˜˜ c ( x , t ) = (cid:88) p ∈ P m p ( t ) δ ( x − x p ( t )) (3)where δ ( · − x p ( t )) represents the Dirac distribution centred at x p ∈ R .Clearly (3) is not a function so we consider a kernel ζ and re-define the mass distribution of the particles { ( x p , m p ) , p ∈ P } as˜ c ( x , t ) = (cid:90) Ω ˜˜ c ( x (cid:48) , t ) ζ ( x − x (cid:48) ) d x (cid:48) (3) = (cid:88) p ∈ P m p ( t ) ζ ( x − x p ( t )) . (4)The function ζ does not have to be smooth and to simplify the rest of thiswork we choose it to be the characteristic function of the rectangle K thatis centred at the origin ∈ R ζ ( x ) = X K ( x ) , x ∈ R . (5)The choice of K (shape, size, and location) is justified in Sections 2.2.1 and3. We understand the particles as isolated cancer cells or cancer-cell aggregatesof similar size and masses. To maintain similar masses, we split and mergethe particles according to their mass and position. In particular, when theparticles represent an isolated cancer cell, we set m ref to be the referencemass of one cell and K its (two-dimensional) size , and proceed as follows: Splitting.
A particle ( x p , m p ) with mass m p > m ref is split into two par-ticles ( x p , m p ), ( x p , m p ) of the same position x p = x p = x p and mass m p = m p = m p . From that moment onwards, these two particles areconsidered different from each other. We consider typically the physiological parameters of the HeLa cells as reference erging. A small particle ( x p , m p ) with mass m p < m ref is merged withanother small particle ( x q , m q ) if they are close to each other i.e. (cid:107) x p − x q (cid:107) < diam( K ) , where (cid:107) · (cid:107) describes the two-dimensional Euclidean norm. The result-ing particle is set to have the cumulative mass of the two particles andto be located at their (combined) centre of mass (cid:18) m p x p + m q x q m p + m q , m p + m q (cid:19) . (6)Given that the distance between the particles is sufficiently small, iterationsof the merging and splitting processes lead to particles with masses m p ∈ [ m ref , m ref ].Besides the merging and splitting procedures, we do not consider other pro-cesses that alter the masses of the particles. Moreover, we do not includeany further interactions between the particles in this work as we try tobe consistent with the dynamics that are usually assumed by macroscopicdeterministic models, such as (1a)–(1c). We assume that the particles perform a biased random motion that is com-prised of two independent processes: a directed-motion part that representsthe haptotactic response of the cells to gradients of the ECM-bound adhesionsites, and a random/stochastic-motion part that describes the undirected ki-nesis of the cells as they sense the surrounding environment. We reproducethis way, at the particle level, the diffusion and -taxis dynamics prescribedby the macroscopic deterministic cancer invasion models, see e.g. [4].We understand this complex phenomenon as a geometric Brownian motion .In the general case, the corresponding stochastic differential equation (SDE)that it follows, would have the form d X pt = µ ( X pt , t ) dt + σ ( X pt , t ) d W pt , for p ∈ P, (7)where X pt represents the position of the particles in physical space (here R ), and W pt a Wiener process. Here, µ and σ are the drift and diffusion coefficients that encode the assumptions made on the directed and random10arts of the motion of the particles. Clearly, if more complex dynamics andinteractions between the particles and/or the environment are assumed, theSDE (7) should be adjusted accordingly.For the needs of this paper, we discretise (7) by an Itˆo-type explicit Euler-Maruyama particle motion scheme : x pt + τ = x pt + A ( x pt ) τ + B ( x pt ) · Z p √ τ , for p ∈ P , (8)cf. [24, 25] and Appendix B. Here, τ > A denotes the advection operator A = v + ∇ · D , (9)that encodes the advection velocity v adjusted by the drift term D . Thisparticular choice of A is made so that (8) converges, in the many-particlelimit N → ∞ , to the desired reaction-less version of the (A-)RD system(1a)–(1c). The more intuitive scheme with A = v would not converge inthe case of a non-constant diffusion tensor D ; we refer to [5, 24, 42, 47, 48]for the proofs of these claims and further discussions.In comparison to the usual macroscopic cancer invasion models, the advec-tion velocity v in (9) corresponds to the advection/-taxis term. The squarematrix B is related to the diffusion tensor D by B · B T = 2 D . (10)The typical Laplace operator L u = d ∆ u , would correspond here to a diago-nal D with inputs d . Moreover, in (8), Z p is a vector of normally distributedvalues of zero mean and unit variance. Modelling reactions
Although the MCs participate in several reaction processes (such as theEMT, MET, the proliferation of the ECs, the production of MMPs, and thedegradation of the ECM), the particle motion scheme (8) does not includeany reaction terms. We account for them in the following way:Some of the MC particles undergo MET to ECs, and subsequently are trans-formed to density via the particle-to-density operator that will be introducedin Section 2.3. These MCs are removed from the system of the MC parti-cles. The new EC density is added to the existing one and participates11ormally in the system (1a)–(1c). Conversely, a part of the EC density un-dergoes EMT towards MC density which is then transformed into particlesvia a density-to-particle operator defined in Section 2.3. These newly formedMCs are then added to the system of the existing MC particles.Moreover, at every (time instance and) timestep of the method, the fulldistribution of MC particles is transformed temporarily to density (withoutundergoing MET to ECs), via the particle-to-density operator. They par-ticipate then in the proliferation of the ECs, the production of the MMPs,and the degradation on the ECM, cf. (1a)–(1c). We give more details onthe combination of the EC and MC phases in Section 2.4.
In this section we describe the forward particle-to-density and the backward density-to-particle phase transition operators.We assume at first, that the domain Ω is regular enough to be uniformlypartitioned in equal rectangles/ partition cells { M i , i ∈ I } Ω = (cid:91) i ∈ I M i , (11)where every M i is an affine translation of the generator cell K . Note that K is the same as the support of the characteristic function in (5). Clearly | M i | = | K | = K > Remark . The partition cells M i , i ∈ I should not be confused with thediscretization cells of the numerical method used to solve (1a)–(1c). Thelatter constitute an instance of a sequence of computational grids of zero-converging step size, whereas the former have a step size that representsphysical properties of biological cells and remain fixed for all computationalgrid resolutions.Using the partitioning of Ω to { M i , i ∈ I } , we represent every measurable c : Ω → R by its simple-function decomposition (cid:88) i ∈ I c i ( t ) X M i ( x ) , (12) We mostly consider the diameter of the HeLa cell. X M i is the characteristic function of the set M i ⊂ Ω, and c i ( t ) themean value of c ( · , t ) over M i c i ( t ) = 1 K (cid:90) M i c ( x , t ) d x . (13)Clearly, this representation conserves the mass of c ( · , t ) (cid:88) i ∈ I Kc i ( · , t ) = (cid:90) Ω c ( x , t ) d x (14)On the other hand a particle, indexed here by p ∈ P , can be representedeither by its position and mass ( particle formulation )( x p ( t ) , m p ( t )) , (15)or by the characteristic function with density value ( density formulation ) m p ( t ) K X K p ( x ) , (16)where K p is the affine translation of the generator cell K centred at x p .Clearly (16) implies that the mass m p of the particle is uniformly distributedover K p .Although the K p , p ∈ P and the M i , i ∈ I in (12) are equivalent up toaffine translations (to the K ), they do not in general coincide. The M i , i ∈ I form a fixed partition of the domain, cf. (11), whereas the K p , p ∈ P “follow” the position of the particles (16).Based on the “dual” description (15) and (16) of the particles, we set forththe transition operators between particles and densities. Let { ( x p ( t ) , m p ( t )) , p ∈ P } be a collection of particles. Using (4), we definethe forward particle-to-density operator F , { ( x p ( t ) , m p ( t )) , p ∈ P } F −→ c ( x , t ) . (17)To define the target function c ( x , t ), we go through all the particles, indexedhere by p ∈ P , and consider their corresponding density formulation (16).13 x y x y Figure 3: Two-dimensional graphic representation of the forward particle-to-density operator F . (left:) We consider a support K p ( p ∈ P ) aroundthe location x p of every particle. The mass of every particle m p , (shownas points) is uniformly distributed over the respective support K p . Thegrid represents the partitioning of the domain. (right:) A view from abovereveals that the supports K p can overlap with several cells of the partition.The corresponding masses are assigned to the partition cells using (19).The support K p of the particles, overlaps with (possibly) several of thepartition cells M i , i ∈ I . In each of these partition cells, we assign thecorresponding portion of the particle mass m p (cid:12)(cid:12) M i = m p K (cid:12)(cid:12) K p ∩ M i (cid:12)(cid:12) . (18)Next, we account for the contribution of all particles p ∈ P at the cell M i by c i ( t ) = (cid:88) p ∈ P K m p (cid:12)(cid:12) M i (18) = (cid:88) p ∈ P m p ( t ) K (cid:12)(cid:12) K p ∩ M i (cid:12)(cid:12) , for i ∈ I . (19)In view now of (12) and (19), we deduce the density function c ( x , t ) (as asimple function) over the full domain Ω as c ( x , t ) = (cid:88) i ∈ I c i ( t ) X M i ( x ) , x ∈ Ω . (20)Refer to Fig. 3 for a graphical representation of the forward particle-to-density operator F in two dimensions. Since the sets K p , p ∈ P and M i , i ∈ I are two-dimensional quadrilaterals of the samedimensions, every K p overlaps at most four M i s. y x Figure 4: Graphical representation of the backward density-to-particle op-erator B . We compute the mass m i of the density function c ( x , t ) (surface),over every partition cell M i , i ∈ I (quadrilateral grid on the xy plane), using(22). We then define the particle as ( x i , m i ) where the location x i is givenby (23). Conversely, we define the backward density-to-particle operator B for a givendensity function c ( x , t ) by { ( x p ( t ) , m p ( t )) , p ∈ P } B ←− c ( x , t ) , (21)in the following way: in every partition cell M i , i ∈ I , we assign one particlewith mass m i ( t ) = (cid:90) M i c ( x , t ) d x . (22)and position x i ( t ) = the (bary)centre of M i . (23)For practical considerations, we set in the numerical simulations a minimumthreshold value on the densities, below which no transition to particles takesplace. This threshold value is quite small and is used to avoid very largenumber of particles of negligible-mass. Refer to Fig. 4 for a graphicalrepresentation of the backward density-to-particles operator. We denote again the two types of cancer cells, EC and MC by the supescripts α and β respectively, and consider for t ≥ w ( x , t ) = ( c α ( x , t ) , m ( x , t ) , v ( x , t )) . At the same physical time t , we write the MC particles as P β ( t ) = (cid:110)(cid:16) x βp ( t ) , m βp (cid:17) , p ∈ P ( t ) (cid:111) , (24)and, accordingly, the overall system is given by the tuple (cid:16) w ( x , t ) , P β ( t ) (cid:17) , x ∈ Ω , t ≥ . (25)In the evolution of the overall system, we consider the EMT and METprocesses separately from the rest of the dynamics of the system (1a)–(1c) . The EMT triggering mechanism is not one of the main foci of this work.Instead, we assume a simplified approach where a randomly chosen part ofthe ECs (in density formulation) c α EMT undergoes EMT to give rise to MCs(still in density formulation) c α EMT EMT −−−→ c β EMT . The newly created MC density c β EMT is transformed to MC particles via thedensity-to-particle operator B given in (21) c β EMT B −→ (cid:110) ( x βp , m βp ) , p ∈ P EMT (cid:111) , (26)where x βp , m βp follow from (22), (23) and P EMT is the corresponding set ofindexes. Subsequently, the family of MC particles is updated as the amalgamof the existing and the newly created particles (cid:110) ( x βp , m βp ) , p ∈ P (cid:111)(cid:124) (cid:123)(cid:122) (cid:125) existing MC particles (cid:93) (cid:110) ( x βp , m βp ) , p ∈ P EMT (cid:111)(cid:124) (cid:123)(cid:122) (cid:125) newly created MC particles = (cid:110) ( x βp , m βp ) , p ∈ P new (cid:111) , (27) To ease the presentation and since the EMT and MET are assumed to be instantaneous and tautochronous , we drop the dependence of the density variables and the particles on x and/or t . P new is a re-enumeration of the multiset P (cid:93) P EMT .Overall, combining the density and particle phases, the EMT operator reads R EMT (cid:16) c α , (cid:110)(cid:16) x βp , m βp (cid:17) , p ∈ P (cid:111)(cid:17) = (cid:16) c α − c α EMT , (cid:110) ( x βp , m βp ) , p ∈ P new (cid:111)(cid:17) . (28) As with the EMT, the triggering mechanism of the MET is not one of thefoci of this paper. We instead assume an approach where each of the MCparticles (cid:110)(cid:16) x βp , m βp (cid:17) , p ∈ P (cid:111) undergoes MET to ECs randomly (cid:110) ( x βp , m βp ) , p ∈ P (cid:111) MET −−−→ (cid:8) ( x αp , m αp ) , p ∈ P MET (cid:9)(cid:124) (cid:123)(cid:122) (cid:125) newly created EC particles . (29)The resulting EC particles are instantaneously transformed to density viathe particle-to-density operator F given in (17): (cid:8) ( x αp , m αp ) , p ∈ P MET (cid:9) F −→ c α MET . In operator form, the MET reads R MET (cid:16) c α , (cid:110)(cid:16) x βp , m βp (cid:17) , p ∈ P (cid:111)(cid:17) = (cid:16) c α + c α MET , (cid:110) ( x βp , m βp ) , p ∈ P new (cid:111)(cid:17) , (30)where P new is a re-enumeration of the set difference P \ P MET . The evolution of the overall system of ECs is controlled by (1a)–(1c) and(7). We study this combined system of PDEs and SDEs numerically andpostpone any analytical investigations to a follow up work. To this end, weconsider the model (2) and first set W n = (cid:110) w n ( i,j ) = (cid:16) c n ( i,j ) , m n ( i,j ) , v n ( i,j ) (cid:17) , ( i, j ) ∈ M x × M y (cid:111) , P β,n = (cid:110)(cid:16) x β,np , m βp (cid:17) , p ∈ P n (cid:111) , to denote numerical approximations of the density and particle variables w ( x , t ) and P β ( t ), respectively, at the instantaneous time t = t n . Here, M x M y denote the resolution of the grid along the xy − directions. We referto Appendix A for further information on the numerical method employedon W ; we focus here on the combination of the two phases by considering an operator splitting approach. In particular, for t ∈ [ t n , t n +1 ], t n +1 = t n + τ n ,we assume that:— During the time period [ t n , t n + τ n ], the system evolves, without theinfluence of the EMT or the MET, as (cid:16) W n , P β,n (cid:17) −→ (cid:16) W n +1 / , P β,n +1 / (cid:17) (31a)with W n +1 / = N [ t n ,t n + τ n ] (cid:16) W n , P β,n (cid:17) , (31b) P β,n +1 / = (cid:110)(cid:16) x β,n +1 / p , m β,n +1 / p (cid:17) , p ∈ P n +1 / (cid:111) , (31c)where N [ t,t + τ ] is the numerical solution operator responsible for thespatiotemporal evolution of the system (1a)–(1c) —without EMT andMET. Here, the x β,n +1 / p , p ∈ P n is given by the Itˆo-type particlemotion scheme (8), re-written here with respect to the local variables x β,n +1 / p = x β,np + A (cid:16) x β,np (cid:17) τ n + B (cid:16) x β,np (cid:17) · Z p (cid:114) τ n . (32)The number of particles, their indices and masses remain unchangedduring this step [ t n , t n + τ n ], i.e. P n +1 / = P n and m β,n +1 / p = m β,np , ∀ p ∈ P n . Altogether, the combined evolution operators of the two phases readsfor this time period as: M τ n (cid:16) W n , P β,n (cid:17) = (cid:16) W n +1 / , P β,n +1 / (cid:17) . (33)— At t = t n + τ n , the EMT and MET processes take place; they areassumed to be instantaneous and tautochrone . They are representedby the R EMT and R MET operators introduced in (28) and (30) respec-tively. For consistency, we scale them by the time step τ n and changetheir notation to R EMT τ n and R MET τ n , respectively.18n effect, the tuple (cid:0) W n +1 / , P β,n +1 / (cid:1) develops as (cid:16) ˜ W n +1 / , ˜ P β,n +1 / (cid:17) = R τ n (cid:16) W n +1 / , P β,n +1 / (cid:17) , (34)where R τ n denotes the parallel application of R EMT τ n and R MET τ n .— During [ t n + τ n , t n +1 ], the two phases evolve again without the influ-ence of the EMT and MET as (cid:16) ˜ W n +1 / , ˜ P β,n +1 / (cid:17) −→ (cid:16) W n +1 , P β,n +1 (cid:17) , where, in a similar way as in [ t n , t n + τ n ], W n +1 = N [ t n + τ n ,t n +1 ] (cid:16) W n +1 / , ˜ P β,n +1 / (cid:17) , (35) P β,n +1 = (cid:110)(cid:16) x β,n +1 p , m β,n +1 p (cid:17) , p ∈ P n +1 (cid:111) . (36)Again, N [ t n + τ n ,t n +1 ] represents the numerical method for the solutionof the system (1a)–(1c), P n +1 = P n +1 / and m β,n +1 p = m β,n +1 / p , ∀ p ∈ P n +1 . In the above, x β,n +1 p , p ∈ P n +1 / , is given by the Ito-type scheme (8) x β,n +1 p = x β,n +1 / p + A (cid:16) x β,n +1 / p (cid:17) τ n B (cid:16) x β,n +1 / p (cid:17) · Z p (cid:114) τ n . (37)We combine the evolution operators of the two phases as: M τ (cid:16) ˜ W n +1 / , ˜ P β,n +1 / (cid:17) = (cid:16) W n +1 , P β,n +1 (cid:17) . (38)Overall, using (33), (34), and (38), we can write the combined evolutionoperator for the time period [ t n , t n +1 ] as a splitting method of the form (cid:16) W n +1 , P β,n +1 (cid:17) = M τn R τ n M τn (cid:16) W n , P β,n (cid:17) . (39) Note that the R EMT τ n acts on the EC density and the R MET τ n on the MC particles description symbol values and unitsEC dens. diff. coef. D α d − EC dens. prol. coef. ρ α − MC part. diff. coef. | B | . d − MC part. hapt. coef. 30 cm mol − d − MC part. ref. mass m ref × − grMC part. ref. diam. | K | × − cmEMT prob. 5 × − EMT rate µ a × MET prob. 0MMP diff. coef. D m λ αv mol − d − ECM MC dens. degr. λ βv mol − d − We perform and present three numerical experiments to exhibit the dy-namics and combination of the two phases. As the focus of this paper ismore technical, we postpone the more biological relevant experiments anddiscussion for a follow-up work.The implementations of numerical schemes and algorithms, and the simu-lations of the experiments included in this paper have been conducted in[33].
Experiment . We set Ω = [ − , × [ − ,
2] andconsider the initial EC density c α ( x ,
0) = (cid:16) e − x + x ) − . (cid:17) + , (40a)with x = ( x , x ) ∈ Ω, where ( · ) + denotes the positive part function. TheECM is non-uniform and exhibits a gradient towards the upper-right partof the domain v ( x ,
0) = 0 .
045 (2 x + 3 x ) + 0 . . (40b)Initially, no MC particles nor MMPs are present.We close the system with no-flux boundary conditions for the EC and MMPdensity and reflective boundary conditions for the MCs particles. As theECM is modelled as an immovable part of the system, it does not translo-cate, hence no boundary conditions are needed. The parameters for this Each particle that escapes the domain is returned to its last position within thedomain. It then “decides” its new direction randomly. description symbol values and unitsEC dens. diff. coef. D α d − EC dens. prol. coef. ρ α − MC part. diff. coef. | B | × − cm d − MC part. hapt. coef. 1 × − cm mol − d − MC part. ref. mass m ref × − g r MC part. ref. diam. | K | × − c m EMT prob. 1EMT rate µ a D m λ αv
20 cm mol − d − ECM MC dens. degr. λ βv
200 cm mol − d − experiment are given in Table 1 and the result of the simulation is shown inFig. 5.The phenomena that is observed in this experiment can be described asfollows:The ECs undergo EMT to MC and new particles appear in the system.The particles “sense” the gradient of the ECM and respond haptotacticallyto it. Their motion incorporates also a random component; the resultingmigration is a biased-random motion. The simplified EMT that we assumein this experiment, takes place in every partition cell K i , i ∈ I with aprobability that is denoted as “EMT prob.” in Table 1. The set union of allthe partitions cells where EMT takes place defines the set A that appears in(1a). The rate µ a at which the EMT occurs in A is a given constant. Thereis no proliferation of the ECs assumed in this experiment, neither diffusionnor MET. Hence, the losses that the ECs suffer due to the EMT appear as“holes” in their density profile and are not replenished with time. Experiment . A typical phenomenon that macro-scopic cancer invasion models exhibit, is the appearance of a propagatingfront that invades the ECM faster than the rest of the tumour, see e.g[7, 8, 45]. This front is followed by an intermediate distribution, whereas thebulk of the tumour stays further behind. This phenomenon is due primarilyto the degradation of the extracellular chemical or the ECM landscape bythe cancer cells. Such phenomena have been observed previously both inmathematical models and in biological experiments, see for example [49, 4].In accordance to these experimental and modelling findings, we exhibit herethe ability of the particle description of our model to reproduce such phe-nomena. In particular, we show that as the cancer cells degrade the ECM,21 a) t = 0 .
01 (b) t = 0 .
06 (c) t = 0 . t = 0 .
40 (e) t = 0 .
85 (f) colorbars
Figure 5: Experiment 3.1 (EMT and flow). The experiment is visualisedover [ − , . (a): An initial circular EC tumour (shown by isolines) residesover an ECM that exhibits a gradient towards the north-east direction. (b):The DDC density undergoes EMT and gives rise to clusters of MC particles(black stars). (c): Due to the diffusion and the haptotaxis, the particlesescape the initial tumour and migrate along the gradient of the ECM. (d):We do not assume proliferation for the ECs, hence the losses of their densitiesthat suffer towards MCs are not replenished. This gives rise to “holes” inthe initial tumour. (e): The phenomenon continues as long as parts of theEC density transform to MC particles. (f): Common colorbars for the ECM(left) and the EC densities (right) in all sub-figures.22 a) t = 1 (b) t = 40 (c) t = 90(d) t = 130 (e) t = 180 (f) t = 230 Figure 6: Experiment 3.2 (Self-generated gradient). (a): Over an initiallyuniform ECM (background) resides a number of MC particles. As the matrixis initially uniform, the motion of the particles is mostly Brownian. (b):The particles degrade the matrix and create a gradient. This self-inducedgradient is sensed by the particles that are closer to the “interface”. In effect,their motion is mostly haptotaxis driven. (c): As the particles invade theECM, they continue to produce MMPs, and to degrade the matrix. Theyfollow the new gradient that they have induced. (d-f): The migration of theparticles in the front is persistent in direction and speed, while the particlesin the rear (where the ECM is depleted) perform mostly a Brownian motion.23hey induce a gradient on it, and subsequently they respond to this gradientby performing a directed and sustainable invasion.For this experiment we consider the domain Ω = [ − . , . × [0 , v ( x ,
0) = 0 . , x ∈ Ω . (41)On the upper part of the domain, an initial EC density is found c α ( x ,
0) = 10 − X S ( x ) , x ∈ Ω , (42)with S = (cid:8) x = ( x , x ) ∈ Ω (cid:12)(cid:12) x > .
01 sin(5 πx ) + 1 . (cid:9) . Before the be-ginning of the simulation, the EC density c α ( x ,
0) is completely transformedto MC particles. The MCs secrete MMPs that participate in the degradationof the ECM. No MET takes place in this experiment. The correspondingmodelling parameters are given in Table 3 and the simulation results in Fig.6.In view of (8), all the particles perform a biased random motion; since theECM is uniform, this motion is initially purely Brownian. As the ECM isdegraded, a gradient in the matrix is formed. The particles that reside closerto this “interface” sense the gradient and respond haptotactically to it. Thedirected part of their motion dominates and drives the particles to highermatrix densities. As the particles continue their invasion of the ECM theykeep on producing MMPs, degrading the ECM, and following the newlycreated gradient. Their motion is persistent in direction and speed.We can now address with our model, particular questions of experimentalinterest: what is the minimum number of cancer cells needed to induceand sustain a consistent in direction and speed invasion of the ECM? Howdoes the remodelling of the matrix affect the self-generated gradient motion?Such questions would among others serve as a bridge between experimentalobservations and mathematical models. Their study has to be the topic ofa follow-up work the relevant experimental data should also be analysed, aswas done, for example, in [55].
Experiment . This experiment is motivated by the organotypic invasion assays where cancer cells are plated over a collagen gelthat contains healthy tissue, and where their invasion is studied over time,see for example [35, 50] and Fig. 7.We employ the complete set of dynamics of the system and consider the do-main Ω = [ − , occupied by an ECM of initial density v ( x ,
0) constructed24igure 7: Timecourse (days 3, 9, and 14) study of the invasion of squamuscell carcinoma cells (black matter) on an organotypic assay with humanfibroblast cells (gray matter). The invasion occurs in the form of cancer cell“islands” formed in front of the main body of the tumour. We reproducethe same phenomenon in the invasion Experiment 3.3 and in Fig. 8. Theseimages are taken from [35] (PERMISSION REQUESTED).Table 3: Parameters, units, and sourses corresponding to Experiment 3.3and Fig. 8. description symbol values and units sourcesEC dens. diff. coef. D α . × − cm d − [8]EC dens. prol. coef. ρ α . − [36]MC part. diff. coef. | B | × − cm d − [46]MC part. hapt. coef. 3 cm mol − d − (our estimate)MC part. ref. mass m ref × − gr B10NUMB3R5 (HeLa cell)MC part. ref. diam. | K | × − cm B10NUMB3R5 (HeLa cell)EMT prob. 1 × − (our estimate)EMT rate µ a × − (our estimate)MET prob. 2 × − D m λ αv × − cm mol − d − [3]ECM MC dens. degr. λ βv × − cm mol − d − [3] & (our estimate) by 64 randomly chosen extremal values per direction that are interpolated ina piecewise linear way. Small perturbations of the form of additive Gaussiannoise are also included.An initial density of ECs is found in the upper part of the domain c α ( x ,
0) = 0 . X S ( x ) , x ∈ Ω , (43)with S ( x ) = (cid:8) x = ( x , x ) ∈ Ω (cid:12)(cid:12) x > .
05 sin(5 πx ) + 0 . x + 1 . (cid:9) . Ini-tially, neither MC particles nor MMPs exist in the system. The parametersfor this experiment can be found in Table 3 and the simulation results arepresented in Fig. 8. 25he ECs proliferate and diffuse, but most notably transform via EMT toMC particles. These MC particles do not proliferate but they are veryaggressive in their motility. As they escape the main body of the tumour,they undergo MET back to ECs. New EC concentrations appear, they growdue to proliferation, and give rise to tumour “islands”. These “islands”merge with each other as well as with the main body of the tumour. Themain characteristic and novelty of our hybrid model, is that it predicts, in anatural way, the appearance of these tumour “islands” outside of the mainbody of the tumour.The growth of the tumour with the combined dynamics of the ECs and MCspossesses several interesting properties. The tumour grows much faster thanit would, if it was comprised only of the ECs. This is so, since the newEC “islands” that arise after the MCs have escaped the main body of thetumour, undergo MET, exploit uninhabited locations, and grow “to all thedirections”. On the contrary, in the main body of the tumour, only the ECsfound in the periphery contribute to the growth of its support.Moreover, the independent and aggressive migration of the MCs providesthem with faster access to the circulatory network and the possibility totranslocate to secondary places within the organism. As they possess theability to give rise to EC “islands” at the new locations, new tumours mightappear, and metastasis will have occurred. Although it is not our aim inthe current paper to reproduce particular experimental scenarios, a directcomparison of the simulation results in Fig. 8 with the organotypic assayimages in Fig. 7 exhibits clearly that this phenomenon is reproduced in avery natural way by our modelling approach.Another sought-for property, in cancer invasion modelling, is that the MCsremain undetected while they invade the ECM. It is not until a new ECstumour has been established that it can grow to be of any detectable size.Again, this property is reproduced by our modelling approach in a verynatural way. We have proposed in this work a new modelling approach to study thecombined invasion of the ECM by two types of cancer cells, the ECs and theMCs. The proposed framework is a multiscale hybrid model that treats the26 a) t = 1 (b) t = 150 (c) t = 200(d) t = 270 (e) t = 330 (f) colorbars Figure 8: Experiment 3.3 (ECM invasion). Shown here the time evolutionof the ECM (landscape) and the EC (isolines) over the domain Ω = [ − , .(a): An initial uniform density of the ECs evolves according to the system(1a)–(1c), and mostly proliferates rather than diffuses. (b): The MC parti-cles that are produced by the EMT (not shown here) escape the main bodyof the tumour, invade the ECM more freely than the EC density, undergoMET and eventually give rise to new EC “islands”. (d)-(e): These “islands”grow mostly due to proliferation and eventually merge with the main bodyof the tumor. (f): The colorbars for the ECM (left) and the EC density(right) are common to all figures. 27Cs in a macroscopic and deterministic manner and the MCs in an atomisticand stochastic way.We assume that the MCs are much fewer than the ECs and that they em-anate from the ECs via a dynamic EMT-like cellular differentiation program.We also assume that the MCs give rise, via the opposite MET-like cellularprogram, to ECs; this is a key property in the metastasis of the tumour.For simplicity, we assume that both types of cancer cells perform a biasedrandom motion, and that the MCs are much more “aggressive” in theirmigration than the ECs.We encode this information through a hybrid approach: the spatiotemporalevolution of the ECs, the ECM, and the rest of the environmental compo-nents are dictated by a macroscopic deterministic model, (1a)–(1c). TheMCs on the other hand are considered as separate particles that evolve ac-cording to a system of SDEs, (7). The transition between the two types ofcancer cells is conducted by the density-to-particles and particles-to-density operators given in (17) and (21).This new modelling approach allows us to reproduce in a very natural way,several biological relevant phenomena encountered in the invasion of cancerthat are not easily addressed with the usual modelling approaches. Ourfocus though in this work lies with the description and the handling of themathematical model and the numerical method; we only present here basicbiological situations and postpone the more elaborate investigations for afollow-up work.With the atomistic component of our model, we are able to reproduce a sus-tainable invasion of the ECM by means of a self-induced haptotaxis gradient as shown in Experiment 3.2. Such behaviour is observed in biological situa-tions and becomes crucial to several biological processes like wound healing .The detailed study of such cases falls beyond the scope of the current paper;here we use this experiment as an indication that our model can reproducebiologically relevant situations.With the full model, we are able to reproduce the spread of the tumourand the invasion of the ECM in the form of invasion “islands”, Experiment3.3. These are well known to appear in many cases of cancer and are quitechallenging to reproduce by either macroscopic or atomistic cancer invasionmodels. With our approach these invasion “islands” appear in a very nat-ural way, and —most notably— they appear outside the main body of thetumour. 28hat is also very natural in our approach, is that the MC particles escapethe main body of the tumour and remain undetected while they invade theECM. It is only after they have established new “islands” in the vicinity ofthe original tumour or in another location within the organism that theycan be detected. This is another sought-after property in the field of cancerinvasion modelling.For the sake of presentation, we have only considered here some of the funda-mental properties of cancer growth that our model can reproduce. Still theysuffice to warrant extensions and investigations of more realistic biologicalsituations and experimental settings. To mention but a few: extension tothree-dimensional space, more realistic EMT and MET transitions, interac-tions between cancer cells of the same and different type including collisions,adhesions, short or long range interactions, and the collective behaviour ofcancer cells. Data Management
All the computational data output is included in the present manuscript.
Acknowledgements
NS was partly funded from the German Science Foundation (DFG) underthe grant SFB 873: “Maintenance and Differentiation of Stem Cells in De-velopment and Disease”. For the work AM is partly supported by fundingfrom the European Union Horizon 2020 research and innovation programmeunder the Marie Sk(cid:32)lodowska-Curie grant agreement No 642866, the Com-mission for Developing Countries, and was partially supported by a grantfrom the Simons Foundation. AM is a Royal Society Wolfson Research MeritAward Holder funded generously by the Wolfson Foundation. The authors(NS, AM, MC) would also like to thank the Isaac Newton Institute forMathematical Sciences for its hospitality during the programme [CouplingGeometric PDEs with Physics for Cell Morphology, Motility and PatternFormation] supported by EPSRC Grant Number EP/K032208/1.29 eferences [1] V. Andasari, A. Gerisch, G. Lolas, A.P. South, and M.A.J. Chaplain. Mathe-matical modelling of cancer cell invasion of tissue: biological insight from math-ematical analysis and computational simulation.
J. Math. Biol. , 63(1):141–71,2011.[2] A.R.A Anderson. A hybrid mathematical model of solid tumour invasion: theimportance of cell adhesion.
Math. Med. Biol. , 22:163–186, 2005.[3] A.R.A Anderson and M.A.J. Chaplain. A mathematical model for capillarynetwork formation in the absence of endothelial cell proliferation.
Appl. Math.Letters , 11:109–114, 1998.[4] A.R.A. Anderson, M.A.J. Chaplain, E.L. Newman, R.J. C. Steele, and A.M.Thompson. Mathematical modelling of tumour invasion and metastasis.
J.Theor. Medicine , 2:129–154, 2000.[5] L. Arnold.
Stochastic Differential Equations: Theory and applications . JohnWiley, New York, 1974.[6] X. Blanc, C. Le Bris, and P.-L. Lions. Atomistic to continuum limits forcomputational materials science.
Math. Mod. Num. Anal. , 41:391–426, 2007.[7] M.H. Byrne, M.A.J. Chaplain, G.J. Pettet, and D.L.S. McElwain. A mathe-matical model of trophoblast invasion.
J. Theor. Medicine , 1:275–286, 1999.[8] M.A.J. Chaplain and G. Lolas. Mathematical modelling of cancer cell invasionof tissue: the role of the urokinase plasminogen activation system.
Math.Modell. Methods. Appl. Sci. , 15:1685–1734, 2005.[9] M.A.J. Chaplain and G. Lolas. Mathematical modelling of cancer invasion oftissue: Dynamic heterogeneity.
Net. Hetero. Med. , 1:399–439, 2006.[10] M.A.J. Chaplain and M.E. Orme. Travelling waves arising in mathematicalmodels of tumour angiogenesis and tumour invasion.
FORMA , 10:147–170,1996.[11] N.E. Deakin and M.A.J. Chaplain. Mathematical modeling of cancer invasion:the role of membrane- bound matrix metalloproteinases.
Front. Oncol. , 3:70,2013.[12] P. Domschke, D. Trucu, A. Gerisch, and M.A.J. Chaplain. Mathematical mod-elling of cancer invasion: Implications of cell adhesion variability for tumourinfiltrative growth patterns.
J. Theor. Biol. , 361:41–60, 2014.[13] C. Engwer, C. Stinner, and C. Surulescu. On a structured multiscale modelfor acid-mediated tumor invasion: The effects of adhesion and proliferation.
Math. Models Methods Appl. Sci. , 27:1355, 2017.[14] P. Friedl and K. Wolf. Tumour-cell invasion and migration: diversity andescape mechanisms.
Nat. Rev. Cancer , 5:362–74, 2003.
15] R.A. Gatenby and E.T. Gawlinski. A reaction-diffusion model of cancer inva-sion.
Cancer Res. , 56:5745–53, 1996.[16] J. Giesselmann, N. Kolbe, M. Lukacova-Medvidova, and N. Sfakianakis. Ex-istence and uniqueness of global classical solutions to a two species can-cer invasion haptotaxis model. (accepted in) Discrete Cont. Dyn.-B , 2017,(arXiv:1704.08208).[17] R.A. Gingold and J.J. Monaghan. Smoothed particle hydrodynamics: theoryand application to non-spherical stars.
Mon. Not. R. Astron. Soc. , 181:375–389.[18] D. Hanahan and R.A. Weinberg. The Hallmarks of Cancer.
Cell , 100:57–70,2000.[19] D. Hanahan and R.A. Weinberg. The Hallmarks of Cancer: The next genera-tion.
Cell , 144:646–671, 2011.[20] F.H Harlow. A machine calculation method for hydrodynamic problems.
LAMS , 1965.[21] H. Hatzikirou, L. Bruscha, C. Schaller, M. Simon, and A. Deutsch. Predictionof traveling front behavior in a lattice-gas cellular automaton model for tumorinvasion.
Comput. Math. Appl. , 59:2326–2339, 2010.[22] S.M. Iacus.
Simulation and Inference for Stochastic Differential Equations .Springer, 2008.[23] C.A. Kennedy and M.H. Carpenter. Additive Runge-Kutta schemes forconvection-diffusion-reaction equations.
Appl. Numer. Math. , 1(44):139–181,2003.[24] P.K. Kitanidis. Particle-tracking equations for the solution of the advection-dispersion equation with variable coefficients.
Water Resour. Res. , 30:3225–3227, 1994.[25] P.E. Kloeden and E. Platen.
Numerical Solution of Stochastic DifferentialEquations . Springer Berlin Heidelberg, Berlin, 1992.[26] N. Kolbe, J. Katuchova, N. Sfakianakis, N. Hellmann, and M. Lukacova-Medvidova. A study on time discretization and adaptive mesh refinementmethods for the simulation of cancer invasion: The urokinase model.
Appl.Math. Comput. , 273:353–376, 2016.[27] A.N. Krylov. On the numerical solution of the equation by which in technicalquestions frequencies of small oscillations of material systems are determined.
Otdel. mat. i estest. nauk. , VII(4):491–539, 1931.[28] O. Lakkis, A. Madzvamuse, and Ch. Venkataraman. Implicit–explicit timestep-ping with finite element approximation of reaction–diffusion systems on evolv-ing domains.
SIAM J. Numer. Anal. , 51:2309–2330, 2012.
29] L.A Liotta, M.G. Saidel, and J. Kleinerman. Diffusion model of tumor vascu-larization and growth.
Bull. Math. Biol. , 39:117–128, 1977.[30] C.D. Madsen and E. Sahai. Cancer dissemination–lessons from leukocytes.
Dev. Cell , 19:13–26, 2010.[31] Ch. Makridakis, D. Mitsoudis, and Ph. Rosakis. On atomistic-to-continuumcouplings without ghost forces in three dimensions.
Appl. Math. Res. Express ,1:87–113, 2014.[32] A. Marciniak-Czochra and M. Ptashnyk. Boundedness of solutions of a hap-totaxis model.
Math. Mod. Meth. Appl. S. , 20:449–476, 2010.[33] MATLAB.
MATLAB version 8.6.0 (R2015b) . The MathWorks Inc., Natick,Massachusetts, 2015.[34] P. Mehlen and A. Puisieux. Metastasis: a question of life or death.
Nat. Rev.Cancer , 6:449–58, 2006.[35] M. Nystr¨om, G.J. Thomas, I.C. Stone, M. Mackenzie, I.R. Hart, and J.F. Mar-shall. Development of a quantitative method to analyse tumour cell invasionin organotypic culture.
J. Pathol. , 205:468–475, 2005.[36] M.E. Orme and M.A.J. Chaplain. Two-dimensional models of tumour an-giogenesis and anti-angiogenesis strategies.
IMA J. Math. Appl. Med. Biol. ,14:189–205, 1997.[37] K.J. Painter and T. Hillen. Mathematical modelling of glioma growth: the useof Diffusion Tensor Imaging (DTI) data to predict the anisotropic pathwaysof cancer invasion.
J. Theor. Biol. , 323:25–39, 2013.[38] L. Peng, D. Trucu, P. Lin, A. Thompson, and M.A.J. Chaplain. A multiscalemathematical model of tumour invasive growth.
Bull. Math. Biol. , 79:389–429,2017.[39] L. Preziosi.
Cancer Modelling and Simulation . Chapman and Hall/CRC, 2003.[40] I. Ramis-Conde, M.A.J. Chaplain, and A.R.A. Anderson. Mathematical mod-elling of cancer cell invasion of tissue.
Math. Comput. Modell. , 47:533–545,2008.[41] I. Ramis-Conde, D. Drasdo, M.A.J. Chaplain, and A.R.A. Anderson. Mod-elling the influence of the e-cadherin - β -catenin pathway in cancer cell invasionand tissue architecture: A multi-scale approach. Biophys. J. , 95:155–165, 2008.[42] P.A. Raviart. Particle numerical models in fluid dynamics. In K.W. Mortonand M.J. Baines, editors,
Numerical Methods for Fluid Dynamics , pages 231–253. Clarendon Press, Oxford, 1986.[43] F. Sabeh, R. Shimizu-Hirota, and S.J. Weiss. Protease-dependent versus -independent cancer cell invasion programs: three-dimensional amoeboid move-ment revisited.
J. Cell Biol. , 185:11–9, 2009.
44] D.K. Schl¨uter, I. Ramis-Conde, and M.A.J. Chaplain. Multi-scale modellingof the dynamics of cell colonies: insights into cell adhesion forces and cancerinvasion from in silico simulations.
J. R. Soc. Interface , 12, 2015.[45] N. Sfakianakis, N. Kolbe, N. Hellmann, and M. Lukacova-Medvidova. A multi-scale approach to the migration of cancer stem cells: Mathematical modellingand simulations.
Bull. Math. Biol. , 79:209–235, 2017.[46] C.L. Stokes and D.A. Lauffenburger. Analysis of the roles of microvessel en-dothelial cell random motility and chemotaxis in angiogenesis.
J. Theor. Biol. ,152:377–403, 1998.[47] R.I. Stratonovich. A new representation for stochastic integrals and equations.
SIAM J. Control , 4(2):362–371, 1966.[48] A.F.B. Thompson and D.E. Dougherty. Patricle-grid methods for reactingflows in porous media with application to Fisher’s equation.
Appl. Math.Modeling , 16(7):374–383, 1992.[49] L. Tweedy, D.A. Knecht, G.M. Mackay, and R.H. Insall. Self-generatedchemoattractant gradients: Attractant depletion extends the range and ro-bustness of chemotaxis.
PLOS Biology , 14(3):1–22, 2016.[50] A. Valster, N.L. Tran, M. Nakada, M.E. Berens, A.Y. Chan, and M. Symons.Cell migration and invasion assays.
Methods , 37:208–215, 2005.[51] H.A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant ofBi-CG for the solution of nonsymmetric linear systems.
SIAM J. Sci. Comput. ,13(2):631–644, 1992.[52] J. Wang, L. Zhang, C. Jing, G. Ye, H. Wu, H. Miao, Y. Wu, and X. Zhou.Multi-scale agent-based modeling on melanoma and its related angiogenesisanalysis.
Theor. Biol. Med. Model. , 10:41, 2013.[53] B. Weigelt, J.L. Peterse, and L.J. van’t Veer. Breast cancer metastasis: mark-ers and models.
Nat. Rev. Cancer , 8:591–602, 2005.[54] M. Winkler and Y. Tao. Energy-type estimates and global solvability in a two-dimensional chemotaxis–haptotaxis model with remodeling of non-diffusibleattractant.
J. Diff. Eq. , 257:784–815, 2014.[55] F.W. Yang, Ch. Venkataraman, V. Styles, V. Kuttenberger, E. Horn, Z. vonGuttenberg, and A. Madzvamuse. A computational framework for particleand whole cell tracking applied to a real biological dataset.
J. Biomech. ,49(8):1290–1304, 2016. (Numerical method for the ARD model (2) ) We use a second order
Implicit-Explicit Runge-Kutta (IMEX-RK)
Finite Volume (FV) numerical method that was previously developed in [26, 45] where we refer formore details, see also [28]. Here we provide some basic description of the method.We consider a generic ARD system of the form w t = A ( w ) + R ( w ) + D ( w ) , (44)where w represents the solution vector, and A , R , and D the advection , reaction ,and diffusion operators respectively.We denote by w h ( t ) the corresponding (semi-)discrete numerical approximation —indexed here by the maximal spatial grid diameter h — that satisfies the system ofODEs ∂ t w h = A ( w h ) + R ( w h ) + D ( w h ) , (45)where the numerical operators A , R , and D are discrete approximations of theoperators A , R , and D in (44) respectively.Our method of choice for solving (45) is an Implicit-Explicit Runge-Kutta (IMEX-RK) method based on a splitting in explicit and implicit terms in the form ∂ t w h = I ( w h ) + E ( w h ) . (46)The actual splitting depends on the particular problem in hand but in a typicalcase, the advection terms A are treated explicitly in time, the diffusion terms D implicitly, and the reaction terms R partly explicit and partly implicit.More precisely, we employ a diagonally implicit RK method for the implicit part,and an explicit RK for the explicit part W ∗ i = w nh + τ n i − (cid:88) j =1 ¯ a i,j E j + τ n ¯ a i,i − E i − , i = 1 . . . s W i = W ∗ i + τ n i − (cid:88) j =1 a i,j I j + τ n a i,i I i , i = 1 . . . s w n +1 h = w nh + τ n s (cid:88) i =1 ¯ b i E i + τ n s (cid:88) i =1 b i I i , (47)where s = 4 are the stages of the IMEX method, E i = E ( W i ), I i = I ( W i ), i = 1 . . . s , { ¯ b, ¯ A } , { b, A } are respectively the coefficients for the explicit and theimplicit part of the scheme, given by the Butcher Tableau in Table 4, [23]. We solvethe linear systems in (47) using the iterative biconjugate gradient stabilised Krylovsubspace method [27, 51]. − − − − − B (An explicit numerical scheme for the SDE (7) ) We consider an Ito process X = { X t , t ≤ t ≤ T } that satisfies the geometricBrownian motion SDE d X t = α X t dt + β X t d W t , (48)where X t denotes the position in space, and where α ∈ R and β > X n +1 = X n + α X n τ + β X n ∆ W t . (49)By setting ∆ W t = Z √ τ with Z ∼ N (0 , X n +1 = X n + α X n τ + β X n Z √ τ (50)which is a simpler version of the Itˆo-type scheme that we employ in (8).For further details on the numerical treatment of (48) and other SDEs we refer to[22, 25].(50)which is a simpler version of the Itˆo-type scheme that we employ in (8).For further details on the numerical treatment of (48) and other SDEs we refer to[22, 25].