Mathematical modelling, analysis and numerical simulations for the influence of heat shock proteins on tumour invasion

J. Math. Anal. Appl. 408 (2013) 597–614

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications journal homepage: www.elsevier.com/locate/jmaa

Mathematical modelling, analysis and numerical simulations for the influence of heat shock proteins on tumour invasion Gülnihal Meral a,∗ , Christina Surulescu b a

Department of Mathematics, Faculty of Arts and Sciences, Bülent Ecevit University, 67100, Zonguldak, Turkey

b

TU Kaiserslautern, Felix-Klein-Zentrum für Mathematik, Paul-Ehrlich-Str. 31, 67663, Kaiserslautern, Germany

article

info

Article history: Received 10 August 2012 Available online 19 June 2013 Submitted by Juan J. Nieto Keywords: Tumour invasion Heat shock proteins Multiscale model Time delay Reaction–diffusion equations

abstract Invasion is a key property of tumour cells for building metastasis. A class of functionally related proteins called heat shock proteins (HSPs), the expression of which is increased when cells are exposed to elevated temperatures or other stress, play an important role in the invasion. In this paper we develop a mathematical model focusing on the effect of HSPs on the tumour cell migration. The resulting multiscale setting accounts both for the microscopic, intracellular level on which these proteins are acting and for the macroscopic level of the cell population. We prove the local existence of a unique positive solution and perform numerical simulations in order to illustrate the behaviour of cancer cells w.r.t. the fibre density of the tissue and the matrix degrading enzymes, along with the effect of HSPs. © 2013 Elsevier Inc. All rights reserved.

1. Introduction Cancer cells feature the ability to invade the surrounding tissue and then disseminate to different parts of the body. The turning point of cancer is the establishment of metastasis, which is defined as the formation of secondary tumours at distal sites. Thereby, the first step is the invasion and penetration into adjacent tissues. The directed movement of cancer cells may be influenced by two mechanisms: chemotaxis and haptotaxis. The former defines the cell motion in response to the concentration of a substance called chemoattractant (or chemorepellent). However, such gradients may lack in the solution; instead, adhesive molecules could be present in increasing amount along the extracellular matrix (ECM). Since cells have to adhere to the ECM fibres in order to be able to move, they will migrate from a region of low concentration of those adhesive molecules to an area with a higher concentration. This is what is called haptotaxis. Cell locomotion is most often influenced by both mechanisms. Thereby, the contact with the surrounding tissue stimulates the production of proteolytic enzymes, which degrade the tissue fibres, hence creating interstices to be occupied during the migration process towards neighbouring blood vessels. The multiscale character of emergence and development of a tumour becomes evident if one considers the mutual interaction between the cells and the tissue, but also the intracellular dynamics. The interesting spatial scales range from the subcellular level to the macroscopic one (tissue and cell populations), while the time scales stretch from seconds (or even shorter) for processes inside the cell up to months for the doubling times of tumours. Tumour cells migrate through the ECM upon using different strategies. One can classify their migration in two types: amoeboid and mesenchymal. In the former case the cells exert a minimal action on the ECM; rather, they drastically modify their form by reorganizing their cytoskeleton and ‘squeeze’ between the fibres of the tissue. On the contrary, the mesenchymal motion causes an active reshaping of the ECM. These types of motion can also alternate during cell migration and are

∗

Corresponding author. E-mail addresses: [email protected], [email protected] (G. Meral), [email protected] (C. Surulescu).

0022-247X/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jmaa.2013.06.017

598

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

influenced among other factors by the amount of heat shock proteins they produce. Heat shock proteins (HSPs), also called stress proteins, are present in all cells in all life forms also under perfectly normal conditions. Their production is enhanced when a cell undergoes various types of environmental stresses like hyper- or hypothermy, or oxygen deprivation [22]. HSPs act like chaperons, i.e. they protect proteins in the process of folding and are also known to be involved in the production of a very important group of proteases called matrix degrading enzymes (MDEs) which are secreted by the cancer cells to degrade the ECM [13,19]. Furthermore the HSPs also seem to influence the restructuring of the cell’s cytoskeleton, e.g., upon regulating actin polymerization [12,24]. Tumour invasion and metastasis have been studied by many authors. Thereby, one can distinguish two main approaches, one of them being the so-called kinetic approach, where a mesoscopic model in the form of an integro-differential PDE is considered for the evolution of the cell density, possibly coupled with integro-differential and/or reaction–diffusion equations for the ECM fibre density and the chemotactic signal. Then with an appropriate scaling the macroscopic limit is deduced, leading usually to a Keller–Segel model or some hyperbolic systems; see e.g., [6,10] and the references therein. Some models for cancer cell migration through tissue networks have been proposed and analysed in [14] and more recently in [15], respectively [16]. The latter is a multiscale setting accommodating the dynamics of integrin binding on the cell surface and aligns to the very general kinetic theory of active particles (KTAP) proposed by Bellomo et al.; see e.g., [4]. The other main approach is to model directly at the macroscopic level, which leads to systems of reaction–diffusion (transport) equations (RD(T)Es), possibly coupled with ordinary differential equations modelling processes inside or on the surface of a cell. We refer the reader e.g., to [2,5,11,26]. There are only few references concerning mathematical models for HSP, mainly focusing on the intracellular dynamics for the HSP concentration in response to external stimuli [8,22,23,27]. Recently a PDE model has been proposed [26] to account for the effect of HSP90 on tumour invasion. Thereby, a system of RDTEs in one dimension has been written for the dynamics of cancer cells, of the ECM density and for the MDE concentration, upon relying on a model by Chaplain and Lolas [5]. The effect of HSP90 has been included in the transport and diffusion terms via sensitivity functions having a simple description with the aid of one or two ODEs. The authors also accounted for (distributed) time delay for the process of cytoskeleton reorganization. However, the models assume either the regulating effect of HSP on the production of MDEs or on the cytoskeleton structure (corresponding to mesenchymal, respectively amoeboid invasion), whereas their joint effect seems to be crucial for cell migration. In this paper we propose a mathematical model for the influence of the HSPs on tumour invasion accounting both for the mesenchymal and amoeboid types of invasion. The proposed model basically consists of a system of reaction– diffusion–transport equations (RDTE) for the density of cancer cells, for the ECM and the MDE concentration. This basic system is then modified in order to include the effects of HSP. Since both the production of proteolytic enzymes (MDEs) and the reorganization of the cytoskeleton involve several biochemical reactions, which need time to happen, a certain time lag is introduced in the model. This leads to a delay differential equation to be coupled with the previous system of RDTE. Moreover, our model is more flexible, since it involves nonconstant diffusion and haptotaxis coefficients and also allows for proliferation of cancer cells and re-establishment of the ECM fibres. All these considerations result into a more realistic setting, but at the same time lead to a more challenging system to analyse and simulate. We prove the existence of a unique weak solution to this system using an iterative procedure employed in different contexts in [17,21,25]. We also perform the corresponding numerical simulations in the one-dimensional case upon using a finite difference method. The difficulty arising from the strongly nonlinear diffusion and haptotaxis coefficients is overcome using a nonstandard finite difference scheme [7]. Our numerical results confirm the expected behaviour of the invasion including the effects of time delay. A comparison of the effects of several types of delay is carried out as well. 2. The mathematical model In this section a PDE model is derived for the cancer cell density c (x, t ), the extracellular matrix density v(x, t ) and the matrix degrading enzyme concentration m(x, t ) involved in the invasion process. 2.1. Governing equations (a) Cancer cells There are two important factors determining the change in cancer cell density due to dispersion: the random motility Jrandom and the directional flow Jcancer . The former characterizes the cell diffusion into the tissue and is given by Jrandom = −Dc (c , v)∇ c with the random motility function Dc (c , v). Since the spread of cancer cells is conditioned by their neighbours (contact inhibition) and surroundings (here the ECM) the diffusion coefficient (random motility function) should reflect these dependencies. On the other hand, Jcancer corresponds to the cancer cell flux due to spatial gradients of stimulating chemotactic and haptotactic responses: Jcancer = Jchemotaxis + Jhaptotaxis = χC (m)c ∇ m + χH (v)c ∇v

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

