Mathematical modelling of cancer invasion : The multiple roles of TGF-β pathway on tumour proliferation and cell adhesion

In this paper, we develop a non-local mathematical model describing cancer cell invasion and movement as a result of integrin-controlled cell–cell adhesion and cell–matrix adhesion, and transforming growth factor-beta (TGF-β) effect on cell proliferation and adhesion, for two cancer cell populations with different levels of mutation. The model consists of partial integro-differential equations describing the dynamics of two cancer cell populations, coupled with ordinary differential equations describing the extracellular matrix (ECM) degradation and the production and decay of integrins, and with a parabolic PDE governing the evolution of TGF-β concentration. We prove the global


Introduction
Cellular adhesion, i.e. cell-cell and cell-matrix adhesion, and cellular proliferation are fundamental features of multicellular organisms, linked to maintenance of order in the organisms, e.g.tissue formation, stability and breakdown. 3These interactions between cells and the extracellular matrix (ECM) are mediated through cell surface receptors, a major group of which is represented by the integrins, 68 and various cytokines and chemokines.Another group of molecules involved in cell-cell adhesion is represented by the cadherin families. 26There are several signalling pathways that control normal cell processes like cell proliferation, division, cellular adhesion and apoptosis, with transforming growth factor β (TGF-β) pathway to be one of the most critical.Belonging to a large family of multifunctional polypeptides, TGF-β regulates the proliferation, differentiation, adhesion, migration and apoptosis of many cell types, including endothelial cells, hematopoietic cells and lymphocytes, 47 and ECM production. 31arious signals, including integrin, Notch, Wnt, TNF-a, and EGF signals, have been reported to cooperate or synergize with TGF-β signalling and stimulate tumour invasion and metastasis. 47Experimental studies 44 showed that the loss of TGF-β responsiveness is one of the events that initiate fibrotic disease and malignant progression of cancer, as well as cancer metastasis. 65TGF-β induces morphological, biochemical and transcriptional changes towards a mesenchymal phenotype, a process called epithelial to mesenchymal transition (EMT) (see Refs. 42 and 51 and many references therein).EMT occurs when epithelial cells lose their epithelial cell characteristics and become mesenchymal.Mesenchymal cells can return to an epithelial phenotype, a process called mesenchymal-epithelial transition (MET).Through these processes, cancer cells become metastatic and form new colonies at distant sites.
Experimental studies 33,41 have shown that tumours consist of heterogeneous populations of cells, which are the result of genetic instability.Intra-tumour heterogeneity appears in almost all phenotypic cell features: from cell morphology, to gene expression, motility, proliferation, immunogenicity and metastatic potential. 45hile both normal cells and cancer cells appear to be heterogeneous for various characteristics (e.g.surface antigens), cellular heterogeneity is shown 54 to be more pronounced in malignant neoplasms.Experimental studies have shown complex interactions between clonal cancer cell sub-populations in heterogeneous tumours: from stable coexistence to competitive exclusion. 38The metastatic and invasive potential of heterogeneous tumours is influenced by the interactions among the of bounded solutions to our system as the vanishing viscosity limit of a classical solution to an associated parabolic problem.In Sec. 4, we undertake numerical simulations to investigate the effect of TGF-β on cancer cell movement, and observe a range of patterns obtained for different values of parameters of the model.Finally, in Sec. 5, we summarise our results and give some concluding remarks.

