Document Detail

A non-local evolution equation model of cell-cell adhesion in higher dimensional space.
Jump to Full Text
MedLine Citation:
PMID:  23289870     Owner:  NLM     Status:  Publisher    
A model for cell-cell adhesion, based on an equation originally proposed by Armstrong et al. [A continuum approach to modelling cell-cell adhesion, J. Theor. Biol. 243 (2006), pp. 98-113], is considered. The model consists of a nonlinear partial differential equation for the cell density in an N-dimensional infinite domain. It has a non-local flux term which models the component of cell motion attributable to cells having formed bonds with other nearby cells. Using the theory of fractional powers of analytic semigroup generators and working in spaces with bounded uniformly continuous derivatives, the local existence of classical solutions is proved. Positivity and boundedness of solutions is then established, leading to global existence of solutions. Finally, the asymptotic behaviour of solutions about the spatially uniform state is considered. The model is illustrated by simulations that can be applied to in vitro wound closure experiments.
Janet Dyson; Stephen A Gourley; Glenn F Webb
Publication Detail:
Type:  JOURNAL ARTICLE     Date:  2013-1-7
Journal Detail:
Title:  Journal of biological dynamics     Volume:  -     ISSN:  1751-3766     ISO Abbreviation:  J Biol Dyn     Publication Date:  2013 Jan 
Date Detail:
Created Date:  2013-1-7     Completed Date:  -     Revised Date:  -    
Medline Journal Info:
Nlm Unique ID:  101299725     Medline TA:  J Biol Dyn     Country:  -    
Other Details:
Languages:  ENG     Pagination:  -     Citation Subset:  -    
a Mansfield College , University of Oxford , Oxford , OX1 3TF , UK.
Export Citation:
APA/MLA Format     Download EndNote     Download BibTex
MeSH Terms

From MEDLINE®/PubMed®, a database of the U.S. National Library of Medicine

Full Text
Journal Information
Journal ID (nlm-ta): J Biol Dyn
Journal ID (iso-abbrev): J Biol Dyn
Journal ID (publisher-id): tjbd
ISSN: 1751-3758
ISSN: 1751-3766
Publisher: Taylor & Francis
Article Information
Download PDF
© 2013 The Author(s). Published by Taylor & Francis.
Received Day: 4 Month: 5 Year: 2012
Revision Received Day: 28 Month: 11 Year: 2012
Electronic publication date: Day: 7 Month: 1 Year: 2013
Print publication date: Month: 12 Year: 2013
Volume: 7 Issue: Suppl 1
First Page: 68 Last Page: 87
PubMed Id: 23289870
ID: 3957472
DOI: 10.1080/17513758.2012.755572

A non-local evolution equation model of cell–cell adhesion in higher dimensional space
Janet Dysona
Stephen A. Gourleyb*
Glenn F. Webbc
a Mansfield College, University of Oxford, Oxford OX1 3TF, UK
b Department of Mathematics, University of Surrey, Guildford, Surrey GU2 7XH, UK
c Department of Mathematics, Vanderbilt University, Nashville, TN 37240, USA
* Corresponding author. Email:
Author Emails:;

1.  Introduction

Cell invasion is important in many biological processes, particularly in the invasive stage of cancer and in embryonic development [23, 32]. There has recently been some interest in how to mathematically model cell–cell and cell–extracellular matrix adhesions in the context of cancer invasion, since the latter is believed to be characterized by a number of processes that include loss of cell–cell adhesion and enhanced cell–matrix adhesion in addition to active migration, cell proliferation and the secretion of matrix degrading enzymes [15]. Accurate modelling of cell adhesion is, therefore, important, but it is challenging to do so using continuous mathematical models because in such an approach one uses continuous variables for cell densities. Individual cells are not recognized as such and, therefore, there is no representation of cell boundaries.

A modelling approach that aims to address these difficulties was proposed by Armstrong et al. [2]. The ideas have been taken further and applied to various contexts by, for example, [3, 12, 15, 23, 31]. For a cell density p(x, t) the basic model in one spatial dimension is

[Formula ID: M1]

for x ∊ (−∞, ∞), t > 0 with initial condition p(x, 0) = ϕ0(x). Equation (1) incorporates a direct representation of cell-cell contact via a non-local flux term, the integral term in Equation (1). Spatially structured models of tumour invasion that incorporate other aspects such as haptotaxis and age structure were considered in [911]; see also [34, 35].

The purpose of the present paper is to formulate and analyse an extension of Equation (1) to the case of N spatial dimensions, considering only cell-cell adhesion. Again letting p(x, t) denote the cell density at (x, t), but now with x ∊ ℝN, t > 0 and Bρ denoting the N-dimensional ball centred at 0 and of radius ρ, we propose the equation

[Formula ID: M2]

with initial condition

[Formula ID: M3]

The function f(p) models cell loss and gain, and for biological reasons we assume that f(0) = 0 and that there is a number P1 > 0 such that f(p) > 0 for p ∊ (0, P1), and f(p) < 0 for p > P1. These assumptions imply that at lower densities, there is cell gain, but at large enough densities cell loss occurs more rapidly than the generation of new cells via division, due to the effects of crowding.

Equation (2) admits the possibility of two components to cell motion. Motion of cells is assumed to be partly due to Fickian diffusion, which gives rise to the Laplacian term. The need for a Fickian term is debatable, and Painter et al. [23] note from simulations of their more complex cell–matrix model that, in the absence of a Fickian term, and unusually for a continuous model, non-invasive growth is possible for some parameter regimes while invasive growth takes place in others. The other component to cell motion is the movement of cells due to adhesion as modelled by the cell adhesion term in Equation (2). It is an advective flux term and models motion of cells that is caused by mechanical forces attributable to adhesive bonds that cells have formed with their closer neighbours. Cells are able to sense their surroundings over some radius ρ, which is believed to be of the order of several cell diameters [23], due to their ability to deform and extend protrusions [8]. These protrusions cause adhesive bonds to form with other cells, and the resulting forces cause cell movement. In fact, the cell adhesion term in Equation (2) is an advective flux term of the form −▽ · (p(x, t)U(x, t)), where U(x, t) is the velocity of the cells at x ∊ ℝN at time t. If we imagine the cells as tiny spheres moving through a viscous fluid, then Stokes' law leads us to suppose that they are subject to a resistive force which is directly proportional to the velocity. This analogy leads us to suppose that the velocity of a cell is proportional to the net adhesive force on it due to the bonds formed with nearby cells. So, the velocity U(x, t) can effectively also be thought of as the adhesive force on a cell due to bonds with these other nearby cells within the sensing radius ρ (the range over which a cell can detect its surroundings). This force is taken to be of the form