599

where χC and χH are the chemotactic and haptotactic functions, respectively. In a realistic setup the former should also depend on the matrix degrading enzyme concentration m and the latter on the extracellular matrix concentration v . However, for simplicity the chemotactic term is often ignored, which will be done in this work, too. Due to mass conservation it follows that

∂c + ∇ · (Jrandom + Jcancer ) = Pc , ∂t

(1)

where Pc is the proliferation term. A natural choice is a logistic growth slowed down by the competition for space between the cancer cells and ECM fibres. Thus the resulting equation for the cancer cell takes the form

  v ∂c c −α , = ∇ · (Dc (c , v)∇ c ) − ∇ · (χH (v)c ∇v) + µc c 1 − ∂t Kc Kv

(2)

where µc is the proliferation rate of the cells, Kc and Kv are the carrying capacities for the cancer cell and the extracellular matrix, respectively. Thereby, 0 < α < 1 is a parameter characterizing growth reduction of the cells due to the competition (e.g., for space) with the tissue fibres (a similar choice has been made in [11] in a related modelling framework). In our model we consider two different random motility functions Dc characterizing the diffusion of cells into the tissue to be one of the forms (D denotes some positive constant).

• Dc (c , v) = D(1 − γ Kcc vKv ), γ < 2 a positive constant. • Dc (c , v) = 1+θD c v , θ a positive constant. Kc Kv

These choices account for the joint effect of the cancer cells and of the ECM fibres on the diffusion. Moreover the term c v illustrates their interaction yielding motion hampering, due to the cells adhering to ECM fibres, then – according to their current polarization – releasing the existing focal adhesions and re-adhering at some other sites on their surface membrane etc. The scaling with the respective carrying capacities of cells and fibres takes into account the limited interaction possibilities, while the choice γ < 2 prevents the diffusion coefficient from becoming degenerate or even negative.1 Another way to avoid such problems is to consider a random motility coefficient like our second choice above, which cannot degenerate or become negative but still accounts for the effects of cell–fibre interactions on reducing the diffusivity. In previous models the diffusion coefficient Dc is taken constant, whereas the haptotaxis coefficient χH is chosen to be K K2

either constant [5,26] or of the form χH (v) = (K K 1+Kv v)2 , with Ki (i = 1, 2, 3) and Kv positive constants [5]. The second 2 v 3 choice involves the sensitivity of cancer cells to the surrounding tissue, however it gives a very small weight to haptotaxis. Instead, we rely on the Michaelis–Menten kinetics and take

• χH (v) =

K1 v , Kv +v

allowing for enhanced haptotactic effects. Thereby, K1 is as above some positive constant. (b) Extracellular matrix The MDEs degrade the ECM upon contact and the ECM components reestablish and remodel themselves while competing with the diffusive cancer cells for space. This is again described by way of a logistic growth. Further, the ECM does not diffuse, but can be only degraded by the cells with the aid of MDEs:

  c v ∂v = −δv mv + µv v 1 − − ∂t Kc Kv

(3)

where δv is the degradation rate and µv describes the ECM production. The restructuring of tissue damaged by proteolytic activity can be assumed e.g., for a slowly evolving cancer. For instance, tissue inhibitors of metalloproteinases (TIMPS) are involved in the remodelling of connective tissue [1,18] in response to degradation by matrix metalloproteinases. Like in the case of cancer cell proliferation we assume logistic growth along with crowding effects of cells competing with fibres for space. (c) Matrix degrading enzymes Matrix degrading enzymes are produced by the tumour cells while in contact with the ECM fibres; they diffuse throughout the ECM and undergo some natural decay:

∂m = Dm ∆m − δm m + µm c v, ∂t where Dm is the constant diffusion coefficient and δm , µm are the degradation and production rates, respectively.

1 In view of the properties for the solution approximation proved in Lemma 3.1.

(4)

600

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

Summarizing, the complete system of partial differential equations is obtained in the form

  ∂c v c −α = ∇ · (Dc (c , v)∇ c ) − ∇ · (χH (v)c ∇v) + µc c 1 − ∂t Kc Kv   ∂v c v = −δv mv + µv v 1 − − ∂t Kc Kv ∂m = Dm ∆m − δm m + µm c v. ∂t

(5)

2.2. Effect of HSP in the model Heat shock proteins have a regulating effect on the production of MDEs and on the cytoskeleton structure, which can be put in correspondence to mesenchymal and amoeboid invasions, respectively. In [26], each of these effects is taken into account separately, however both are crucial for cell migration and we will jointly consider them in our model. The concentration level of HSP influences the ability of cancer cells to secrete MDEs in order to degrade components of the surrounding tissue [13,19] and is therefore expected to modulate the invasion. In the model we account for these effects upon modifying the degradation term for the ECM concentration and the production term in the equation for MDEs. Then system (5) takes the following form:

  c v ∂c = ∇ · (Dc (c , v)∇ c ) − ∇ · (χH (v)c ∇v) + µc c 1 − −α ∂t Kc Kv   ∂v c v = −δv mv h + µv v 1 − − (6) ∂t Kc Kv ∂m = Dm ∆m − δm m + µm c v h, ∂t where h(t ) is a given time dependent and smooth enough function of the HSP concentration. Here and everywhere in the following we use the HSP70 function, the shape of which has been approximated by a Gamma distribution (see Fig. 1) to reproduce the trend of empirical HSP70 data from [23,27]. We did not employ the HSP90 profile in [26], since we could not find any databased evidence for that particular shape and moreover that function was not explicitly specified there. Furthermore, the level of HSP concentration also affects the cell ability to reorganize its cytoskeleton [24], which in turn can influence the cell’s diffusion and haptotaxis. Hence both the random motility function Dc and the haptotactic function χH should depend on the HSP concentration and we propose that Dc (c , v, I ) = Dc (c , v)I (h),

χH (v, I ) = χH (v)I (h),

(7)

where I (h) describes the functional relationship between the HSP concentration and the flexibility of the cell. Since the reorganization of the cytoskeleton involves several biochemical reactions, which need time to happen, the cell’s flexibility depends on some time delay. A reasonable description of time lag effects in biological systems is to assume that past influences the present state over a certain time interval. The most popular way to account for a delay is to assume a constant one, denoted here by τ and describing the time lag for the production of the HSP function h, which modifies the production term for MDEs in (6) into µm c v h(t − τ ). Moreover, the change in flexibility leads to an equation of the form dI

= −qI + h(t − τ ), q > 0 (8) dt to be coupled with the system (6). As a second attempt we use the classical distributed delay (see, e.g. [20]), a special case (q = 1) of which has also been proposed for this particular problem in [26]: I (h(t )) =

∞



e−qτ h(t − τ )dτ

(9)

0

with the rate parameter q being a positive constant. Next the usual trick is to differentiate the integral expression (9) to obtain an ODE for the cell flexibility: dI dt

= −qI + h(t ),

q > 0.

(10)

2.3. Boundary and initial conditions Since the invasion takes place within an isolated system, according to the in vitro experimental protocol, we assume that there is no flux of tumour cells and MDEs across the (smooth enough) boundary of the domain Ω ⊂ Rn . Thus, the boundary

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

601

10

9

8

7

6

5

4

3

2

1

0

10

20

30