The Mathematical Model of TGF-β Regulatory
Network in Cancer TGF-β plays a crucial role in embryonic development, wound healing and cancer.Moreover, TGF-β signalling stimulates EMT in certain epithelial cells (see Ref. 51  and many references therein) and consequently induces various diseases, including cancer.The way that TGF-β interacts with cancer cells varies between early and late stages of cancer (see Fig. 1), making its behaviour difficult to analyse.We consider a two-population model describing the behaviour of an early stage cancer population and a late stage cancer (descendant clone) population, which interact with each other, as well as with the ECM, via long-range integrin-controlled adhesive and repulsive forces 18,23 on bounded spatial domain Ω ⊂ R n with smooth boundary ∂Ω.For T > 0 let Ω T = (0, T ) × Ω. Denote by u 1 (t, x) and u 2 (t, x) the density of early and late stage cancer cells, respectively, at position x and time t, by f (t, x) the ECM density, and by c(t, x) the density of integrin receptors on the surface of cancer cells (receptors involved in cell-cell and cell-matrix interactions).Finally, we denote by b(t, x) the TGF-β concentration.For compact notation, we define the vectors u(t, x) = (u 1 (t, x), u 2 (t, x)) and υ(t, x) = (u(t, x), f(t, x)) .The multiple roles of TGF-β pathway on tumour 5 Cancer cells dynamics.Cancer cells can switch from a homogeneous type of invasion to a heterogeneous type of invasion described by (directionally moving) invading chains. 12Therefore, we assume that the movement of the early stage cancer cell population u 1 is governed by random motility (which underlines a homogeneous type of invasion), as well as directed motility in response to cell-cell and cell-matrix adhesive forces (which underlines the heterogeneous type of invasion). 9et D u describe the random motility coefficient and F 1 [u, f, c, b] describe the nonlocal directed motility.In contrast, the late stage cancer cell population, u 2 , moves only in a directed manner (hence exhibiting a heterogeneous type of invasion) in response to cell-cell and cell-matrix adhesion forces (described by a non-local term ).Moreover, the u 1 cells can mutate into u 2 cells at a constant rate M .TGF-β has been found to have bidirectional functions in the progression of cancer.
In early stages of cancer, TGF-β is an antiproliferative and proapoptotic signal, while in late stages of cancer it acts as a tumour promoter. 57Thus, we have the following equations describing the dynamics of the two cancer cell populations: ) Taking into account the effect of TGF-β on cancer cell proliferation and assuming that both u 1 and u 2 cells can proliferate in a logistic manner (to describe the observed slow-down in tumour growth following the loss of nutrients 37 ), we choose the growth functions to be given by where r 1 and r 2 are the growth rates of the u 1 and u 2 populations, respectively, k u is the carrying capacity, b m is the maximum TGF-β concentration, and c b is a coefficient related to the effect of TGF-β on cancer cell proliferation/decay.In particular, the term (−1) i models the anti-tumour effect of TGF-β on early tumours (i = 1), and the pro-tumour effect on late tumours (i = 2).Note that these growth functions incorporate also the principle of competition between clonal sub-populations in heterogeneous tumours. 38he non-local cell-cell and cell-matrix adhesion and repulsion forces for cancer cell populations u 1 and u 2 , are described by a function that depends on cell densities, ECM and integrin densities, and concentrations of TGF-β molecules given by the following relation where , describe the nature of the cell-cell and cell-matrix adhesive forces.These functions increase when the cell density and ECM density increase, and accordingly they decrease when the cell density and ECM density decrease.The functions g i , i = 1, 2, are given by Integrins are molecules known to have a regulative role in cell-cell and cellmatrix adhesion, 4 while the role of TGF-β in cellular adhesion is dual: (i) Promotes cell-matrix adhesion by inducing the synthesis and the secretion of ECM-adhesion molecules laminin and fibronectin and the upregulation of integrin expression for these matrix-adhesion molecules 30,66 ; (ii) Decreases cell-cell adhesion. 52,67Thus, to define these adhesion strength functions we consider the integrin density, c, and TGF-β concentration, b.Since cell mutation could lead to more integrins, 35 we consider strength functions with different integrin levels for each of the two populations.The more integrins a cell has, the stronger its adhesion force.Therefore, biologically realistic choices for these adhesion strength functions are the increasing, bounded and positive functions given by: where a i , a bi , d, d b , e i , e bi and s i * , s * , c i * , i = 1, 2, are positive real numbers.
ECM dynamics.The ECM is considered as non-motile matter, with changes to its density due to degradation by u 1 and u 2 cell populations upon contact at rates α > 0 and β > 0, respectively, and ECM density remodelling back to normal levels, at a constant rate of δ > 0.Moreover, TGF-β induces the synthesis of ECM adhesion molecules, 66,67 at a rate of θ β > 0. Thus the dynamics of ECM, f (t, x), is described by: where f m is the maximum ECM density at which the ECM fills up all available physical space.

