Non-local Parabolic and Hyperbolic Models for Cell Polarisation in Heterogeneous Cancer Cell Populations

Tumours consist of heterogeneous populations of cells. The sub-populations can have different features, including cell motility, proliferation and metastatic potential. The interactions between clonal sub-populations are complex, from stable coexistence to dominant behaviours. The cell–cell interactions, i.e. attraction, repulsion and alignment, processes critical in cancer invasion and metastasis, can be influenced by the mutation of cancer cells. In this study, we develop a mathematical model describing cancer cell invasion and movement for two polarised cancer cell populations with different levels of mutation. We consider a system of non-local hyperbolic equations that incorporate cell–cell interactions in the speed and the turning behaviour of cancer cells, and take a formal parabolic limit to transform this model into a non-local parabolic model. We then investigate the possibility of aggregations to form, and perform numerical simulations for both hyperbolic and parabolic models, comparing the patterns obtained for these models.


Introduction
Collective cell movement can be observed in many types of cells and plays an important role in many physiological processes, including wound healing, embryonic development and metastasis of cancer cells (Friedl and Wolf 2003;Rørth 2009). Cancer cells' movement and aggregations are influenced by external factors (e.g. concentration of nutrients), as well as internal factors, such as the social forces among cells (i.e. attrac-B Vasiliki Bitsouni vbitsouni@gmail.com; vbitsouni@maths.dundee.ac.uk 1 tion, repulsion and polarisation). These internal factors lead to self-organised cell aggregation and the formation of a wide variety of patterns, playing a crucial role in cell movement. Experimental studies (Omelchenko et al. 2003;Rørth 2012) have shown that alignment (polarisation) has been reported as the initial cellular response in wound healing and cancer invasion. By alignment, we mean the process where cells turn to adapt their orientation to that of their neighbours, which leads to a polarised group of cells having the same orientation in space and travelling large distances together. In contrast, there are non-polarised groups in which all cells move individually, while the group as a whole can remain stationary or drift slowly (Lutscher 2002;Firtel and Meili 2000). Cells interact with their neighbours and change their shape and direction of movement as a result of this collective movement and a process known as contact inhibition of locomotion (CIL) (Vicente-Manzanares and Sánchez-Madrid 2000), which plays a crucial role in cancer invasion and metastasis. During this process, cells alter their direction of movement when contact other cells in order to avoid collision.
Although the exact mechanism that makes the cells cooperate with each other and migrate collectively in one direction is not fully clear, it has been observed that "leader" cells at the front of outgrowths (e.g. an epithelial cell at the edge of an epithelial sheet that adopts a fibroblast-like morphology extending a wide lamellipodium) are accompanied by many "follower" cells along the sides, both migrating to distant sites (Haga et al. 2005;Omelchenko et al. 2003). As cells move in a collective manner, only the cells in the free edge will produce lamellipodia, while cells inside the group will form smaller protrusions or no protrusion. Although CIL process will lead to a change in the direction of movement of the cells in the edges, the whole group of cells will follow this movement as a result of cell-cell interactions, ending up in the realignment of the cell populations (Mayor and Carmona-Fontaine 2010). Moreover, recent experimental results showed that the polarisation and migration of cells within an epithelial monolayer are coordinated over spatial distances greater than ten cell diameters (Angelini et al. 2010;Petitjean et al. 2010;Das et al. 2015).
In the mathematical literature, there are various models that consider the effect of non-local social interactions on the collective movement of cells and animals. A large number of models for the collective movement of animals consider the interplay between all three social interactions: repulsion, attraction and alignment (Canizo et al. 2010;Cavagna et al. 2010;Gautrais et al. 2012;Huth and Wissel 1992;Kunz and Hemelrijk 2003;Lukeman et al. 2010). Some of these models consider non-local turning rates and constant speeds (see, e.g. Buono and Eftimie 2015;Fetecau 2011). Other models investigate the effect of social interactions also on animals speed (Fetecau and Eftimie 2010;Topaz et al. 2006). In regard to the models for the collective movement of cells, the majority of these models focus on attractive-repulsive interactions (Armstrong et al. 2006;Bitsouni et al. 2017Bitsouni et al. , 2018Domschke et al. 2014;Painter et al. 2015;Sherratt et al. 2009). Very few non-local models incorporate cell alignment [see, for instance, Mogilner and Edelstein-Keshet (1995)]. Therefore, it is very important to develop non-local models that consider cell polarisation and describe the way that all three social forces affect the velocity and the turning behaviour of cells.
In this paper, we introduce a new model describing the interplay between cell polarisation and cell repulsive-attractive interactions. In contrast to the models mentioned in the previous paragraph, here we consider both non-local speed and turning rates. To this end, we derive a model of nonlinear non-local first-order hyperbolic equations describing the dynamics of polarised early-and late-stage cancer cell populations. In addition to cell movement and cell turning behaviours (which depend on repulsive, attractive and polarising forces), we also consider mutation and proliferation. We investigate numerically the patterns generated by this hyperbolic model-by focusing on the effect of the following parameters: (i) magnitude of repulsive/attractive/polarisation (alignment) interactions; (ii) turning rates; (iii) proliferation rates; and (iv) baseline speed. Since the majority of papers describing collective movement of cells are of parabolic type, in this paper we also take a parabolic limit to investigate the preservation of patterns in this limit.
To keep the model as simple as possible, we focus only on the effect of attraction/repulsion/alignment on cell-cell interactions [while ignoring the interactions between cells and the extracellular matrix (ECM)]. This is consistent with other mathematical approaches in the literature of collective movement of cells [see, e.g. the individual-based models in Arboleda-Estudillo et al. (2010), Chang et al. (2013), Hirashima et al. (2013), Méhes and Vicsek (2014) and Woods et al. (2014)]. We emphasise here that in contrast to these studies that develop discrete models (which are difficult to be investigated analytically), here we present a continuous model where we apply analytical and computational techniques to understand pattern formation in cellular aggregations. This paper is organised as follows. In Sect. 2, we present a model of non-local nonlinear hyperbolic equations describing the dynamics of two sub-populations of polarised cancer cells, with different levels of mutation. In Sect. 3, we derive the parabolic limit of this hyperbolic model. In Sect. 4, we perform linear stability analysis of both hyperbolic and limiting parabolic models to investigate the ability of these models to form cell aggregation. In Sect. 5, we investigate numerically the spatiotemporal patterns obtained by the hyperbolic model and compare the results with the patterns obtained by the limiting parabolic model. We conclude in Sect. 6 with a discussion of the results.