[Formula ID: M4]

We think of expression (4) as a force, but it is really a velocity. Expression (4) is assumed to contain dimensional constants such as the viscosity of the medium (perhaps embodied within the functions g and h), to ensure that it has the units of velocity. Viewing expression (4) as a force, we are asserting that the total force on a cell at position x is the sum, over all cells within the sensing radius ρ, of the local forces attributable to bonds formed with the cells at the nearby points x + ξ, ξBρ. The function h : ℝ → ℝ describes how the magnitude of the force depends on |ξ|, and is a positive function since adhesive forces are always directed towards cell centres. The vector ξ in front of h(|ξ|), which is not present in Equation (1), gives direction to the force on cells at x due to bonds with cells at x + ξ. Another formulation would be to make this direction vector a unit vector and write (ξ/|ξ|) (|ξ|) rather than ξh(|ξ|), as some other authors have done, see for example Painter et al. [23]. There is no real distinction; we are simply taking h(|ξ|) = (|ξ|)/|ξ|.

The magnitude of the adhesive forces from cells at the nearby location x + ξ will depend on the number of adhesive attachments made by cells at x + ξ to a cell at x, and hence on the cell density at x + ξ. This explains the term g(p(x + ξ, t)) in the integrand of expression (4), and the function g(p) describes how the forces depend on the local cell density. We anticipate that g(p) should increase linearly with p if cell density is not too great, because in this situation a doubling of the number of cells at x + ξ should roughly double the number of adhesive attachments made to cells at x. However, at higher densities, particularly as close cell packing is approached, the tendency to form attachments drops off and so g(p) in fact decreases once p has risen above a certain threshold, and is zero for all p above a critical cell density P2 > 0 corresponding to close cell packing. The shift in balance between cell division and cell loss should occur at an achievable cell density. Thus, for biological realism, P1 < P2.

There is unquestionably a need for analytic study of equations of the form (2) in spatial domains of dimension N ≥ 2. Recent interest in this kind of equation has been largely in the areas of cancer cell invasion and wound healing. Gerisch and Chaplain [15] and Painter et al. [23] both consider systems containing a term of the same form as the non-local term in Equation (2), and with the same interpretation – cell velocity due to adhesive bonds with other cells – with an emphasis on two dimensions so that Bρ becomes a disk of radius ρ. These studies have been largely numerical, which is time-consuming in the 2D case. In addition to their numerical simulations, Gerisch and Chaplain [15] examine the possibility of expanding the non-local term, using Taylor series in the case when ρ is small. This analytical approach does yield some useful insight, although the approximated equations can allow unrealistic phenomena such as blow up. Two- and three-dimensional studies are extremely important in the cancer modelling context because of the manner in which cancer cells invade tissue and the possibility of phenomena such as fingering at the invasive front (see Figure 5 in [23]). The shape of a tumour-host boundary is an important diagnostic indicator ([23], and the references therein). Straight and/or sharp boundaries tend to imply non-invasive tumours, whereas if a tumour has a diffuse and/or wavy boundary it is likely to be invasive.

Another important reason for considering equations of the form (2) in higher space dimensions is that the inclusion of another space dimension can be an important test of the stability and robustness of a one-dimensional pattern, and therefore of the plausibility of a type of mathematical model in a situation where many details are unknown. This was an important issue in Armstrong et al. [3] who considered a system with the same type of adhesive flux term as in Equation (2) as a model for somite formation. Segmentation in a number of organisms proceeds through the formation of somites, which will eventually become repeated structures such as the vertebrae. Somitogenesis is, therefore, seen as a one-dimensional process (along the anterior-posterior axis of the presomitic mesoderm) and an important step in the model validation process can be to include a second spatial dimension and investigate whether the model is still capable of forming somites, as was done in [3].

Using fractional powers of the diffusion operator, we establish results about the existence and uniqueness of solutions and their positivity and boundedness in spaces of uniformly continuous functions. We also look at solutions in Lp spaces and show that under suitable conditions the non-zero uniform state is asymptotically stable in L2 and that under more stringent conditions all solutions converge to this uniform state. Finally, we look at simulations which apply to wound closure experiments.

2.  Preliminaries

Initially, we will work in the spaces

In either case, take the norm ||u||k = Σ|Λ|≤k ||DΛu||∞, where for i = 1, …, N, Di : 𝒟(Di)

[Formula ID: M5]

so Di is a closed linear operator. Set D = (D1, …, DN), with domain 𝒟(D) = 𝒟(D1) × … × 𝒟(DN). If Λ = (Λ1, …, ΛN), DΛ = D1Λ1DNΛN. We define |Λ| = Λ1 + … + ΛN and say Λ ≤ μ for N-tuples Λ and μ if the inequality holds component-wise. Note that from [22] page 33, uBUC0(k)(ℝN) implies that, for |Λ| ≤ k, DΛu → 0 as |x| → ∞.

We write k when results apply to either Xk or X0k and set X = X0.

We define à : 𝒟(Ã) ⊂ XX, by

For ω > 0 set A = Ã + ωI, 𝒟(A) = 𝒟(Ã) so that the spectrum of A is contained in the open right half plane.

We rewrite the problem (2)–(3) abstractly as

[Formula ID: M6]

where p(t) : [0, ∞) → X, t ≥ 0, f : ℝ → ℝ, g : ℝ → ℝ, and

[Formula ID: M7]

with ith component Ki(ϕ)(x).

Let T(t) be the analytic semigroup in X generated by −Ã, so

[Formula ID: M8]

[Formula ID: M9]

and T(t) : X0X0. Let Tk(t) be the analytic semigroup which is T(t) restricted to the space k (k an integer); the corresponding generators −Ãk are given by the operator −Â but restricted to different domains [22].

