## Abstract

Cancer is a disease characterized by abnormal division of cells combined with the malignant behaviour of these cells. Malignant cancer cells tend to spread, either by direct growth into adjacent tissue through invasion and dispersal or by implantation into distant sites by metastasis. In this paper, we will concentrate on the first type of tumour growth and consider a reaction–diffusion model for competition cells with *nonlinear diffusion term* modelling contact inhibition between normal and tumour cell populations for which wave propagation is usually observed in clinical data. Mathematically, based on a combination of perturbation methods, the Fredholm theory and the Banach fixed point theorem, we theoretically justify the existence of the travelling wave solution. Numerical simulations are finally illustrated to confirm our rigorous results.

## 1. Introduction

Cancer is a disease developed, with few exceptions, from mutations on single somatic cells that may divide uncontrollably, invade adjacent normal tissues and give rise to secondary tumours (metastasis) on sites different from their primary origin (Hanahan & Weinberg 2000). Although tumour progression involves a complex network of interactions among cancer cells and their host microenvironment (Israel 1996), it is observed that all neoplasms evolve according to a universal scheme of progression (Clark 1991; Evan & Vousden 2001). For tumour growth, viewed usually as a competition process between tumour cells and surrounding normal tissue cells, numerous mathematical models have recently been investigated. These include some examples based on ordinary differential equations (ODEs) modelling tumour growth and tumour–host interaction as competing populations (Michelson & Leith 1991; Gatenby 1995), and reaction–diffusion systems modelling the dispersal behaviour of tumour cell growth (Sherratt 1993; Gatenby & Gawlinski 1996; Sherratt 2000) and also some particularly novel models that reflect the cancer evolution and its interaction with the immune system (Sherratt 1993; Bellomo & Preziosi 2000).

The ODE model of Gatenby (1995) is an important reference to the study of tumour growth. Gatenby (1995) modelled the interaction of the tumour and normal cells as populations competing for space and other resources in some small volume of tissue within an organ. Using the Lotka–Volterra population competing principle (see Murray 2002), Gatenby established some critical mathematical parameters that can be used to control the outcome of different stages of neoplasm–host competition. He also proposed novel modes of therapies based on tumour classification and treatment strategies, and his results predict that therapies directed only at the tumour population will generally be inadequate. Successful tumour therapy requires enhancement of the competitive state of adjacent normal cells. Treatment strategies must therefore be developed to increase the carrying capacity of the environment for normal cells.

With regard to interaction between the normal and tumour cell populations, diffusion has also been used very successfully in models for spatial spread (Sherratt 1993; Gatenby & Gawlinski 1996; Sherratt 2000). Diffusion terms in reaction–diffusion equations of tumour growth are broadly divided into two categories. One is the linear diffusion term (e.g. Sherratt 1993) and the other is the nonlinear one (Gatenby & Gawlinski 1996; Sherratt 2000). Sherratt (1993) developed a reaction–diffusion model with the linear diffusion term for the initial growth of a tumour. This model also incorporated the immune response to the cancer cells. The author analysed the partial differential equations system for travelling wave solutions and obtained a lower bound on the wave speeds. Under biologically relevant approximations, a necessary and sufficient condition was derived for the existence of a travelling wave solution. Moreover, this model predicted that there is a critical level of immune response above which the immune system prevents the initial growth of the tumour. For the nonlinear diffusion model, it first goes back to Gatenby & Gawlinski (1996), in which the authors presented a nonlinear reaction–diffusion system with three equations describing the spatial distribution and temporal development of tumour tissue, normal tissue and excess H^{+} ion concentration. In view of the fact that the neoplastic tissue is unable to spread without the surrounding healthy tissue, the authors assumed that the flux of neoplastic tissue is dependent on the normal cells' movement. They found that the tumour spreading speeds determined via the marginal stability analysis of the model were consistent with the tumour growth rates *in vivo*.