A Non-local Hyperbolic Model for Cancer Cell Polarisation
In this section, we introduce a new non-local model that incorporates the tendency of cancer cells to align with other cells that are within a range (alignment range). The model describes the movement of two cancer cell populations, an early-and a latestage population. Here, we assume that the movement of cancer cells is governed by directed motility in response to cell-cell interactions, choosing to ignore the cell-ECM interactions.
Let Ω ⊂ R be a bounded interval. Let I T = [0, ∞) be the time interval. We denote by u + 1 (t, x) (u − 1 (t, x)) the density of early-stage cancer cells at (t, x) that move to the right (left), and, respectively, by u + 2 (t, x) (u − 2 (t, x)) the density of late-stage cancer cells at (t, x) that move to the right (left). The total cancer cell population density is given by the relation u 1 = u + 1 + u − 1 for the early-stage cancer cell population and, respectively, by u 2 = u + 2 + u − 2 for the late-stage cancer cell population. For compact notation, we define the vector u (t, x) = (u 1 (t, x) , u 2 (t, x)) T . We also define the cell population flows by υ i = u + i − u − i , i = 1, 2, for early-stage (i = 1) and late-stage (i = 2) cancer cells. Thus, we derive the following hyperbolic system of conservation laws that describe the evolution of densities of left-moving and right-moving earlyand late-stage cancer cells: where Γ ± u are the density-dependent speeds and λ + u i (λ − u i ) are the density-dependent turning rates for the cancer cells initially moving to the right (left) which then turn to the left (right). We denote by M the mutation rate of cancer cells and by r i , i = 1, 2, the proliferation rate of population u i . Note that we consider a non-dimensionalised model, where the cancer cell densities u i , i = 1, 2, are non-dimensionalised by the carrying capacity for the cells, k u , leading to logistic growth functions with unit-valued carrying capacity for the cells. A non-dimensionalisation of the model and two tables with the model variables and parameters are presented in "Appendix A".
The turning rates are functions of the cell-cell interactions, y ± u + 1 , u − 1 , u + 2 , u − 2 , described as in Eftimie et al. (2007): where the constants λ r i and λ b i , i = 1, 2, represent a baseline random turning rate and a biased turning rate, respectively. The dimensionless functionals y ± u + 1 , u − 1 , u + 2 , u − 2 (see "Appendix A.1") of the densities of right-moving, u + i , and left-moving, u − i , cancer cells incorporate non-local interactions between the two sub-populations of polarised cells and can be described by the following relation: where y ± j u + 1 , u − 1 , u + 2 , u − 2 , j = a, r , al, denote the attraction, repulsion and alignment functionals, respectively, which influence the likelihood of a cancer cell to turn to the left (+) or to the right (-). We note here that stronger interaction forces lead to higher turning rates. Let us define R s > 0 to be the cells sensing radius, 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). The attraction and repulsion interactions are described by the following non-local terms (Buono and Eftimie 2015;Colombi et al. 2015Colombi et al. , 2017: with q a and q r describing the magnitudes of attractive and repulsive interactions, respectively, and K a (x) and K r (x) describe the spatial ranges over which these interactions take place. We denote by K (x) := q a K a (x)−q r K r (x) the attraction-repulsion kernel, assuming that it is attractive at medium/long ranges (i.e. at the edges of the cell, and over the neighbouring cells) and repulsive at very short ranges (i.e. over the cell surface). The non-local alignment term is given by the relation (Buono and Eftimie 2015): with q al describing the magnitude of alignment and K al (x) describing the spatial range over which alignment takes place. Let us now focus on the density-dependent speeds Γ ± u . Here, we choose the non-local speeds of the two-population cancer cells to be described by non-negative, bounded and increasing functionals of the attractive-repulsive cell-cell interactions. Thus Γ ± u are given by the following relations where γ is a constant baseline speed describing the behaviour of the cancer cell populations in the absence of cell-cell interactions (see Fetecau and Eftimie 2010). We denote by g u := tanh y + u . Since function tanh (·) is an odd function, then for q al = 0 relation (6) becomes: For the non-local terms, we choose translated Gaussian kernels with s j representing half the length of the interaction ranges and m j = s j /8 representing the widths of the interaction kernels [the constants m j , j = a, r , al, 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)]. Note that other studies used discontinuous Morse-type repulsion-attraction kernels (Fetecau and Eftimie 2010), which have a more realistic shape (with highest repulsion at x = 0), but which can cause density blow-up [a different class of repulsion-attraction kernels in higher dimensions, which are also discontinuous at the origin where they have the highest density, but which are always positive (in contrast to the more classical Morse kernels that can be positive and/or negative depending on parameter values), was recently discussed by Carrillo et al. (2016)]. To avoid this type of unrealistic aggregation behaviour, we have chosen translated Gaussian kernels (8).
We study the hyperbolic model (1) on a finite domain of length L, that is, x ∈ [0, L], with wrap-around boundary conditions for the non-local social interactions. Thus, we have a problem with a discrete spectrum, where for L large we can approximate the process of pattern formation on an unbounded domain. To complete the model, we have to impose boundary conditions. Note that since system (1) is hyperbolic, we have to follow the characteristics of the system when imposing these boundary conditions. For this reason, u + i , i = 1, 2, are prescribed only at x = 0, while u − i , i = 1, 2 are prescribed only at x = L. For this model, we choose periodic boundary conditions, where the cancer cells move on a circular domain, leaving the domain at one end and entering it again at the other end. The boundary conditions are described by:

