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. attraction, 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. 2017, 2018; Domschke 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 late-stage 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 be a bounded interval. Let be the time interval. We denote by () the density of early-stage cancer cells at that move to the right (left), and, respectively, by () the density of late-stage cancer cells at that move to the right (left). The total cancer cell population density is given by the relation for the early-stage cancer cell population and, respectively, by for the late-stage cancer cell population. For compact notation, we define the vector . We also define the cell population flows by , for early-stage () and late-stage () cancer cells. Thus, we derive the following hyperbolic system of conservation laws that describe the evolution of densities of left-moving and right-moving early- and late-stage cancer cells: where are the density-dependent speeds and () 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 , the proliferation rate of population . Note that we consider a non-dimensionalised model, where the cancer cell densities , are non-dimensionalised by the carrying capacity for the cells, , 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, , described as in Eftimie et al. (2007):where the constants and , represent a baseline random turning rate and a biased turning rate, respectively. The dimensionless functionals (see “Appendix A.1”) of the densities of right-moving, , and left-moving, , cancer cells incorporate non-local interactions between the two sub-populations of polarised cells and can be described by the following relation:where , 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 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. 2015, 2017):with and describing the magnitudes of attractive and repulsive interactions, respectively, and and describe the spatial ranges over which these interactions take place. We denote by 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 describing the magnitude of alignment and describing the spatial range over which alignment takes place.

Let us now focus on the density-dependent speeds . 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 are given by the following relationswhere 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 . Since function is an odd function, then for relation (6) becomes:For the non-local terms, we choose translated Gaussian kernelswith representing half the length of the interaction ranges and representing the widths of the interaction kernels [the constants , are chosen such that the support of more than of the mass of the kernels is inside the interval (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 ), 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, , 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, , are prescribed only at , while are prescribed only at . 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. . 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 (), and the alignment terms (for ) incorporate left- and right-moving cancer cells (which would lead to terms depending on the flow ).

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 ( and ), we are left with two equations for the total densities of cancer cells ()—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 ] and of the turning rates (i.e. the term that appears in ) [see also the approach in Buono and Eftimie (2015)].

With these assumptions, we introduce a small dimensionless parameter and use it for the following rescalingWe denote by and . 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 , 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 in this simplified system, we obtain the following parabolic equations where , are the diffusion coefficients, and , 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 , the non-local terms and 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 , assuming that cancer cells are spread evenly over the domain. We denote the constant total density of the populations by . From the right-hand side of Eqs. (1a)–(1d), we have the following system: which has the solutions and . 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 are given byIf 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. ), then the steady states are given byfor .

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 and with , 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 . Let , be the Fourier transform of the interaction kernel , given by the following relationWe denote by the Fourier sine transform of kernel K, and, respectively, by the Fourier cosine transform of kernel . Throughout this study, we will consider translated Gaussian kernels given by relation (8). Then, the Fourier transform of these kernels is given byand the Fourier sine and cosine transforms are given byTo simplify the results of this section (using the fact that ), we set the following parameter values:Substituting now the expressions , into the system (1) and using the above relations, we obtain the following dispersion relations:

For the steady state , we have: with and

For the steady state , we have: with and

Table 2

A list of model parameters with their units and their non-dimensional values, which we used during numerical simulations

Param.	Description	Dimensional units	Non-dim. value (p~)	Reference
R_s	Sensing radius	length	1	Armstrong et al. (2006)
γ	Constant baseline speed of cancer cells	length/time	0.1	Estimated
q_a	Magnitude of attraction	length^2/cell	1.2	Estimated
q_r	Magnitude of repulsion	length^2/cell	0.1	Estimated
q_al	Magnitude of alignment	length^2/cell	3	Estimated
s_a	Attraction range	length	1	Estimated
s_r	Repulsion range	length	0.25	Estimated
s_al	Alignment range	length	0.5	Estimated
m_a	Width of attraction kernel	length	1/8	Estimated
m_r	Width of repulsion kernel	length	0.25/8	Estimated
m_al	Width of alignment kernel	length	0.5/8	Estimated
M	Mutation rate	1/time	0.0002	Calabrese and Shibata (2010) and Mareel et al. (1991)
λ1r	Baseline random turning rate of cancer cell population u_1	1/time	0.2	Estimated
λ2r	Baseline random turning rate of cancer cell population u_2	1/time	0.2	Estimated
λ1b	Biased turning rate of cancer cell population u_1	1/time	0.8	Estimated
λ2b	Biased turning rate of cancer cell population u_2	1/time	1	Estimated
r_1	Growth rate of u_1	1/time	0.1	Morani et al. (2014)
r_2	Growth rate of u_2	1/time	0.2	Morani et al. (2014)

Equations (21) and (26) show that the steady states are unstable, i.e. , when or . Examples of such dispersion relations are shown in Figs. 1a and 2a. There is a range of k-values for which is positive, and thus, aggregation can arise from spatial perturbations of the steady states (see Fig. 1a) and (see Fig. 2a). Note that similar results (not shown here) are obtained for any steady state , with .

