1. Introduction

Math. Model. Nat. Phenom.


2013

Non-homogeneous random walks, subdiffusive migration of cells

and anomalous chemotaxis

S. Fedotov111Corresponding author. E-mail: sergei.fedotov@manchester.ac.uk, A. O. Ivanov and A. Y. Zubarev

School of Mathematics, The University of Manchester, Manchester, M13 9PL, UK

Department of Mathematical Physics, Ural Federal University, Ekaterinburg, 620083, Russia

Abstract. This paper is concerned with a non-homogeneous in space and non-local in time random walk model for anomalous subdiffusive transport of cells. Starting with a Markov model involving a structured probability density function, we derive the non-local in time master equation and fractional equation for the probability of cell position. We derive the fractional Fokker-Planck equation for the density of cells and apply this equation to the anomalous chemotaxis problem. We show the structural instability of fractional subdiffusive equation with respect to the partial variations of anomalous exponent. We find the criteria under which the anomalous aggregation of cells takes place in the semi-infinite domain.

Key words: anomalous random walks, cell migration, aggregation

AMS subject classification:

1. Introduction

Random cell movement plays a very important role in embryonic morphogenesis, wound healing, cancer cells proliferation, and many other physiological and pathological processes [34]. The microscopic theory of the migration of cells and bacteria towards a favorable environment (chemotaxis) is based on random walk models [8, 19, 31, 32]. The “velocity-jump” models concern with self-propelled motion involving the runs and tumbles, while “space-jump”models deal with the cells making jumps in space. Much of the literature on the theoretical studies of cells motility has been concerned with Markov random walk models (see, for example, [4, 8]). However, the analysis of random movement of wild-type and mutated epithelial cells shows the anomalous dynamics of cell migration [7] (see also [27]). Over the past few years there have been several attempts to model non-Markovian anomalous cell transport involving subdiffusion and superdiffusion [7, 9, 10, 12, 13, 17]. In this paper we shall deal with a non-Markovian “space-jump”model that describes the non-homogeneous in space subdiffusive transport of cells.

1.1. Markov random walk model.

First let us consider a Markov model for random cell movement along one-dimensional lattice such that all steps are of equal length . We define the probability

(1.1)

that the position of cell  is at point  at time We introduce at each point  the rate of jump to the left  and the rate of jump to the right . This random walk is called a generalized birth-death process [6]. The master equation for  can be written as

(1.2)

This model corresponds to the case when intervals between jumps at point are exponentially distributed with parameter When the cell makes a jump from the position  it jumps to the right with probability and to the left with the probability [6].

The dependence of and on space can be introduced in different ways depending on how cells sense the surrounding environment. For the local chemotaxis models, the rates  and are the functions of the local concentration of the chemotactic substance There exist several non-local and barrier models that are different in terms of the dependence of rate functions on the chemotactic substance [4, 32]. For example, the rates and can depend on the concentration of the chemotactic substance at neighbouring positions and as in (3.6). In the continuous limit, the master equation (1.2) can be reduced to the classical advection-diffusion equation in which the cell flux involves the standard diffusion term and the advection term due to chemotaxis.

If we consider only positive values of we need to implement boundary conditions at the point . Here we assume that if cell hits the wall on the boundary, it is reflected with the probability and absorbed by wall with the probability Then one can write In the limit we obtain

(1.3)

where

Non-uniform stationary solution of master equation (1.2) can be interpreted as cell aggregation phenomenon [32]. In particular, if there is no absorption on the boundary (), the stationary solution can be easily found from (1.2) and (1.3). We obtain

(1.4)

where

provided the series is convergent.

1.2. Anomalous random walks

It is tempting to generalize the master equation (1.2) for the anomalous case by replacing the time derivative with the Caputo derivative [2, 24, 26]

(1.5)

as it is done in [33] for a fractional linear birth–death process. Here is the anomalous exponent: Although this generalization is very attractive from a mathematical point of view, it is not appropriate for a non-homogeneous medium for which the exponent depends on . The non-homogeneous fractional equation for can be written as

(1.6)

where

is the Riemann-Liouville fractional derivative with varying order

(1.7)

Here is the anomalous exponent corresponding to the site and the anomalous rate coefficients and have to be determined, see (3.12). The crucial point here is that the anomalous exponent depends on the site . The fractional equation (1.6) cannot be rewritten in terms of Caputo derivative (1.5). It turns out that even small non-homogeneous variations of the exponent lead to a drastic change of in the limit  [14]. It means that the subdiffusive fractional equations with constant anomalous exponent are not structurally stable. If, for example, the point has the property that for all and one can find that

(1.8)

as This result has been interpreted as anomalous aggregation of cells at the point [12]. In this paper we shall find the conditions for anomalous aggregation for the semi-infinite interval It should also be noted that non-homogeneous variations of the exponent destroy the Gibbs-Boltzmann distribution as a long time limit of the fractional Fokker-Planck equation [14]. Of course, for the constant value of the formulation in terms of Caputo or Riemann-Liouville operators are equivalent, as long as proper care is taken of the initial values [2, 24, 26].