Parabolic Limit for Non-local Interactions
In this section, we take a formal parabolic limit to investigate the connection between the hyperbolic model (1) and other published non-local parabolic models for collective cell dynamics, which have density-dependent speed (see, e.g. Domschke et al. 2014;Painter et al. 2015 and the references therein). To study the parabolic limit of our hyperbolic model, we assume that there is no alignment, i.e. q al = 0. The main reason for choosing to ignore alignment is that we aim to obtain closed-form parabolic equations for the total densities of cancer cells (u 1,2 ), and the alignment terms (for q al = 0) incorporate left-and right-moving cancer cells (which would lead to terms depending on the flow v 1,2 ).
Recent experimental results on cancer cell motility have shown that there are differences in terms of speed and directional persistence between 2D and 3D cell migrations (Wu et al. 2014). Moreover, while cells have been observed to move persistently in 2D at short time scales (i.e. their behaviour corresponding to hyperbolic dynamics), they displayed correlated random walk at long time scales (i.e. their behaviour corresponding to parabolic dynamics). We can explain the experimental behaviour observed on short and long time scales by taking the parabolic limit of the hyperbolic model (1). As discussed in Hillen and Painter (2013), there are two approaches for this parabolic limit: (i) an appropriate scaling of space and time and (ii) large turning rates and large speeds. Since the two approaches are equivalent (Hillen and Painter 2013), throughout this study we chose to focus on the scaling of the turning rates and speeds.
To transform model (1) into a parabolic model, we follow the classical approach in (Kac 1974;Hillen and Stevens 2000) and differentiate with respect to t and x the sum and difference of Eqs. (1a)-(1b) and also Eqs. (1c)-(1d). After eliminating the equations for the cell fluxes ( , we are left with two equations for the total densities of cancer cells (u 1,2 = u + 1,2 + u − 1,2 )-see Eqs. (51) and (52) in "Appendix B" (as well as the details of the calculations shown in "Appendix B").
Next, we rescale the turning rates and speeds, by assuming that the cancer cells move very fast and change direction even faster [with respect to other normal cells in the tissue; see also the experimental study in Tzvetkova-Chevolleau et al. (2008)]. Moreover, since it makes sense to assume that high speeds and high turning rates lead also to a reduced sensitivity to the environment (and implicitly to other neighbouring cells), we consider also a scaling of the density-dependent components of the speed [i.e. the term g [u] that appears in Γ ± [u] = γ (1 ± g [u])] and of the turning rates (i.e. the term p(y ± [u]) that appears in λ ± u 1,2 ) [see also the approach in Buono and Eftimie (2015)].
With these assumptions, we introduce a small dimensionless parameter > 0 and use it for the following rescaling We denote by This reduction in the sensitivity in the environment leads to the following rescaling We then substitute these rescaled parameters and functions into a reduced system for the total cell densities u 1,2 = u + 1,2 + u − 1,2 , which was obtained from model (1); see Eqs. (51)-(52) in "Appendix B". This approach has been previously described in detail in Hillen and Levine (2003) and Hillen and Stevens (2000). As we then take the limit → 0 in this simplified system, we obtain the following parabolic equations , the growth functions. Here, the initial conditions are given by the functions To fully define the parabolic model (12), we need to impose boundary conditions. To be consistent with the hyperbolic model (1), we impose again periodic boundary conditions on a finite domain of length L: We note that, since we assumed q al = 0, the non-local terms f u and g u now depend only on the repulsive and attractive interactions.

Linear Stability Analysis
In this section, we investigate the possibility of pattern formation for models (1) and (12) via linear stability analysis. To this end, we focus on model parameters, including the magnitudes of social forces (i.e. attraction, repulsion, alignment) between cancer cells, and their role on pattern formation.

Linear Stability Analysis of the Hyperbolic Model
We start with the linear stability analysis of the hyperbolic model (1). First, we look for the spatially homogeneous steady states u ±, * i , i = 1, 2, assuming that cancer cells are spread evenly over the domain. We denote the constant total density of the populations by u * i , i = 1, 2. From the right-hand side of Eqs. (1a)-(1d), we have the following system: which has the solutions u * 1 , u * 2 = (0, 0) and u * 1 , u * 2 = (0, 1). Note that, for biological realism, we consider only non-negative solutions. If we consider the states where both cell populations are evenly spread in both directions over the domain, then these states u +, * 1 , u −, * 1 , u +, * 2 , u −, * 2 are given by (0, 0, 0, 0) and (0, 0, 0.5, 0.5) .
If we consider populations that are evenly spread over the domain, but where more individuals are facing one direction compared to the other direction (i.e. u +, * 2 = u −, * 2 ), then the steady states are given by for 0 ≤ u +, * 2 ≤ 1. Now that we know the steady states, we proceed with the study of the local stability of these solutions under small perturbations caused by spatially non-homogeneous terms. We let u where k and λ are the wave number and frequency, respectively. Due to the finite domain (with wrap-around boundary conditions), we have that the wave number, k, takes only discrete values k j = 2 pj/L, j = 1, 2, 3, . . . . Let K ± j , j = a, r , al, be the Fourier transform of the interaction kernel K j , given by the following relation We denote byK s (k) =K + (k) −K − (k) = q aK s a (k) − q rK s r (k) the Fourier sine transform of kernel K , and, respectively, byK c al (k) =K + al (k) +K − al (k) the Fourier cosine transform of kernel K al . Throughout this study, we will consider translated Gaussian kernels given by relation (8). Then, the Fourier transform of these kernels is given byK and the Fourier sine and cosine transforms are given bŷ To simplify the results of this section (using the fact that u +, * 1 = u −, * 1 = 0), we set the following parameter values: Substituting now the expressions u ± j = u ±, * j + A ± u j e ikx+λt , j = 1, 2, into the system (1) and using the above relations, we obtain the following dispersion relations: -For the steady state (0, 0, 0, 0), we have: and -For the steady state 0, 0, u +, * 2 , 1 − u +, * 2 , we have: and as and increase λ λ  (21) for D 1 , E 1 (blue) and D 2 , E 2 (red); b the effect of γ on the graph of Re (λ 2 (k)); c the effect of r 2 on the graph of Re (λ 2 (k)); d the effect of λ r 2 on the graph of Re (λ 2 (k)). The continuous curves represent the Re(λ(k)), while the dotted curves represent the I m(λ(k)). The model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers k j = 2π j/L, j = 1, 2, . . . (Color figure online) Equations (21) and (26) show that the steady states are unstable, i.e. Re (λ (k)) > 0, when D l (k) < 0 or E l (k) < 0, l = 1, . . . , 4. Examples of such dispersion relations are shown in Figs. 1a and 2a. There is a range of k-values for which Re (λ (k)) is positive, and thus, aggregation can arise from spatial perturbations of the steady states (0, 0, 0, 0) (see Fig. 1a) and (0, 0, 0.5, 0.5) (see Fig. 2a). Note that similar results (not shown here) are obtained for any steady state 0, 0, u +, * 2 , 1 − u +, * 2 , with 0 ≤ u +, * 2 ≤ 1.
The effect of the parameters on the stability of the steady states We now use the dispersion relations (21) and (26) to study the effect of the key parameters on pattern formation. We investigate the stability of the spatially homogeneous steady states (0, 0, 0, 0) and (0, 0, 0.5, 0.5) by increasing (or decreasing) the parameters connected to the dispersion relations. Precisely, we show the effect of the parameters on the graph of the eigenvalue with the maximum real part, i.e. Re (λ 2 (k)) and Re (λ 4 (k)) of the dispersion relations (21) and (26), respectively. However, we note that the effect of  (26) for D 3 , E 3 (blue) and D 4 , E 4 (red); b the effect of γ on the graph of Re (λ 4 (k)); c the effect of r 2 on the graph of Re (λ 4 (k)); d the effect of λ r 2 on the graph of Re (λ 4 (k)); e the effect of q a on the graph of Re (λ 4 (k)); f the effect of q r on the graph of Re (λ 4 (k)). The continuous curves represent the Re(λ(k)), while the dotted curves represent the I m(λ(k)). The model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers k j = 2π j/L, j = 1, 2, . . . (Color figure online) the parameters (e.g. λ r 1 and r 1 , and the rest parameters) is similar on the graphs of Re (λ 1 (k)) and Re (λ 3 (k)).
For the tumour-free steady state (0, 0, 0, 0), we can see from the dispersion relation (21) that its stability does not depend on the magnitudes of cell interactions, although these magnitudes are crucial in the nonlinear interactions that control cell aggregation patterns. In Fig. 1b, c, d, we see that the stability of this steady state depends mainly on the baseline speed, the proliferation rates and the baseline random turning rates, with the latter not having a significant impact on the instability of the steady state. Precisely, for zero proliferation rates (r 1 = r 2 = 0) the steady state is stable (see Fig.  1c), while an increase in them or in the baseline random turning rates (λ r 1,2 ) results in a shift to the right of the wave number that will emerge. Note also that an increase in the baseline speed (γ ) has an inverse result in the dispersion relation, with a shift to the left of the wave number that will emerge.
For the steady state (0, 0, 0.5, 0.5), there are more parameters that affect the stability properties and as it can be seen by relations (29) and (30) the stability of this steady state depends on the magnitudes of cell interactions as well. We can see in Fig. 2b, c that the baseline speed and the proliferation rates have an opposite effect on the stability changes of (0, 0, 0.5, 0.5), compared to that on the cancer-free steady state. In Fig. 2e, f we see that a decrease in the magnitude of attraction leads to the change in the stability of this steady state, while an increase in the magnitude of repulsion leads to the shift of the critical wave numbers to the right. Note that although the magnitude of alignment does not seem to affect the stability wave number, it is crucial though in the pattern formation, as we will see in the following section.

Remark 1
We should mention here that the effect of the mutation rate, M, on the dispersion relation is not significant. Precisely, as M appears only on the functions D 1 , E 1 and D 3 , E 3 , any changes in the values of M will not lead to stability change, but only to a reduction on the eigenvalues λ 1 and λ 3 (up to below zero) as M increases. Note that we always refer to the greater eigenvalues of Eqs. (21) and (26) To investigate the effect of the scaling parameter on the dispersion relation, in Fig. 3a, b we show the stability of the steady states (0, 0, 0, 0) and (0, 0, 0.5, 0.5), respectively, after applying the rescaling given in Sect. 3 and taking q al = 0. Although the wave numbers k j = 2π j/L, j = 1, 2, . . . are discrete (represented by diamondshaped points in x-axis of Figs. 1 and 2), here they are plotted as a continuous axis to show clearly the effect of on the imaginary part of the dispersion relation described by dotted curves. We see in Fig. 3b that for = 1 there are wave numbers for which we can have Re λ k j > 0 and I m λ k j > 0, giving rise to travelling patterns. As decreases (e.g. = 0.5), we note that I m λ k j = 0 at some wave numbers where Re λ k j > 0, and thus stationary pulses are expected to be obtained. As → 0, the imaginary part of the eigenvalues will be always zero, and numerically we expect to observe stationary patterns.

Linear Stability Analysis of the Parabolic Model
Next, we investigate the conditions under which aggregations can arise for the limiting parabolic model (12). We first calculate the spatially homogeneous steady states of the parabolic model. We see from Eqs. (12a)-(12b) that the ODE model associated with system (12) is described by system (14), which has the solutions u * 1 , u * 2 = (0, 0) and (0, 1) (considering again only non-negative solutions).  Table 2. For graphical purposes, the discrete wave numbers has been omitted (Color figure online) Proceeding with the linear stability analysis of the spatial system (12), we apply small spatial perturbations to the homogeneous steady states: u 1 = u * 1 + A u 1 e ikx+λt and u 2 = u * 2 + A u 2 e ikx+λt with |A u 1 |, |A u 2 | 1. Substituting these terms into system (12), using the parameter values (20) and replacing u * 1 = 0, yields the following dispersion relation: Therefore, for the steady state (0, 0) we have the solutions: Fig. 4 Plot of the eigenvalues obtained by dispersion relation (31). a λ 1 (k) = −k 2 D u 1 − M + r 1 (blue) and λ 2 (k) = −k 2 D u 2 + r 2 (red) for the steady state (0, 0)

(b) (a)
, for the steady state (0, 1). The model parameters are given in Table 2  The dispersion relation (31) for the steady state (0, 1). a The effect of λ r 2 on the graph of Re (λ 2 (k)); b the effect of q r on the graph of Re (λ 2 (k)); the rest of the model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers k j = 2π j/L, j = 1, 2, . . . (Color figure online) and for the steady state (0, 1), the solutions: As in the case of the hyperbolic model, we see in Fig. 4 that there is a range of k-values for which Re (λ (k)) > 0, and thus, aggregation can arise from spatial perturbations of the steady states (0, 0) and (0, 1) of the parabolic model. However, taking into account that we have q al = 0 and all eigenvalues have zero imaginary part, we expect that the numerical simulations will show stationary patterns. Figures 1a and 4a, as well as Figs. 2a and 4b, are similar, and the effect of the key parameters on the dispersion relation is the same. In Fig. 5, we only include the parameters which seem to have a slightly different effect on the dispersion relation for the steady state (0, 1) of the parabolic model, compared to the one of the hyperbolic model for the steady state (0, 0, 0.5, 0.5). Although in Fig. 2d the effect of λ r 2 on the stability of the steady state is not very clear, in Fig. 5a we can see that a decrease in λ r 2 leads to a shift of the graph below and a change in stability for very low λ r 2 . Moreover, when we increase q r (see q r = 1.1 in Fig. 5b), the steady state becomes stable and we do not expect to have aggregation.
Remark 2 Note that k = 0 is unstable for the trivial steady states [see Figs. 1 and 4a; a similar behaviour being observed in other non-local hyperbolic models ] and stable for the non-trivial states (see Figs. 2 and 4b; this behaviour is similar to the classical Turing-type instability).

Numerical Results
To understand the behaviour of systems (1) and (12), we investigate them numerically. The aim of this section is to study the effect of the cell-cell interactions, baseline speed, proliferation and turning rates on the pattern formation for both models. The numerical results presented in this section are based on the investigation of the parameter sets that lead to pattern formation, i.e. predicted unstable wave numbers, as shown in the context of the linear stability analysis presented in the previous section. Note that similar patterns for different parameter sets or different initial conditions are not shown here for a better flow of this paper.
We use a time-splitting approach to discretise our model. We discretise the spacetime plane choosing a time step t = 0.001 and a space step x = 0.01. We use a Crank-Nicolson scheme to propagate the solution of the diffusion terms for the parabolic equations (12), obtained with the formal parabolic limit of Eqs. (1). For the time propagation of the advection terms in both models (1) and (12), we use the Nessyahu-Tadmor scheme (Nessyahu and Tadmor 1990). Finally, for the time propagation of the reaction terms in (1) and (12), we use a fourth-order Runge-Kutta algorithm, where the integrals are further discretised using the Simpson's rule. All simulations are performed on a domain of length L = 10 with periodic boundary conditions (introduced to approximate the dynamics on an infinite domain). To deal with the integrals at the boundaries of the domain, we wrap them around the domain. The simulations ran for times up to t = 2000, but we show the dynamics for time that the patterns are more clear. The parameters used in the numerical simulations are listed in Table 2 in "Appendix A".

Pattern Formation for the Non-local Hyperbolic Model
Let us focus first on the numerical simulations for the non-local hyperbolic model (1). 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 To begin we first run numerical simulations for small random perturbations of the steady states (0, 0, 0, 0) and (0, 0, 0.5, 0.5). As stated in the context of linear stability analysis in the previous section, the mutation rate, M, does not have any significant effect on the stability of the steady states. The numerical simulations obtained for M = 0.05 show similar patterns (not shown here) with those for M = 0.0002, but with u 1 population vanishing much faster and population u 2 reaching greater density values.
The effect of proliferation rate on cancer cell movement and aggregation As one of the key parameters of the stability of the steady state (0, 0, 0, 0) is the proliferation rate, in Fig. 6 we see that as the proliferation rates increase, e.g. r 1 = 0.3 and r 2 = 0.4 (see Fig. 6a', b') the u 1 and u 2 exhibit larger number of smaller rotating waves, as it was expected from the increased number of critical wave number shown in linear stability analysis (see Fig. 1c). Increasing further the proliferation rates r 1,2 (e.g. up to 0.7) leads to more aggregation waves. Choosing now the same proliferation rate for the two populations, e.g. r 1 = r 2 = 0.1, we see in Fig. 6a", b" an interesting effect of clonal competition, with population u 1 dominating the dynamics. This could be explained by the very low mutation rate and an equal competitive effect (i.e. −r 1 u 1 u 2 = −r 2 u 1 u 2 ). Note that similar behaviours (where u 1 dominates the dynamics) have been observed for other sets of parameters with equal cell proliferation rates: e.g. r 1 = r 2 = 0.6 or 0.7 (not shown here).
The effect of turning rates on cancer cell movement and aggregation If we increase the baseline random turning rates, from λ r i = 0.2 to λ r i = 0.4, i = 1, 2, we see in Fig. 7 that the cancer cell populations change their movement from rotating waves (panels a, b) to a combination of rotating waves and stationary pulses that are connected through splitting/merging behaviours that result from high random cell turning rates (panels a', b'). However, if we reduce the baseline random turning rates, from λ r 1,2 = 0.2 to λ r 1,2 = 0.1, the rotating waves persist (not shown here). This behaviour was obtained for both steady states.

The effect of baseline speed on cancer cell movement and aggregation
To investigate the effect of the baseline speed on pattern formation, we focus on the random initial conditions (34) and any of the two steady states (since perturbations of both states give rise to similar patterns). First we observe in Fig. 7a, b, a", b" that an increase in cells' speed (from γ = 0.1 to γ = 1) leads to a change in cells' movement from rotating waves (panels a, b) to standing waves (panels a", b" shown only for 495 < t < 500 to improve the clarity of the figures). Moreover, we note that a reduction in the values of the baseline speed, γ , leads to the spread of populations over the whole domain. An example of such behaviour [for γ = 0.01 and pulse-like initial conditions (35)] is shown in Fig. 10.  (34)]. a, b Total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for r 1 = 0.1 and r 2 = 0.2; a', b' total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for r 1 = 0.3 and r 2 = 0.4; a", b" total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for r 1 = r 2 = 0.1. The rest of model parameters are given in Table 2 (Color figure online) The effect of attraction and repulsion on cancer cell movement and aggregation Consider again the random initial conditions (34) applied to the steady state (u +, * 1 , u −, * 1 , u +, * 2 , u −, * 2 ) = (0, 0, 0.5, 0.5) [note that similar results have been obtained also for the steady state (0, 0, 0, 0)-not shown here]. To investigate the effect of attraction and repulsion on pattern formation, we focus on two different cases: (i) q al = 3 and q a = 6 q r = 0.1 and (ii) q al = 3 and q a = 1.2 q r = 6.5. We see in Fig. 8a, b that the increase in cell-cell attraction leads to a (relatively) small number of large stationary cell aggregations. In contrast, the increase in cell-cell repulsion leads to the  (34)]. a, b Total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for λ r 1,2 = 0.2 and γ = 0.1; a', b' total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for λ r 1,2 = 0.4 and γ = 0.1; a", b" total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for λ r 1,2 = 0.2 and γ = 1. The rest of model parameters are given in Table 2 (Color figure online) formation of a much larger number of smaller stationary cell aggregations (see Fig.  8a', b').
The effect of alignment on cancer cell movement and aggregation Next we investigate the effect of alignment on collective cell movement and aggregation. In Fig. 9, we show the numerical results when q al = 0.5 and λ r 1,2 = 0.1, for initial conditions consisting of a rectangular pulse [see (35)]. The left-moving and right-moving cancer cell populations seem to exhibit a semi-zigzag pattern, characterised by a periodic  (34)]. a, b Total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for q a = 6, q r = 0.1 and q al = 3; a', b' total density of u 1 = u + 1 + u − 1 and u 2 = u + 2 + u − 2 for q a = 1.2, q r = 6.5 and q al = 3. The rest of model parameters are given in Table 2 (Color figure online) transition between two different types of sub-patterns: stationary aggregations and travelling aggregations-see the first two columns in Fig. 9. This periodic transition seems to be similar to a heteroclinic connection. However, due to the complexity of the theory behind these heteroclinic connections in infinite-dimensional dynamical systems, it is beyond the purpose of this study to investigate them further. This will be the subject of future research.
We note here that increasing the magnitude of alignment to q al = 3 leads again to rotating waves, while decreasing it to q al = 0 leads to standing waves or stationary aggregations. In both cases, the patterns are similar to those shown in previous figures, and thus not included.
Finally, we investigate the combined effect of two parameters: cell alignment (q al ) and cell baseline speed (γ ). In contrast to the previous simulations, we now consider small speeds (i.e. γ = 0.01) and slightly larger turning rates (i.e. λ r 1,2 = 0.2). We see in Fig. 10a, b that when alignment is absent (q al = 0), the u 1 cells form stationary aggregations (which eventually vanish due to the mutation term) and the u 2 cells spread throughout the whole domain.
When alignment is present (q al > 0), we see in Fig. 10a', b' that some subpopulations of u 1 and u 2 cells move quickly to the left and to the right, reaching the domain boundaries. As in the previous case, the u 1 cells are eventually eliminated, while the u 2 cells spread throughout the whole domain. Note that similar results were obtained also with the random initial conditions (34)-not shown here.