It is easy to see that for ϕXk, k > 0, |Λ| ≤ k

[Formula ID: M10]

so we also have T(t) : XkXk,

[Formula ID: M11]


Set Ak = Ãk + ω. We will use the fractional powers Akα of Ak, and exploit the result that, for suitable α, the operator DiAk−α : XkXk is bounded, see [24]. We will write 𝒟α,k = 𝒟(Akα) in Xk, 𝒟0α,k = 𝒟(Akα) in the space X0k and α,k = 𝒟(Akα) in the space k. We give these spaces the graph norm.

Thus the mild form of (2)–(3) in Xk is

[Formula ID: M12]

[Formula ID: M13]


Observe that for ϕ𝒟(D), gC1(ℝ), hL1(0, ρ), we have K(ϕ) ∊ 𝒟(𝒟),

[Formula ID: M14]


[Formula ID: M15]

Note that in our proofs of existence, we will take f, g : ℝ → ℝ but we will then prove that if ϕ0 ≥ 0 then p(t) ≥ 0. Thus, provided there is enough smoothness at 0 it is only the properties of f and g on [0, ∞) that are relevant as we can then define f and g appropriately for p < 0.

We also have [17, 18, 22, 24]

  • (i)  if α > 0, then for every t > 0 the operator AkαTk(t) is a bounded linear operator in k and for any > 0 there exists C1 > 0 such that
    [Formula ID: M16]
  • (ii)  if 0 ≤ α ≤ 1, then Ak−α : kk is a bounded linear operator so there exists a constant C2 ≥ 1 such that for all u ∊ Xk,
    [Formula ID: M17]
  • (iii)  if ½ < α < 1, then DiAk–α : kk is a bounded linear operator so there exists a constant C3 ≥ 1 such that for all u ∊ Xk,
    [Formula ID: M18]

The following lemma is proved in very much the same way as Lemma 1 in [12]

Lemma 1

Suppose that for some k ≥ 0, fCk+1(ℝ), gCk+2(ℝ) and hL1(0, ρ), ½ < α < 1 and f(0) = 0. Then G : kk and

  • (a)  if ||Ak–αϕ||kR1there exists a constant K1(R1) such that
    [Formula ID: M19]
  • (b)  if ||ϕi||k ≤ R2, i = 1, 2, then there exists a constant K2(R2) such that
    [Formula ID: M20]
  • (c)  if k ≥ 2, and ||ϕ||k–2R3
    [Formula ID: M21]

3.  Existence, positivity and boundedness

We can now prove local existence, and regularity, as in Theorem 1 [12]. The continuous dependence on the initial data follows immediately from [17] Theorem 3.4.1.

Theorem 1

Suppose that for some k ≥ 0, fCk+1(ℝ), gCk+2(ℝ), f(0) = 0, that h ∊ L1(0, ρ), and that ϕ0Q ∊ 𝒟0α,k, where ½ < α < 1, and either Q = 0 or there exists P1 > 0 such that f(P1) = 0 and Q = P1. Then the problem (2) and (3) has a unique mild solution p(t) such that p(t) − Q ∊ C([0, T0]; 𝒟0α,k), for T0 > 0 small enough.

In addition p(t) is the classical solution on [0, T0] of the problem (2) and (3) in the sense that

[Formula ID: M22]

and for 0 < tT0, p(t) ∊ 𝒟(Ak), and

[Formula ID: M23]

Also, under the above conditions, if k ≥ 2m,

[Formula ID: M24]

In particular if k = 2, p(t) − QC1([0, T0]; X0) and Equation (23) also holds at t = 0.

Furthermore, there is continuous dependence on the initial data. Let [0, 0) be the maximal interval of existence of the solution p(t). Now let pn(t) be the solution of Equation (2) with initial data pn(x, 0) = ϕn(x) where ϕnQ ∊ 𝒟0α,k. Suppose that ||Akα(ϕn − ϕ0)||k → 0 as n → ∞. Then, given any t1 ∊ [0, 0), pn(t) is defined on [0, t1] for large enough n and ||Akα(pn(t) − p(t))||k → 0 as n → ∞, uniformly on [0, t1].


To deal with the case Q = P1 we set p(x, t) = w(x, t) + P1 in Equations (2) and (3) to get

[Formula ID: M25]

where f(w) = f(w + P1), and g(w) = g(w + P1). If we set

[Formula ID: M26]

the abstract form of Equation (25) is

[Formula ID: M27]

If G(ϕ) = D · ((A−αϕ + P1)K(A−αϕ)) − f(A−αϕ), then it is easy to see that G satisfies Lemma 1 and we can apply the same methods to the resulting equation in w.

We can now obtain positivity and boundedness.

Proposition 1

Suppose that hC1 [0, ρ], f ∊ C3(ℝ), gC4(ℝ), f(0) = 0, and that ϕ0QD0α,2, where ½ < α < 1 and Q = 0 or Q = P1 > 0 where f(P1) = 0. Let ϕ0(x) ≥ 0 for all x ∊ ℝN. If p(t) is the unique classical solution of Equations (2) and (3) on [0, T0], then p(x, t) ≥ 0 for all 0 ≤ tT0, x ∊ ℝN.

Proof Take M1 > 0 and M2 ≥ 0 such that

[Formula ID: M28]

and, using the mean value theorem,

[Formula ID: M29]

Choose C sufficiently large that

[Formula ID: M30]

where ωN is the surface area of ∂B1, and define u(x, t) = p(x, t) eCt, so

[Formula ID: M31]

We now prove that p(x, t) ≥ 0 on ℝN × [0, T0]. Suppose not, then u(x, t) will be strictly negative somewhere in this set and hence attains a global minimum at some (x0, t0) ∊ ℝN × (0, T0]. Thus for all i

[Formula ID: M32]

In the ith integral in the second term in the right-hand side of Equation (31), we may replace ∂/∂xi by ∂/∂ξi and therefore an integration by parts shows that the sum at the point (x0, t0) equals

The above quantity is bounded in absolute value by

Evaluating Equation (31) at (x0, t0), and making use of (32) and the above bound, we obtain

[Formula ID: M33]

