Aggregation and travelling wave dynamics in a two-population model of cancer cell growth and invasion.

Cells adhere to each other and to the extracellular matrix (ECM) through protein molecules on the surface of the cells. The breaking and forming of adhesive bonds, a process critical in cancer invasion and metastasis, can be influenced by the mutation of cancer cells. In this paper, we develop a nonlocal mathematical model describing cancer cell invasion and movement as a result of integrin-controlled cell-cell adhesion and cell-matrix adhesion, for two cancer cell populations with different levels of mutation. The partial differential equations for cell dynamics are coupled with ordinary differential equations describing the ECM degradation and the production and decay of integrins. We use this model to investigate the role of cancer mutation on the possibility of cancer clonal competition with alternating dominance, or even competitive exclusion (phenomena observed experimentally). We discuss different possible cell aggregation patterns, as well as travelling wave patterns. In regard to the travelling waves, we investigate the effect of cancer mutation rate on the speed of cancer invasion.


Introduction
Normal cells proliferate, divide and die in a highly controlled manner. This proliferation requires mitogenic signals, which are transmitted into the cell by the transmembrane receptors and bind signalling molecules, i.e. diffusible growth factors, extracellular matrix (ECM) components and cellcell adhesion molecules (Hanahan & Weinberg, 2000). When a cell divides, its DNA is copied by the two new cells. However, if the DNA is damaged or not copied correctly, the new damaged cells will either die or start to proliferate in an uncontrolled manner, creating a signalling of oncogenes that act by mimicking growth signalling (Hanahan & Weinberg, 2000), eventually leading to an abnormal mass of tissues from cells that differ in clinically important phenotypic features (Marusyk et al., 2012). More precisely, tumour formation is a result of clonal expansion driven by somatic mutation, developed by a single precursor (monoclonal) that undergoes genetic and biological changes (Khalique et al., 2007;Nowell, 1976). This transformation of normal cells to cancer cells is a multistep process described via 4 V. BITSOUNI ET AL. sub-populations have been shown to be distinct not only in their adhesion capabilities but also in their motility and metastatic potential (Marusyk & Polyak, 2010). Moreover, we are interested in investigating the effect of mutation on the spreading speed of the mutated cancer cells, since the invasion levels of the mutant cancer cells are higher compared to those of 'normal' cancer cells (Chapman et al., 2014;Wojciechowska & Patton, 2015). Also, while previous nonlocal models have assumed constant adhesive interactions, here we will assume that these adhesive interactions depend on the level of integrins. Note that instead of assuming a kinetic (mesoscale) approach for cell-dependence on integrins as in Lorenz & Surulescu (2014), Engwer et al. (2015), Perthame et al. (2016), Engwer et al. (2017), Hunt & Surulescu (2017), here we follow a microscale-macroscale approach similar to Stinner et al. (2014Stinner et al. ( , 2015Stinner et al. ( , 2016 and Meral et al. (2015a,b) and couple the parabolic-hyperbolic model for spatial cell dynamics with an Ordinary Differential Equation (ODE) for integrins dynamics. We then investigate this complex multiscale model in terms of pattern formation, with particular attention being payed to the role of mutation rate on the coexistence (or not) of cell sub-populations that form stationary aggregations or travelling wave patterns.
The layout of this paper is as follows. In Section 2 we formulate a model of partial integrodifferential equations for the dynamics of two cancer cell populations, coupled with ordinary differential equations describing ECM and integrins dynamics. In Section 3 we use linear stability analysis to investigate the ability of the model to form cell aggregations. In Section 4 we investigate numerically some types of spatio-temporal patterns exhibited by this nonlocal model. In Section 5 we investigate numerically and analytically travelling wave behaviours. We conclude in Section 6 with a summary and discussion of the results.

Derivation of the model
Cancer invasion is a complex process involving cell-cell interactions, and interactions between cells and non-cellular components (Calvo & Sahai, 2011). To develop our model we consider two populations of cancer cells, a 'normal' (original clone) cancer population and a 'mutant' (descendant clone) population, which interact with each other as well as with the ECM via long-range adhesive and repulsive forces (Deman et al., 1976;Geiger, 1991)

5
Let Ω ⊂ R d denote a bounded spatial domain (here we consider only d = 1; for a 2D version of the model see Appendix A), with periodic boundary conditions. Let I T = [0, ∞) be the time interval. Denote by u 1 (x, t) the density of 'normal' cancer cells at position x and time t, and by u 2 (x, t) the density of 'mutant' cancer cells at position x and time t. We also denote by f (x, t) the ECM density, and by c (x, t) the density of integrin receptors on the surface of cancer cells (receptors involved in cell-cell and cell-matrix interactions). For compact notation, we define the vectors u (x, t) = (u 1 (x, t) , u 2 (x, t)) and υ (x, t) = u (x, t) , f (x, t) .
Dynamics of cancer cells. Experimental studies (Chapman et al., 2014) have shown that cancer cells can switch from a homogeneous type of invasion to a heterogeneous type of invasion described by invading chains (Chapman et al., 2014;Friedl & Wolf, 2003;Wojciechowska & Patton, 2015). Here, we assume that u 1 cells can mutate into u 2 cells at a constant rate M. We can derive a model for heterogeneous cancer cell populations by considering the equations where J i , i = 1, 2, are the cell fluxes and G i u , i = 1, 2, are the growth functions of populations u i , i = 1, 2. The fluxes describe the factors that define cell movement. In this study we assume that the movement of the 'normal' cancer cell population u 1 is governed by random motility (which underlines a homogeneous type of invasion), as well as directed motility in response to cell-cell and cell-matrix adhesive forces (which underlines the heterogeneous type of invasion) (Calvo & Sahai, 2011). In contrast, the 'mutant' u 2 cell population moves only in a directed manner (hence exhibiting an heterogeneous type of invasion) in response to cell-cell and cell-matrix adhesion forces. Biologically, this directed movement of 'mutant' (more invasive) cancer cells can be explained by the increase of macrophage density near highly mutated cancer cells (Lin et al., 2006), which promotes the directed movement and invasion of these cancer cells, and decreases their random movement (Goswami et al., 2005;Hagemann et al., 2005). Therefore, in this study we assume that the movement of u 2 population due to random walk can be considered negligible. In this case, the total fluxes will be where J D is the flux due to Fickian diffusion, given by J D = −D∂u 1 /∂x, with D to be the diffusion coefficient, and J a i , i = 1, 2, are the adhesive fluxes for cell-cell and cell-matrix adhesion. In the case of u 2 the diffusion is considered very small and thus we ignore it. Cell-cell/cell-matrix adhesion-mediated directed cancer cell migration occurs as a result of the social forces, i.e. attraction and repulsion, between cells and between cells and ECM components, when adhesive bonds are formed and broken. Cell-cell adhesion is described as the adhesion between cells of the same population, as well as cross-adhesion between cells of the two sub-populations. The cellcell/cell-matrix adhesion forces are created through the binding of adhesion molecules, e.g. integrins, at cell surface. To model these forces, let us define R > 0 to be the sensing radius of the cells, i.e. the maximum range over which cells can detect other surrounding cells, which biologically represent the extent of the cell protrusions (e.g. filopodia) (Armstrong et al., 2006). Let g i υ (x + r, t) , c (x, t) i = 1, 2, describe the nature of the cell-cell and cell-matrix adhesive forces created at x due to signalling 6 V. BITSOUNI ET AL.
with cell/ECM components at x + r. These functions increase when the cell density and ECM density increase, and accordingly they decrease when the cell density and ECM density decrease. The functions g i , i = 1, 2, are chosen as in Bitsouni et al. (2017) to be given by and is the cell-cell cross-adhesion strength function between the two populations, and C i (c (x, t)) is the adhesion strength function between population u i and ECM. We should mention here that a similar term was considered before by Chaplain et al. (2011). Other studies (Armstrong et al., 2006;Gerisch & Chaplain, 2008) consider volume filling effect to avoid unbounded aggregations. While unbounded aggregations are possible in this type of nonlocal models, these have not been observed in our system (at least not for the parameter ranges used in this study). We define these adhesion strength functions in terms of the integrin density c: the more integrins a cell has, the stronger its adhesion force (Gallant et al., 2005). However, the adhesion strength reaches a plateau for a large integrin density (Maheshwari et al., 2000). Thus, we require an increasing, bounded function, so that we avoid extreme phenomena for large values of integrin density. A biologically realistic choice is a sigmoid function. This choice is supported by various other experimental studies in the literature which have shown that different biological processes can be described by sigmoid functions: from cell adhesion via integrins (Cutler & García, 2003;Michael et al., 2009), to oxygen saturation (Shiao & Ou, 2007). Note also that since cell mutation could lead to more integrins (Kidera et al., 2010), we consider strength functions with different integrin levels for each of the two populations. To this end, we choose the adhesion strength functions to be given by: where d, a i , b i , i = 1, 2, and s*, s i *, c i *, i = 1, 2, are positive real numbers. Within this sensing radius, let us now define K cc and K cm , the kernels for cell-cell and cell-matrix adhesion ranges, respectively: These interaction kernels describe the dependency of the force magnitude on the distance r from x.
Here K ij , i, j = 1, 2, are the interaction kernels between population i and population j (i.e. cell-cell 7 interactions), and K i , i = 1, 2, are the kernels for cell-matrix interactions (for the two populations). For simplicity, throughout this study we will consider (see also Domschke et al., 2014): Moreover, we assume that these kernels are attractive at medium/long ranges (i.e. at the edges of the cell) and repulsive at very short ranges (i.e. over cell surface), and thus can be defined as with q a and q r describing the magnitudes of attractive and repulsive interactions, and K 1,2 a (r) and K 1,2 r (r) describing the spatial ranges over which these interactions take place. We will discuss various examples of attractive (K a ) and repulsive (K r ) kernels in Section 3, in the context of linear stability analysis.
Therefore, the local cell-cell and cell-matrix adhesion forces will be described by the product of functions g i and K i (r) , i = 1, 2. To describe the nonlocal cell-cell and cell-matrix adhesion and repulsion forces (Deman et al., 1976;Geiger, 1991), we sum all such forces by integrating over space to obtain the total forces: where η (k) = (−1) k , k = 0, 1, is the direction of the forces. Thus, the full nonlocal interaction terms are described by and We assume that the adhesive fluxes are proportional to the density of the cells and the nonlocal adhesion forces, F i , i = 1, 2. Thus, we obtain the following expressions for the two adhesive fluxes: (2.12) Substituting (2.12) into the general mass conservation equations, we have the following equations describing the dynamics of the two cancer cell populations: We assume that both u 1 and u 2 cells can proliferate in a logistic manner (to describe the observed slowdown in tumour growth following the loss of nutrients (Laird, 1964)). Thus, the growth functions are given by where r 1 and r 2 are the growth rates of the u 1 and u 2 populations, respectively, and k u is the carrying capacity. Note that these growth functions incorporate also the principle of competition between clonal sub-populations in heterogeneous tumours (Leith et al., 1989).
ECM dynamics. Cancer cell populations degrade the ECM upon contact. Moreover, ECM density is remodelled back to normal levels by cancer and healthy cells. Thus, the dynamics of ECM, f (x, t), can be described by: where α and β are the positive rate constants of ECM degradation by u 1 and u 2 cell populations, respectively, k u is the carrying capacity of the cancer cell populations and f m is the maximum ECM density at which the ECM fills up all available physical space. Following the approach of Anderson et al. (2000); Gerisch & Chaplain (2008), we choose to ignore the ECM remodelling terms for this model, i.e. we always use γ = 0.
Integrin dynamics. We assume that the level of integrins depends on cancer cell density. Moreover, we assume that cell mutation changes the density of receptors (since in highly metastatic cancers, the expression of integrins is up-regulated (Kidera et al., 2010)). We follow an approach similar to (Stinner et al., 2014(Stinner et al., , 2015(Stinner et al., , 2016Meral et al., 2015a,b) and assume that the dynamics of integrins c(x, t) can be described by an ODE: where q is the decay rate of c, and p 1 and p 2 are the production rates of integrins by population u 1 and u 2 , respectively. To model the increase in receptors on mutant cancer cells, we assume that p 2 > p 1 (see Table 2). Note that integrins are produced at the positions where these cells exist, and as cells move to different spatial positions the integrins will be produced at these new positions, too (this can be seen more clearly in Figs 4-10). Moreover, the production of integrins at various positions in space is counterbalanced by the decay term (−qc), which ensures that integrins will be removed from regions where cells do not exist. For the initial conditions, we start with only population u 1 (x, 0) = u 0 1 (x). The second population u 2 , arises from mutations in the population u 1 after a period of time, and thus u 2 (x, 0) = 0. For the ECM, we assume that the tumour has already degraded some of its surrounding tissues: (2.17) Finally, the integrin density, c, is proportional to the initial tumour cell density To complete our model we need to impose boundary conditions. Throughout this study we assume that the two cancer cell populations move on a circular finite domain of length L, i.e. x ∈ [0, L]. The periodic boundary conditions that model this movement are Note that with this choice of boundary condition, the nonlocal terms are wrapped around the domain. We will return to this aspect in Section 4.
Remark 2.1 Note that equations (2.13)-(2.16) describe a 'microscale-macroscale' model, where the spatially-averaged microscale integrin dynamics (described by an ODE) is coupled with the spatial dynamics of cancer cells (described by nonlocal PDEs). As mentioned above, this modelling approach follows various other recent studies on cancer cell-matrix adhesion via integrins (Stinner et al., 2014(Stinner et al., , 2015(Stinner et al., , 2016Meral et al., 2015a,b). Moreover, the coupling of PDEs for cell movement with ODEs for the dynamics of cell receptors has been considered not only for cancer cells but also for other types of cells, such as 'Dictyostelium' cells (Höfer et al., 1995a,b). However, the past 10-15 years have seen the development of structured cell population models (described by 'mesoscale' models) where the macroscopic cell dynamics is dependent on a microscopic variable that characterizes the cells at molecular level (Erban & Othmer, 2004Xue et al., 2011;Lorenz & Surulescu, 2014;Engwer et al., 2015;Perthame et al., 2016;Engwer et al., 2017;Hunt & Surulescu, 2017). In Appendix B we present a modified mesoscale version of model (2.13)-(2.16) where we assume that the two cancer cell populations are structured by an internal variable that describes the integrin level. (Here we use the terminology of 'mesoscale' and 'microscale-macroscale' models in the same spirit as Stinner et al.

Non-dimensionalization of the model
To non-dimensionalize system (2.13) and equations (2.15)-(2.16), we define the following quantities: Here, L 0 is a length scale, defined as the maximum invasion distance of the cancer cells at the 'normal' of invasion (Anderson et al., 2000). Usually L 0 is in the range of 0.1-1 cm. We rescale time with τ := L 2 0 /D τ , where D τ is the characteristic diffusion coefficient (∼ 10 −6 cm 2 s −1 ). In addition, we rescale the cancer cells, the ECM and the integrins with k u , f m and c m , respectively. Here k u is the carrying capacity of the cancer cell populations and it is taken to be ∼ 6.7 · 10 7 cell/volume, and f m is the maximum ECM density at which the ECM fills up all available physical space and it is taken to be equal to 4 mg/volume, as in Domschke et al. (2014). Finally, c m is the maximum integrin density and it is taken to be 5 · 10 4 integrins per cell (as in Benedetto et al. (2006)). For the kernels K i (r) , i = 1, 2, we choose (as in Domschke et al. (2014)) the dimensionless functionsK i (r) , i = 1, 2, given bỹ Therefore, we have for the nonlocal terms Finally, we obtain the dimensionless parameters: After dropping the tildes for notational convenience, we obtain the following non-dimensionalized system From now on, we will always refer to the non-dimensional quantities. In the next two sections, we discuss the dynamics of model (2.22).

Linear stability analysis of the model
Next, we investigate the conditions under which the two populations of cancer cells form aggregations. We first calculate the spatially homogeneous steady states of the model, and then apply standard linear stability analysis to investigate the conditions (for parameter values) under which the two cancer cell populations can form aggregations.
To start, we look for the homogeneous steady states of the ODE model associated to system (2.22), that describes the growth and mutation of the two cancer cell populations and the temporal dynamics of ECM and integrins (i.e. no spatial movement): which has the following steady state solutions u * 1 , u * 2 , f * , c * : Here we consider only the non-negative solutions, since we require biological realism. Note that all these homogeneous steady states have u 1 = 0 (so the more invasive population u 2 persists longer). However, as we will see in Section 4, the model can exhibit non-homogeneous steady states with u 1 (x, t) = 0.
We now investigate the stability of the steady states for the spatially homogeneous and spatially heterogeneous systems, to see if the introduction of spatial dynamics (i.e. cell movement and cellcell/cell-matrix adhesion) can lead to instabilities. First, we substitute the steady states (3.2) into 12 V. BITSOUNI ET AL.
the Jacobian matrix of the spatially homogeneous system (3.1) and calculate the four eigenvalues corresponding to each steady state. Thus, for the steady state (0, 0, f * , 0) we have the eigenvalues Since λ 3 > 0, this state is unstable. For the steady state (0, 1, 0, p 2 /q) we obtain the eigenvalues Since all eigenvalues are negative, the state (0, 1, 0, p 2 /q) is always linearly stable to homogeneous perturbations.
Next, we investigate the linear stability analysis of the spatial system (2.22). To this end, we apply small spatial perturbations to the previous homogeneous steady states: , whereū 1 ,ū 2 ,f andc denote the small perturbations. (Note that, to avoid negative solutions when we perturb the zero steady states, we considerū 1 ( 0.) Substituting these into the system (2.22) (after the adhesion strength functions S i (c), S(c), C i (c), i = 1, 2, have been expanded in Taylor series), and ignoring non-linear terms, leads to the following linearized system: We look for solutions of the form A u 1 e ikx+λt , A u 2 e ikx+λt , A f e ikx+λt and A c e ikx+λt with |A u 1 |, |A u 2 |, |A f |, |A c | 1, where k and λ are the wave number and frequency, respectively. Then system (3.5) reduces to dr are the Fourier sine transforms of the kernels K 1,2 (r). For simplicity, throughout the rest of this study, we will assume that K 1 (r) = K 2 (r) =: K(r) and thusK s Cellular aggregations form when the steady states (u * 1 , u * 2 , f * , c * ) are unstable to spatial perturbations, and for this to happen we require the maximum real part of eigenvalues of the Jacobian matrix of system (3.6) to be positive, for at least one discrete value of k > 0.
• For the state u * 1 , u * 2 , f * , c * = (0, 0, f * , 0) we obtain the eigenvalues Thus, the steady state u * 1 , u * 2 , f * , c * = (0, 0, f * , 0) is unstable. • For the state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q) we obtain the eigenvalues Recall that in the absence of diffusion and advection, the steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q) is linearly stable. Thus, for aggregation to form, we require this state to be unstable to spatially inhomogeneous perturbations (Benson et al., 1993), i.e. the maximum real part of the eigenvalues to be greater than zero. Since λ 1 , λ 2 , λ 4 < 0, the steady state (0, 1, 0, p 2 /q) is unstable if the following dispersion relation holds: We notice that this dispersion relation, and therefore the possibility of cell aggregations to form, depends on the choice of the interaction kernel. Mogilner & Edelstein-Keshet (1999) showed that the type of kernel affects the possibility of movement and/or aggregation. Throughout the rest of this study, 14 V. BITSOUNI ET AL.
we will consider translated Gaussian attraction and repulsion kernels (as in Eftimie et al. (2007); Carrillo de la Plata et al. (2015)): where s a and s r represent half of the length of attraction and repulsion ranges, respectively, with s r < s a . Also, m j = s j /8, j = a, r, represent the width of the attractive and the repulsive interaction ranges. (The constants m j , j = a, r, are chosen such that the support of more than 98% of the mass of the kernels is inside the interval [0, ∞) (Eftimie et al., 2007).) In Section 3.1 we will investigate how different types of kernels lead to different dispersion relations. The Fourier sine transform of kernel (3.10) is given bŷ (3.11) Therefore, relation (3.9) becomes where the imaginary part of the eigenvalue λ 3 (k) will be zero.
An example of such dispersion relation is shown in Fig. 2. Note the steady-state/steady-state mode interaction that occurs when we increase the parameters q a and q r (i.e. the strength of attractive and repulsive cell-cell interactions). There is a range of k −values for which Re (λ 3 (k)) is positive, and thus aggregations can arise from spatial perturbations of the steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q). This means that the steady state (0, 1, 0, p 2 /q), which is linearly stable for the homogeneous system, is destabilized due to the spatial dynamics (diffusion and advection). Setting diffusion coefficient zero we observe that there is still a range of k −values for which Re (λ 3 (k)) > 0, implying the adhesion-driven instability.

Aggregation with different types of kernel
In Mogilner & Edelstein-Keshet (1999) it was shown that even kernels give rise to group drift, while odd kernels have a greater effect in regions where the distribution of the density is uneven (e.g. the edges of the population), leading to stationary groups. Considering the importance of the symmetry of the interaction kernels, we next discuss the cases of even and odd kernels, and compare these results with our previous results where we used a translated Gaussian kernel (see equation (3.10)), to investigate the possibility of cell aggregations to form from spatial perturbation of the steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q).  Table 2. The diamonds on the x-axis represent the discrete wave numbers k j = 2π j/L, j = 1, 2, . . . . The two critical wave numbers that become unstable at the same time (giving rise to steady-state/steadystate mode interactions) are k 12 and k 13 .

• Aggregation with even kernel
We consider the situation that the interaction kernel is an even function with respect to the y-axis. We choose the following Gaussian kernel (Mogilner & Edelstein-Keshet, 1999): where q a , q r , m a and m r have the same meaning as in equation (3.10). The Fourier sine transform of the above even kernel is zero:K s (k) = 0. We know (see equation (3.8)) that the Jacobian matrix of system (3.6) has three negative eigenvalues λ 1,2,4 < 0. For the even kernel considered here we also have λ 3 (k) = −r 2 < 0. (3.14) Therefore, the steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q) is stable, and no aggregations will form.

• Aggregation with an odd kernel
We investigate now the stability of u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q) when we consider an odd kernel with respect to the origin. To this end, we choose   Table 2. The solid curve represents the real part of λ 3 , while the diamonds represent the discrete wave numbers k j = 2π j L , j = 0, 1, 2, . . . .
where q a , q r , m a and m r have the same meaning as in equation (3.10). Then the Fourier sine transform of the kernel (3.15) iŝ As before, the Jacobian matrix of system (3.6) has three negative eigenvalues: λ 1,2,4 < 0. The stability of the state (0, 1, 0, p 2 /q) is given by the sign of Re (λ 3 (k)) = −r 2 + 2k 2 s * 2 R 1 + tanh a 2 p 2 q q a m 2 a e − (kma) 2 2 − q r m 2 r e − (kmr ) 2 2 . (3.17) In Fig. 3 we plot Re(λ 3 (k)) against the wave number k for this steady state. We observe that for values of q a greater than those in the translated Gaussian case (see Fig. 2), there is a range of k-values for which Re (λ 3 (k)) is positive, thus the steady state (0, 1, 0, p 2 /q) can become unstable, and spatial aggregations could arise.

Numerical results
In this Section we investigate numerically the stability of model (2.22). To discretize our model we use a time-splitting approach. First, we use a Crank-Nicolson scheme to propagate the solution of the diffusion term. Then, we use the Nessyahu-Tadmor scheme (Nessyahu & Tadmor, 1990) for the timepropagation of the advection terms. Finally, for the time-propagation of the reaction terms we use a forth order Runge-Kutta algorithm, where the integrals are further discretized using the Simpson's rule. All simulations are performed on a domain of length L = 10 with periodic boundary conditions, and for time up to t = 500. Because of the periodic boundary conditions, we wrap the integrals around the domain. The initial conditions for the cancer cell populations are either small random perturbations of spatially homogeneous steady states or small random perturbations of rectangular-shaped aggregations located in the middle of the domain with u c 1 = 0 and u c 2 as specified in the caption of the figures. The initial conditions for the ECM and integrin density depend on cancer cell density, and are given by relations (2.17) and (2.18), for all numerical simulations.

Stationary aggregation patterns.
To check the validity of our results obtained via linear stability analysis (see Fig. 2), we first run simulations for small random perturbations of the spatially homogeneous steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q), and for q a = 0.09 and q r = 0.01. The results, presented in Fig. 4, show 12-13 stationary pulses (corresponding to critical wave numbers k 12 − k 13 ). The spread of cancer cells over the whole domain (panels (a), (b)) leads to the degradation of ECM (panel (c)). The integrin density (panel (d)) follows the patterns in the density of cancer cells. We observe here a coexistence between a low u 1 population and a high u 2 population.
To check the effect of initial data on model dynamics, in Fig. 5 we investigate the behaviour of model (2.22) when we change the initial conditions to a rectangular pulse describing an initial aggregation of tumour cells. The magnitudes of attractive-repulsive interactions (and all other parameter values) are the same as in Fig. 4. The short-term results in Fig. 5(a)-(d) show a travelling pulse behaviour exhibited by the u 2 population that moves away from population u 1 (which is mostly stationary-although even this population spreads a bit). This faster spreading of the u 2 cells compared to the u 1 cells is consistent with experimental observations (Chapman et al., 2014;Wojciechowska & Patton, 2015). The movement of u 2 population eventually stops near the boundary (due to the periodic boundary conditions, the leftmoving and right-moving subgroups can sense each other across the boundary). The fast-moving u 2 cells degrade the ECM (panel (c)). We also note the high density of integrins associated with the u 1 population (panel (d)). The long-term results in Fig. 5(a )-(d ) show the formation of new aggregations of cells at distant positions in space. For t ∈ (100, 250), the two tumour populations coexist. For t > 250, population u 1 is slowly eliminated and population u 2 dominates the dynamics of model (2.22). We also note that despite some chaotic-like dynamics exhibited by the u 1 and u 2 populations for t < 300, in the long term the system approaches stationary pulses defined by 13 peaks (corresponding to unstable wave number k 13 ).
In Figs 6 and 7 we investigate the effect of various mutation rates on the dynamics of the u 1 and u 2 populations. In Fig. 6 we choose M = 0.001 (see Table 2), and observe that population u 2 vanishes for t > 400, while population u 1 persists and forms small high-density aggregations of cells. This unexpected behaviour might be explained by the combined effect of three factors: (i) the mutation rate (M = 0.001) that is much smaller than the growth rates (r 1,2 = 0.1) of the two cancer cell populations, 18 V. BITSOUNI ET AL. (ii) the lack of diffusion for u 2 and (iii) the competition for nutrients between the two cell populations (which is modelled implicitly by the logistic terms G 1,2 ; see equations (2.14)). This eventually leads to an overall increase in the first population, u 1 , and a decrease in the second (highly mutated) population, u 2 . In Fig. 7 we increase the mutation rate to M = 0.05, while chosing the growth rates fixed at r 1 = r 2 = 0.1. We observe that faster mutation rates lead to a decrease and eventual elimination of the u 1 population. The u 2 population increases and dominates the long-term dynamics of model (2.22).
Travelling wave patterns. Choosing again rectangular pulse initial conditions, we reduce the magnitudes of cell-cell attractive and repulsive interactions to q a = 0.00025 and q r = 0.0005. The numerical simulations in Fig. 8 show travelling waves that propagate in opposite directions at a constant speed. Similar behaviour is obtained for q a > q r . The evolution of the travelling waves is shown in Fig. 9 for t = 25j, j = 1, . . . , 9. We observe that the waves connect the unstable steady states (0, 1, 0, p 2 /q), with p 2 /q = 0.5, and (0, 0, f * , 0), with f * = 1.
Finally, we check the effect that an increase in the mutation rate M has on the travelling wave behaviour exhibited by model (2.22). Figure 10 shows that for high mutation rates (M = 0.05) the first population still exhibits a travelling wave before vanishing.  (2.22). The initial conditions for the two cancer cell populations are described by a rectangular pulse (see (4.1)) with u c 2 = 1.0. We use the interaction kernel given by (3.10), and the model parameters given in Table 2     Remark 4.1 To understand the long-term behaviour of model (2.22) with odd kernels, we also ran numerical simulations using the interaction kernel given by (3.15). The spatial patterns obtained in this case were similar to the patterns obtained for the translated Gaussian kernel (see relation (3.10)).

Summary of model variables and parameters
Here we present two tables with the model variables and parameters. In Table 1 we list the model variables with their units. In Table 2 we list the parameters of our model and their corresponding units and non-dimensional values used in the simulations.
• The sensing radius was based on the range of values given in (Armstrong et al., 2006;Gerisch & Chaplain, 2008). In this study, we chooseR = 0.99.
• Attraction and repulsion ranges were chosen to be smaller or equal to sensing radius, with the repulsion range to be smaller than the attraction range (Green et al., 2010).

22
V. BITSOUNI ET AL. Fig. 9. The evolution of the travelling waves under small perturbations for q r = 0.0005, q a = 0.00025 and the rest of the model parameters given in Table 2, of (a) cancer cell density u 1 of the first ('normal') population; (b) cancer cell density u 2 of the second ('mutant') population; (c) ECM density f ; (d) integrin density c.
• Experimental studies (Cillo et al., 1987;Hill et al., 1984;Mareel et al., 1991) have shown that the mutation rate ranges between M = 10 −3 /day and M = 0.1/day. Thus the non-dimensional value of the mutation rate is in the range betweenM = 0.001 andM = 0.1 (for highly aggressive tumours).
• The parameters a i , d, b i , s * i , s * , c * i , i = 1, 2, were based on the range of the adhesion strength parameters used in Armstrong et al. (2006).

Speed of travelling waves
The numerical simulations showed that the model exhibits travelling wave solutions (see Figs 8 and 9). In this section we estimate the minimum speed at which these waves could propagate (assuming that the travelling waves do exist). We assume that close to the invading front, the nonlinear differential equations describing the spread of a population have similar speeds as their linear approximation (see the approaches in Medlock & Kot (2003); Mollison (1991)). Therefore, we first calculate the linearized system (3.5) at the steady state u * 1 , u * 2 , f * , c * = (0, 1, 0, p 2 /q). In travelling wave coordinates Let s 1 , s 2 , s 3 , s 4 > 0 be the steepness (i.e. exponent of an exponential decay/increase profile) of u 1 , u 2 , f , c. The minimum of w s j for s j > 0, j = 1, . . . , 4, gives an upper bound on the rightward travelling wave speed of the nonlinear system. Thus, we make the general exponential ansatz First we investigate the s j , j = 1, . . . , 4, numerically. In Fig. 11 we fit these exponential ansatz profiles (given by (5.2)) to the numerical solutions u 1 (x, t), u 2 (x, t), f (x, t) and c(x, t) at time t = 175 (dotted curve). We obtain the following ansatz profiles (continuous curve in Fig. 11): Hence,we deduce that the travelling wave profiles of U 1 (z) , U 2 (z) and C (z) have the same steepness s := s 1 = s 2 = s 4 = 3.17, while the profile of F (z) has steepness s 3 = 6. Therefore, we replace the exponential ansatz given by (5.2) with U 1 (z) = A 1 e −sz , U 2 (z) = A 2 e −sz , F (z) = A 3 e s 3 z and C (z) = A 4 e −sz . Let us assume that K (r) has a moment generating function,M (s), defined as   M (s) := +∞ −∞ K (r) e −sr − e sr dr (and such a function exists for the translated Gaussian kernels (3.10)). Then we obtain from system (5.1): Taking now the determinant of the system equal to zero we have the characteristic relation Since we are looking for the minimum positive speed with respect to s > 0 and s 3 > 0, in Fig. 12(a) we plot implicitly equation (5.8) in the (s, s 3 , w)-space. Note that the steepness coefficient s 3 (for the ECM degradation profile) does not have any significant effect on the minimum positive wave speed w (of the invading u 1 and u 2 populations). For this reason, in Fig. 12(b) we plot the relation (5.8) in the (s, w)-plane for fixed s 3 > 0. We see that indeed we cannot have a travelling wave with positive speed unless the wave has a steepness s > 3.17. Moreover, faster waves have higher steepness.
Finally, we are interested in investigating how the speed of the invading waves correlates with their steepness, as we vary different model parameters. In Fig. 13(a) we see that a decrease in the diffusion coefficient D by one order of magnitude leads to a reduction in the velocity w > 0. In Fig. 13(b) we see that an increase in the mutation rate M (from 0.001 to 0.05) leads to a slightly lower velocity w > 0 and Fig. 11. Plot of the numerical simulations of the travelling wave profile, shown in Figs 8 and 9, for (a) population u 1 and the exponential ansatz function E 1 (x) given by relation (5.3) with A 1 = 0.063 · exp (24.5675); (b) population u 2 and the exponential ansatz function E 2 (x) given by relation (5.4) with A 2 = 0.8 · exp (24.2505); (c) ECM density, f , and the exponential ansatz function E 3 (x) given by relation (5.5) with A 3 = 0.82 · exp (−52.5); (d) integrin density, c, and the exponential ansatz function E 4 (x) given by relation (5.6) with A 4 = 0.4 · exp (24.092). Here, q r = 0.0005, q a = 0.00025 and the rest of the model parameters are given in Table 2. a higher steepness. (Note that for s > 100, the invading speed obtained for M = 0.001 and M = 0.05 is the same-not shown here.)

Conclusions and discussion
In this paper we introduced a model of integro-differential equations describing the dynamics of two cancer cell populations: a 'normal' cell population that exhibits both random and directed movement, and a 'mutant' cell population that exhibits only directed movement. The model incorporated both nonlocal cell-cell interactions and cell-matrix interactions. Moreover, unlike other models in the literature, in our model these interactions were not constant but depended on the cellular level of integrins.
Linear stability analysis of the non-dimensional model showed that aggregations could arise only via real bifurcations (the system could not exhibit Hopf bifurcations). Numerical results showed that these aggregations were described by a large number of stationary pulses (e.g. 13 pulses corresponding to the 28 V. BITSOUNI ET AL. critical unstable wave number k 13 ). Moreover, numerics also showed the existence of travelling waves. We investigated the speed of these waves, which seemed to be affected by the diffusion coefficient and the mutation rate of cells.
The rate at which cancer cells mutate seemed to play a critical role in our model. In Figs 6 and 7 (as well as Figs 8 and 10) we showed that depending on the magnitude of the mutation rate, either the u 1 or the u 2 cell populations can be eliminated (or, in some cases they can co-exist-see Figs 4 and 8). The existence of these dominant behaviours exhibited by the u 1 or u 2 populations are consistent with the principle of competitive exclusion of clonal sub-populations in heterogeneous tumours Keats et al., 2012;Leith et al., 1989). In Fisher et al. (2013) the authors interpreted the experimental data, which showed suppression and reappearance of cancer clones in myeloma patients (Keats et al., 2012) and chronic lymphocytic leukaemia patients (Schuh et al., 2012), by suggesting that two subclones can exist in a dynamic equilibrium. While all these experimental studies record the survival/suppression of tumour clones in various cancers, they do not offer a mechanistic explanation for the factors that could lead to these behaviours. In contrast, our numerical results (see Figs 6 and 7, 8 and 10) offer such a mechanistic explanation by identifying the magnitude of the mutation rate as a factor that could explain the experimentally observed suppression and reappearance of cancer clones.
This issue of cancer heterogeneity has significant implications on cancer drug therapies, since it can lead to drug resistance. For example, clinical studies have shown the emergence of imatinib-resistant mutation clones in patients with chronic myeloid leukaemia, which can co-exist with subclones that carry different imatinib-resistant mutations in treatment-naïve patients (Shah et al., 2002). Overall, these clinical observations suggest that a single drug might not be successful in treating a genetically heterogeneous tumour, since sub-populations of cancer cells with drug-resistant mutations could become dominant, thus leading to therapeutic failure.
We note that the spatial distribution of the two cancer sub-populations, with the original u 1 population in the centre of the aggregation and the mutated u 2 clones towards the edges of the aggregation (see Fig. 5 for t ∈ [0, 10] or Fig. 6 for t ∈ [100, 300]), is consistent with experimental studies on the spatial relationship between clonal sub-populations of hepatocellular carcinoma (HCC) tumour (Ling et al., 2015). In this experimental study, the authors investigated the clonal diversity of a HCC-15 tumour (and the spatial distribution of these clones), and showed the ancestral clones being positioned in the middle of the tumour, with the descendant clones radiating outward.
The numerical results of the model presented in this paper show that cell-cell/cell-matrix adhesion combined with cell proliferation (in the presence of cell competition) and cell mutation, can impact which tumour clones survive (see Figs 4-8 and 10). Therefore, our model puts forward the idea that in order to predict the emergence of a particular cancer clone, we should first quantify cell-cell and cell-matrix adhesive interactions for the various tumour clones, and their proliferative capacities.
In this paper we focused on a 1D model. However, real life cell dynamics occurs in 2D or 3D. In Appendix A we extended model (2.22) to two spatial dimensions. Moreover, to verify if the linear stability results for the 1D model (2.22) were valid also in 2D, we applied linear stability analysis. We showed that for similar kernels, we obtained similar dispersion relations with no imaginary parts (and hence no Hopf bifurcations). Therefore, we expect that the 2D model (A1) would exhibit stationary pulses similar to the ones exhibited by the 1D model (2.22). Future work will consider extending the numerical results in 2D. The nonlocal terms F 1 [ u, f , c] and F 2 [ u, f , c] are described by equations (2.10) and (2.11), respectively. The microscale evolution of the internal state variable c is described by differential equation (2.16): The incorporation of this internal cell structure c into the mesoscopic model (B1) follows the same approach as in Erban & Othmer (2004 One can eliminate the derivative ∂/∂c from equations (B1a)-(B1b) by integrating these equations over c, and obtaining new macroscopic equations defined in terms of the averaged tumour cell variables U 1,2 (x, t) = u 1,2 (x, t, c) dc (Erban & Othmer, 2005;Engwer et al., 2015;Domschke et al., 2107). Other studies obtained different macroscopic limits of the mesoscopic structured-population models by considering the assumption that the internal (structure) variable evolves on a much faster time scale that is described by a small parameter ε (Perthame et al., 2016). If we would take the same approach and assume that the integrin dynamics is much faster compared to the dynamics of cancer cells and ECM (note that in Table 2 the production rates of integrins by u 1 and u 2 cells are estimated, due to a lack of available data), then the time-evolution of integrins would be given by ε ∂c ∂t = p 1 u 1 + p 2 u 2 − qc, ( B . 3 ) which, to leading order, gives rise to the following steady-state equation: c = (1/q)(p 1 u 1 (x) + p 2 u 2 (x)). If we substitute this state into (B1) we recover model (2.22), with the exception of equation (2.22d) which is now a steady-state equation. We note here that while different 'microscale-macroscale' models (Stinner et al., 2014(Stinner et al., , 2015(Stinner et al., , 2016Meral et al., 2015a,b) and 'mesoscale' models (Lorenz & Surulescu, 2014;Engwer et al., 2015;Perthame et al., 2016;Engwer et al., 2017;Hunt & Surulescu, 2017) for tumour-ECM interactions via integrins have been developed over the last two decades, to our knowledge there are not many studies that compare and contrast the outcomes of different ways the integrins are incorporated into these different models.
A comparison between the effects of integrins on the dynamics of the structured-population model (B1) and the original model (2.22) (both assuming heterogeneous cell populations), including an investigation into the possibility of having travelling wave solutions for the structured model (B1), will be the subject of a future study. (Such an investigation is beyond the scope of this study due to its complexity, which would increase considerably the length of this paper.)