Remark 3
To ensure that the cell aggregation patterns shown above did not depend on the boundary conditions, or on the numerical scheme used for discretisation, we also ran numerical simulations for the cases shown in Figs. 7a", b" and 8a', b' when we doubled the domain size and refined the grid mesh. In all of these cases, the results showed no significant differences.

Pattern Formation for the Limiting Parabolic Model
In this section, we run simulations for the limiting parabolic model given by (12).
As in the hyperbolic model (1), we choose the initial conditions for the cancer cell populations to be small random perturbations of the spatially homogeneous steady states (0, 0) and (0, 1) (see Sect. 4.2)  Fig. 11 The spatiotemporal patterns obtained with the hyperbolic model (1) for q al = 0 after scaling, and the parabolic model (12). a, b Standing waves obtained by (1) after scaling for = 1; a', b' stationary pulses obtained by (1) after scaling for = 0.5; a", b" stationary pulses obtained by (12) when → 0. The initial conditions for the two cancer cell populations are described by small random perturbation of the steady state (0, 0, 0.5, 0.5), for the rescaled hyperbolic model, and (0, 1), for the parabolic model [see (34) and (36)]. The rest of model parameters are given in Table 2 (Color figure online) First, we assume that q al = 0 and investigate how the travelling and stationary patterns predicted by the unstable wave numbers in the linear stability analysis of the hyperbolic model (1) (see Fig. 3), are preserved in the limit to the macroscopic parabolic model. If we consider small random perturbations of the hyperbolic steady state (0, 0, 0.5, 0.5), we see in Fig. 11a, b that for = 1 we obtain standing waves, and as decreases to 0.5, the cancer cells exhibit stationary pulses (Fig. 11a', b'), as  Table 2 (Color figure online) expected by the linear stability analysis results. Assume now → 0, and focus on the parabolic model (12). For initial conditions that are small random perturbations of the steady state (0, 1), we obtain stationary pulses (Fig. 11a", b") that are similar to the ones for the rescaled hyperbolic model when = 0.5 (see Fig. 11a', b'). Similar pattern formation results are obtained for the hyperbolic steady state (0, 0, 0, 0) (and for the state (0, 0) corresponding to the parabolic model) for the case λ r 1,2 = 0.1. If we increase the random turning rates to λ r 1,2 = 0.2, we obtain stationary pulses for every 0 < ≤ 1 in the rescaled hyperbolic model and for → 0 in the parabolic model, which was expected also from the linear stability analysis (see Fig. 3, where I m λ k j = 0 at the wave numbers where Re λ k j > 0).
Finally, since we expect that the repulsive-attractive interactions could affect also the dynamics of the parabolic model, in Fig. 12 we investigate the effect of increasing the attraction strength to q a = 6 (same as in Fig. 8a, b, but here q al = 0). We see that the parabolic model exhibits a large number of small stationary pulses, similar to those obtained by the hyperbolic model in Fig. 8.
Reducing now the speed (and assuming pulse-like initial conditions (37), as for the hyperbolic model), we notice that the patterns displayed in Fig. 13a, b are similar to those obtained for the hyperbolic model presented in Fig. 10a, b.
From the numerical simulations of both models, we conclude that the hyperbolic model can exhibit both moving and stationary behaviours, in contrast to the parabolic model that can exhibit only stationary behaviours.