using inequality (29), the fact that u(x0, t0) < 0 and that C satisfies inequality (30). This is a contradiction, hence p(t) ≥ 0.

The following two results are proved using similar arguments to that used above. (See [12] for the proof when N = 1.)

Proposition 2

Suppose that hC1[0, ρ], h ≥ 0, fC3(ℝ),gC4(ℝ), f(0) = 0 and that ϕ0QD0α,2where ½ < α < 1, and Q = 0 or Q = P1 > 0 where f(P1) = 0. Suppose there exists P2 > P1such that g(P) ≥ 0 for P ∊ [0, P2], f(P2) < 0 and that

[Formula ID: M34]

Let 0 ≤ ϕ0(x) ≤ P2for all x ∊ ℝN. If p(t) is the unique classical solution of Equations (2) and (3) on [0, T0], then 0 ≤ p(x, t) ≤ P2for all x ∊ ℝN, 0 ≤ tT0.

Proposition 3

Suppose that h ∊ C1 [0, ρ], h ≥ 0, g bounded, fC3(ℝ), gC4(ℝ), f(0) = 0, and that ϕ0Q ∊ 𝒟0α,2where ½ < α < 1, and Q = 0 or Q = P1 > 0 where f(P1) = 0. Suppose there exists P3 > P1such that for P > P3, f(P) < 0 and

[Formula ID: M35]

Let 0 ≤ ϕ0(x) ≤ P3for all x ∊ ℝN. If p(t) is the unique classical solution of Equations (2) and (3) on [0, T0], then 0 ≤ p(x, t) ≤ P3for all x ∊ ℝN, 0 ≤ tT0.

The following global existence result is proved similarly to Theorem 2 in [12]. The basic idea runs as follows. There will be global existence if ||A2αp(t)||2 is bounded on bounded intervals of time. The method uses Proposition 2 or 3 to show first that p(t) is bounded on bounded intervals of time. Thus by Lemma 1(a) ||G(Aαp(t))|| is bounded linearly, so we can use the Gronwall-type inequality in [17] 1.2.1 to show that ||Aαp(t)|| is suitably bounded. Then we can use Lemma 1(c) to get a linear bound on ||G(A2αp(t))||2, so, using the Gronwall-type inequality again, ||A2αp(t))||2 is also bounded as required

Theorem 2

Suppose that hC1[0, ρ], h ≥ 0, fC3(ℝ), gC4(ℝ). Let ϕ0Q ∊ 𝒟0α,2where ½ < α < 1, and Q = 0 or Q = P1 > 0 where f(P1) = 0 and ϕ0(x) ≥ 0 for all x ∊ ℝN. Suppose that either

  • (a)  There exists P2 > P1such that ϕ0(x) ≤ P2for all x ∊ ℝN, g(P) ≥ 0 for P ∊ [0, P2], f(P2) < 0 and inequality (34) holds. Or
  • (b)  g is bounded and there exists P3 > P1such that ϕ0(x) ≤ P3for all x ∊ ℝN, and for P > P3, f(P) < 0 and inequality (35) holds.

Then the solution p(t) in Theorem 1 is global.

4.  Initial data in Lp(ℝN)

We can also set the problem in Lp(ℝN). For any 1 < p < ∞ define the operator, ÃLp : 𝒟(ÃLp) ⊂ Lp(ℝN) → Lp(ℝ) by

[Formula ID: M36]

Similarly, we can also define ÃL1 : 𝒟(ÃL1) ⊂ L1(ℝN) → L1(ℝN), for a suitable domain 𝒟(ÃL1) [16]. For 1 ≤ p < ∞, ω > 0, set ALp = ÃLp + ωI, 𝒟(ALp) = 𝒟(ÃLp). It is well known [19] that −ÃLp generates a positive analytic semigroup {TLp(t)}t≥0, with ||TLp(t)|| ≤ 1. Thus inequalities (16) and (17) hold also in Lp(ℝN). Further, in the case 1 < p < ∞, from [20] Propositions 4.1.7 and 1.1.4 and Example 1.3.9, if ½ < β < α < 1,

[Formula ID: M37]

Thus, if we define Di : 𝒟(Di) ∊ Lp(ℝN) → Lp(ℝN),

[Formula ID: M38]

then inequality (18) holds in Lp(ℝ). To show that inequality (18) also holds in L1(ℝN), it is sufficient to show that 𝒟(AL1α) ↪ W1,1(ℝN). To do this note first that it follows from [20] Example 1.3.9, Example 1.3.11 and Remark 1.3.7 that if 1 < s < 2 and 0 < β < α < 1, then

Now, by [20] Proposition 4.1.7, [16] Theorem 1.7 and Prop 4.8, and [20] Prop 1.1.4, for 1 < s < 2

(The definition and properties of the Besov spaces Bsp,q, and their relationship to the Sobolev-Slobodeckii spaces Ws,p, can be found in [1, 16, 20].)

As before, from [19], TLp(t) satisfies Equation (8), hence if ϕXkLp(ℝN), TLp(t)ϕ = Tk(t)ϕ. Equally ϕ ∊ 𝒟(ALp) ∩ 𝒟(ALp) implies ALpϕ = Akϕ.

We now have the following local existence theorem in L2(ℝN):

Theorem 3

Suppose that fC1(ℝ), gC2(ℝ), that hL1(0, ρ), f(0) = 0 and that ϕ0Q ∊ 𝒟(AL2α), where ¾ < α < 1, and either Q = 0 or there exists P1 > 0 such that f(P1) = 0 and Q = P1. Then, for N ≤ 3, the problem (2) and (3) has a unique mild solution p(t) such that p(t) − QC([0, T0]; 𝒟(AL2α)), for T0 > 0 small enough.

In addition p(t) is the classical solution on [0, T0] of the problem (2) and (3) in the sense that

[Formula ID: M39]

and for 0 < tT0, p(t) ∊ 𝒟(AL2), and

Also there is continuous dependence on the initial data. Let [0, 0) be the maximal interval of existence of the solution p(t). Now let pn(t) be the solution of Equation (2) with initial data pn(x, 0) = ϕn(x) where ϕnQ ∊ 𝒟(AL2α). Suppose that ||AL2α(ϕnϕ0)||L2 → 0 as n → ∞. Then, given any t1 ∊ [0, 0), pn(t) is defined on [0, t1] for large enough n and ||AL2α(pn(t) − p(t))||L2 → 0 as n → ∞, uniformly on [0, t1].

