1. Introduction
Tumor cell invasion is a complex process that involves both the proliferation and migration of cancer cells. Proliferation refers to the ability of cancer cells to divide and grow, while migration refers to their ability to move and invade surrounding tissues. Researchers have devoted large efforts in investigating the development of malignant tumors, particularly those that are characterized by dysregulated cell migration and uncontrolled proliferation. Gliomas, a type of brain tumor, are known for their high rates of cell proliferation and motility, as well as their ability to invade normal tissue away from the multicellular tumor. Clinical studies and investigations have indeed shown that the proliferation rate of migratory cells in gliomas is often lower in the invading tumor tissue compared to the tumor core [
1,
2]. This relationship between cell migration and proliferation, known as the migration–proliferation dichotomy, has been shown to play a critical role in the invasiveness of malignant gliomas and is an important factor that must be taken into consideration when developing treatments for these tumors. Despite advances in medical technology and research, the switching between two phenotypes, as well as the invasion process, remain not fully understood. This presents a challenge in developing effective treatments and accurate prognoses for patients with this type of tumor. The mathematical modeling of the migration and proliferation dichotomy of cancer cells has been an active area of research in recent years, as it provides a valuable tool for understanding the complex interplay between these two processes in cancer progression and invasion. The migration–proliferation dichotomy in gliomas has been studied using a variety of approaches. These include a stochastic two-state switching model, where the migration–proliferation of cancer cells was formulated in terms of a continuous time random walk (CTRW) [
3,
4,
5,
6,
7]; multiparametric modeling of the phenotype switching [
8]; the simulation of multiscale glioma growth and invasion using an agent-based approach [
9,
10]; and numerical and analytical approaches, where the phenotype switching depends on oxygen in a threshold manner, which were developed in the framework of reaction–diffusion equations [
11]. It was proposed by Khain et al. [
12,
13] that the motility of cancer cells is a function of their density. The phenomenon of phenotypic switching caused by the density effect was also discussed in [
14,
15]. In recent years, there has been a significant amount of research focused on the “Go-or-Grow” mechanism and its effect on glioma invasion. The role of the “Go-or-Grow” mechanism in glioma invasion, which is thought to be linked to oxygen availability in the environment of a growing tumor, was proposed by Hatzikirou et al. [
16]. Multiscale modeling for acid- and vasculature-mediated glioma invasion and the effects of go-or-grow dichotomy and tissue anisotropyin were discussed in [
17]. Mansury and colleagues [
18] introduced the use of evolutionary game theory in an agent-based brain tumor model as a means of studying the evolution and dynamics of brain tumors. An alternative approach to studying the growth and progression of tumors was proposed by Böttger et al. [
19], in which a model that specifically incorporates tumor cell dynamics was developed to analyze the interplay of proliferation and migration processes in detail. Recently, Godlweski et al. [
20] identified a single microRNA that plays a key role in the regulation of glioma cell proliferation, migration, and responsiveness to glucose deprivation. New approaches to modeling glioma invasion have emerged in recent years, with sources [
21,
22,
23] providing examples.
Previous research by Fedotov et al. [
24] demonstrated that the use of the standard diffusion approximation for transport processes, in conjunction with logistic growth, results in an overestimation of the overall propagation rate in glioma. This finding was verified in [
5,
6]. The presence of anomalous diffusion in the stochastic movement of cells is frequently observed, thereby highlighting the significance of fractional diffusion equations and fractional derivative in modeling some biological phenomena, such as the cancer cell motility [
25] and tumor growth [
26]. The persistent random walk model constitutes a well-established framework that can be easily augmented to incorporate extensions, such as cell reactions [
27,
28]. However, most existing methods and techniques for addressing the issue of persistent random walks with reactions focus on Markovian switching between two states [
29,
30,
31]. According to a study described in [
32], cell motility can be explained by a persistent random walk model, which demonstrates that cells have a memory of their past velocities.
In [
33], the authors proposed a model that incorporates a three-state switching mechanism and generates anomalous subdiffusion. The resulting governing equations for this formulation are
where
and
are the probability density functions (PDFs) of particles moving with constant speed
in the positive or negative direction for exponentially distributed running times with rate
, and
is the PDF of a particle that has zero velocity for Mittag–Leffler distributed resting times. Here,
,
, and
describe three possible transitions that the particle can make from the rest, and
is the integral escape rate from the rest state and is defined as
where
denotes the Riemann–Liouville derivative and is given by
It is important to highlight that for non-Markovian switching states in random walk theory, a general formula for the escape rate can be represented through convolution of the memory kernel and the density, which will be illustrated later.
Although the model (
1)–(
3) involves a non-Markovian anomalous rest state, the switching running state is Markovian. The main objective of this paper is to expand upon the work of the anomalous stochastic transport model of particles outlined in [
33] by incorporating non-Markovian running states. Specifically, we will assume that the running time is arbitrarily distributed [
34], rather than solely following an exponential distribution as described in [
33]. A second extension of [
33] is the implementation of a reaction term into the non-Markovian model. Another main objective of this paper is to extend the analysis of the migration–proliferation dichotomy model presented in [
6]. For a biologically motivated reaction–diffusion model, it may be worthwhile to enhance the model by incorporating additional states for the cells or the particles. Our approach addresses this by by including an additional state for migrating cells, in addition to the existing migratory and proliferative states. We are interested in investigating the impact of considering the movement of cells in both positive and negative directions with finite velocities
, as opposed to the movement of cells in a single direction with a fixed velocity, as previously considered in [
5,
6]. Our goal is to create a more realistic theory for the phenomenon of phenotype cell switching. Our model provides a description of the complex cell transport in the framework of a persistent random walk with memory effects, taking into account anomalous transport and the resulting slow mobility of cancer cells.
2. Three-State Proliferation–Migration Model Involving a Residence Time Variable
In this paper, we provide an analysis of the migration and proliferation of tumor cells using a three-state persistent random walk model on a one-dimensional space with continuous sample paths (no jumps). The model considers transitions between moving states with constant velocity (migratory phenotype) via an intermediate resting state with zero velocity (proliferating phenotype). More specifically, we suggest a model in which the process of tumor cell invasion consists of three possible states: two moving states, in which the tumor cells randomly move (migrate) with a constant velocity in the positive (state +) or negative (state −) direction but are unable to proliferate and can switch phenotypes; and a resting state (state 0), in which the tumor cells are allowed to proliferate but do not migrate and can also switch phenotypes.
One way to characterize the migration–proliferation dichotomy is through the use of a commonly used phenomenological model that employs reaction–diffusion equations. The basic structure of our model is as follows: The migratory phenotype cell remains in state + during a running time (state +’s residence time) or in state − during a running time (state −’s residence time). In these moving states, the cell starts to move with constant speed at time in the positive (state +) or negative (state −) direction, and after a random time (running time) it switches to a proliferating phenotype cell. Following a waiting time (state 0’s residence time) spent in state 0, the cell makes the choice to switch to some next state. Specifically, it either returns to a migratory phenotype (state + or state −) or remains at its resting state (state 0). Clearly, the fundamental aspect of this process can be represented by a three-state random process, where the cell exhibits one of three distinct states: positive directional movement of tumor cells, negative directional movement of tumor cells, or a state of rest for tumor cells. We assume that the residence times are random variables.
To describe the transition process between these states, we introduce the switching rate,
, as a function of residence time
such that [
35]
where
and
are PDFs of the residence times and the corresponding survival functions, respectively.
are defined as
In this case, and represent the conditional probability of transition from migrating states to a proliferating state in the small interval , given that there is no transition up to time . The product has a similar meaning for the transition from a proliferating state to one of migrating states. Thus, is the switching rate from state + to state 0, is the switching rate from state − to state 0, and determines the switching rate from state 0 to either state + or state −. Since the switching rate is dependent on the duration of time the cell has stayed in a specific state, our approach has the advantage of being adaptable to various forms of the residence time probability distribution. This is noteworthy due to recent evidence indicating that the PDF of the residence time may not be exponential. In this study, we introduce three PDFs: , and . The functions and represent the probability densities that the position of cells with a migratory phenotype (moving state) is within the interval at time t and moving with positive and negative velocities, respectively. The function represents the probability density that the position of cells with a proliferating phenotype (resting state) is within the interval at time t and has zero velocity.
In this section, we aim to develop a migration and proliferation model for cancer cells and derive non-Markovian master equations to describe the dynamics of the cells transport process. Transition probabilities
,
and
can be introduced that depend on the residence time
spent in the migrating and proliferating states, respectively [
35]. It is useful to define the structured densities of cancer cells as dependent on the residence time
[
36]. Let
and
be the densities of migrating cells at point
x at time
t, whose residence time in a migratory phenotype state (state + or state −, respectively) lies in the interval
. The density of cells in a proliferating state that corresponds to the aforementioned migrating cells is
.
The balance equations for structured densities , and can be written as follows.
For a proliferating state,
The nonlinear function
is the proliferation rate, where
. For example, we can use the logistic growth for cell proliferation involving the cell proliferation rate
U and the carrying capacity
K such that [
6]
We assume that the initial conditions for the system are given by
where
is the initial densities of cancer cells
. It corresponds to the case when the residence time of all cells at
equals to zero.
The boundary conditions at zero running/waiting time (
) are
where
,
and
describe the conditional transition probabilities from proliferation state to a state of positive velocity, negative velocity, or remaining at its resting state, respectively. Specifically,
represents the probability of transitioning to a state of positive velocity,
, while
represents the probability of transitioning to a state of negative velocity,
. On the other hand,
represents the transition probability that the resting cell remains at rest again, where
. Generally, the probabilities
,
and
can facilitate the formulation of the effect of the external force.
The mean densities at point
x at time
t,
, can be obtained by integrating the structured densities
over residence time variable
3. Non-Markovian Three-State Model
The aim of this section is to set up a non-Markovian model for the migration and proliferation of cells by eliminating the residence time variable
and find equations for
,
and
by solving the partial differential Equations (
8)–(
10) together with the boundary conditions (
13)–(
15) at
and initial condition (
12) at
.
The product
gives the phenotype switching rate corresponding to a particular residence time
. If we denote the fluxes between migrating and proliferating cells by
,
and
, then the switching term
can be obtained by integrating
over variable
from 0 to
t:
It follows from Equations (
13)–(
15) and (
17) that
Applying the method of characteristics, we obtain the following solutions to Equations (
8)–(
10):
For migrating states
For a proliferating state
It can be observed that all of the derived solutions (
19)–(
21) contain an exponential factor
that can be understood to represent the survival function
:
Considering the boundary conditions specified in Equations (
13)–(
15), the solutions represented by Equations (
19)–(
21) can be expressed in terms of the survival function from (
22) and the fluxes between the migrating and proliferating states (switching terms) from (
17) as follows:
It should be noted that the residence time PDF,
, can be expressed in terms of the switching rate,
, as follows [
35]:
The nonlinear master equations for unstructured densities
,
and
can be obtained by differentiating (
16) with respect to time
t together with (
8)–(
10):
By using (
13)–(
17), we obtain a system of integro-differential equations for the PDFs
,
and
in terms of the switching terms
,
and
:
Note that the above Equations (
30)–(
32) can be obtained by using the Fourier–Laplace transform technique. The switching terms
and
describe the average flux of cells from migrating states to proliferating state and
describes the average flux of cells from the proliferating state to migrating states.
The terms accounting for changes in state
,
and
are expressed in terms of
,
and
, respectively, as follows (see
Appendix A for the details of the derivation):
where
is the memory kernel defined by its Laplace transform [
37]:
Here,
and
are the Laplace transforms of the residence time PDF
and the survival function
, respectively. The nonlinear master Equations (
30)–(
32) together with interaction terms (
33)–(
35) are a generalization of a linear system of equations obtained in [
33,
38].
The Fourier–Laplace transform of the total density
, where the proliferating rate
, can be written as (see
Appendix B)
where
,
.
In what follows, we will examine various distributions for residence time PDFs, , including power law distributions.
5. Anomalous Three-State Model with Mittag–Leffler Distributed Rest Times
In this section, we consider an anomalous switching case and determine an average position of a cancer cell,
, and the mean squared displacement (MSD),
, when there is no proliferating process in state 0:
. Given that the proliferation of cancer cells is influenced by numerous conditions, it is assumed that a characteristic scale for the proliferating residence time is not present [
6]. As a result, the PDF for the residence time in the proliferating state exhibits a power law behavior. Thus, we introduce a switching term,
, that incorporates the Riemann–Liouville fractional derivative, which is equivalent to the waiting time for a cell in the rest state, with a corresponding PDF derived from a Mittag–Leffler function. Its survival function is expressed as
then
where
is a parameter with units of time and
is the one-parameter Mittag–Leffler function. The parameter
is the measure of the strength of the resting state [
6]. A decrease in
increases the probability of a longer residence time in the proliferating state. As
approaches 1, the Markovian case with exponentially distributed resting times,
, is recovered.
The residence times in the proliferating state will be distributed approximately as
for large values of
(as
), leading to the formation of “heavy” or power law tail distributions. Power law distributions have been widely observed in various empirical studies of stochastic processes with intricate underlying mechanisms. Additionally, numerous examples have demonstrated the significance of power law waiting times in real-world phenomena.
The Laplace transforms of residence time PDF
and its corresponding survival function
are
The Laplace transform
corresponding to Equation (
52) can be approximated by
for small
s[
40]. The mean waiting time
is infinite in this case. The memory kernel
in terms of its Laplace transform is
In this situation, it is more convenient to find the switching term
using the Laplace transform version of Equation (
35):
Inserting
into the above equation and then taking the inverse Laplace transform gives
where
is the Riemann–Liouville fractional derivative defined in (
5). Contrary to the proliferating process, an average transport time is finite. Therefore, we assume that the residence time PDFs for migrating states are exponentially distributed with rate
where
is constant. The Laplace transform of
and
are
and the Laplace transform of their corresponding survival functions
and
are
In order to comprehend the implications of incorporating a heavy-tailed waiting time distribution for the rest state, two basic scenarios are analyzed: a single velocity model and a symmetric dual velocity model, both of which possess a non-Markovian rest state. The first and second moments at the long-term limit are performed for each of these models.
5.1. Single Migration State Model
The aim of this subsection is to show that if the migrating cancer cells move with a constant velocity
, then the mean cell position
increases as
for
. Consider the particular case when all cells at
start to move to the positive direction with the velocity
from the point
. To achieve this, we set
,
and
, for which the initial density of cells in migrating states are
and
, and the initial density of cells in proliferating state,
, is zero. Then,
and
. The system of integro-differential equations for the PDFs
and
in terms of the switching term
is
It follows from Equation (
37) that
In the limit
, we find from (
54), (
59) and (
60) that
To find the average position of cancer cells for large time asymptotic
, we use the Laplace transform of
defined by
where
is the Fourier–Laplace transform of the total density of cells
.
By using Equation (
65), we find
Finally, taking the inverse Laplace transform, the average position of cells is
which is sublinear. The average position of cancer cells was estimated in [
6] using a different idea, yet the same anomalous behavior was obtained. By employing a similar approach, the MSD,
, of this model can be determined and show that
. The anomalous advection reflects the memory effect associated with the slow motion of cells due to the power law of the residence time distribution for the proliferating state with an infinite mean residence time (
52).
5.2. Symmetric Active States Model
The objective of this subsection is to demonstrate by calculation the MSD that, upon expanding the previously considered model by incorporating an additional active state with a velocity of , results in subdiffusion in the long-time asymptotic limit.
Equations (
30)–(
32) can be combined and reduced to a system of governing equations for the total density function
and
.
By adding (
30)–(
32), we obtain
Now multiplying (
30) and (
31) by
and subtracting the resulting equations, and by setting
,
and
, we obtain
Then, by differentiating (
69) with respect to
t and (
70) with respect to
x and eliminating the flux
, and combining it with (
32), we arrive at the governing equations
We consider now the symmetrical initial conditions for which the cells start to move from the point
at
as follows:
and
, and the initial density of cells in proliferating state,
, is zero. Then
and
. From (
71) and (
72), and by using
, we obtain a reduced system of integro-differential equations for the PDFs
and
:
Performing the Fourier–Laplace transform to (
73) and (
74), we obtain
Mathematically, the MSD is defined by
which in the Fourier–Laplace space reads
Calculating the MSD at the long time limit using the above formula, we find
and subsequent Laplace inversion, we obtain
Once more, it is observed that the power law PDF for the residence time in the proliferating state results in subdiffusive behavior. The calculation of the second moment is performed utilizing the Laplace technique in [
33]; however, the same anomalous behavior is observed.
6. Conclusions
We developed a three-state non-Markovian model to characterize the migration–proliferation dichotomy of tumor cells with general residence time distributions. Specifically, the model differentiates between a migratory state, where cells move without undergoing proliferation, and a proliferative state, where they do not migrate but instead undergo division. The objective of this study was to provide a mesoscopic description of the anomalous transport and reactions associated with the migration–proliferation dichotomy mechanism. This work further analyzed the two-state non-Markovian model for the migration–proliferation dichotomy of cancer cells introduced in [
6] by incorporating an additional state for migrating cells. It also extended the stochastic transport of particles outlined in [
33]. In comparison to the work of Han et al. (2021), we enhanced the model by introducing non-Markovian running states. In particular, we assumed that the running time follows an arbitrary distribution instead of being limited to an exponential distribution as stated in [
33]. Additionally, we incorporated a reaction term into the non-Markovian model as another extension of [
33]. The developed probabilistic approach for the transport of tumor cells with directional in migration speed and a rest state was described using a persistent random walk with an arbitrary residence time distribution, while the proliferation rate was modeled using a nonlinear function of the densities of all three cell types. A set of balance equations for cancer cells with two phenotypes, which can switch randomly between cell proliferation and migration, was derived with a non-local switching term involving the Riemann–Liouville fractional derivative. To achieve this, the Markovian model was employed as the initial approach, under the premise that the transition probabilities depend on the residence time variable.
To investigate the nature of our model in the anomalous case involving power law residence time distributions, we calculated the average position of cancer cells in a single migration state model that corresponds to sublinear growth in time: , for . The MSD for symmetric bi-directional migration model was estimated and showed that , for . We demonstrated through analytical methods that if cells move in a migratory state with constant velocity , using power law residence time distributions for the proliferative state (the Mittag–Leffler distributed waiting times for rests) results in subdiffusive behavior in the long time asymptotic limit. This describes the impact of an anomalous rest state on persistent random walks of cells with finite velocity: the longer a cell survives in a proliferative state, the smaller the probability of switching to a migratory state.
Our proposed model offers a general framework for analyzing reaction–transport systems characterized by non-Markovian and anomalous transitions between active and inactive states. This model is suitable for bi-directional intracellular transport [
41] that involves a resting state for power law distributed times. Recently, it was experimentally found that the residence time of active movement follows a power law distribution [
42]. The presence of subdiffusive behavior in a stochastic transport system due to anomalous switching constitutes a beneficial aspect of the proposed model, enhancing its potential applications, particularly in describing biological or ecological movement, including the fields of population theory and cellular biology. Our intention is to apply the model to intracellular transport.