1.3. Anomalous diffusion with reaction.

Another extension of traditional Markov random walks models is non-Markovian theory of anomalous transport with reaction dynamics [11, 28, 30, 35, 36, 37]. In particular, this theory has been used for the analysis of the proliferation and migration dichotomy of cancer cells [9, 10, 13]. In this paper we consider the inhibition of cell growth by anticancer therapeutic agents. To model this inhibition we introduce the random death process with non-uniform death rate parameter. We assume that during time interval at point  each cell has a chance of dying, where is the death rate [20]. It is easy to take into account this process for Markov models. We just add the term to the right hand side of the master equation (1.2). On the contrary, the anomalous master equation involves a non-trivial combination of transport and death kinetic terms because of memory effects [18, 28, 1]. In this paper we shall derive the following fractional equation

(1.9)

1.4. Mean field master equation for the density of cells.

Instead of the probability for an individual cell one can consider the mean density of cells as a function of space and time . The master equation (1.2) can be rewritten as the equation for the density by changing the variables as and :

(1.10)

where is the jump size, is the death rate. The advantage of this equation is that one can easily take into account various non-linear effects by assuming the dependence of the rate functions and on the average density

In the anomalous subdiffusive case, the master equation for mean field can be obtained from (1.9). It can be written as a mass balance equation

(1.11)

where is the total flow of cells from the point to

(1.12)

and is the total flow of cells from the point to

(1.13)

Here and are the anomalous rate functions, see (3.12). One can see that the flow of cells depends on the death rate . It means that in the anomalous case one cannot separate the transport of cells from the death process [18]. This phenomenon does not exist in the Markovian case. For the Markov model (1.10) the flux is independent from

When the density is conserved (), the master equation (1.11) can be approximated by the fractional Fokker-Planck equation [2, 25, 26]

(1.14)

This is an example of the fractional equation with varying anomalous exponent [5]. Note that as , see (3.18).

The purpose of the next section is to set up a non-Markovian discrete-space random walk model describing cell motility involving memory effects, the death process and subdiffusive transport.

2. Non-Markovian discrete-space random walk model

2.1. Random cell motility

There exist numerous mechanisms that facilitate random cell movement [34]. In this paper we adopt the following random model of cell motility. When the cell makes a jump to position , the time the cell spends here before it makes a jump to point or is random. It is called the residence time or waiting (holding) time. We define the residence time at position as

(2.1)