Recently, Sherratt (2000) introduced a new nonlinear diffusion term that reflects the phenomenon known as contact inhibition of migration between different cell populations. Based on the Lotka–Volterra competition (Murray 2002), Sherratt considered the following reaction–diffusion model for the early stages of solid tumour growth:(1.1)

As stated in Sherratt (2000), *u*(*x*, *t*) and *v*(*x*, *t*) represent the normal cell and tumour densities, respectively. Mathematically, both variables are rescaled to ensure that the normal cell density *u*≡1 in normal tissue. When the parameter *γ* is greater than 1, the tumour cell population assumes a better growth rate than that of the normal cells. In the kinetic items, −(*u*+*v*) serves the decline effect in cell division rate owing to overcrowding. Sherratt (2000) introduced the new nonlinear diffusion term based on the consideration of contact inhibition which reflects the phenomenon that the movement of one population is affected by that of the other one. In accordance with the classical diffusion model for the spatial spread of population, the product of the fraction *u*/(*u*+*v*) and the overall flux −*∇*(*u*+*v*) indicates the flux for population *u*. The same situation happens with the population *v*(*x*, *t*) in the second equation. When the ratio *v*/*u* is a constant, the term is equal to (∂/∂*x*)*u*, and this can be seen from the following identity (see Sherratt 2000 for more details):By neglecting the highest derivatives in the travelling wave equations corresponding to system (1.1), Sherratt (2000) solved two coupled first-order equations and obtained the leading-order wavefront solutions.

To the best of our knowledge, no rigorous work has previously been done for the existence of travelling wavefronts when *the nonlinear diffusion terms are present*, as shown in system (1.1). The purpose of this paper is to revisit the competition equation with a new motility term modelling contact inhibition between cell populations, theoretically justifying the existence of a travelling wave solution with large wave speed in system (1.1) by adopting the asymptotic method developed by Faria *et al*. (2006) and Ou & Wu (2007). Although here we use the same idea as that in Faria *et al*. (2006) and Ou & Wu (2007), there are crucial differences in that the singularly perturbed terms are nonlinear and thus result in particular difficulties in transforming the differential system into an integral system so that the fixed point theorem can be applicable.

The remainder of this paper is organized as follows: in §2, we prove the existence of a travelling wave solution. Simulations are provided in §3 to not only confirm our result but also compare the difference between the leading-term solution and the real solution.

## 2. Travelling wave solution with large wave speed

In this section, we prove the existence of a spreading wave in our system with large wave speeds. Our first result, theorem 2.1, is essentially from Sherratt (2000). But for technical reasons, we need to reanalyse the result since we have changed the wave variable into *z*=*x*+*ct* instead of *z*=*x*−*ct* in Sherratt (2000). As can be seen below, the main result in this section is our theorem 2.4, which provides us with the justification for the existence of wavefronts.

To study the existence of travelling wave solutions with large wave speeds in system (1.1), we look for solutions that are functions of the travelling wave variable *z*=*x*+*ct*, with *u*(*x*, *t*)=*U*(*z*) and *v*(*x*, *t*)=*V*(*z*). In view of the symmetry of *u* and *v* in (1.1), we will rewrite the equations in terms of *V*(*z*) and *N*(*z*)≡*U*(*z*)+*V*(*z*)−1, which gives(2.1)

It is easy to see that (*N*, *V*)=(0, 0) and (*N*, *V*)=(*γ*−1, *γ*) are two equilibria solving the above system. Thus, we take the boundary conditions for (*N*, *V*) as

For wavefronts with large wave speeds, we want to develop the idea, originated from Canosa (1973) for the Fisher equation, to show the existence. We thus rescale the travelling wave coordinate by writing *ζ*=*z*/*c* for equations (2.1) to have(2.2)