Fig. 1

The dispersion relation (21) for the steady state . a Plot of the larger eigenvalues , obtained by dispersion relations (21) for (blue) and (red); b the effect of on the graph of ; c the effect of on the graph of ; d the effect of on the graph of . The continuous curves represent the , while the dotted curves represent the . The model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers (Color figure online)

Fig. 2

The dispersion relation (26) for the steady state . a Plot of the larger eigenvalues , obtained by dispersion relations (26) for (blue) and (red); b the effect of on the graph of ; c the effect of on the graph of ; d the effect of on the graph of ; e the effect of on the graph of ; f the effect of on the graph of . The continuous curves represent the , while the dotted curves represent the . The model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers (Color figure online)

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 and 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. and of the dispersion relations (21) and (26), respectively. However, we note that the effect of the parameters (e.g. and , and the rest parameters) is similar on the graphs of and .

For the tumour-free steady state , 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 () the steady state is stable (see Fig. 1c), while an increase in them or in the baseline random turning rates () 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 , 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 , 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 and , any changes in the values of M will not lead to stability change, but only to a reduction on the eigenvalues and (up to below zero) as M increases. Note that we always refer to the greater eigenvalues of Eqs. (21) and (26), given by the relation .

To investigate the effect of the scaling parameter on the dispersion relation, in Fig. 3a, b we show the stability of the steady states and , respectively, after applying the rescaling given in Sect. 3 and taking . Although the wave numbers are discrete (represented by diamond-shaped 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 there are wave numbers for which we can have and , giving rise to travelling patterns. As decreases (e.g. ), we note that at some wave numbers where , and thus stationary pulses are expected to be obtained. As , the imaginary part of the eigenvalues will be always zero, and numerically we expect to observe stationary patterns.

Fig. 3

Plot of the eigenvalue with the maximum real part, after parabolic scaling, of a relation (21) for the s.s. for ; b relation (26) for the s.s. for . The continuous curves represent the , while the dotted curves represent the . The rest of the model parameters are given in Table 2. For graphical purposes, the discrete wave numbers has been omitted (Color figure online)

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 and (considering again only non-negative solutions).

Proceeding with the linear stability analysis of the spatial system (12), we apply small spatial perturbations to the homogeneous steady states: and with . Substituting these terms into system (12), using the parameter values (20) and replacing , yields the following dispersion relation:

Therefore, for the steady state we have the solutions:and for the steady state , 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 , and thus, aggregation can arise from spatial perturbations of the steady states and of the parabolic model. However, taking into account that we have and all eigenvalues have zero imaginary part, we expect that the numerical simulations will show stationary patterns.

Fig. 4

Plot of the eigenvalues obtained by dispersion relation (31). a (blue) and (red) for the steady state ; b (blue) and (red), for the steady state . The model parameters are given in Table 2. The continuous curves represent the , as the imaginary part of the eigenvalues is zero (represented by dotted lines) in the case of the parabolic model [see relations (31)–(33)]. The diamonds on the x-axis represent the discrete wave numbers (Color figure online)

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 of the parabolic model, compared to the one of the hyperbolic model for the steady state . Although in Fig. 2d the effect of on the stability of the steady state is not very clear, in Fig. 5a we can see that a decrease in leads to a shift of the graph below and a change in stability for very low . Moreover, when we increase (see in Fig. 5b), the steady state becomes stable and we do not expect to have aggregation.

Fig. 5

The dispersion relation (31) for the steady state . a The effect of on the graph of ; b the effect of on the graph of ; the rest of the model parameters are given in Table 2. The diamonds on the x-axis represent the discrete wave numbers (Color figure online)

Remark 2

Note that is unstable for the trivial steady states [see Figs. 1 and 4a; a similar behaviour being observed in other non-local hyperbolic models (Eftimie et al. 2017)] 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 space–time plane choosing a time step and a space step . 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 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 , 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 statesor small random perturbations of rectangular-shaped aggregations located in the middle of the domainTo begin we first run numerical simulations for small random perturbations of the steady states and . 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 show similar patterns (not shown here) with those for , but with population vanishing much faster and population 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 is the proliferation rate, in Fig. 6 we see that as the proliferation rates increase, e.g. and (see Fig. 6a’, b’) the and 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 (e.g. up to 0.7) leads to more aggregation waves. Choosing now the same proliferation rate for the two populations, e.g. , we see in Fig. 6a”, b” an interesting effect of clonal competition, with population dominating the dynamics. This could be explained by the very low mutation rate and an equal competitive effect (i.e. ). Note that similar behaviours (where dominates the dynamics) have been observed for other sets of parameters with equal cell proliferation rates: e.g. or 0.7 (not shown here).

Fig. 6

Patterns exhibited by the hyperbolic model (1). The initial conditions for the two cancer cell populations are described by small random perturbations of the steady state [see (34)]. a, b Total density of and for and ; a’, b’ total density of and for and ; a”, b” total density of and for . The rest of model parameters are given in Table 2 (Color figure online)

The effect of turning rates on cancer cell movement and aggregation If we increase the baseline random turning rates, from to , 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 to , the rotating waves persist (not shown here). This behaviour was obtained for both steady states.

Fig. 7

Patterns exhibited by the hyperbolic model (1). The initial conditions for the two cancer cell populations are described by small random perturbations of the steady state [see (34)]. a, b Total density of and for and ; a’, b’ total density of and for and ; a”, b” total density of and for and . The rest of model parameters are given in Table 2 (Color figure online)

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 to ) leads to a change in cells’ movement from rotating waves (panels a, b) to standing waves (panels a”, b” shown only for 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 and pulse-like initial conditions (35)] is shown in Fig. 10.

