Time delay in necrotic core formation.

A simple model of avascular solid tumour dynamics is studied in the paper. The model is derived on the basis of reaction-di(cid:11)usion dynamics and mass conservation law. We introduce time delay in a cell proliferation process. In the case studied in this paper the model reduces to one ordinary functional-di(cid:11)erential equation of the form that depends on the existence of necrotic core. We focus on the process of this necrotic core formation and the possible in(cid:13)uence of delay on it. Basic mathematical properties of the model are studied. The existence, uniqueness and stability of steady state are discussed. Results of numerical simulations are presented.


INTRODUCTION
The process of tumour growth and its dynamics is one of the most intensively studied processes in the recent years.There have appeared many papers devoted to it (cf.[1,2,6,8,9,12] and references therein).This process can be divided into several different stages, starting from the very early stage of solid tumour without necrotic core inside (cf.[6]).In the present paper we focus on the next stage, i.e. the process of necrotic core formation.In this stage there are three main cellular processes -proliferation, apoptosis and necrosis.It should be marked that solid tumour growth leads to the limited size, which is shown theoretically (compare [18,19]) and experimentally ( [2, 3, 10, 12, 20]) as well.
Following [6,12] we introduce time delay to the model presented in [8].In [12] the model of necrotic core formation without time delay was considered.On the other hand, the model studied in [6] does not include a process of necrotic core formation.The aim of this paper is to derive and analyse the model of necrotic core formation with the presence of time delay in the proliferation process, as proposed in [8].
The model we study is based on the idea of symmetric growth of avascular multicellular spheroids (MCS), which was described in [8].It is assumed that the tumour growth depends on the nutrient (like glucose or oxygen) concentration.However, the mean time of the nutrient diffusion is much shorter then the mean time of tumour doubling that leads to a quasi-steady state approximation.The following notation is used in the paper: • R(t) and R nec (t) denote the external and necrotic radius of MCS at time t, respectively; • σ ∞ denotes the external concentration of nutrients, which is assumed to be constant; σ N is the minimal nutrient concentration needed for proliferation; σ = σ ∞ − σ N and we assume σ > 0; • Γ, a, b are positive coefficients of proliferation, apoptosis and necrosis, respectively and s is a scaling constant.
Following [8] we introduce time delay to the model, namely in proliferation process.A heuristic argument for introducing time delay is the duration of the mitosis process that could be important.Time delays can be introduced also for another cellular processes (compare [7,8,13,14]) or for next stages of tumour growth, e.g. for angiogenesis as in [4], where delays explain oscillations that appear in vascular tumour growth ( [15]).Finally, in the case of the presence of necrotic part, the tumour evolution is governed by the following system of equations (see [8] for the detailed derivation) and R nec is described by the implicit function of R in Eq. (1a).

DERIVATION AND ANALYSIS OF THE MODEL
At the beginning of this Section we notice that Eqs. ( 1) are well posed only if there exist positive solutions to (1a) for R = R(t) and R = R(t − τ ) which are smaller then R(t) and R(t − τ ), respectively.Now, we derive a final model and study the existence and properties of solutions to it.Eq. (1a) define the implicit function which describes the connection between R and R nec .
Following the analysis presented in [12] we obtain that if R is small, i.e.R < R = 6σ Γ , then there is no positive solution to Eq. (1a) and then it is reasonable to assume that R nec = 0. On the other hand, if then Eq. (1a) has exactly one positive solution.
We consider four cases.If R(t − τ ) > R, and R(t) > R, then the tumour growth is described by Eqs.(1).
The right-hand side of Eq. ( 6) is a continuous Lipschitz function of R(t) and R(t − τ ).This implies that the model is well posed and its solutions exist and are unique (see [17]) .
is the steady state for Eq. ( 4).We focus on the process of necrotic core formation and hence, R 0 (t) < R for all t ∈ [−τ, 0], i.e. there is no necrotic part at the beginning.We have two cases.
Consequently, Th 2.1 from [6] implies that the steady state R is stable independently on τ .Following the proof of Th. 3.1 from [6] we prove that if the initial data satisfies R 0 (t) < R, then the steady state R is globally stable.Hence, no necrotic core is formed in this case.
2. If R > R, then following the proof of Th 3.1 from [6] we show that R(t) reaches the level R for some t > 0 and the necrotic core is formed.
For the necrotic core formation, the only interesting case is the second one.Hence, combining the inequality which defines the second case with the assumption σ > 0 one gets As in [12] the asymptotic behaviour of R nec (R) can be approximated as follows: On the other hand, if R nec → 0, i.e.R → R, the asymptotic is the following: Using the formula for R nec (R) that was calculated in [12] we have This yields, that the size of the proliferation ring decreases as the tumour radius increases and give the following estimate for The same analysis as in [12] shows that there exists at least one steady state R. In the next section we focus on the problem of uniqueness and stability of R.