40 t

50

60

70

80

Fig. 1. Concentration of HSP70.

conditions are

∇ c · n = 0, at (0, T ) × ∂ Ω ∇ m · n = 0 at (0, T ) × ∂ Ω ,

(11)

where n denotes the outward unit normal vector to ∂ Ω . The initial conditions are given by c (x, 0) = c0 ,

v(x, 0) = v0 ,

m(x, 0) = m0

(12)

for x ∈ Ω . Here c0 , v0 and m0 are strictly positive functions; c0 and m0 satisfy the boundary conditions (11). 3. Local existence and uniqueness of solutions In this section, we prove the existence and uniqueness of a weak solution to system (5) with boundary and initial conditions given by (11) and (12). Thereby, we use the random motility function Dc (c , v) = D(1 − γ Kc vK ), which is c v

c Kv , since the former one can become negative for certain values of the interaction more challenging than Dc (c , v) = K DK c Kv +θ c v parameter γ , leading to nonphysical oscillations and possibly blowup (see Section 5). In the mathematical analysis below this is avoided by a proper selection of γ and of the initial data. However we use the latter diffusion coefficient throughout the numerical simulations in order to elude the mentioned problems with parameter values leading to oscillations.2 The existence proof relies on an iterative procedure.3 Consider the function spaces:

X := L∞ (0, T ; H 1 (Ω )), Y := u ∈ L2 (0, T ; H 2 (Ω )) : ut ∈ L2 (0, T ; L2 (Ω )) ,





(13)

Z := L∞ (0, T ; L2 (Ω )), W := L∞ (0, T ] × (Ω ).

Definition 3.1. A weak solution of the system (6) with the boundary conditions (11) and initial data (12) is a triple (m, v, c ) of functions in X × Y × Z , such that for all φ ∈ H 1 (Ω ) the following three equations are satisfied:

Ω

∂m φ dx + ∂t

Ω

∂v φ dx + h(t ) ∂t

 



 Ω

Dm ∇ m∇φ dx +

 Ω

δv mvφ dx =

Ω

δm mφ dx = h(t )

 Ω

 Ω

µm c vφ dx a.e. in [0, T ]

  c v µv v 1 − − φ dx Kc

Kv

2 For more comments on this we refer the reader to Section 5. 3 The same proof also goes through (and is actually even easier to do) in the case where the alternative diffusion coefficient is chosen.

(14) (15)

602

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

 Ω

∂c φ dx + ∂t

 Ω

Dc (c , v, I )∇ c ∇φ dx −

 Ω

χH (v, I )c ∇v∇φ dx =

 Ω

  c v µc c 1 − −α φ dx. Kc

Kv

(16)

Theorem 3.1. There exists T > 0 such that the system (6) with the boundary conditions (11) and initial data (12) satisfying

v0 ∈ W 2,∞ (Ω ),

m0 ∈ L∞ (Ω ) ∩ H 1 (Ω ) ∩ C (Ω ),

c0 ∈ H 1 (Ω )

(17)

and m0 ≥ Cm > 0,

0 < v0 ≤

Kv 2

,

0 < c0 ≤

Kc 2

,

(18)

with Cm some positive constant has a unique solution (m, c ) ∈ (X × X ) ∩ (Y × Y ) and v ∈ Z . Proof. The proof will consist of three steps. In the first one, we will construct a sequence with appropriate properties. Then we will show that this is a Cauchy sequence. Due to the completeness of the spaces we work on, this will lead to the convergence of the sequence to some limit function which is the solution to system (6). Finally, we will show that this solution is unique. a. Construction of the iterative sequence Let (m0 , c 0 ) ∈ (X × X ) ∩ (Y × Y ) and v 0 ∈ Z be the weak solutions to

∂ m0 − Dm ∆m0 + δm m0 = 0 ∂t ∂v 0 + δv m0 v 0 h(t ) = 0 ∂t       ∂ c0 v0 K1 v 0 0 0 −∇ · D 1−γ I (h)∇ c 0 + ∇ · I ( h ) c ∇v = 0, ∂t Kv Kv + v 0

(19) (20) (21)

and for k ∈ N0 let (mk , c k ) ∈ (X × X ) ∩ (Y × Y ) and v k ∈ Z be the weak solution to

∂ mk+1 − Dm ∆mk+1 + δm mk+1 = µm c k v k h(t ) ∂t   ∂v k+1 ck vk + δv mk+1 v k+1 h(t ) = µv v k 1 − − ∂t Kc Kv    k k+1       k+1 ∂c c v K1 v k+1 ck vk k+1 k+1 k +1 k −∇ D 1−γ I (h)∇ c +∇ I (h)c ∇v = µc c 1 − −α ∂t Kc Kv Kv + v k+1 Kc Kv

(22) (23)

(24)

with the corresponding initial and boundary conditions (11) and (12). Lemma 3.1 (Properties of the Iterative Sequence). Under assumptions (17) and (18) there exists T > 0 such that: (i) for every k ∈ N0 , there exists a unique weak solution to the system (19)–(21) and (22)–(24) with conditions (11) and (12) such that mk , c k ∈ L2 (0, T ; H 2 (Ω )) ∩ L∞ (0, T ; H 1 (Ω ))

,

(25)

∈ L (0, T ; L (Ω ))

(26)

v , ∇v , v ∈ L (0, T ] × (Ω ),

(27)

mkt k

ctk

2

k

k t

2

∞

whereby the subscript ‘t’ is used for the time derivative. (ii) Furthermore, for all k ∈ N0 , the functions mk , v k , c k are positive and they satisfy the following inequalities: mk (x, t ) ≥ Cm e−δm t ,

v k (x, t ) ≤

Kv 2

,

c k (x, t ) ≤ Kc

for a.e. x ∈ Ω , t ∈ [0, T ]

(28)

and the estimates

∥mk ∥X + ∥mk ∥L2 (0,T ;H 2 (Ω )) ≤ C (Ω , T )(∥c0 ∥H 1 (Ω ) + ∥m0 ∥H 1 (Ω ) )

(29)

∥v ∥X ≤ C (Ω , T )∥v0 ∥H 1 (Ω )

(30)

∥c k ∥X + ∥c k ∥L2 (0,T ;H 2 (Ω )) ≤ 2C (Ω , T )∥c0 ∥H 1 (Ω ) ,

(31)

k

with C (Ω , T ) being a generic constant.

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

603

Proof of Lemma 3.1. The proof relies on mathematical induction. (i) Induction start: k = 0

• With the substitution m0 (x, t ) = M 0 (x, t )e−δm t (19) becomes the homogeneous heat equation and using the assumption (17) that m0 ∈ H 1 (Ω ), this leads to the weak solution (see e.g. Chapter 7 in [9]) to m0 ∈ L2 (0, T ; H 2 (Ω )) ∩ L∞ (0, T ; H 1 (Ω )) m0t ∈ L2 (0, T ; L2 (Ω )) with

∥m0 ∥X + ∥m0 ∥L2 (0,T ;H 2 (Ω )) ≤ C (Ω , T )∥m0 ∥H 1 (Ω ) . Moreover, from the expression of the exact solution of (19) written in terms of the initial condition m0 and the heat kernel, it obviously follows that it is positive.

• Eq. (20) has the solution    t m0 (x, s)h(s)ds v 0 (x, t ) = v0 (x) exp −δv

(32)

0

where h(s) is the HSP function which in our simulations is considered to be the Gamma distribution, thus a smooth enough function. One can easily see that the solution v 0 in (32) is positive and it satisfies

    t   0  ≤ ∥v0 ∥H (Ω ) ∥v ∥X =  v exp −δ m ( s ) h ( s ) ds 0 v 1   0

0

(33)

X