where and are the independent random times of jump to the left and right respectively. The idea here is that there exist internal cellular signals involving two ”hidden” independent random alarm clocks. If one of the clocks goes off first, say , the cell moves to the right to the point . The other clock ”tells” the cell to move left to the point if it goes off first (. Note that migration of cells is a highly complicated dynamic process which is regulated by both intercellular signals and the surrounding environment. Since we do not know the exact mechanism of cell motility we use a stochastic approach involving two random times and for jumping to the left and right. Note that if the random times and are exponentially distributed with the rates and respectively, we have a classical Markov model with the master equation (1.2). If the random variables and are not exponentially distributed, the standard Markov approach does not work. In this section we consider the non-Markovian case when and are independent positive random variables with general survival functions

(2.2)

The Markov model (1.2) corresponds to the following choice

(2.3)

It is convenient to introduce the rate of escape (hazard function) from the point as

(2.4)

If we denote the survival function at the point as

and the residence time probability density function as

then [6]

(2.5)

Now we determine this rate function in terms of statistical characteristics of random residence times and It follows from the definition of the residence time at position (2.1) that the survival function can be written as a product

where and are defined by (2.2). Differentiation of this equation with respect to gives

(2.6)

where the transition densities and are defined as

(2.7)

The formula (2.6) is the particular case of the general expression for the residence time PDF in terms of the transition densities (see formula (5) in the classical paper by van Kampen [22]). These transition densities have a clear probabilistic meaning. For example, is the probability that the cell’s jump to the left occurs in the time interval since the cell arrived at point . If we divide both sides of (2.6) by the survival function and use the formula (2.5), we obtain

(2.8)

where the rate of jump to the right and the rate of jump to the left are defined as

(2.9)

Note that the transition rates and can be introduced from the very beginning as it is done in [22]. By using (2.5) and (2.8), we write the survival function as

(2.10)

The residence time probability density function

takes the form

For the Markov case for which and are independent of the residence time variable we obtain from (2.10) the standard exponential survival function

corresponding to the Markov master equation (1.2).

2.2. Structured probability density function

If the residence time probability density function is not exponential, the random walk is non-Markovian. The standard method to deal with non-Markovian stochastic processes is to add auxiliary variables to the definition of the random walk to make the process Markovian [6]. Here we introduce the structured probability density function involving residence time as auxiliary variable. The structural density gives the probability that the cell position at time is at the point and its residence time at point is in the interval This is a standard way to deal with non-Markovian random walks [6, 28]. Suppose that cells die at random at rate that depends on The density obeys the balance equation

(2.11)

We consider only the case when the residence time of random walker at is equal to zero, so the initial condition is

(2.12)

where . The boundary condition in terms of residence time variable ( can be written as [6]

(2.13)

In what follows we consider only positive values of . In which case, we have to specify the boundary condition for . We write

(2.14)

This equation tells us that when cells escape from the point and move to the left with the rate , they are adsorbed by the wall with probability , and reflected back to the position with the probability Note that this boundary condition can be written in many different ways, for example, the cells can be reflected to state . One can also introduce a residence time PDF for a wall such that the reflection is not instantaneous.

We solve (2.11) by the method of characteristics

(2.15)

The structural density can be rewritten in terms of the survival function (2.10) and the integral arrival rate

as

(2.16)

Our purpose now is to derive the master equation for the probability

(2.17)

Let us introduce the integral escape rate to the right and the integral escape rate to the left as

(2.18)

Then the boundary conditions (2.13) and (2.14) can be written in a very simple form:

(2.19)

and

(2.20)

It follows from (2.9), (2.12), (2.16) and (2.18) that

(2.21)
(2.22)

Substitution of (2.12) and (2.16) to (2.17), gives

(2.23)

It is convenient to introduce the integral escape rate as the sum of the escape rate to the right and the escape rate to the left

(2.24)

The balance equation for can be written as

(2.25)

To obtain a closed equation for we need to express and in terms of By applying the Laplace transform to (2.21), (2.22) and (2.23), we obtain

and

In the Laplace space we have the following expressions for escape rates

(2.26)

Using inverse Laplace transform and shift theorem we find

(2.27)

where and are the memory kernels defined by Laplace transforms

(2.28)

It follows from (2.19), (2.24), (2.25) and (2.27) that the master equation for the probability is

(2.29)

for The balance equation for is

or

(2.30)

where The master equation for can be rewritten in terms of the probability flux from the point to

(2.31)

as

(2.32)

In the next section we shall derive fractional master equation for

3.  Anomalous subdiffusion in heterogeneous media

We now turn to the anomalous subdiffusive case. We assume that the longer cell survives at point , the smaller the transition probability from becomes. It means that the transition rates and are decreasing functions of residence time We assume that

(3.1)

where is a parameter with units of time. Both and play a very important role in what follows. From (2.10) and (3.1) we find that the survival function has a power-law dependence

where the exponent

(3.2)

depends on the state Residence time probability density function has the Pareto form

(3.3)

The anomalous subdiffusive case [2, 28] corresponds to

We can notice from (3.1) that the ratios and to are independent of the residence time variable that is

In this case it is convenient to introduce the probabilities of jumping to the right

(3.4)

and to the left

(3.5)

Note that these jump probabilities are completely determined by the anomalous exponents and In the standard CTRW theory, these jump probabilities are given independently [2, 28].

Let us consider the non-local model for which the jump probabilities and depend on the chemotactic substance as follows

(3.6)

where the parameter is determined from These jump probabilities describe the bias of cells with respect to the spatial gradient [4, 32]. One can obtain [17]

(3.7)

The transition PDF’s and can be rewritten in terms of and as

(3.8)

The asymptotic approximation for the Laplace transform of the waiting time density of the Pareto form (3.3) can be found from the Tauberian theorem [15]

with

(3.9)

We obtain from (2.28) and (3.8) the Laplace transforms of the memory kernels

(3.10)

Therefore, the integral escape rates to the right and to the left in the subdiffusive case are

(3.11)

Here is the Riemann-Liouville fractional derivative with varying order defined by (1.7). The anomalous rate functions  and are

(3.12)

with the anomalous exponent defined in (3.2).The master equation (2.29) takes the form of non-homogeneous fractional equation

(3.13)

for

For with , we obtain

(3.14)

The fractional probability flux from the point to is

(3.15)

The equation (3.13) can be rewritten in terms of the probability flux as

(3.16)

3.1. Fractional Fokker-Planck equation for cells density and chemotaxis

In this subsection we consider the continuous case () and find the drift together with diffusion coefficient in the fractional Fokker-Planck equation (1.14). It follows from (3.12) that the drift is proportional to the difference in the anomalous exponents and , since

(3.17)

The difference can be approximated in the different ways. In the case of chemotaxis this difference is proportional to the gradient of the local concentration of the chemotactic substance Using (3.7), we obtain

In the limit , we have the standard chemotaxis model

(3.18)

where the case   corresponds to the negative taxis: the drift is in the direction of the decrease in the value of the chemotactic substance . The fractional Fokker-Planck equation (1.14) takes the form