Integrin dynamics.
We assume that the level of integrins depends on cancer cell density, such that cell mutation changes the density of receptors (since in highly metastatic cancers, the expression of integrins is up-regulated 35 ).Moreover, TGF-β signalling up-regulates the integrin expression, 4,42 at a rate of p 3 .Therefore, the dynamics of integrins, c(t, x), can be described by: where q is the decay rate of c, and p 1 and p 2 are the production rates of integrins by u 1 and u 2 cancer cell populations, respectively.To model the increase in receptors on highly mutated cancer cells, we assume that p 2 > p 1 (see Table A.2).
TGF-β dynamics.Finally, TGF-β is assumed to diffuse freely in the spatial domain, after being released by u 1 and u 2 cells, and decay at a rate of q b > 0. Therefore, the dynamics of TGF-β, b(t, x), is described by: where D b is the TGF-β diffusion coefficient and λ(u) is the TGF-β production term.Here, we choose λ(u) = µ 1 u 1 + µ 2 u 2 , with µ 1 and µ 2 to be the production rates of TGF-β by u 1 and u 2 , respectively.The relations (2.1)-(2.9)are summarised in the following system: ) ) We impose the following initial conditions: (2.11) Finally, we assume that there is no-flux of both cancer cells and TGF-β proteins on the boundary of the domain, and where ν is the outward unit normal vector to ∂Ω.

Existence of Solution
To prove the existence of solution for system (2.10), we use the theory of semigroups combined with the vanishing viscosity method (to transform Eq. (2.10b) into a parabolic equation).Then we show that in the vanishing viscosity limit, we obtain weak solutions for (2.10).We note that the steps in this proof of existence The link between these two proofs is based on the extraction of the appropriate estimates for the vanishing viscosity limit, for our more complex system of non-local parabolic-hyperbolic equations coupled with ODEs.