which is just the estimation (30) for k = 0. Using Eq. (32) and ideas similar to those in [21] we obtain the estimations for (27)

v 0 ∈ L∞ ((0, T ] × Ω ) and vt0 ∈ L∞ ((0, T ] × Ω ). Moreover, again by (32) we can easily deduce that

∥∇v 0 ∥L∞ (Ω ) ≤ C (Ω , T )

(34)

as well, where we made use of the fact that m0 is the solution of the heat equation, with the corresponding regularity implications. In the following we consider for simplicity the case with distributed delay. Hence I (h) is the solution of Eq. (9) and it is given (via (10)) in the form I (h(t )) = I (0)e−qt +

t



h(s)e−q(t −s) ds,

(35)

0

which is positive and in L∞ ([0, T ]). The proof of (25), (26) and (31) follows the ideas in [9,21]. However, we need to elude a difficulty arising from the haptotaxis term, i.e. from the last term on the left hand side of Eq. (21). In order to apply the Galerkin method we first let cl0 (t , x) :=

l 

dil (t )wi (x),

(36)

i=1 i 1 2 where {wi }∞ i=1 is an orthogonal basis of H (Ω ) and an orthonormal basis of L (Ω ). In (36) dl (t ) (0 ≤ t ≤ T , i = 1, 2, . . . , l) are chosen such that

dil (0) = (c0 , wi )

 and

d dt

cl0 , wi



+ B[cl0 , wi ; t ] = 0,

(37)

with (·, ·) denoting the inner product in L2 (Ω ). In Eq. (37), for c , v ∈ H 1 (Ω ) and a.e. 0 ≤ t ≤ T we use the time dependent bilinear form B[c , v; t ] :=

  n Ω i,j=1

aij (x, t )cxi vxj +

n  i=1

bi (x, t )cxi v + d(x, t )c v dx,

(38)

604

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

where

v 0 (x, t )



aij (x, t ) = D 1 − γ

Kv

K1 v 0 (x, t )

bi (x, t ) =



I (h)δij ,

i, j = 1, . . . , n

∂xi v 0 (x, t )I (h),

Kv + v 0 (x, t )

d(x, t ) = ∇ ·



K1 v 0 (x, t ) Kv + v 0 (x, t )

i = 1, . . . , n

 ∇v 0 (x, t )I (h) .

We first multiply Eq. (37) by dil (t ) and do the summation for i = 1, 2, . . . , l to obtain



d dt

cl0 , cl0



  + B cl0 , cl0 ; t = 0.

In order to apply Theorem 7.1.2 in [9], we have to show that aij , bi , d ∈ L∞ (Ω )

(39)

and that there exist some constants ρ > 0 and κ ≥ 0 such that

  ρ∥cl0 ∥2H 1 (Ω ) ≤ B cl0 , cl0 ; t + κ∥cl0 ∥2L2 (Ω )

(40)

for all 0 ≤ t ≤ T and l = 1, 2, . . . . (39) follows from the fact that I (h) ∈ L∞ ([0, T ]) and v0 ∈ W 2,∞ (Ω ). In order to show (40), we first note that ∂∂t + L is a parabolic operator, where L denotes the operator from which the bilinear form B in (38) stems. This means that for a.e. x ∈ Ω and all ξ ∈ Rn there exists a constant ϑ > 0 such that n 

aij (x, t )ξi ξj ≥ ϑ|ξ |2

(41)

i,j=1

for all (x, t ) ∈ Ω × [0, T ] and where ϑ is less than or equal to the smallest eigenvalue of the matrix A(x, t ) = (aij )i,j=1,...,n =

(D(1 − γ v

0 (x,t )

)δij )i,j . Kv Thus in view of (41) we can write ϑ



cl0 2 dx

|∇ | Ω

≤B

cl0



,

cl0



;t +

n 



i

∥b ∥L∞(Ω )

i=1

Ω

cl0

cl0

|∇ || |dx + ∥d∥L∞(Ω )

 Ω

(cl0 )2 dx.

Now, if we use Young’s inequality for the second term on the right hand side for an ϵ small enough such that ϵ < ϑ2 , then applying Poincaré’s inequality for the integral on the left hand side gives

ϑ 2(C˜ 2 + 1)

∥cl0 ∥2H 1 (Ω ) ≤ B[cl0 , cl0 ; t ] + ϑ

which is (40) with ρ = 1 4ϵ

2(C˜ 2 +1)

+ ∥d∥L∞ (Ω ) .



2 C˜ 2 + 1

n

i=1

∥bi ∥L∞ (Ω )

 + C ∥cl0 ∥2L2 (Ω ) ,

and κ = ˜ 22 + C where C˜ (n, Ω ) is the constant in Poincaré’s inequality and C = C +1

Since (39) and (40) are satisfied, Theorem 7.1.2 in [9] provides the energy estimates

 ∂

 

0 max ∥cl0 ∥L2 (Ω ) + ∥cl0 ∥L2 (0,T ;H 1 (Ω )) +  ≤ C ∥c0 ∥L2 (Ω ) .  ∂ t cl  2 0≤t ≤T L (0,T ;H 1 (Ω ))

(42)

Next consider the symmetric bilinear form A[ , cl0

cl0

 ] := Ω

a(x, t )(∇ cl0 )2 dx



0 with a(x, t ) = D 1 − γ Kv

v





I (h)In . Then multiply (37) by

d i d dt l

(t ), sum up for i = 1, . . . , l and integrate from 0 to T to obtain

 0 2    T  ∂ cl  d 1 0 0   dt + A[cl , cl ] dt  ∂t  2 0 0 dt L2 (Ω )       1 T Dγ I (h) ∂v 0 D T v0 (∇ cl0 )2 dxdt + (−qI (t ) + h(t ))(∇ cl0 )2 dxdt =− 1−γ 2 0 Ω Kv ∂t 2 0 Ω Kv T

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614 T 

 +

Ω

0 T



1

≤

2

K1 v 0 I (h)

Dγ I (h) ∂v 0

 Ω

0

cl0 ∇v 0 ∇

Kv + v 0

+ K1 ∥I (h)∥L∞ ([0,T ])

0

∂ cl ∂t

Ω

cl0

605

 0 dxdt

(∇ cl0 )2 dxdt +

∂t  T

Kv



0

D 2



∇v ∇

T

 0

∂ cl0 ∂t



 

1−γ

Ω

v0 Kv



(−qI (t ) + h(t ))(∇ cl0 )2 dxdt

dxdt ,

(43)

for which the first term on the right hand side can be written as (see [21]) 1

T



2

Dγ I (h) ∂v 0



∂t

Kv

Ω

0

(∇ cl0 )2 dxdt ≤

Dγ

 0  ∂v   ∥I (h)∥L∞ ([0,T ])  ∥cl0 ∥2L2 (0,T ;H 1 (Ω ))  ∂t  ∞ Kv L ((0,T ]×Ω )

= C1 ∥cl0 ∥2L2 (0,T ;H 1 (Ω )) Dγ

(44)

∂v 0

with C1 := K ∥ ∂ t ∥L∞ ((0,T ]×Ω ) ∥I (h)∥L∞ ([0,T ]) . v For the second term on the right hand side of Eq. (43) we get the estimate D

T



2

  Ω

0

1−γ

v0



Kv

(−qI (t ) + h(t ))(∇ cl0 )2 dxdt ≤ C2 ∥cl0 ∥2L2 (0,T ;H 1 (Ω ))

(45)

with a constant C2 depending on D, γ , Kv , q, ∥v 0 ∥L∞ ((0,T ]×Ω ) and on supt h(t ). For the third integral on the right hand side of Eq. (43), we use integration by parts with respect to t to obtain T

 0

 Ω

cl0 ∇v 0 ∇