Let *ϵ*=1/*c*^{2}, then *ϵ* is small when *c* is large. In particular, when *ϵ* is zero, system (2.2) reduces to(2.3)

Obviously, system (2.3) has two uniform steady-state solutions, (0, 0) and (*γ*−1, *γ*). By solving system (2.3), we have the following result, which is actually due to Sherratt (2000).

*Assume γ*>1. *Then system* (*2.3*) *has a heteroclinic orbit* (*N*_{0}(*ζ*), *V*_{0}(*ζ*)) *connecting the two uniform steady-state solutions* (0, 0) *and* (*γ*−1, *γ*).

From system (2.3), we can deduce that(2.4)which upon the use of separation of variable gives(2.5)where *A*≥0 is a constant of integration. From the second equation of (2.3), we have(2.6)

Let

Substituting into (2.6) and using (2.5) gives(2.7)Solving (2.7) by separation of variables gives

Therefore, we obtain the solutions

Since *A*≥0, the only case giving a positive solution for *N* is *A*=0. Thus, the wavefront solutions of (2.3) are(2.8)where *B*≥0 is a constant of integration.

Obviously, *N*_{0}(*ζ*) and *V*_{0}(*ζ*) satisfyThis proof is complete. ▪

It is easy to see that the origin (0, 0) is a saddle and the equilibrium (*γ*−1, *γ*) is a stable node.

As we said earlier, our primary aim in this paper is to rigorously establish the existence of a travelling wave solution to system (1.1) for large wave speeds, and we now proceed to achieve this. We will prove that the travelling wavefront to (2.2) can be approximated by the corresponding wavefront (*N*_{0}(*ζ*), *V*_{0}(*ζ*)) of (2.3) when *ϵ* is small. For later use, we now introduce some notations. Denote *C*=*C*(*R*, *R*) as the Banach space of continuous and bounded functions from *R* to *R* equipped with the standard norm . Let ,andwhere all the norms are defined byand

We observe that can be approximated by and hence assume that(2.9)where *N*_{0} and *V*_{0} satisfy system (2.3), i.e.(2.10)and *W*_{1} and *W*_{2} are subject to the boundary condition *W*_{i}(±∞)=0, *i*=1, 2. Substituting (2.9) into (2.2) and letting a prime be d/d*ζ*, using the first equation of (2.10), we have(2.11)where(2.12)From the second equation of (2.2), we have(2.13)

Now we prove that there exists a satisfying (2.11) and (2.13) when *ϵ* is small.

Since the equation *ϵλ*^{2}−*λ*−1=0 has two real roots *λ*_{1} and *λ*_{2} with(2.14)equation (2.11) is equivalent to the following integral equation:(2.15)For later use, by taking the derivative we can also obtain(2.16)Inserting formula *H* of (2.12) into (2.14) and using the fact thatthis gives(2.17)Let(2.18)then, by (2.10), we have from (2.13)Multiplying e^{ζ} to both sides of the above equation yieldsIntegrating it from −*∞* to *ζ*, we have(2.19)To work out the first integral in (2.19), we insert (2.18) and use the integration by parts, and make use of (2.10) and (2.16) to have(2.20)where(2.21)and(2.22)It is easy to show from (2.14) thatThus, we have from (2.17) and (2.20) that(2.23)where(2.24)and(2.25)

We will aim to show the existence of the solution to system (2.23). To this end, we define a linear operator *L* from the l.h.s. of (2.23) as *L*: *C*_{0}×*C*_{0}→*C*_{0}×*C*_{0} and

It is readily seen that *L*(*W*)∈*C*_{0}×*C*_{0} when *W*∈*C*_{0}×*C*_{0}. In order to use a fixed point theorem to verify the existence of the solution *W*∈*C*_{0}×*C*_{0} to system (2.23), we want to establish some estimates for the terms in the r.h.s. of (2.23) when *W*∈*C*_{0}×*C*_{0}. We have the following lemmas.

*For any given ϵ and* , *there exists a constant M*_{0} *independent of ϵ such that**uniformly for all ζ*∈(−*∞*, *∞*), *where B*(1/2) *is the ball in C*_{0}×*C*_{0} *with radius* 1/2 *and centre at the origin*.

Here we prove only the second estimate. The inequality of *A*_{1} can be done in a similar manner. From (2.25) we know thatUsing (2.21) we can get the first term(2.26)For the second term, when , we recall formula (2.22) to have(2.27)From (2.12), we have(2.28)Since *N*_{0} and *V*_{0} are bounded, a combination of (2.26)–(2.28) gives the desired result.

*For each δ*>0, *there is a σ*>0 *such that for any two elements ϕ and* *and* *we have that*(2.29)*and**where B*(*σ*) *is the ball in C*_{0}×*C*_{0} *with radius σ and centre at the origin*.

As in lemma 2.2, here we prove the estimate only for *A*_{1} and leave the estimate for *A*_{2} to interested readers. For any two elements and from (2.24) we haveDirect computation on the above formula givesSince the integral and 1/*λ*_{2}=*O*(*ϵ*) (see (2.14)), choosing *ϵ* and *σ* sufficiently small, we can obtain the estimate (2.29). The proof of this lemma is complete. ▪

Now we are ready to prove our main result.

*Assume γ*>1. *Then there is a constant c*^{*}>0 *such that, for any c*>*c*^{*}, *system* (*2.2*) *has a travelling wave solution* (*N*(*ζ*), *V*(*ζ*)) *connecting the two equilibria* (0, 0) *and* (*γ*−1, *γ*); *moreover, the wave profile* (*N*(*ζ*), *V*(*ζ*)) *converges to* (*N*_{0}, *V*_{0}) *when the wave speed c*→∞.

We first define a linear operator (2.30)whereand

The formal adjoint equation of *TΨ*=0 is defined byi.e.(2.31)

From now onwards, we will use the same argument as that in Ou & Wu (2007) to proceed with our proof. It contains five steps as follows.

*Step 1*. This step is about the property of the linear operator *T*. We want to prove that if is a solution of (2.31), then *Φ*=0. Moreover, we have *T*=*C*×*C*, where *T* is the range space of *T*.

Indeed, *N*_{0}→γ−1 and *V*_{0}→γ hold when *ζ*→∞. Thus, when *ζ* is large, the coefficient matrix satisfiesThis means system (2.31) tends asymptotically to the following system with constant coefficient matrix:(2.32)The eigenvalues *λ* of the coefficient matrix in equation (2.32) satisfy(2.33)The two roots of (2.33) have positive real parts as *γ*>1 and we thus know that any bounded solution to (2.32) must be the zero solution, and that any solution to (2.31) other than the zero solution must grow exponentially as *ζ*→∞. Then any bounded solution to (2.30) should be the zero solution. By the Fredholm theory (see lemma 4.2 in Palmer 1984) we have that *T*=*C*×*C*.

*Step 2*. We want to show that if is a bounded solution of *TΨ*=*Θ*, where then we have

Indeed when *ζ*→+∞, the system(2.34)tends asymptotically to(2.35)

For (2.35), it is easy to know that the *ω*− limit set of every bounded solution is just the critical point *Ψ*=0. Using theorem 1.8 from Mischaikow *et al*. (1995), we conclude that every bounded solution component of (2.34) also satisfies

Similarly, when *ζ*→−*∞*, system (2.34) tends asymptotically to(2.36)

Obviously, the coefficient matrix of (2.36) has two eigenvalues *λ*_{1}=−1<0 and *λ*_{2}=*γ*−1>0, then every bounded solution must satisfySince the main result in Mischaikow *et al*. (1995) is only valid for the *ω*− limit set, we thus invert the time from *z* to −*z* and use the result in Mischaikow *et al*. (1995) again to have that any bounded solution to (2.34) must satisfyHence the claim of step 2 holds.

*Step 3*. For the linear operatordefined bywe will prove that *L*, the range space of *L*, is equal to *C*_{0}×*C*_{0}, that is, for each we have a so thatTo this end, we suppose that *ξ*=*W*−*u* and deduce a system for as follows:

Differentiating both sides gives(2.37)In view of the results that *T*=*C*×*C* in step 1 and *u∈C*_{0}×*C*_{0}, it follows from step 2 that there exists: satisfying (2.37) and *ξ*(±∞)=0. Going back to the variable *W*, we have *W*=*ξ*+*u∈C*_{0}×*C*_{0}.

*Step 4*. Denote *N*(*L*) as the null space of operator *L*. It follows that (see lemma 5.1 in Faria *et al*. 2006) there is a subspace in *C*_{0}×*C*_{0} so thatObviously, is a Banach space. If we let be the restriction of *L* to , then is one-to-one and onto. Using the well-known Banach inverse operator theorem, we conclude that is a bounded linear operator.

*Step 5*. When the domain of *L* is restricted to , system (2.23) can be written asFrom lemmas 2.2 and 2.3, it follows that there exists *σ*>0 and 0<*ρ*<1 such that for andwhere

It is easy to know that for any , we haveTherefore, *F*(*ζ*, *Φ*) is a uniformly contractive mapping for . The Banach contraction principle gives that for *ϵ*∈[0,*σ*) (or ), system (2.23) has unique solution Returning to the original variable, we obtain that is a heteroclinic connection between the two equilibria (*γ*−1,*γ*) and (0,0). The convergence of to as *ϵ*→0 is easily seen. This completes the proof. ▪

## 3. Simulation

In this paper, we study the reaction–diffusion model with the nonlinear diffusion term modelling contact inhibition between the normal and tumour cell populations. By adopting the asymptotic method developed by Faria *et al*. (2006) and Ou & Wu (2007), which is a combination of perturbation methods, the Fredholm theory and the Banach fixed point theorem, we have theoretically justified the existence of a travelling wave solution with large wave speed in system (1.1).

In order to illustrate the validity of the theoretical result obtained in §2, we perform numerical calculations using the software Matlab.

Consider the following system:(3.1)with initial conditions(3.2)where *n*=*u*+*v*−1 and *γ* is taken to be 2, the constant *ξ* is the decay rate of the initial wavefront.

As shown below, the simulations here are not only for the purpose of confirming our theoretical results but also for comparing the leading-term travelling wave and the real wave. The contrast between the nonlinear diffusion and linear diffusion effects is also given, because we may have used some new numerical schemes that are different from that in Sherratt (2000). Figure 1 is about the numerical solution to (3.1) subject to free boundary conditions and initial function (3.2) with *ξ*=0.1. We find from figure 1*b* that the solution stabilizes to a wavefront with a speed which can be computed by the formulasee formula (2.5) in Sherratt (2000).

Now we compare the leading-term wavefront with the true solution. Graphs of solutions by the leading-order terms (2.3) and solutions computed from (3.1) are presented (figure 2). When the wave speed is large, say *c*=10, the two solutions match perfectly. Even in the small speed case when *c*=1, the absolute error is within 0.018.

Next we contrast the spread of wavefronts of the nonlinear diffusion system (system (3.1)) and that of the linear system(3.3)with the same initial condition (3.2). Our result is that when the wave speed is large (or the initial decaying rate *ξ* is small), the difference is not apparent because both of them are dominated by the same leading-term contribution *N*_{0} and *V*_{0} (figure 3*a*). But if the initial decaying rate is big, there is a significant difference (figure 3*b*). In figure 3*a*,*b*, the spreading speeds in the linear diffusion case (system (3.3)) are faster than those in the nonlinear diffusion case (system (3.1)). This agrees with our expectation owing to the contact inhibition effect in (3.1).

Finally, we should mention that even though numerical simulation can produce wavefronts with wave speeds that are not large, theoretical proof of the existence of such wavefronts seems extraordinarily difficult and this remains as our future work.

## Acknowledgments

Work by H.Z. was supported by the China Scholarship Council and the Natural Science Foundation of Human Province. Work by C.O. was supported by an NSERC grant of Canada.

## Footnotes

- Received June 20, 2007.
- Accepted January 22, 2008.

- © 2008 The Royal Society