Finally, if in addition ϕ0Q ∊ 𝒟(AL1α), then also p(t) − QC([0, T0]; 𝒟(AL1α)).


The proof of the first part follows from Theorems 3.3.3 and 3.4.1 in [17] provided we can show that G : L2(ℝN) → L2(ℝN) and G : L2(ℝN) → L2(ℝN) are locally Lipschitz.

If we consider G then it consists of a sum of products. Each product has just one term of the form DiAL2−αϕ. The rest are functions of AL2−αϕ ∊ D(AL2α). But, for ¾ < α < β < 1, and N ≤ 3

(see [1] 7.34) and the local Lipschitzness follows. Similarly for G.

Note further that, also using inequality (18) in L1 (ℝN), G, G : L1 (ℝN) ∩ L2(ℝN) → L1 (ℝN) ∩ L2(ℝN) and are locally Lipschitz continuous in L1(ℝN). For ϕ0 ∊ 𝒟(AL1α) ∩ 𝒟(AL2−α), we set up the iteration: ν0 = AL1α ϕ0,

then the iteration converges in L1(ℝN) and L2(ℝN) and by considering the restrictions to bounded subsets of ℝN, it can be seen that the limits are equal almost everywhere. Similarly for G, and the required invariance follows.

5.  Stability of the uniform state

We now consider the local asymptotic stability of the uniform state p = P1. For simplicity, we will consider the case where f and g are logistic. Each of f and g in Equation (2) can be of the form rp(1 − p/K) with different r and K. After suitable non-dimensionalization, without loss of generality, they can be taken as f(p) = p(1 − p) and g(p) = p(Λ − p), with Λ > 1. The uniform steady state is then p = 1. To examine its stability, set p = w + 1, ϕ0 = w0 + 1 to obtain the following equation for w:

[Formula ID: M40]

Theorem 4

Take N ≤ 3. Suppose that hL2(0, ρ), and that ϕ0 − 1 ∊ 𝒟(AL2α), where ¾ < α < 1. For η ∊ ℝN define

[Formula ID: M41]

Then if k(η) < 0 for all η ∊ ℝN the zero solution of Equation (40) is uniformly asymptotically stable in 𝒟(AL2α). Precisely, if γ > 0 is such that k(η) < −γ for all η ∊ ℝN, then there exist ζ > 0, M ≥ 1 such that if ||AL2αw0||L2 ≤ ζ, then

Moreover, k(η) < 0 for all η ∊ ℝN holds if

[Formula ID: M42]


[Formula ID: M43]

Conversely if there exists η0such that k0) > 0, then the stationary solution is unstable.


We apply the results of [17] Section 5.1. The linearization of Equation (40) is

[Formula ID: M44]

Define the bounded linear operator B : 𝒟(AL2α) → L2(ℝN), such that

and the operator H : L2(ℝN) → L2(ℝN) such that H(ϕ) = G(ϕ) − BAL2−αϕ.

Using the same reasoning as in Theorem 3 it can be seen that if ||ϕ||L2R, then there exist constants K(R) and K(R) such that

[Formula ID: M45]

In the case of stability, we note that this implies that H(ϕ) = o(||ϕ||L2) as ||ϕ||L2 → 0. We have already seen that G : L2(ℝN) → L2(ℝN) is locally Lipschitz continuous. Thus Theorem 5.1.1 of [17] holds so that if the solutions of the linearized problem are uniformly asymptotically stable then so too are those of the nonlinear problem.

The right-hand side of Equation (44) is linear and generates an analytic semigroup R(t), say. So, for all u0L2(ℝN), Equation (44) has a global classical solution u(t) = R(t)u0. First take u0L1(ℝN) ∩ L2(N), so that, as in Theorem 3, u(t) ∊ L1(ℝN) ∩ L2(N), and we may take Fourier transforms. If we denote the Fourier transform of u by



[Formula ID: M46]

Thus, using the Plancherel identity,

as required. Thus by density ||R(t)u0||L2 ≤ ||u0||L2 exp(−γt) for all u0L2(ℝN), and the result follows. The inequality (42) follows from completing the square in k(η), while inequality (43) follows from the fact that |sin(η · x)| ≤ |η · x| ≤ |η||x|.

For instability, we apply Corollary 5.1.6 of [17]. From inequality (45) H(ϕ) = O(||ϕ||2L2) as ||ϕ||L2 → 0, so now we require that there exists Λ in the spectrum of the generator of R(t) such that Re Λ > 0. Suppose not, so that there exists M such that, for all ϕL2(ℝN),

[Formula ID: M47]

But, using the expression for û(η, t) in (46), there exists ∊ > 0 such that

If we now choose u0 such that ∫B|û0(ξ + η0)|2 dξ > 0, then using the Plancherel identity gives a contradiction to inequality (47).

6.  Convergence to the uniform state

The next result follows in a similar fashion to Theorem 3 and [12] Theorem 3.

Theorem 5

Suppose that fC3(ℝ), gC4(ℝ), hC1[0, ρ]. Let ϕ0Q ∊ 𝒟0α,2, where ½ < α < 1, and either Q = 0 or Q = P1 > 0 with f(P1) = 0. Suppose also that the hypotheses of Proposition 1 and of either Proposition 2 or Proposition 3 on f, g, h, and ϕ0hold. Take also ϕ0Q ∊ 𝒟(AL2α). Let p(t) with p(t) − QC([0, ∞); 𝒟0α,2) be the unique classical solution of Equations (2) and (3). Then, for each t ≥ 0, p(t) − Q ∊ 𝒟 (AL2) = W2,2(ℝN). Furthermore p(t) − QC1((0, ∞); L2(ℝN)).

This is enough to justify calculations in the proof of the following theorem.

Theorem 6