Existence of approximate solution
We will approximate system (2.10) by the following system: ) ) for 0 < ≤ 1, where The ICs are given by Finally, the BCs corresponding to (2.12)-(2.13)are given by the relations and Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
For the full non-local interaction terms (2.3)-(2.4)we make the following assumptions: where , is a continuous function, which satisfies and 2, (given by (2.6)) are bounded, we assume that there is a constant L N , which depends on the bound for uniformly with respect to (x, y) ∈ Ω 2 .We assume that h 2 , h 3 and h 4 exist, and that Moreover, based on relation (3.2) we can assume that there are constants Λ j > 0, j = 1, . . ., 4, such that for We now consider the sectorial operators in the space X = L p (Ω), with common domain of definition Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
and Re(σ(A j )) > 0, j = 1, 2, 3. Then the fractional powers are well defined with the graph norm Then from Ref. 27, we have the following embeddings: where C 0,r ( Ω) is the space of all Hölder continuous functions with exponent r in Ω.Notice that for (3.18) and (3.19) are satisfied.Moreover, since A 1 , A 2 and A 3 are sectorial operators, then each of −A 1 , −A 2 and −A 3 is the infinitesimal generator of an analytic semigroup {e −tAj } t≥0 , j = 1, 2, 3. Therefore, there exists a positive constant C γ such that the following inequality holds 27 : where 0 < ξ j < Re(σ(A j )), j = 1, 2, 3, and where k γ positive constant.
) and (3.20) are satisfied, then for any T > 0 there exists a unique global-in-time solution , which remains bounded and the bounds are -independent.Moreover, the solution satisfies Proof.We will prove the existence of a local-in-time solution using the Banach contraction theorem.We first focus on the ODEs (3.1c) and (3.1d).We notice that relation (3.18) implies that u 1 , u 2 , b ∈ W 1,p (Ω), and since the functions h 2 : R × R × R → R, h 3 : R × R → R and h 4 : R → R are locally Lipschitz, we have by the property of superposition operator that the value of the functions h 2 , h 3 and h 4 is also in W 1,p (Ω). 58The space W 1,p (Ω) for p > n is an algebra with pointwise multiplication, and thus it follows that the functions (u 1 , u 2 , f , b ) → For a fixed T > 0 we note that functions: belong to the space C([0, T ]; W 1,p (Ω)).Then, the system of the PDEs (3.1a), (3.1b) and (3.1e) with (3.3) can be rewritten as: with and or, equivalently, we write that 27 and the mapping H : ) is defined as the mapping of the right-hand side of Eqs.(3.1a), (3.1b) and (3.1e).Using similar arguments as before for P , and assumption (3.9) we deduce that H : For a fixed T > 0, we define the space where We define the mapping J : where k γ satisfies relation (3.22), and max{ We show that J maps B R into itself and that J is a strict contraction.By using relations (3.21)-(3.22)and (3.29) we obtain, for T small enough, Similarly we obtain and Moreover, for J 3 and J 4 we have Hence, we can choose T sufficiently small such that {M R Cγ as it can be easily proved using inequality (3.22).Thus J maps B R into itself.
Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
If y 1 , y 2 ∈ B R then for t ∈ [0, T ] we have where L R is the Lipschitz constant of H. Similar estimates can be obtained for the differences of the rest of the arguments.Therefore, it follows that Let us now prove the uniqueness of solution.Let ) be two solutions of system (3.1), with the same initial conditions.By linearity y = y 1 − y 2 is a solution of (3.1) with zero initial conditions.Then, since all nonlinear terms are Lipschitz continuous and the components of the solution are L ∞ -bounded functions on bounded time intervals, we have for a constant k 0 , and since y(0) L 2 (Ω) = 0, Gronwall's inequality implies that y(t) L 2 (Ω) = 0 for all t ≥ 0, so y = 0.The equation for the ECM density given by (3.1c) does not involve any spatial derivatives and x behaves as a parameter.Thus, it is an ordinary differential equation in which the right-hand side is zero when f (t, x) = 0 and for which local Lipschitz conditions hold.Therefore, from Picard-Lindelöf theorem we obtain a local unique solution for the initial value problem (3.1c), with f (0, x) = 0.In the same way, we obtain a local unique solution for the initial value problem (3.1c), with f (0, x) ≥ 0. Therefore, since f (0, x) ≥ 0, from uniqueness of solutions we have f (t, x) ≥ 0 for all t > 0, x ∈ Ω.Then from maximum principle arguments it follows from system (3.27 Finally, the equation for the integrin density given by (3.1d) can be treated in a similar manner as Eq.(3.1c) for the ECM density.Again we have an ordinary differential equation in which the right-hand side is greater than or equal to zero when Similarly, we have hence from relations (3.12)-(3.13)and Gronwall's inequality again we have sup Thus, from relations (3.6)-(3.9)we have for all t ∈ [0, T max ): From system (3.27),we rewrite the elliptic operators in the form: where we denote by and and where M γi , M γ3 , i = 1, 2, are constants depending on M ∞ .
We show now the global existence of solution by contradiction.Let us suppose that for T max < ∞ we have where C Tmax := C Cγ Mγ i ,γ,Tmax , i = 1, 2, is a continuous function increasing with respect to T max , while for the function b we have from relation (3.48) that sup Note that k γ and C γ are -independent constants since 0 < ≤ 1 (see Theorem By Eq. (3.1c) and direct calculations, we obtain where . By Gronwall's inequality and previous estimates, we have Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
The multiple roles of TGF-β pathway on tumour 17 This, together with Gronwall's inequality and previous estimates, yields

Vanishing viscosity limit
Now we are ready to take the vanishing viscosity limit → 0 and prove the existence of solution for system (2.10).First we introduce the notion of a weak solution to problem (2.10)-(2.11)with (2.12)-(2.13).
) is called a weak solution of the problem (2.10)-(2.11)with (2.12)-(2.13)if: Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.We show now the L 1 -estimates with respect to time, which will be used in the proof of existence of solution to model (2.10).Theorem 3.2.Let the assumptions of Theorem 3.24 hold.Then for each ρ > 0 there exist nondecreasing functions , such that for any ∈ (0, 1] and for any t, t + ∆t ∈ [0, T ], ∆t ≥ 0 we have for a ball B ρ = {|x| ≤ ρ} that :

The multiple roles of TGF-β pathway on tumour 19
Proof.Let us consider a function g ∈ C 2 0 (Ω) with supp(g) ⊂ B ρ .Then from the estimates obtained by Theorem 3.24 we have Similarly we can obtain estimates for u 1 , f, c and b.Then, as in Ref. 29, the L 1 -estimates of u 1 , u 2 , f , c and b with respect to time t follow.
. This implies that the convergence is even pointwise a.e. for a suitable subsequence.It remains to show the initial conditions (3.63)-(3.66).Let us first define the set E ⊂ [0, T ] such that for all t 0 ∈ [0, T ]\E we have for almost all x ∈ Ω that (t 0 , x) is Lebesgue point of u 1 , u 2 , f, c and b.The set E has Lebesgue measure zero.For any fixed t 0 ∈ [0, T ]\E and ρ > 0 it follows from Theorem 3.2 that The pointwise convergence of The properties of ω ui ρ (t 0 ), i = 1, 2, give (3.63) since u i has compact support.Similarly we obtain relations (3.64)-(3.66).

Non-dimensionalisation of the model
In this section, we investigate numerically the type of patterns exhibited by model (2.10) in the one-dimensional case.Let R s > 0 be the cells sensing radius (i.e. the Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only. The multiple roles of TGF-β pathway on tumour 21 maximum range over which cells can detect other surrounding cells).We consider a bounded domain Ω = [0, R s ], and following the approach in Ref. 24, we choose the non-local terms F i [u, f, c, b], i = 1, 2, to be given by where η(k) = (−1) k , k = 0, 1 and g i , i = 1, 2, as described in Sec. 2 (see relation (2.5)).
Let us define the kernel K, assuming that it is attractive at medium/long ranges (i.e. at the edges of the cell) and repulsive at very short ranges (i.e. over cell surface), and thus can be defined as with q a and q r describing the magnitudes of attractive and repulsive interactions, respectively, and K a (x) and K r (x) describe the spatial ranges over which these interactions take place.We consider translated Gaussian attraction and repulsion kernels (as in Ref. 20): where s a and s r represent half of the length of attraction and repulsion ranges, respectively, with s r < s a .Also, m j = s j /8, j = a, r, represent the width of the attractive and the repulsive interaction ranges.To perform numerical simulations, we first non-dimensionalise system (2.10) by using 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. 2The time scale is defined as τ := L 2 0 /D τ , where D τ is the characteristic diffusion coefficient (∼10 −6 cm 2 s −1 ).Furthermore, we rescale the cancer cells, the ECM, the integrins and the TGF-β with k u , f m , c m and b m , respectively.Here, k u is the carrying capacity of the cancer cell populations and it is taken to be ∼6.7 • 10 7 cell/volume, and f m is the maximum ECM density at which the ECM fills up all available physical space and it is taken to be equal to 4 mg/volume, as in Ref. 19 is the maximum integrin density and it is taken to be 5 • 10 4 integrins per cell (as in Ref. 7), while b m is the maximum TGF-β concentration taken to be equal to 141.59 ng/volume (as in Ref. 34).We choose the dimensionless functions K(r) : Therefore, the non-local terms are given by Fi Finally, we obtain the dimensionless parameters: After dropping the tildes for notational convenience, we obtain the following nondimensionalised system: ) ) )