Conclusion and Discussion
In this paper, we introduced a one-dimensional non-local hyperbolic model describing the interactions between heterogeneous cancer cells. We developed a model where nonlocal turning rates are included and incorporate all three social interactions: attraction, repulsion and alignment, that play a crucial role in cell movement and aggregation. We assumed that a cell changes its movement direction only after weighing the information received from left and right, speeding up and slowing down to catch up with the surrounding cells, or to avoid collisions. The mutation terms and the proliferation terms are chosen to take into account the movement of cells in opposite directions. We should emphasise that we assumed that cancer cells can detect cells that are in front and behind them.
Then, we reduced our non-local hyperbolic model to a non-local parabolic model for cellular aggregations, since parabolic models have been commonly used in the study of the formation and movement of cell and animal aggregations. To this end, we considered a formal parabolic limit which assumed very large turning rates and very large speeds.
Linear stability analysis of both hyperbolic and parabolic model was used to examine the possibility of cell aggregations to form (as a result of spatial perturbations of the steady states). The results showed that for the hyperbolic model, the dispersion relation had nonzero imaginary part (hence it was possible to have Hopf bifurcations, in addition to real bifurcations). In contrast, the dispersion relation for the parabolic model had zero imaginary part (hence aggregations could arise only via real bifurcations). Finally, we ran simulations for the hyperbolic and the related parabolic model and compared the results: the hyperbolic model was more rich in patterns, showing moving cell aggregations, while the parabolic model exhibited mainly stationary cell aggregations. These numerical results are consistent with the linear results obtained via stability analysis, emphasising the more complex behaviour of the hyperbolic model.
The model presented in this paper focussed mainly on the interactions between cells, excluding the important role of ECM on cellular adhesion, movement and aggregation. A straightforward future research direction is to include ECM density in this non-local model that incorporates cell alignment. and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A.1 Non-dimensionalisation of the Model
To obtain the non-dimensional system (1), we use the following quantities: The length scale, L 0 , is in the range of 0.1-1 cm and is defined as the maximum invasion distance of the cancer cells at the early stage of invasion (Anderson et al. 2000). The time scale is defined as τ := L 2 0 /D τ , where D τ is a reference chemical diffusion coefficient, e.g. ∼ 10 −6 cm 2 s −1 (Bray 2001). Furthermore, we rescale the cancer cells with k u . Here, k u is the carrying capacity of the cancer cell populations and it is taken to be ∼ 6.7 · 10 7 cell volume −1 (Domschke et al. 2014). For the kernels K j (s) , j = a, r , al, we define (as in Domschke et al. 2014) the dimensionless functionsK j (s) , i = 1, 2, given bỹ The non-dimensionalisation of the model is straightforward, and we obtain the following dimensionless parameters: After dropping the tildes for notation convenience, we obtain the non-dimensionalised hyperbolic model (1). Note that we assume that both populations u 1 and u 2 proliferate in a logistic manner to describe the observed slow down in tumour growth following the loss of nutrients (Laird 1964). The growth functions G 1 and G 2 for populations u 1 and u 2 , respectively, are then given by the following relations which after the non-dimensionalisation lead to the logistic growth functions with unitvalued carrying capacity for the cells that appear in model (1).