Let hC1[0, ρ] be such that h ≥ 0, and let f(p) = p(1 − p) and g(p) = p(Λ − p) with Λ > 1. Let ϕ0 − 1 ∊ 𝒟0α,2 ∩ 𝒟(AL2α(ℝN)), where ½ < α < 1, and let p(x, t) be the solution of the problem (2) and (3).

  • (a)  If 0 ≤ ϕ0(x) ≤ Λ for all x ∊ ℝN, where Λ satisfies
    [Formula ID: M48]

    then 0 ≤ p(x, t) ≤ Λ for all x ∊ ℝN, t ≥ 0.
  • (b)  If η ≤ ϕ0(x) ≤ Λ for all x ∊ ℝN, where 0 < η < Λ, Λ satisfies inequality (48), and
    [Formula ID: M49]

    then η ≤ p(x, t) ≤ Λ for all x ∊ ℝN, t ≥ 0.
  • (c)  If η ≤ ϕ0(x) < Λ for all x ∊ ℝN, where η > 0, Λ satisfies inequalities (48), (49), and
    [Formula ID: M50]

    then p(x, t) → 1 exponentially as t → ∞ in the sense that for all t > 0
    [Formula ID: M51]


First, to prove (a), note that inequality (48) implies inequality (34) with P2 = Λ, f(P) = P(1 − P), g(P) = P(Λ − P), and thus, solutions remain in the closed interval [0, Λ] for all positive times.

Next, to prove (b), suppose that there exists t* > 0 and a corresponding x* with p(x*, t*) = η, ∂p(x*, t*)/∂xi = 0, ∂2p(x*, t*)/∂xi2 ≥ 0, then

Note that 0 ≤ g(p) ≤ Λ2/4 when p ∊ [0, Λ]. Therefore

by inequality (49). Hence p(x, t) ≥ η for t > 0 and (b) is proved.

Last, to prove (c), define w(x, t) by p(x, t) = 1 + w(x, t), so that w(x, t) satisfies Equation (40) and η − 1 ≤ w ≤ Λ − 1. Multiplying Equation (40) by w(x, t) and integrating with respect to x over ℝN, and then integrating by parts on the Laplacian term and the ∂/∂xi term (this is justified by Theorem 5),

[Formula ID: M52]

The integrand of the middle term in the right-hand side will now be estimated by using the inequality A · B ≤ δ|A|2 + |B|2/4δ with A and B identified by the underbraces. This bounds the integral by a sum of two integrals, one of which cancels with the first term on the right-hand side to give

[Formula ID: M53]

Since Λ > 1, it is easy to show that

[Formula ID: M54]

Recalling that 0 ≤ 1 + w(x, t) ≤ Λ, so using inequality (54) and Theorem IV.15 from [4],

[Formula ID: M55]

Therefore inequality (53) becomes

[Formula ID: M56]

Note that −w3(x, t) ≤ (1 − η)w2(x, t), so that inequality (56) becomes

and (c) follows from inequality (50).

7.  Simulations

We illustrate in Figures 14 the results of the previous two sections with simulations that can be applied to in vitro wound closure experiments. In these experiments, cell cultures are scored in a thin line that closes as the cell population proliferates. A review of mathematical models of wound healing is given in [30] and other references are found in [57, 13, 14, 21, 2529, 33, 36]. In previous work [12], we modelled in vitro wound healing experiments for the model (2) with N = 1, and examined the dependence of solutions on the diffusion parameter δ and the sensing radius ρ. In [12] numerical simulations demonstrated that for larger values of the diffusion parameter δ and smaller values of the sensing radius ρ the wound closed completely across the wound opening, but for smaller values of δ and larger values of ρ complex patterns arose across the wound opening with incomplete closing. Either of the alternative conditions (42) or (43) is sufficient for the stability of the uniform steady state p = 1, and note that both of these conditions are largeness conditions on δ and smallness conditions on both h(·) and on the sensing radius ρ. Similarly conditions (48), (49) and (50) are sufficient for convergence to the uniform steady state and conditions (48) and (49) are smallness conditions on h(·) and ρ while (50) is in addition a largeness condition on δ. In our simulations here, we examine the dependence of the wound healing model with respect to the parameter Λ, which is the non-dimensionalized cell density corresponding to close cell packing, in Equation (2) with g(p) = p(Λ − p).

In Figure 1, N = 1, f(p) = p(1.0 − p), g(p) = p(Λ − p), δ = 0.1, η = 0.2, ρ = 1.0, h(x) = arctan(100.0x)/(10.0x arctan 100.0), and ϕ0(x) = 1.0 − 0.8 exp(−(0.1x)10). In this case, Theorem 6 (48) is satisfied if Λ < 18.94, (49) is satisfied if Λ < 1.457, and (50) is satisfied if Λ < 2.823. We take Λ = 19.0, so that none of the conditions (48), (49), (50) of Theorem 6 are satisfied. The simulation in Figure 1 shows that for this value of Λ, the solution does not remain bounded below by η and does not converge to the equilibrium 1.0. On the other hand, if condition (48) were satisfied, then the solution p would be bounded (0 ≤ p ≤ Λ), and Theorem 2 would give global existence of the solution. Indeed, if we take Λ = 18, then condition (48) holds, and we find that the solution behaves qualitatively the same as in Figure 1. This provides an example where the theory gives global existence of the solution, but the solution does not converge to the equilibrium but instead forms a series of peaks.

In Figures 2 and 3, N = 1, f(p) = p(1.0 − p), g(p) = p(Λ − p), δ = 1.0, η = 0.2, h(x) = arctan(100.0x)/(9.0x arctan 100.0), and ϕ0(x) = 1.0 − 0.8 exp(−(0.1x)10). We consider different values of Λ to illustrate the sensitivity of this parameter for the convergence of the solutions to the equilibrium ≡ 1.0. The hypotheses of Theorem 6 require that Λ < 16.94 (48), Λ < 2.683 (49), and Λ < 2.597 (50). For Λ < 2.597 all three conditions are satisfied, and Theorem 6 implies the solution remains in the interval [0.2, Λ] and converges to 1.0 for all initial data such that 0.2 ≤ ϕ0(x) ≤ Λ. On the other hand if one looks at the stability condition from Theorem 4, then the equilibrium ≡ 1.0 is asymptotically stable for Λ < 17.62 with instability if Λ > 17.62. In Figure 2 (a) Λ = 1.5 and there is convergence to 1.0 (the solution appears to converge monotonically to 1.0 from below, never rising above 1.0). In Figure 2(b)–(d) Λ = 17.0, 18.0, 19.0, respectively, and the solutions develop a series of peaks, symmetric to the right and left of 0.0, and rising above and falling below 1.0. Further examination of this wave behaviour is given in Figure 3 for the cases Λ = 17.0 (a), Λ = 18.0 (b), and Λ = 19.0 (c). For Λ = 17, the amplitude of the waves decreases as t increases and the solution appears to converge slowly to 1.0. For Λ = 18 and Λ = 19, the wave propagates with amplitude dependent on Λ (the greater the value of Λ, the greater the amplitude). In general, for the roles of the model parameters, we postulate that cells concentrate in interior and exterior regions of the wound in a regular pattern depending on the strength of adhesion Λ, the sensing radius ρ, and the diffusion coefficient δ.