Fig. 10

Patterns exhibited by the hyperbolic model (1) showing the cancer cell density for . The initial conditions for the two cancer cell populations consist of a rectangular pulse [see (35)]. a, b Total density of and for ; a’, b’ total density of and for . 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 , [note that similar results have been obtained also for the steady state —not shown here]. To investigate the effect of attraction and repulsion on pattern formation, we focus on two different cases: (i) and and (ii) and . 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 formation of a much larger number of smaller stationary cell aggregations (see Fig. 8a’, b’).

Fig. 8

Patterns exhibited by the hyperbolic model (1). The initial conditions for the two cancer cell populations are described by small random perturbations of the steady state [see (34)]. a, b Total density of and for and ; a’, b’ total density of and for and . The rest of model parameters are given in Table 2 (Color figure online)

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 and , 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 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.

Fig. 9

Patterns exhibited by the hyperbolic model (1) for and . The rest of model parameters as given in Table 2. The initial conditions for the two cancer cell populations consist of a rectangular pulse [see (35)]. a–d show the density of right-moving cancer cells (a, b) and (c, d). a’–d’ show the density of left-moving cancer cells (a’, b’) and (c’, d’). a”–d” show the total density of cancer cells (a”, b”) and (c”, d”) (Color figure online)

We note here that increasing the magnitude of alignment to leads again to rotating waves, while decreasing it to 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 () and cell baseline speed (). In contrast to the previous simulations, we now consider small speeds (i.e. ) and slightly larger turning rates (i.e. ). We see in Fig. 10a, b that when alignment is absent (), the cells form stationary aggregations (which eventually vanish due to the mutation term) and the cells spread throughout the whole domain.

When alignment is present (), we see in Fig. 10a’, b’ that some sub-populations of and cells move quickly to the left and to the right, reaching the domain boundaries. As in the previous case, the cells are eventually eliminated, while the 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 and (see Sect. 4.2)or small random perturbations of rectangular-shaped aggregations located in the middle of the domain

First, we assume that 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 , we see in Fig. 11a, b that for we obtain standing waves, and as decreases to 0.5, the cancer cells exhibit stationary pulses (Fig. 11a’, b’), as expected by the linear stability analysis results. Assume now , and focus on the parabolic model (12). For initial conditions that are small random perturbations of the steady state , we obtain stationary pulses (Fig. 11a”, b”) that are similar to the ones for the rescaled hyperbolic model when (see Fig. 11a’, b’). Similar pattern formation results are obtained for the hyperbolic steady state (and for the state corresponding to the parabolic model) for the case . If we increase the random turning rates to , we obtain stationary pulses for every in the rescaled hyperbolic model and for in the parabolic model, which was expected also from the linear stability analysis (see Fig. 3, where at the wave numbers where ).

Fig. 11

The spatiotemporal patterns obtained with the hyperbolic model (1) for after scaling, and the parabolic model (12). a, b Standing waves obtained by (1) after scaling for ; a’, b’ stationary pulses obtained by (1) after scaling for ; a”, b” stationary pulses obtained by (12) when . The initial conditions for the two cancer cell populations are described by small random perturbation of the steady state , for the rescaled hyperbolic model, and , for the parabolic model [see (34) and (36)]. The rest of model parameters are given in Table 2 (Color figure online)

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 (same as in Fig. 8a, b, but here ). 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.

Fig. 12

Patterns exhibited by the parabolic model (12). The initial conditions for the two cancer cell populations are described by small random perturbations of the steady state [see (36)]. Total density of and for and . The rest of the model parameters are given in Table 2 (Color figure online)

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.

Fig. 13

Patterns exhibited by the parabolic model (12) showing the cancer cell density for . The initial conditions for the two cancer cell populations consist of a rectangular pulse [see (37)]. a, b Total density of and for and . The rest of model parameters are given in Table 2 (Color figure online)

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 non-local 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.

Table 1

A list of model variables with their units

Variable	Description	Dimensional units
u1+	Right-moving early-stage cancer cell density	cell/length
u1-	Left-moving early-stage cancer cell density	cell/length
u_1	Total early-stage cancer cell density	cell/length
u2+	Right-moving late-stage cancer cell density	cell/length
u2-	Left-moving late-stage cancer cell density	cell/length
u_2	Total late-stage cancer cell density	cell/length

Since we are in 1D, length and volume coincide and we express the variables in terms of domain length