A.2 Summary of Model Variables and Parameters
We present here 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.

Parameter estimation
-The sensing radius was based on the range of values given in Armstrong et al. (2006) and Gerisch and Chaplain (2008). In this study, we chooseR s = 1. -Attraction, repulsion and alignment ranges, s j , j = a, r , al, were chosen to be smaller or equal to sensing radius, with s r < s al < s a (Green et al. 2010). -Experimental studies (Cillo et al. 1987;Hill et al. 1984;Mareel et al. 1991) have shown that the mutation rate ranges between 10 −8 / cell/generation and 10 −5 / cell/generation. We assume an average growing cell population of ≈ 10 4 cells/generation, a 1-day generation of cells (since the doubling time is about 1.2 days) (den Breems and Eftimie 2016). Thus, the mutation rate will range between M = 10 −4 /day and M = 0.1/day, which for τ = 10 5 s corresponds to a non-dimensional value of the mutation rate ranging betweenM = 0.000116 and M = 0.116 (for highly aggressive tumours). In this study, we chooseM = 0.0002. -Various experimental studies (Cunningham and You 2015;Morani et al. 2014) have shown that doubling times for tumour cells range from 1 to 10 days. This corresponds to growth rates between ln (2) /10 = 0.07 (days) −1 and ln (2) /1 = 0.7 (days) −1 . We choose τ = 10 5 s; thus,r i ∈ (0.081, 0.81) , i = 1, 2. In this study, we assume thatr 1 = 0.1 andr 2 = 0.2.   Morani et al. (2014) and the rescaled functionals f i u and h i u described in relations (10) and (11), respectively. Substituting these rescaled parameters and functions into Eqs. (51)-(52) and multiplying with 2 , taking the limit as → 0 and dropping the bars (for simplicity), leads to the following parabolic equations where D u i = (γ ) 2 2λ r i , i = 1, 2, are the diffusion coefficients.