Proof:
Notice, that if a solution to our model is negative, then there exists time t 0 such that R(t 0 ) = 0.If R(t 0 − τ ) < R, then the solution fulfills Eq. ( 4) and analysis presented in [6] yields that it cannot be negative.On the other hand, if R(t 0 − τ ) > R, then inequality holds and theorems from [5] yield that solutions to our model are nonnegative.
Lemma 2 (Global existence) Let an initial datum R 0 (t) ≥ 0 for all t ∈ [−τ, 0].Then for any t > 0 the following inequality holds and the solution is globally defined.

Proof:
Assume that Ineq.( 9) does not hold for some time t > 0. Due to continuity of R there exists time t 0 such that R(t 0 ) = R max , R(t) ≤ R max for any t < t 0 and R (t 0 ) ≥ 0. The estimates ( 8) yields which contradicts the assumption (9).The upper estimates together with nonnegativity of the solution yields that it is globally defined.

Steady state
It is obvious that a trivial steady state always exists.For its stability we refer to [6].In this Section we focus on the problem of uniqueness of the positive steady state for (6) in the case when a formation of necrotic core is possible.In order to reduce number of coefficients we denote Thus, we can rewrite the right-hand side of Eq. (1b) in the form and Eq.(1a) as First we state the following Let y(x) be a solution to (11) having the property 0 ≤ y(x) ≤ x.Then there is a unique solution to equation F (x, y(x)) = 0 for x ≥ η.
Using an implicit function theorem we obtain: We would like to show, that y − y 0 grows.This yields that y − y 0 has at most one extrema and this give that there exists only one solution to F (x, y(x)) = 0 since the number of solution to F (x, y(x)) = 0 is odd.We substract and then differentiate y (x) − y 0 (x) with respect to x.Since the denominator is positive we consider only the nominator which has the form V (x, y) =25(x − y) 5 (x + 2y)(2x 4 + 2x 3 y + 9x 2 y 2 + 4xy 3 + y 4 ) (13a) +30(x − y) 3 β(x + y)(3x 4 + x 3 y + 9x 2 y 2 + 3xy 3 + 2y 4 ) (13b) Notice, that for β > 0 lines (13a)-(13c) are nonnegative.Hence, we consider the second part of V , presented in the lines (13d)-(13f).Subtracting a common part 9y(x − y)(x + 2y) which is positive we consider We show, that V (x, y) > 0, assuming that η Thus, using the assumption η It can be readily calculated, that Assumption of Lemma give g(η) > 0. Calculating g (x) it is easy to find, that if g (x 0 ) = 0, then x 0 < η, due to Assumption of Lemma.
Analogous calculations lead to the following estimate.There exists a unique positive steady state if the following condition 27γ

Proof:
First notice, that ( 7) is equivalent to (12).Hence, Lemma 3 yields that the steady state exists and is unique.
In order to study stability of the steady state we substitute x(t) = R 3 (t) and y(t) = R 3 nec (t).Hence, Eq. (1b) takes the form Without lost of generality we may assume that sΓ 30 = 1 (we can rescale time to eliminate this coefficient).Thus, we have Calculating the characteristic quasi-polynomial at the steady state x and denoting ȳ = y(x) one obtains where For τ = 0 the uniqueness of steady state, and the condition F (η, 0) > 0 implies that the steady state is stable and A + B < 0. It is easy to see that B − A < 0 for every positive parameters that implies stability of R independently on τ -for details see [11], compare also [16].

