Math. Model. Nat. Phenom.
Non-homogeneous random walks, subdiffusive migration of cells
and anomalous chemotaxis
S. Fedotov111Corresponding author. E-mail: firstname.lastname@example.org, 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:
Random cell movement plays a very important role in embryonic morphogenesis, wound healing, cancer cells proliferation, and many other physiological and pathological processes . 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  (see also ). 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
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 . The master equation for can be written as
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 .
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.2. Anomalous random walks
as it is done in  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
is the Riemann-Liouville fractional derivative with varying order
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 . 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
as This result has been interpreted as anomalous aggregation of cells at the point . 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 . 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 . 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.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 :
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
where is the total flow of cells from the point to
and is the total flow of cells from the point to
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 . This phenomenon does not exist in the Markovian case. For the Markov model (1.10) the flux is independent from
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 . 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
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
The Markov model (1.2) corresponds to the following choice
It is convenient to introduce the rate of escape (hazard function) from the point as
If we denote the survival function at the point as
and the residence time probability density function as
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
where the transition densities and are defined as
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 ). 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
where the rate of jump to the right and the rate of jump to the left are defined as
The residence time probability density function
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 . 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
We consider only the case when the residence time of random walker at is equal to zero, so the initial condition is
where . The boundary condition in terms of residence time variable ( can be written as 
In what follows we consider only positive values of . In which case, we have to specify the boundary condition for . We write
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
The structural density can be rewritten in terms of the survival function (2.10) and the integral arrival rate
Our purpose now is to derive the master equation for the probability
Let us introduce the integral escape rate to the right and the integral escape rate to the left as
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
The balance equation for can be written as
In the Laplace space we have the following expressions for escape rates
Using inverse Laplace transform and shift theorem we find
where and are the memory kernels defined by Laplace transforms
for The balance equation for is
where The master equation for can be rewritten in terms of the probability flux from the point to
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
where the exponent
depends on the state Residence time probability density function has the Pareto form
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
and to the left
Let us consider the non-local model for which the jump probabilities and depend on the chemotactic substance as follows
The transition PDF’s and can be rewritten in terms of and as
Therefore, the integral escape rates to the right and to the left in the subdiffusive case are
Here is the Riemann-Liouville fractional derivative with varying order defined by (1.7). The anomalous rate functions and are
For with , we obtain
The fractional probability flux from the point to is
The equation (3.13) can be rewritten in terms of the probability flux as
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
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
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