In Figure 4, we show that the convergence to the equilibrium ≡ 1.0 corresponding to the healed wound depends on the initial conditions. In these simulations, the parameters are Λ = 2.0, δ = 1.0, ρ = 3.0, f(p) = p(1.0 − p), g(p) = p(Λ − p), and h(x) = 9.0 arctan(100.0x)/(4.0x arctan 100.0). In this figure, the only change in the simulations is the initial data – the parameters are the same for all. For initial conditions corresponding to a shallow, wide wound, and to a narrow, deep wound, closure occurs, but for an initial condition corresponding to a wound with the same width as the first and the same depth as the second, closure does not occur. A wider, shallow wound also closes. In general, the behaviour of the solutions is highly sensitive to initial conditions. We note that the function k(η) defined in Equation (41) is negative for all η ∊ ℝ, which by Theorem 4 implies that the equilibrium ≡ 1.0 is stable in this example.

We have chosen to illustrate the behaviour of solutions for the model using the simple case of an equilibrium corresponding to wound healing experiments, in which wound closure corresponds to convergence to a constant function in the x-direction, perpendicular to the direction of the scoring. Numerical simulations corresponding to higher dimensional cases will be conducted in future work.

[1]. Adams R.A.,Fournier J.J.F.. Sobolev SpacesYear: 20032nd ed.AmsterdamAcademic Press
[2]. Armstrong N.J.,Painter K.J.,Sherratt J.A.. A continuum approach to modelling cell–cell adhesionJ. Theoret. Biol.Year: 20062439811316860344
[3]. Armstrong N.J.,Painter K.J.,Sherratt J.A.. Adding adhesion to a chemical signalling model for somite formationBull. Math. Biol.Year: 20097112418766407
[4]. Brezis H.. Analyse fonctionnelleYear: 1983ParisMasson
[5]. Byrne H.M.,Chaplain M.A.J.,Evans D.L.,Hopkinson I.. Mathematical modelling of angiogenesis in wound healing: Comparison of theory and experimentJ. Theoret. Med.Year: 20002175197
[6]. Chen X.,Friedman A.. A free boundary problem arising in a model of wound healingSIAM J. Math. Anal.Year: 2000324778800
[7]. Diegelmann R.F.,Evans M.C.. Wound healing: An overview of acute, fibrotic and delayed healingFront. Biosci.Year: 2004928328914766366
[8]. Doherty G.J.,McMahon H.T.. Mediation, modulation and consequences of membrane-cytoskeleton interactionsAnn. Rev. Biophys.Year: 200837659518573073
[9]. Dyson J.,Sánchez E.,Villella-Bressan R.,Webb G.. An age and spatially structured model of tumor invasion with haptotaxisDiscrete Contin. Dyn. Syst. BYear: 200784560
[10]. Dyson J.,Villella-Bressan R.,Webb G.. A spatially structured model of tumor growth with cell age, cell size and mutation of cell phenotypesMath. Model. Nat. Phenom.Year: 20072369100
[11]. Dyson J.,Villella-Bressan R.,Webb G.. An age and spatially structured model of tumor invasion with haptotaxis IIMath. Popul. Stud.Year: 2008157395
[12]. Dyson J.,Gourley S.A.,Villella-Bressan R.,Webb G.. Existence and asymptotic properties of solutions of a non-local evolution equation modelling cell-cell adhesionSIAM J. Math. Anal.Year: 201042417841804
[13]. Friedman A.,Xue C.. A mathematical model for chronic woundsMath. Biosci. Eng.Year: 2011825326121631128
[14]. Friedman A.,Hu B.,Xue C.. Analysis of a mathematical model of ischemic cutaneous woundsSIAM J. Math. Anal.Year: 20104220132040
[15]. Gerisch A.,Chaplain M.A.J.. Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesionJ. Theoret. Biol.Year: 200825068470418068728
[16]. Guidetti D.. On elliptic systems in L1Osaka J. Math.Year: 199330397429
[17]. Henry D.. Geometric Theory of Semilinear Parabolic EquationsLecture Notes in Mathematics,Year: 1981840New YorkSpringer-Verlag
[18]. Lieb E.,Loss M.. AnalysisGraduate Studies in MathematicsYear: 2001142nd ed.Providence, RIAmerican Mathematical Society
[19]. Lorenzi L.,Lunardi A.,Metafune G.,Pallara D.. Analytic semigroups and reaction diffusion problemsYear: 2005 Internet Seminar 2004–2005, February 16, Available at∼lunardi/LectureNotes.html..
[20]. Lunardi A.. An introduction to interpolation theoryYear: 2007 Internet Seminar, February Available at∼lunardi/LectureNotes.html..
[21]. Menke N.B.,Ward K.R.,Witten T.M.,Bonchev D.G.,Diegelmann R.F.. Impaired wound healingClin. Dermatol.Year: 200725192517276197
[22]. Mora X.. Semilinear parabolic problems define semiflows on Ck spacesTrans. Amer. Math. Soc.Year: 19832782155
[23]. Painter K.J.,Armstrong N.J.,Sherratt J.A.. The impact of adhesion on cellular invasion processes in cancer and developmentJ. Theoret. Biol.Year: 20102641057106720346958
[24]. Pazy A.. Semigroups of Linear Operators and Applications to Partial Differential EquationsYear: 1983BerlinSpringer-Verlag
[25]. Pettet G.,Chaplain M.A.J.,McElwain D.L.S.,Byrne H.M.. On the role of angiogenesis in wound healingProc. R. Soc. Lond. Ser. BYear: 199626314871493
[26]. Pettet G.J.,Byrne H.M.,McElwain D.L.S.,Norbury J.. A model of wound-healing angiogenesis in soft tissueMath. Biosci.Year: 199613635638755336
[27]. Ross J.S.,Fletcher J.A.,Linette G.P.,Stec J.,Clark E.,Ayers M.,Symmans W.F.,Pusztal L.,Bloom K.J.. The Her-2/neu gene and protein in breast cancer 2003: Biomarker and target of therapyOncologistYear: 20038430732512897328
[28]. Roy S.,Biswas S.,Khanna S.,Gordillo G.,Bergdall V.,Green J.,Marsh C.B.,Gould L.J.,Sen C.K.. Characterization of a preclinical model of chronic ischemic woundPhysiol. GenomicsYear: 20093721122419293328
[29]. Schugart R.C.,Friedman A.,Zhao R.,Sen C.K.. Wound angiogenesis as a function of tissue oxygen tension: A mathematical modelPNASYear: 20081052628263318272493
[30]. Sherratt J.A.,Dallon J.C.. Theoretical models of wound healing: Past successes and future challengesC. R. Biol.Year: 2002325555756412187641
[31]. Sherratt J.A.,Gourley S.A.,Armstrong N.J.,Painter K.J.. Boundedness of solutions of a nonlocal reaction-diffusion model for adhesion in cell aggregation and cancer invasionEur. J. Appl. Math.Year: 200920123144
[32]. Turner S.,Sherratt J.A.. Intercellular adhesion and cancer invasion: A discrete simulation using the extended Potts modelJ. Theoret. Biol.Year: 20022168510012076130
[33]. Wang S.E.,Shin I.,Wu F.Y.,Friedman D.B.,Arteaga C.L.. HER2/Neu (ErbB2) signaling to Rac1-Pak1 is temporally and spatially modulated by transforming growth factor βCancer Res.Year: 200666199591960017018616
[34]. Webb G.F.. Theory of Nonlinear Age-Dependent Population DynamicsMonographs and Textbooks in Pure and Applied MathematicsYear: 1985New YorkMarcel Dekker
[35]. Webb G.F.. Magal P.,Ruan S.Population models structured by age, size and spatial positionStructured Population Models in Biology and EpidemiologyLecture Notes in Mathematics, No. 1936, Mathematical Biosciences SeriesYear: 2008BerlinSpringer149
[36]. Xue C.,Friedman A.,Sen C.K.. A mathematical model of ischemic cutaneous woundsPNASYear: 2009106167821678719805373