NUMERICAL SIMULATIONS AND DISCUSSION
In this Section we present results of numerical simulations.The aim of these simulations is to compare the behaviour of solutions to the model presented in this paper with the model which do not consider the necrotic core formation as well as illustrate some possible behaviour of the solution to Model (6).We study also dependence of solutions on the model parameters, particularly on time delay.Th. 4 implies that for a wide range values of coefficients there should not be a quantitative change in the behaviour of the solutions.At the beginning, we use the following values of parameters σ ∞ = 12.0 , σ N = 1.0 , Γ = 30.0, a = 2.0 , b = 2.0 , s = 1.0 , τ = 0.6 .
In this case, the coefficients used in Section 3 are the following The constant function R 0 (t) = R 0 = 0.5 is taken as an initial one.The notation used on the pictures is the following: the solid line (R 1 ) denotes the tumour radius in the model with necrotic core formation; the dashed line (R 1 nec ) describes the radius of the necrotic core; the dotted line (R 2 ) describes the radius of tumour without necrotic core formation.The stars denote points at which the tumour radius reaches the level R, i.e. the necrotic core is formed (or disappeared).First notice, that periodic solutions which appear in the model without necrotic core formation are not present if we consider this process (see Figs. 1).The coefficients used in the cases presented on Figs. 1 do not fulfill the assumption of Th. 4, since β < 0. The case where β > 0 but Ineq.(12) do not hold is presented on Figs. 2. It turns out that for a wide range of coefficients (for which) the dependence of solutions on time delay is not relevant.
Figs. 1 show also the dependence of the solution on parameters Γ and σ N .If we put Γ small, then (as it could be expected) the tumour growth is slower at the begging.However, the behaviour of solutions for greater values of time t is similar to those with larger Γ.In fact, notice that changes of Γ do not change the condition (12).
In Figs. 3 and on the right-hand of Figs. 2 the solutions form "stairs".Notice, that solutions stabilize very fast at some level.Then, after some time when delayed functions reach the higher values the solution grows very rapidly to the new quasi-steady state.This behaviour suggests that the convergence of solutions assuming given, constant delayed function, is very fast.Notice that the time the solution needs to approximate to the steady state is greater for larger values of delay parameter.These are the only observed influence of time delay.We would like to point out, that in the cases presented in Figs. 3, solutions to the model without necrotic core formation become negative vary fast (for t around 4 at the case presented on the left-hand figure and for t around 12 for the case presented on the right-hand one).
The conclusion is that the process of necrotic core formation is very important.The behaviour of solutions is more stable in this case.Theorem 4 yields that for a wide range of coefficients steady state is asymptotically stable.On the other hand, computer simulations suggest that Assumption of Th. 4 can be weakened.
The analysis shows that if R < R, then the steady state is globally stable.Simulations suggest that it remains stable when the necrotic core is formed.Next, we study the dependence on the minimal nutrient concentration σ N .Although the qualitative behaviour is similar for all positive values of σ N , the surprising numerical result is that the steady state is not a decreasing function of σ N (see Fig. 4).For σ N = 0 the tumour radius stabilises around 4.2, then this level increases and the maximal value is achieved for σ N between 6 and 7 and then the value of steady state decreases.
On the other hand, the width of proliferation ring decreases as σ N increases (i.e.σ decreases) which could be expected.Finally, we want to present that delay may be important.Figs. 5 show that undumping oscillations may arise.However, to obtain it negative values of parameter σ N was used, which implies that Assumptions of Th. 4 are not valid for this case.Hence, there is no sense to consider it from the biological point of view.
For the values of coefficients used in the simulations the qualitative behaviour of the model (6) with different values of delay parameter is similar to the model without delay presented in [12].However, for large values of time delay the influence of the delay is noticeable (see Figs. 3) On the other hand, Figs. 5 shows that delay could be important for some values of parameters.However, we have only shown that the steady state is locally stable for some range of parameters, the results of the simulations give a hint that there might exists only one globally stable steady state for any nonnegative coefficients.The solutions to Model (6) are more stable than the solutions to the model presented in [6], in which the process of necrotic core formation was not considered.

ACKNOWLEDGMENT
This paper was presented in the framework of 6-th EU Programme MRTN-CT-2004-503661.

Figure 4 :
Figure 4: The dependence of the value of steady state (on the y-axes) on the minimal nutrient coefficient σ N (on x-axes);

Figure 5 :
Figure 5: Comparison of solutions to the model with necrotic core formation depending on parameter σ N .On the left-hand figure σ N = −12, and −14 on the right-hand one.