∂ cl0 ∂t



 dxdt = Ω

 0  ∇v ∇ cl |T0 dx −

T

 0

 Ω

∇ cl0 ν(x, s)dxds

where

ν(x, s) =

∂ cl0 0 ∇v + cl0 (−δv m0 v 0 h). ∂t

This provides the bound also for this integral by using the energy estimates (42) and the estimates (27) and (29). Thus following the same steps (starting with c0 ∈ H 1 (Ω )) as in Theorem 7.1.5 in [9] proves that there exists a unique weak solution c 0 (x, t ) of Eq. (21) together with the corresponding boundary conditions (11) such that c 0 ∈ L2 (0, T ; H 2 (Ω )) ∩ L∞ (0, T ; H 1 (Ω )) ct0 ∈ L2 (0, T ; L2 (Ω )) and

∥c 0 ∥X + ∥c 0 ∥L2 (0,T ;H 2 (Ω )) ≤ C (Ω , T )∥c0 ∥H 1 (Ω ) . Furthermore, the positive initial condition c0 guarantees the positivity of c 0 (x, t ) by the weak maximum principle. (ii) Induction assumption: the lemma holds for an arbitrary k ∈ N0 . (iii) Induction step: the lemma holds for the term k + 1 of the iterative sequence.

• Indeed, using the induction hypothesis (25), (27) and (31), with an adequate embedding constant c1 (Ω , T ) we have  T  T ∥c k v k ∥2L2 (Ω ) h2 (t )dt ≤ ∥v k ∥2L∞ (Ω ) ∥h∥2L∞ ([0,T ]) c1 ∥c k ∥2H 1 (Ω ) dt 0

0

≤ 4c1 ∥v k ∥2L∞ (Ω ) ∥h∥2L∞ ([0,T ]) C 2 (Ω , T )T ∥c0 ∥2H 1 (Ω ) < ∞. By the theory of linear partial differential equations, there exists a unique weak solution to (22) together with the corresponding boundary and initial conditions and it satisfies mk+1 ∈ L2 (0, T ; H 2 (Ω )) ∩ L∞ (0, T ; H 1 (Ω )) mtk+1 ∈ L2 (0, T ; L2 (Ω )) and

∥mk+1 ∥X + ∥mk+1 ∥L2 (0,T ;H 2 (Ω )) ≤ C (Ω , T )(∥c0 ∥H 1 (Ω ) + ∥m0 ∥H 1 (Ω ) ), with an adequate constant C (Ω , T ).

606

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

The lower bound in (28) for mk+1 is obtained by using the auxiliary function θ k+1 (x, t ) := mk+1 (x, t ) − Cm e−δm t as it is done in [21]. Further, (23) is an inhomogeneous linear differential equation with solution

   t β(x, s)eα(x,s) ds , v k+1 (x, t ) = e−α(x,t ) v0 (x) +

(46)

0

where

  t   mk+1 (x, s)h(s)ds, α(x, t ) = δv 0   c k (x, s) v k (x, s)   − . β(x, s) = µv v k (x, s) 1 − Kc

Kv

By the induction hypothesis, we can obviously conclude that v k+1 ∈ W . In order to show that vtk+1 ∈ W we use Eq. (30) to write

∥vtk+1 ∥W ≤ µv ∥v k ∥W ∥1 −

ck Kc

−

vk Kv

∥W + δv ∥mk+1 ∥W ∥v k+1 ∥W ∥h∥L∞ ([0,T ]) < ∞,

where we also made use of the induction hypothesis and the fact ∥mk+1 ∥W < ∞ (see [9]). Thus we proved (27). Furthermore, in the same way as in the case with k = 0, we can show that ∇v k+1 ∈ L∞ (Ω ). In order to prove the positivity of v k+1 (x, t ), using (28) one can observe that

|β(x, s)| ≤ |µv v k (x, s)(1 − v k (x, s))| ≤ µv |v k (x, s)| ≤ µv

Kv 2

.

Since mk+1 (x, t ) ≥ Cm e−δm t , there exists a positive constant C˜ m with H =

v k+1 (x, t ) ≤

t 0

h(ν)dν leading to

Kv −δv C˜ m H µv K v Kv e + T ≤ 2 2 2

for an adequately chosen T (depending on the bounds of the HSP function h), which implies the positivity of v k+1 . By (46) and the induction hypothesis, we can write

∥v

k+1

∥H 1 (Ω )

2   t   −α(x,t ) −α(x,t ) α(x,s)   v0 (x) + e β(x, s)e ds = e 0 H 1 (Ω )  t   2 k   c (s) v k (s) 2 k ≤ 2∥v0 ∥H 1 (Ω ) + 2µv  v (s) 1 − − ds  Kc

0

Kv

H 1 (Ω )

 t 2  (v k (s))2  k  ( s ) − ≤ 2∥v0 ∥2H 1 (Ω ) + 2µ2v  v ds   1 Kv 0 H (Ω )   2  2 t    1  t k k   (v (s))2 ds ≤ 2∥v0 ∥2H 1 (Ω ) + 4µ2v  v ( s ) ds +   1   1 2 K H (Ω )

0

≤ (2 + 4µv C (Ω , T ))T ∥v ∥ 2

2

2 0 H 1 (Ω )

v

0



H (Ω )

,

which proves (30). In order to show the positivity of c k+1 an auxiliary function

ξ k+1 (x, t ) := ηt +

Kc 2

− c k+1 (x, t )

is used (see [21]) and substituted in the corresponding Eq. (24), leading to

        ∂ξ k+1 c k v k+1 K1 v k+1 I (h) η t + Kc k+1 k +1 k+1 −∇ D 1−γ I (h)∇ξ −∇ −ξ + ∇v ∂t Kc Kv K2 Kv v k+1 2   ck vk ≥ η − µc c k 1 − −α ≥0 Kc

Kv

upon choosing η > µc K2c and µc < 1.

(47)

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

607

By construction ξ k+1 (x, 0) ≥ 0 and it follows that c k+1 < ηt +

Kc 2

≤ Kc ,

(48)

for t ≤ T3 with an adequate T3 depending on Kc and the choice of η (T3 is in any case less than µ1 ). This shows that the c right hand side of (24) is positive. Therefore, the weak maximum principle together with the hypothesis c0 > 0 ensures the positivity of c k+1 . In order to apply Theorem 7.1.5 in [9], we have to show that the coefficients D(1 − γ are in the space L∞ ((0, T ] × Ω ) and

c k v k+1 Kc Kv

)I (h) and

K1 v k+1 I Kv +v k+1

(h)∇v k+1

  vk ck −α ∈ L2 (0, T ; L2 (Ω )). µc c k 1 − Kc

Kv

To this end, we start with

    k k+1   D 1 − γ c v I (h)  ∞ K K c

v

L

((0,T ]×Ω )

      c k v k+1  = ess sup0≤t ≤T ess supx∈Ω D 1 − γ I (h) Kc Kv    (48) v k+1  ≤ ess sup0≤t ≤T ess supx∈Ω D∥I (h)∥L∞ ([0,T ])  1 + γ  Kv

< ∞. The boundedness of

K1 v k+1 I Kv +v k+1

(h)∇v k+1 in W follows from the fact that v k and ∇v k are in W , which we already proved.

For the reaction term modelling logistic growth with crowding we have T

 0

   2  T   k 2 k  k ck vk k c 1 − c − α v  dt = c 1 − − α dxdt  Kc Kv L2 (Ω ) Kc Kv 0 Ω   T  (ck )4 (αv k c k )2 k 2 ≤ 3 (c ) + 2 + dxdt 2 0

(28)

≤

27 4

Ω

Kc

Kv

Kc2 |Ω |T .