[Figure ID: F1]
Figure 1. 

Simulation with ρ = 1.0, Λ = 19.0, δ = 0.1, η = 0.2, f(p) = p(1.0 − p), g(p) = p(Λ − p), h(x) = arctan (100.0x)/(10.0x arctan 100.0) and initial data ϕ0(x) = 1.0 − 0.8 exp(−(0.1x)10). The red plane and lines correspond to η = 0.2 and the green plane and lines correspond to 0.0. (a) Graph of p(x, t), 0 ≤ t < 2.0. The wound does not close, and instead the solution exhibits an evolving pattern of peaks and valleys which do not converge to 1.0. (b) Graphs of p(x, .0), p(x, .5), p(x, 1.0), p(x, 2.0). The solution stays above 0.0, but not above η = 0.2.

[Figure ID: F2]
Figure 2. 

Simulation with ρ = 1.0, Λ = 1.5 (a), Λ = 17.0 (b), Λ = 18.0 (c), Λ = 19.0 (d), δ = 1.0, η = 0.2, f(p) = p(1.0 − p), g(p) = p(Λ − p), h(x) = arctan(100.0x)/(9.0xarctan 100.0), and initial data ϕ0(x) = 1.0 − 0.8 exp(−(0.1x)10). For the higher values of Λ, complex patterns of oscillations develop across the wound.

[Figure ID: F3]
Figure 3. 

The wound healing simulation wave patterns as in Figure 2 with adhesion strength parameters (a) Λ = 17.0, (b) Λ = 18.0, and (c) Λ = 19.0. For Λ = 17.0, the solution is shown at times t = 100.0 (blue), t = 95.0 (red) and t = 90.0 (green). The amplitude envelope of the oscillation packet decreases to a minute value for the times shown. For Λ = 18.0 the solution is shown at times t = 200.0 (blue), t = 150.0 (red) and t = 100.0 (green). The envelope increases minutely for the times shown. For Λ = 19.0, the solution is shown at times t = 150.0 (blue), t = 125.0 (red), t = 100.0 (green), t = 75.0 (yellow) and t = 50.0 (cyan). The envelope widens uniformly for the times shown. The simulations indicate convergence to the uniform equilibrium ≡ 1.0 in (a) Λ < 17.62, but not in (b) and (c), Λ > 17.62.

[Figure ID: F4]
Figure 4. 

Wound healing simulations with parameters Λ = 2.0, δ = 1.0, ρ = 3.0, f(p) = p(1.0 − p), g(p) = p(Λ − p), h(x) = 9.0 arctan(100.0x)/(4.0x arctan 100.0), and initial data ϕ0(x) = 1.0 − 0.4 exp(−(0.1x)10) (top left), ϕ0(x) = 1.0 − 0.95 exp(−(3x)12) (top right), ϕ0(x) = 1.0 − 0.5 exp(−(0.2x)12) (bottom left) and ϕ0(x) = 1.0 − 0.95 exp(−(0.1x)10) (bottom right). The wounds top left and top right both close but the wound bottom right which has the same width as the first and the same depth as the second fails to close. The widest wound, bottom left, also closes.

Article Categories:
  • Research Article

Keywords: cell adhesion, non-local, advection, existence, boundedness, positivity.

Previous Document:  An Autologous Platelet Rich Plasma Stimulates Periodontal Ligament Regeneration.
Next Document:  Inhibition of COX-2 expression by topical diclofenac enhanced radiation sensitivity via enhancement ...