Pattern formation
To discretize our model we use a time-splitting approach.We use a Crank-Nicolson scheme to propagate the solution of the diffusion term.Then, we use the Nessyahu-Tadmor scheme 53 for the time-propagation of the advection terms.Finally, for the time-propagation of the reaction terms 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).For this reason, the integrals are wrapped-up at the boundaries.The simulations ran for times up to t = 1000, but for clarity in Figs.The initial conditions for the cancer cell populations are small random perturbations of rectangular-shaped aggregations located in the middle of the domain u i (0, x) = u c i + rand(0, 10 −4 ), x ∈ (L/2 − 1, L/2 + 1), 0, everywhere else, ( with u c 1 = 0 and u c 2 = 0.1.For the ECM density, f , we assume that the tumour has already degraded some of its surrounding tissues: Finally, the integrin density and TGF-β concentration, c and b, respectively, are proportional to the initial tumour cell density c(0, x) = 0.5u 1 (0, x) + 0.5u 2 (0, x) (4.9) and b(0, x) = 0.05u 1 (0, x) + 0.05u 2 (0, x).(4.10) To investigate the effect of TGF-β signalling on cell proliferation, movement and aggregation (the last two aspects being controlled by cell adhesion), we focus on three possible cases for the magnitudes of cell-cell and cell-matrix adhesion.For each of these three cases, we investigate the dynamics of u 1 and u 2 populations  Reducing now the magnitudes of cellular adhesion forces to s * i = s * = c * i = 0.1, i = 1, 2, we see in Fig. 5 that irrespective of the absence/presence of TGF-β, population u 1 forms a stationary aggregation that eventually vanishes for large times, while population u 2 exhibits a travelling wave.This behaviour might be explained by the combined effect of high mutation rate and clonal competition (see the logistic growth terms in (2.2)), since adhesive forces are very small and lead to the spread of population u 2 .In contrast to the dynamics in Figs.3(b) and 3(b ), and 4(b) and 4(b ), where the u 2 cells seem to travel slower towards the boundaries in the presence of TGF-β (compared to the absence of TGF-β), in Fig. 5 the u 2 cells travel faster to the boundaries in the presence of TGF-β.We deduce from here that the spread of tumour cells depends both on the magnitude of adhesive forces as well as on the presence of TGF-β molecules.

The multiple roles of TGF-β pathway on tumour 27
We note here that we investigated numerically the case when TGFβ is present, but does not influence cell proliferation or cell adhesion (i.e. for θ β , p 3 , µ 1 , µ 2 = 0, and c b = a bi = d b = e bi = 0, i = 1, 2).The patterns (not shown here) that we obtained were similar to those presented in Figs.2(a) and 2(b), 3(a) and 3(b), 4(a) and 4(b), 5(a) and 5(b).This suggests that the effect of TGF-β on the cancer cell density is greater than the effect on the ECM and on the integrins density.

Conclusion and Discussion
In this paper, we introduced a model of integro-differential equations describing the dynamics of early stage and late stage cancer cell populations, under the effect of TGF-β signalling.The model was then used to investigate the role of TGF-β on cellular adhesion and proliferation.
We first proved the global existence of bounded solutions to our non-local model by taking a vanishing viscosity approach and approximating our model with a nonlocal parabolic PDE.The proof used the Banach contraction mapping theorem, the Möser-Alikakos method and the vanishing viscosity method.
We then investigated numerically the solution of this non-local model, paying particular attention to the effect of TGF-β on cell-cell and cell-matrix interactions.We showed that: (i) In the absence of TGF-β, the magnitudes of cell-cell and cellmatrix interactions influenced the formation of cancer cell aggregation at specific position in space (see Figs.While the numerical investigation of cancer spread uncovered an interesting combined effect of cell adhesion and presence/absence of TGF-β, in the future we plan to investigate analytically the travelling waves and study the effect of parameters related to TGF-β and cell adhesion on the speed of these waves.• Various experimental studies 15,50 shown that doubling times for tumour cells range from 1−10 days.This corresponds to growth rates between (ln(2)/10, ln(2)/1) = (0.07, 0.7).In this study, we assume that r1 = 0.1 and r2 = 0.2.• Experimental studies 14,28,43 have shown that the mutation rate ranges between M = 10 −3 /day and M = 0.1/day.Thus the non-dimensional value of the mutation rate is in the range between M = 0.001 and M = 0.1 (for highly aggressive tumours).In this study, we choose M = 0.05.• The parameters a i , d, e i , s * i , s * , c * i , i = 1, 2, were based on the range of the adhesion strength parameters used in Ref. 3.
• Experimental studies 35 have shown greater production of integrins for mutated cancer cells.Thus, we choose p 1 < p 2 .• Experimental studies 16,17,39,40 have shown that the half-lifes of the integrins range from 0.04−4 days.This corresponds to a decay rate between (ln(2)/4, ln(2)/0.04)= (0.17, 17.3).In this study, we assume that q = 0.3.• The remodelling rate was chosen to be greater than cell proliferation rate, as considered also in Ref. 11.

Fig. 1 .
Fig. 1.A caricature summarising the dual role of TGF-β in cancer progression.Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.

. 5 )
where S i (c, b) is the cell-cell self-adhesion strength function for populations u i , S(c, b) is the cell-cell cross-adhesion strength function between the two populations, and C i (c, b) is the adhesion strength function between population u i and ECM.

. 56 )Remark 3 . 1 .
Bounds(3.50),(3.54) and (3.56) contradict(3.49),therefore the solution exists globally.It is easy to see how the proof of Möser-Alikakos method in Ref. 13 can be used in our case to show the -independent L ∞ -estimates given in relation(3.44).Following along the same lines with the proof of Lemma 9.3.1 in Ref.13, we obtain a constant c (as described in relation (9.3.11) of the proof in Ref. 13) such that c := max{const.a 0 , 1, D|Ω|, K 2 }, (3.57)where a 0 will be each of the diffusion coefficients D u , and D b .Using the fact that 0 < D u , , D b ≤ 1 it follows that c is -independent.By this result, it follows that the rest of the estimates given by relations (3.45)-(3.48),(3.50)-(3.51),(3.54) and (3.56) are -independent.
Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.when TGF-β is absent and does not influence cell proliferation or cell adhesion (i.e. for c b = a bi = d b = e bi = θ β = p 3 = µ 1 = µ 2 = 0, i = 1, 2), and when TGF-β is present and influences both cell proliferation and cell adhesion (i.e. for c b = 20 and a bi , d b , e bi , θ β , p 3 , µ 1 , µ 2 , = 0, i = 1, 2).(i) Cell-cell adhesion < cell-matrix adhesion.To investigate the effect of greater cell-matrix adhesion, we choose s * 1 = 1.8, s * 2 = 0.6, s * = 1, c * 1 = 1.9 and c * 2 = 2.5 and the rest of model parameters as given in Table A.2.We see in Figs.2(a) and 2(b) that in the absence of TGF-β, the population of early stage cancer cells (u 1 ) decreases, while the population of late stage cancer cells (u 2 ) increases and dominates the long-term dynamics.This behaviour is expected due to the mutation term "−Mu 1 ", and due to large cell-matrix adhesion, which impedes cells to move and thus leads to the formation of stationary pulses for t > 50.Considering now the effect of TGF-β, we see in Figs.2(a ) and 2(b ) that population u 1 vanishes faster, due to the presence of antiproliferative and proapoptotic signals from TGF-β (described by c b > 0 in Eq. (4.6a)).Population u 2 persists and increases significantly, due to the promoting effects of TGF-β on the late stages of cancer, which also induces the movement of the cancer cells (via EMT) thus leading to their spread over the domain until they reach the boundaries.