Thus (see [9]), there exists a unique weak solution c k+1 (x, t ) for Eq. (24) with c k+1 ∈ L2 (0, T ; H 2 (Ω )) ∩ L∞ (0, T ; H 1 (Ω )) ctk+1 ∈ L2 (0, T ; L2 (Ω )). b. Existence In this part of the proof we show that the iterative sequence constructed above is a Cauchy sequence, which will lead to the existence of the solution (m, v, c ) as its limit. For k ∈ N0 we can write (see [9])

 k+1  m − mk 2 ≤ C (Ω , T ) X

T



  µm (c k v k h(t ) − c k−1 v k−1 h(t ))22

L (Ω )

0

dt

   2  2  K2 ≤ 2C (Ω , T )µ2m C32 max Kc2 , v T ∥h∥2L∞ (0,T ) v k − v k−1 X + c k − c k−1 X 4

≤

1  4

    v k − v k−1 2 + c k − c k−1 2 , X X

where C3 = C3 (Ω , T ) is an embedding constant, T ≤ [8C (Ω , T )µ2m C32 max{Kc2 ,

(49) Kv2 4

}∥h∥2L∞ (0,T ) ]−1 .

Replacing v k and respectively v k+1 in Eq. (23) and subtracting we have

    ∂ k+1 ck vk c k−1 v k−1 k k+1 k+1 k k k k−1 (v − v ) + δv (m v h(t ) − m v h(t )) = µv v 1 − − − µv v 1− − . ∂t Kc Kv Kc Kv

(50)

608

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

If we multiply both sides of (50) by v k+1 − v k and integrate with respect to x we obtain

2 d  k+1 v − v k L2 (Ω ) ≤ dt



 

k k−1 2µv  − v − v



(v k )2

−

Kv

(v k−1 )2



Kv

 −

vk c k Kc

−

 v k−1 c k−1   2 K c

L (Ω )



    + 2δv h(t ) v k (mk+1 − mk )L2 (Ω ) v k+1 − v k L2 (Ω )  2  2 1 ≤ µv (1 + h(t )) v k+1 − v k L2 (Ω ) + µv Kv2 h(t ) mk+1 − mk H 1 (Ω ) 4

1

+ 3µv ∥v k − v k−1 ∥2L2 (Ω ) + µv Kv2 ∥c k − c k−1 ∥2L2 (Ω ) .

(51)

4

After applying Gronwall’s inequality to Eq. (51) we get

 2  2   k+1 2 2  v − v k Z ≤ D1 (Ω , T ) v k − v k−1 Z + mk+1 − mk X + c k − c k−1 X T ,

(52)

where



T



D1 (Ω , T ) = exp µv

0

and D1 (Ω , T ) T ≤ k +1

k

1 . 4

   1 1 (1 + h(t ))dt max C4 µv Kv2 ∥h∥L∞ , 3µv , C4 µv Kv2 4

4

Here C4 denotes again an embedding constant.

− c it follows (compare with [9, Section 7.1])    2  T k k k−1   k+1 2 v k−1  µc c k 1 − c − α v − µc c k−1 1 − c  c − α − c k X ≤ C (Ω , T )  

For c

Kc

0

Kv

Kc

Kv

dt

L2 (Ω )

≤ D2 (Ω , T )T [∥c k − c k−1 ∥2X + ∥v k − v k−1 ∥2Z ], with D2 (Ω , T ) = C (Ω , T ) max{21µ2c C5 , 4µ2c By (49), (51) and (52) we can write

Kc2 Kv2

(53)

} and D2 (Ω , T )T ≤ 14 , where C5 denotes an embedding constant.

2 2  2  2   k+1 2  c − c k X + v k+1 − v k Z + mk+1 − mk X ≤ c k − c k−1 X + v k − v k−1 Z   which shows that c k , v k , mk is a Cauchy sequence in X × Z × X . c. Uniqueness This obviously follows from the fact that the constructed sequence is a Cauchy one.

(54)



4. Nondimensionalization In order to solve the system (5), the problem should be written in terms of dimensionless variables. To this aim rescale c˜ :=

c Kc

,

v˜ :=

v Kv

,

˜ := m

mDm

µm K c

,

 x˜ =

µv Dm

x,

t˜ :=

t T

,

where T is the reference time unit. The transformations (55) lead to the dimensionless parameters:

˜ = D

D µv T Dm

µ ˜ v = µv T ,

,

K˜ 1 =

Kv µv TK1 Dm

˜ m = µv T , D

,

µ ˜ c = µc T ,

δ˜m = δm T ,

δ˜v =

µm K c δ v T Dm

µ ˜ m = D m Kv T

and to the following nondimensionalized system (omit the tildes for ease of notation):

    ∂ c˜ = ∇. D˜ 1 − γ (˜c v˜ ) ∇ c˜ − ∇ · ∂ t˜ ∂ v˜ ˜ v˜ + µ = −δ˜v m ˜ v v˜ (1 − c˜ − v˜ ) ∂ t˜ ˜ ∂m ˜ − δ˜m m ˜ +µ = D˜ m ∆m ˜ m c˜ v˜ . ˜ ∂t



K˜ 1 v˜ 1 + v˜

 c˜ ∇ v˜

+µ ˜ c c˜ (1 − c˜ − α v˜ )

(55)

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

609

5. Numerical results In this section, we present some numerical results in 1D for the solution of our final model

∂c = ∇ · (Dc (c , v)I (h)∇ c ) − ∇ · (χH (v)I (h)c ∇v) + µc c (1 − c − αv) ∂t ∂v = −δv mv h(t ) + µv v(1 − c − v) ∂t ∂m ∂m = Dm ∆m − δm m + µm c v h(t ) or = Dm ∆m − δm m + µm c v h(t − τ ) ∂t ∂t

(56)

to be coupled with dI dt

= −qI + h(t ),

q>0

or dI dt

= −qI + h(t − τ ),

q>0

for the cases with distributed and constant delays, respectively. For the discretization of the model we use the finite difference method. We divide the interval [0, 1] into k parts with k + 1 nodes. We first obtain I n+1 for the (n + 1)-st time level applying forward difference in the corresponding equation for the flexibility. For the equation characterizing the dynamics of the MDE concentration we use forward differences for the time discretization and central differences for the discretization of the space derivative: mni +1 − mni

∆t

 = Dm

mni−+11 − 2mni +1 + mni++11

∆x2

 − δm mni +1 + µm cin vin h(t n+1 ),

i = 0, 1, . . . , k

(57)

with ∆t and ∆x being the increments in the time and space directions, respectively. Since the HSP function h is known at any time level, it is evaluated implicitly, whereas the c and v values are computed explicitly for the evaluation of m at the unknown time level (n + 1) leading to the (k + 1) × (k + 1) linear system of equations Am mn+1 = mn + ϑm n˜

(58)

where n stands for the known time level, mn+1 is the vector containing the values of m for the k + 1 space nodes at the (n + 1)-st time level, Am is the tridiagonal matrix coming from the FDM discretization and ϑm n˜ is the vector with the entries µm cin vin hi (t n+1 ) for i = 0, 1, 2, . . . , k. After solving (58) for mn+1 , these values are used to find the values for the ECM concentration v at the (n + 1)-st level at discretization points solving

vin+1 =

1 1 + ∆t δv mni +1 hi (t n+1 )

[vin + µv vin (1 − cin − vin )]

(59)

for each i = 0, 1, 2, . . . , k. In order to discretize the first two terms on the right hand side of the first equation of the PDE system (6) for the cancer cell density, we use a finite difference scheme as in [3,7].4

∇(D(c )∇ c )|xi =

1



2∆x2

k∈Ni

∇(χH (v)c ∇v)|xi =

1 2∆x2

(D(ckn , vkn+1 , I n+1 ) + D(cin , vin+1 , I n+1 ))(ckn+1 − cin+1 )

 (χH (vkn+1 , I n+1 )ckn + χH (vin+1 , I n+1 )cin )(vkn+1 − vin+1 ),

(60)

k∈Ni

where Ni ⊂ I is the index set pointing at the direct neighbours of xi on the grid. The matrix–vector form of the scheme reads Ac cn+1 = cn + ϑc n˜ ,

(61)

where Ac is the (k + 1) × (k + 1) tridiagonal matrix coming from the nonstandard finite difference discretization, ϑc is the vector coming from the proliferation term, cn+1 and cn are the vectors containing the c values at the discretized points for the time levels n + 1 and n, respectively. n˜

4 This scheme is called there ‘nonstandard’, as it handles the diffusion part explicitly, while the reaction terms are implicitly discretized.

610

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

t=0

1

v 0.9

c 0.8 0.7 0.6 0.5

m

0.4 0.3 0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Fig. 2. Initial data (62) with ε = 0.01.

Table 1 Parameter ranges in the model. Parameters

Range

Source

D (multiple of random motility coefficient Dc ) µc (proliferation of cancer cells) K1 (haptotaxis coefficient) δv (rate of degradation of ECM) µv (proliferation of ECM) Dm (MDE diffusion coefficient) δm (decay of MDEs) µm (production rate of MDEs)

10−5 –10−3 0.05–2 10−3 –1 1–20 0.15–2.5 0.001–10 0.1–10 0.01–5

[5,26] [5] [5,26] [5,26] [5] [5,26] [26] [26]

The boundary conditions for c and m are taken as given in (11). For the initial conditions we assume that there already exist cancer cells having penetrated a short distance into the ECM, while the rest of the space is occupied by the matrix itself. The density of the MDE is proportional to the initial cancer cell density and combining these assumptions leads to initial conditions of the form (compare [26])

 −x 2 , x ∈ [0, 1] and ε > 0 c (x, 0) = exp ε  2 −x v(x, 0) = 1 − exp , x ∈ [0, 1] and ε > 0 ε m(x, 0) = ζ c (x, 0), x ∈ [0, 1] and ζ ∈ [0, 1). 

(62)

In our simulations we fixed the following parameters: D = 0.0005,

Dm = 0.005,

K1 = 0.09,

δm = 0.1,

µm = 0.3,

δv = 10.

The parameter ranges employed in the simulations can be found in Table 1 (taken from [5,26]). In the figures we illustrate the variations of the tumour cell density, ECM density and the MDE concentration in space. The initial conditions are plotted in Fig. 2 and the competition parameter α in the proliferation term for the cancer cells is taken to be 0.5. As mentioned before, the invasion is influenced on the one part by the tissue fibres (a too high density of which would impede the motion) and on the other part by the density of available tumour cells themselves. Actually the interaction of these two populations is the main factor conditioning the spread of the cells. Therefore, we tried to capture this influence not only through the couplings in the reaction terms, but also including it in the form of the random motility coefficients. A first choice, Dc (c , v) = 1 − γ ( Kc + Kv ), has been proposed e.g., by Gatenby & Gawlinski [11] (in a slightly different c v context), however this did not account for a genuine interaction, but only for the joint effect on inhibiting the spread. The form Dc (c , v) = 1 − γ c v (after the nondimensionalization) involves the interaction, but this can lead for large enough values of γ to unrealistic oscillations and blowup, due to the random motility coefficient becoming negative. An example

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

611

t=40 1

v

0.8

0.6

c 0.4

0.2

m 0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Fig. 3. Oscillations with D(c , v) = D(1 − 4c v). Parameter values: Dm = 0.005, K1 = 0.09, D = 0.0005, δm = 0.1, µm = 0.3, δv = 10, µc = 0, µv = 0.

a

b

t=5

v

1

t=15

v

1

0.8

0.8

c

c 0.6

0.6

0.4

0.4

m 0.2

0.2

0

0

0.1

0.2

0.3

0.4

c

0.5 x

0.6

0.7

0.8

0.9

0 0

1

0.1

0.2

0.3

0.4

d

t=50

1

0.8

0.8

c

0.5 x

0.6

0.7

0.8

0.9

1

0.9

1

t=55

1

0.6

m

0.6

c

0.4

0.4

v 0.2

0.2

m 0 0

0.1

v

m 0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

Fig. 4. Comparison of the case with a constant diffusion coefficient Dc (c , v) = D (solid line) with the case D(c , v) = 1+Dc v (dashed line) µc = 0.06, µv = 0.03.

can be seen in Fig. 3 for γ = 4. While keeping all other parameters fixed, there seems to be a critical value of γ (above 2) for which the solution starts to oscillate. However, finding an analytical expression for this oscillation threshold seems to be a nontrivial task, which is even further hindered by the fact that there are no reliable sources for selecting the rest of the D parameters. This motivates another choice of the diffusion coefficient, Dc = 1+θ , still involving the interaction of the cells cv with the ECM, but never becoming negative. The latter will be used throughout the rest of the numerical simulations.

612

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

a

b

t=5

t=10

v

1

v

1 0.8

c

0.8

0.6

c

0.6

0.4

0.4

m 0.2

0

0

0.1

0.2

0.3

0.4

c

0.5 x

0.6

0.7

0.8

0.9

0

1

v

1

0

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

1

t=30 v

1

0.8

0.8

c

c

0.6

0.6

0.4

0.4

0.2

0.2

m

m 0

0

0.1

0.2

0.3

0.4

e

0.5 x

0.6

0.7

0.8

0.9

0

0.2

0.3

0.4

f

t=50

0.5 x

0.6

0.7

0.8

0.9

0.6

0.7

0.8

0.9

t=60 1

0.8

0.8

c

c

0.6

0.1

1

1

0.6

0.4

0.4

0.2

v

0.2

m

m 0

0.1

d

t=20

0

m

0.2

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0 0

v 0.1

0.2

0.3

0.4

0.5 x

1

Fig. 5. Effect of the reestablishment of fibres: µc = 0.15, µv = 0.75 (dashed line), µc = 0.15, µv = 0.0 (solid line). Dc (c , v) = 1/(1 + c v).

In Fig. 4 we illustrate the comparison between the case with constant diffusion Dc = D and the one where Dc = 1+Dc v . As expected, in the former case the cell spread and thus also the decay of the ECM is faster. Fig. 5 shows the effect of including the re-establishment term for the fibres. In that case the cells are able to invade a smaller area compared to the situation where there is no tissue re-establishment. As time passes the tumour cells use the produced MDEs to degrade the tissue fibres and to enhance their invasion, eventually penetrating the entire region under consideration. In the last set of graphics shown in Fig. 6 we compare the cases with delay (both distributed and constant, see (8), (9), (10)) and the case without delay. For the latter we considered the flexibility function of the form I (h) = h. The computations are carried out for q = 2 for both kinds of delay and the constant time lag τ is taken to be 1. These figures confirm a faster invasion in the case without delay. On the other hand, taking a constant time lag shows a slower invasion compared to the case with distributed delay, which is expected, since the former includes the delay directly.

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

a

b

t=10

v

1

613

t=20

v

1

0.8

0.8

c 0.6

0.6

0.4

0.4

m

0.2

0

c

0

0.2

0.1

0.2

0.3

0.4

c

0.5 x

0.6

0.7

0.8

0.9

0

1

m 0

0.2

0.3

0.4

d

t=30

v

1

0.5 x

0.6

0.7

0.8

0.9

1

0.8

0.9

1

t=40 1

v

0.8

0.8

0.6

0.6

c

c 0.4

0.4

0.2

0.2