Fig. 3 .
Fig. 3. (Colour online) Patterns exhibited by model (4.6) showing the effect of TGF-β on cancer cell density for cell-cell adhesion greater than cell-matrix adhesion, i.e. s * 1 = 2.4, s * = 2.1, s * 2 = 2, c * 1 = 1.1 and c * 2 = 0.9.The rest of model parameters are given in Table A.2. (a), (b) Density of u 1 and u 2 populations in the absence of TGF-β.The inset in panel (b) shows the long-term dynamics of u 2 (t, x) (for t ≤ 1000); (a ), (b ) Density of u 1 and u 2 populations in the presence of TGF-β.Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.

Fig. 5 .
Fig. 5. (Colour online) Patterns exhibited by model (4.6) showing the effect of TGF-β on cancer cell density for the same cell-cell and cell-matrix adhesion, i.e. s * i = s * = c * i = 0.1 and a i = d = e i = 0.5, i = 1, 2. The rest of model parameters are given in Table A.2. (a), (b) Density of u 1 and u 2 populations in the absence of TGF-β; (a ), (b ) Density of u 1 and u 2 populations in the presence of TGF-β.Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
2(b), 3(b) and 4(b)); (ii) The consideration of TGF-β leads to the spread of mutated (i.e.u 2 ) cancer cells over the whole domain mainly in a travelling-wave manner (with no cell aggregations; see Figs. 3(b ), 4(b ) and 5(b )).We also emphasise that the speed at which cells spread depended on the presence/absence of TGF-β and on the magnitudes of cell adhesion forces (see Figs. 4(b) and 4(b ) versus Figs.5(b) and 5(b )). .8) Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
. Finally, c m Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.
2-5 we show mainly the dynamics for t ≤ 400.If the patterns do not reach a steady state before t = 400, we add inset figures showing the dynamics for t ≤ 1000. .Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.

Table A .
1.A list of variables with their units.Since we are in 1D, length and volume coincide and we express the variables in terms of domain length.Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.Math.Models Methods Appl.Sci.Downloaded from www.worldscientific.comby UNIVERSITY OF ST ANDREWS on 07/06/17.For personal use only.