m 0

0.1

m 0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

Fig. 6. Cases with distributed delay (solid line), with constant delay (dashed line) and without delay (stars). µc = µv = 0.

6. Conclusion In the present note we proposed a mathematical model for the influence of HSPs on tumour invasion. Our model has a multiscale character, as it couples the dynamics on the subcellular microscale with the macroscopic dynamics of cancer and ECM fibre populations along with those for the concentration field of matrix degrading enzymes. It is obvious that processes on the microscale influence the behaviour on the macrolevel, and multiscale settings of this type – though not being very detailed – provide a way of characterizing such influences. They can be further developed in order to accommodate more level of detail, according to the particular subcellular effect one would like to assess on the macroscale. The multiscality achieved by coupling the micro and macroscopic dynamics is a rather new feature of models in mathematical biology. Even multiscale models involving biological processes on the subcellular level, individual cell behaviour on the mesoscopic scale and macroscopic population dynamics are less wide spread than their single scale counterparts. In spite of the simplicity of coupling, the present micro–macro model provides a first step towards more involved models and offers at the same time a means to test (via numerical simulations) some qualitative assertions arising from experimental observations (e.g., the effect of intracellular factors like HSPs on the tactic behaviour of a cancer cell population and on the surrounding tissue). The model presented in this work is different from the previous ones, since it includes the effects of HSP on both amoeboid and mesenchymal motions and considers nonlinear diffusion and nonconstant, solution dependent haptotaxis coefficients. These settings are more realistic, however lead to a more complicated model, for which we prove the local existence of a unique weak solution by way of a fixed point argument. Furthermore, our numerical simulations show a good agreement with the expected behaviour for the invasion and for the effect of HSP as well as for different choices of the time delay, as commented in Section 5. Acknowledgments We thank Prof. M. Mimura (Meiji University) for discussions concerning the choice of the diffusion coefficients for the tumour cells.

614

G. Meral, C. Surulescu / J. Math. Anal. Appl. 408 (2013) 597–614

G.M. acknowledges the support by the German Academic Exchange Service (DAAD) and the Turkish Council of Higher Education (YÖK) during her visits in 06-09.2010 at the University of Stuttgart and in 07-09.2011 at the University of Münster. C.S. acknowledges the support of the Baden-Württemberg Foundation. References [1] D. Abraham, M. Ponticos, H. Nagase, Connective tissue remodeling: cross-talk between endothelins and matrix metalloproteinases, Curr. Vasc. Pharmacol. 4 (2005) 369–379. [2] A.R.A. Anderson, M.A.J. Chaplain, E.L. Newman, R.J.C. Steele, A.M. Thompson, Mathematical modelling of tumour invasion and metastasis, J. Theoretical Med. 2 (2000) 129–154. [3] R. Anguelov, J.M.S. Lubuma, Contributions to the mathematics of the nonstandard finite difference method and its applications, Num. Meth. PDE 17 (2001) 518–543. [4] N. Bellomo, A. Bellouquid, J. Nieto, J. Soler, Complexity and mathematical tools toward the modelling of multicellular growing systems, Math. Comput. modelling 51 (2010) 441–451. [5] M.A. Chaplain, G. Lolas, Mathematical modelling of cancer invasion of tissue: dynamic heterogeneity, Networks and Heterogeneous Media 1 (2006) 399–439. [6] A. Chauviere, T. Hillen, L. Preziosi, Modelling cell movement in anisotropic and heterogeneous tissue: dynamic heterogeneity, Networks and Heterogeneous Media 2 (2007) 333–357. [7] H.J. Eberl, L. Demaret, A finite difference scheme for a degenerated diffusion equation arising in microbial ecology, Electron. J. Differential Equations Conf. 15 (2007) 77–95. [8] H. El-Samad, H. Kurata, J.C. Doyle, C.A. Gross, M. Khammash, Surviving heat shock: control strategies for robustness and performance, Proc. Natl. Acad. Sci. USA 102 (2005) 2736–2741. [9] L.C. Evans, Partial Differential Equations, AMS, 1998. [10] F. Filbet, P. Laurençot, B. Perthame, Derivation of hyperbolic models for chemosensitive movement, J. Math. Biol. 50 (2004) 189–207. [11] R.A. Gatenby, E.T. Gawlinski, A reaction–diffusion model of cancer invasion, Cancer Res. 56 (1996) 5745–5753. [12] J. Guay, H. Lambert, G. Gingras-Breton, J.N. Lavoie, J. Huot, Regulation of actin filament dynamics by p38 map kinase-mediated phosophorylation of heat shock protein 27, J. Cell Sci. 110 (1997) 357–368. [13] R.K. Hansen, I. Parra, S.G. Hilsenbeck, B. Himelstein, S.A. Fuqua, Hsp 27-induced MMP-9 expression is influenced by the Src tyrosine protein kinase yes, Biochem. Biophys. Res. Commun. 282 (2001) 186–193. [14] T. Hillen, M 5 , mesoscopic and macroscopic models for mesenchymal motion, J. Math. Biol. 53 (2006) 585–616. [15] J. Kelkel, C. Surulescu, On some models for cancer cell migration through tissue, MBE 8 (2011) 575–589. [16] J. Kelkel, C. Surulescu, A multiscale approach to cell migration in tissue networks, Math. Models Methods Appl. Sci. 22 (3) (2012) 1150017-1–115001725. [17] J. Kelkel, C. Surulescu, A weak solution approach to a reaction–diffusion system modelling pattern formation on seashells, Math. Meth. Appl. Sci. 32 (2009) 2267–2286. [18] E. Lambert, E. Dasse, B. Haye, E. Petitfrere, TIMPs as multifacial proteins, Crit. Rev. Oncol. Hematol. 3 (2004) 187–198. [19] P. Lemieux, S. Oesterreich, J.A. Lawrence, P.S. Steeg, S.G. Hilsenbeck, J.M. Harvey, S.A. Fuqua, The small shock protein hsp 27 increases invasiveness but decreases motility of breast cancer cells, Invasion Metastasis 17 (3) (1997) 113–123. [20] N. MacDonald, Biological Delay Systems: Linear Stability Theory, Cambridge University Press, 1989. [21] C. Märkl, G. Meral, C. Surulescu, Mathematical analyis and numerical simulation for a system modelling acid-mediated tumour cell invasion, Int. J. Analysis. 2013 (2013) Article ID 878051, 15 pages. [22] A. Peper, C.A. Grimbergen, J.A. Spaan, J.E. Souren, R. van Wijk, A mathematical model of the hsp70 regulation in the cell, Int. J. Hyperthermia 14 (1) (1998) 97–124. [23] T.R. Rieger, R.I. Morimoto, V. Hatzimanikatis, Mathematical modelling of the eukaryotic heat-shock response: dynamics of the hsp70 promoter, Biophys. J. 88 (3) (2005) 1646–1658. [24] C. Shäfer, P. Clapp, M.J. Welsh, R. Benndorf, J.A. Williams, HSP 27 expression regulates CCK-induced changes of the actin cytoskeleton in CHO-CCK-A cells, Am. J. Physiol. Cell Physiol. 227 (6) (1999) C1032–C1043. [25] C. Surulescu, On the stationary interaction of a Navier–Stokes fluid with an elastic tube wall, Appl. Anal. 86 (2007) 149–165. [26] Z. Szymanska, J. Urbanski, A. Marciniak-Czochra, Mathematical modelling of the influence of heat shock proteins on cancer invasion of tissue, J. Math. Biol. 58 (4–5) (2009) 819–844. [27] Z. Szymanska, M. Zylicz, Mathematical modelling of heat shock protein synthesis in response to temperature change, J. Theoret. Biol. 259 (2009) 562–569.