Existence of traveling wave solutions to data-driven glioblastoma multiforme growth models with density-dependent diffusion

1.   Introduction

World Health Organization (WHO) reported in ^[1] that cancer is the second leading cause of death. Worldwide, nearly 1 in 6 deaths is due to cancer and about 9.6 million deaths registered in 2018. There are more than 100 types of cancer and in all the types the abnormal cells are dividing constantly and form the growths called tumours ^[2].

Cancerous tumours are malignant and they invade the surrounding tissue, moreover, the tumours can spread to other parts of human body, producing the new, metastatic, tumours. The International Agency for Research on Cancer (IARC) classifies the malignancy of tumours by grade I-IV. The brain tumour is considered to have the most severe health consequences. Diagnosis of the brain tumours is very difficult and almost impossible for being detected in the early stages. According to ^[3], the brain cancer and cancer of the central nervous system (CNS) appear in 1.6% of new cases and 2.5% of deaths out of all cancers. According to statistics published in ^[4] in January 2020, approximately 18 000 people may die from brain and CNS tumours within a year.

Glioblastoma Multiforme (GBM) is an agammaessive malignant brain tumour rated by WHO as a grade IV astrocytoma, i.e. of the highest grade and most malignant of all gliomas ^[5]. The patients with GBM have a very low chance to overcome the disease, the median survival time is approximately 14.6 months. GBM is highly invasive and the tumour cells easily dissipate into the normal brain tissue. Consequently, researchers around the world are intensively investigating the questions about the brain neoplasms such as formation, growth, invasion and spread of the tumour. Mathematical modeling of the tumour dynamics and treatment becomes increasingly popular since it allows researchers to study those problems by formulating and testing various clinical scenarios based on plausible assumptions. Earlier mathematical models on the GBM growth can be found in ^[6,7,8] while modeling efforts on GBM treatments can be found in ^[9,10]. Readers interested in modeling GBM growth and treatments may benefit from the review papers on the subject in ^[11,12] and an excellent chapter in ^[13]. More interesting and complicated mathematical models for GBM take into account the natural behavior of the tumour e.g., hypoxic and normoxic cell migration ^[14], phenotypic switching or so called "go-or-grow" cells strategy ^[15,16]. Most of existing modeling efforts focused on qualitative comparison of model simulation results to clinical observations. One of the rare exceptions is the recent GBM growth modeling effort presented in ^[17] which contains careful model formulation, data validation based on the data published in ^[18] and a rigorous mathematical study of the existence of a traveling wave solution in the model describing the growth of the GBM.

A well formulated and data-validated mathematical model may produce tumour growth pattern mimics the experiment's result and predict future tumour growth by simulation. The logistic, Gompertz, Malthusian, von Bertalanffy and Bernoulli growth functions are widely used in mathematical models of tumour growth ^[19,20]. To accurately describe tumour spatiotemporal dynamics, appropriately selected spatial terms are needed in a serious mathematical model. As a result, most existing GBM growth and/or treatment models take the form of reaction-diffusion equations (RDEs) or systems of PDEs involving RDEs, see e.g. ^{[12,16,17,19,21,22,23]}.

The growth of GBM can be described by the reaction-diffusion equations or proliferative-invasive models. Many of such invasive phenomena in biology, medicine, ecology and the problems in physics and chemistry with a traveling wave solution were observed since 1937 when two important papers were published ^[30]. R. A. Fisher in his work "The wave of advance of advantageous genes" ^[24] considered the following nonlinear equation

where $g(u) = ku(1-u)$. He found that the traveling wave solution exists for all $c\geq c_0 = \sqrt{Dk}$ with $c$ a wave speed.

At the same time A.N.Kolmogorov et al.^[25] published the proof that solution under initial condition

converges to the travelling wave $w$ of minimal speed $c_0$ in sense that $u(t, x+m(t))\to w(x)$ for $t\to\infty$ and $m(t)$ taken appropriately, and $m(t)\to c_0.$ Contributions to the study of solution and its convergence were made by the different authors e.g. ^[26,27] but the papers ^[28,29] by Aronson and Weinberger in 1975 resolved the speed discussion ^[30]. They introduced the notion of asymptotic speed of propagation of disturbances and proved that even in higher space dimension, this speed is the minimal velocity of the travelling plane wave. All these influential papers gave rise to publications in application of the reaction-diffusion equations in biology. The equation itself is now referred in the literature as a Fisher-KPP equation (FKPP), Fisher’s equation, Kolmogorov–Petrovsky–Piscounov equation or KPP equation.

The travelling wave $u(x, t)$ represents a solution of FKPP equation and is assumed to have the same shape at all time and the speed of this wave $c$ appears stationary, i.e. it can be written in the form of ^[31] $u(x, t) = u(x-ct) = u(z), \, z = x-ct.$ The dependent variable $z$ is called the wave variable. After substitution it into the PDE, the equation in $x, t$ becomes an ordinary equation in terms of $z.$ The Fisher's equation is invariant in $x$ and thus $u(x, t) = \phi(x+ct)$ is also a wave solution with the boundary conditions $\phi(-\infty) = 0, \, \phi(\infty) = 1.$ If we consider the radius instead of $x$ then the Fisher's equation in radially asymmetric disk case looks as follows

The term $\frac{1}{r}\frac{\partial u}{\partial r}$ depends on $r$ and therefore the substitution $z = r-ct$ does not lead to ODE. Wave solution takes time to form and emerges in locations distant from the origin. For large values of $r$, the solution of Eq (1.2) tends asymptotically to a traveling wave solution with a certain speed $c$ of the classical Fisher's equation.

The first models of GBM growth appeared in ^[32,33] and improved in ^[8]. In ^[8], GBM tumor growth is described by the following Fisher equation ^[20,34],

Here $u(r, t)$ describes GBM cell density at time $t$ in position $r$, $\rho$ is a constant rate of GBM cell proliferation, $u_{M}$ is GBM cell density upper limit (carrying capacity) and $D$ is the GBM cell diffusion rate.

The model presented in ^[18] takes into account the different behavior of proliferating and invasive cells. The tumor core is modeled as a sphere increasing in radius at constant rate $v_c$ and shedding invasive cells at rate $s$. The invasive cells $u_i(r, t)$ diffuse and proliferate as in Fisher equation but move away the tumor spheroid at speed $v_i$. The model of Stein et al. ^[18] is as follows

where the forces such as random diffusion, logistic growth, taxis, cells shedding from the tumor core acting upon invasive cells $u_i.$ Again, $D$ is the diffusion constant, $u_{M}$ is the carrying capacity for the GBM cell density, $\delta$ is a Dirac delta function, $R(t)$ is the radius of tumor core, $R = R_0+v_c t$, $R_0$ is the initial radius of core. This model is motivated by experimental observations and the images taken during the experiment were used for model validation and numerical analysis. Notice that the tumor boundary represented by $R(t)$ is artificially described by a linear function of $t$ instead of as the wave front of a traveling wave of a RDE model in most other existing models. To overcome this major shortcoming of the above model, a modified version of the Stein's model for GBM is presented by Stepien et al. in ^[17] as follows

which is a density-dependent convective-reaction-diffusion equation. Here, the logistic growth term for the number of cells, the taxis term for GBM cells are kept, but the term for the cell shedding is neglected. The diffusion is not constant but depends on density. It is shown that the model generates traveling wave solutions accurately match the experiment images reported in ^[18] without the need of assuming the tumor size growing linearly.

While logistic equation is mathematically preferred by theoretical modeling efforts due to the simplicity of its quadratic function form, other growth functions such as the Bernoulli growth function may fit experimental data equally well or even better. To this end, we consider one-dimensional convective-reaction diffusion equation of GBM proliferation and migration with a generic growth function $g(u)$. Our choice for one dimensional model is motivated by the fact that it is simpler in analysis and fast in numerical computations. Indeed, the previous study done by Stepien et al. ^[17] suggests that one-dimensional density-dependent GBM models can better fit the real data than two-population GBM model of Stein et al. ^[18].

The main purpose of this paper is to give a criteria for the existence of the traveling wave solution for density-dependent reaction-diffusion equation with a generic growth function. Our main result, Theorem 3.1, proves the existence of the traveling wave solution under the condition of

where $D(0)$ is the diffusion level when there is no cell density, $\rho$ is the intrinsic growth rate and $\nu_i$ is cell migration rate. Finally, for the numerical simulations we choose the Bernoulli's growth function, i.e., $g(u) = \rho u\left(1 -\left(\frac{u}{u_M} \right)^{\mu-1}\right), \, \, \mu > 1$, as it generalizes the logistic growth. We performed parameter sensitivity analysis for the chosen diffusion and growth functions. The results are compared with the experimental in vitro data provided in ^[18]. Numerical simulations reveal that the model is able to fit the experimental data better when $\mu$ is close to $3/2$ compared to the logistic equation, i.e., when $\mu = 2$.

2.   Model formulation

Motivated by the work of Stepien et al. ^[17], we consider the following density-dependent reaction-diffusion equation with a general reaction function.

Here $u$ represents the invasive cells of the tumor at a position $x$ and time $t, \; u_M$ is the carrying capacity, $\rho > 0$ is the intrinsic growth rate, $\nu_i > 0$ is the degree at which cells migrate. As in ^[17], inspired by the experimental study of Stein et al. ^[18], we impose following conditions on the density-dependent diffusion function:

● (A1) $D(u^*)$ is a continuous and differentiable function;

● (A2) $D(u^*)$ is a positive and decreasing function for $u^*\ge0$, i.e., $D(u^*) > 0$ and $D'(u^*) < 0$ for $u^*\ge0,$

where $u^* = \frac{u}{u_M}.$ Examples include the following diffusion functions.

where $D_1\ge D_2$, $a > 0$, $n\ge1$ so that $D(u)$ is positive.

For the growth function we know that tumor cells grow fast if $u$ is small and grow slowly as $u$ approaches its carrying capacity $u_M.$ Thus, we impose the following conditions:

● (C1) $f(0) = 1$ and $f(1) = 0;$

● (C2) $f(u^*)$ is a decreasing function for $u^*\ge0$, i.e., $f'(u^*) < 0$ for $u^*\ge0$.

There are many growth functions that satisfy these conditions. Let us define $g(u) = \rho uf\left(\frac{u}{u_M} \right)$. If we choose $f(u^*) = 1-u^*$ that is $f\left(\frac{u}{u_M} \right) = 1-\frac{u}{u_M}$ then $g(u^*) = \rho u^*\left(1-u^*\right)$ is the logistic growth function and the model (2.1) is the same as in ^[17]. We can choose $f(u^*) = 1-(u^*)^{\mu-1}$, $\mu > 1$, then $g(u^*) = \rho u_M u^*\left(1-(u^*)^{\mu-1}\right)$ is the Bernoulli growth function. One can show that all of these choices satisfy above conditions on the function $f$.

3.   Traveling wave solutions

This section is devoted to the analysis of traveling wave solutions with a focus on rigorously establishing the existence of positive traveling wave solutions.

In one-dimensional Cartesian coordinates, the Eq (2.1) has the following form.

Let us substitute

and define a new parameter

to obtain the following equivalent form of the Eq (3.1).

We are seeking to find a traveling wave solution of the Eq (3.2) of the form

where $c\ge 0$ is the speed of the traveling wave and the function $h(y)$ is defined on the whole real line $\mathbb{R}$ and satisfies the following boundary conditions:

By substituting (3.3) into (3.2) and applying the chain rule, we arrive at the following second-order ODE with boundary conditions.

By the choice of the function $D$ we have $D(h(y)) > 0$. Therefore, by dividing the Eq (3.5) by $D(h(y))$ and setting $k = dh/dy$ we reduce (3.5) to the following system of first-order ODEs.

Since $f$ is monotone and $f(1) = 0$, the system (3.6) has two biologically meaningful steady-state solutions: $(h, k) = (0, 0)$ and $(h, k) = (1, 0)$. Linearizing the system (3.6) around these steady-state solution we compute the Jacobian matrices as follows.

and

One can easily see that $\det J_2 = \frac{f'(1)}{D(1)} < 0$ since $D(1) > 0$ and due to the condition (C2) $f$ is a decreasing function. Therefore, we conclude that $(1, 0)$ is a saddle point.

Similarly, we compute $\det J_1 = \frac{1}{D(0)} > 0$ and $\text{tr}J_1 = \frac{p-c}{D(0)} < 0$ if we assume $p < c.$ In this case, $(0, 0)$ is either a stable node or a stable spiral. We must rule out the second case because oscillations are not possible in the tumor cell dynamics. This implies we must require roots from the characteristic equation of the Jacobian at $(0, 0)$ be negative. An elementary algebraic derivation leads to the following condition.

which in terms of the original parameters is equivalent to the following

Next, we aim to prove the existence of a traveling wave solution under the condition (3.7) that satisfies $\lim\limits_{y\to -\infty}(h, k) = (1, 0)$ and $\lim\limits_{y\to \infty}(h, k) = (0, 0)$. Thus, we need to construct a heteroclinic orbit connecting the two steady-state solutions. To this end, we use phase plane analysis and construct a positively invariant region to trap the unstable orbit emanating from the saddle point $(1, 0).$

Consider the line

Note that the horizontal line $k = 0$ is $h-$nullcline of the system (3.6). Moreover, on $L_1$ we have

This implies that $k' < 0$ on $L_1$ since $f(h) > 0$ for $h\in (0, 1).$ Thus, the direction fields on the line $L_1$ is in the negative $k$ direction.

Next, let us consider the line

where $\lambda_1$ is the more negative eigenvalue of $J_1,$ i.e.,

One can verify that $\lambda_1$ satisfies the following equation

The normal vector to the line $L_2$ in the positive $k$ direction is given by $(-\lambda_1, 1).$ Along $L_2$ the system (3.6) is equivalent to

Thus, the inner product between the normal vector to $L_2$ and the vector field gives us

Due to the relation (3.11) we have

Note that by monotonicity of the diffusion function we have $0 < D(h) < D(0)$ and $D'(h) < 0$ for $h\in (0, 1].$ Thus, we have $\frac{h}{D(h)}\ge0$ and $\lambda_1^2[D(h)-D(0)+D'(h)h]\le0.$ Further, since $f(0) = 1$ and $f$ is a decreasing function we have $f(h)-1\le0$ for $h\in (0, 1].$ Therefore, the inner product satisfies

The last inequality means that the direction fields and the normal vector make an acute angle.

Finally, we define the line

It is easy to see that the direction fields across $L_3$ is in the negative $h$ direction. Thus, the triangular region $L$ bounded by $L_1, \; L_2$ and $L_3$ is positively invariant.

Thus, we summarize our discussion in the following theorem.

Theorem 3.1. There exists a traveling wave solution in the form (3.3) of the convention-reaction-diffusion Eq (3.2) which satisfy the boundary conditions (3.4) i.e., $u(x, t)\to 1$ as $x\to -\infty$ and $u(x, t)\to 0$ as $x\to \infty$ with $0 < u(x, t) < 1,$ whose orbit connects the steady states $u = 0$ and $u = 1$ if and only if (3.7) is satisfied.

Proof. Let us show that the unstable separatrix that is leaving the saddle point $(1, 0)$ has intersection with $L.$ To this end, consider $k-$nullclines which satisfy

Since the Eq (3.14) is a quadratic equation the solutions can be computed as

By monotonicity of $f$ we have $f'(1) < 0.$ Thus, the slope of the nullcline $k_1 (h)$ at $h = 1, \, k = 0$ is given by

On the other hand, we can show the eigenvector of $J_2$ corresponding to positive eigenvalue can be chosen as

It is easy to see that

Thus,

Therefore, we conclude that the slope of the eigenvector $v$ is less than $k_{1}'(1).$ Note that the trajectory (unstable manifold) that leaves $(1, 0)$ in the negative $k$ direction has the tangent vector $v$ at $(1, 0).$ Thus, these trajectories leave $(1, 0)$ above $k_1.$ Further, we know that near $h = 1, \; k_1(h)$ is contained in $L$ and the direction fields on the curve $k_1$ is horizontal in the negative $h$ direction. Hence, we have that the unstable separatrix that is leaving the saddle point $(1, 0)$ remains in the region $L.$ Moreover, the $\omega-$limit set of this orbit is also in $L.$ Since $h' = k \le 0$ in the region $L$ and there is no steady state solutions in the interior of $L$, Poincaré$-$Bendixson theorem implies that there is no periodic solution in the interior of $L$. Thus, the $\omega-$limit set cannot be in the interior of $L$ and it must be on the boundary of $L.$ So, one can see that the $\omega-$limit set is $(0, 0)$. Thus, under the condition (3.7) we conclude that there exists a heteroclinic orbit connecting $(0, 0)$ and $(1, 0)$. This finalized the proof as it proves the existence of a traveling wave solution.

Remark 3.1. Since Eq (2.1) is more general then the one considered in ^[17], Theorem 3.1 is an extension of the main result (Theorem 3.1) in ^[17].

4.   Numerical results

In this section, we show numerical results of the model, compare it with experimental data obtained in ^[18] and provide optimal model parameters. We focus on the U87WT cell line for our model validation. The experimental data was obtained using GRABIT ^[35] from the experimental work ^[18].

Following the paper ^[17], we have used the following diffusion function,

where the constants are chosen so that $D_1 \ge D_2, \; n\ge1,$ and $a > 0$ to satisfy the conditions $(A1)-(A2).$

The partial differential Eq (3.2) describes normalized concentration $u(x, t)$ which simulates migration of the invasive cell. Since solution vanishes at the boundary (front of the traveling wave solution) and tumor is symmetric with respect to tumor center, we need only to consider half of domain from $x = 0$ to $x = 1$. In other words, $u(x, t) = 0$ for $x = 1$ cm and $D(u)\frac{\partial u(x, t)}{\partial x} = 0$ for $x = 0$ cm. Time step is $0.01$ and $\Delta x = 1/1000$. Equation (3.2) with $f(u) = 1-u^{\mu-1}$ was solved numerically by using pdepe tool in Matlab. We note that $\rho$ can be neglected in sensitivity analysis since it vanishes in the normalized PDE (3.2).

We have used the relative error estimation from ^[17].

where $N$ is total number of days ($N = 7$), $M$ is total number of cell density data points at day 3 ($M = 17$) and $q$ is number of optimized parameters ($q = 6$).

4.1.   Parameter estimation

We have conducted numerical experiments on finding optimal parameters $(D_1, D_2, a, n, \nu_i, \mu)$ using fminsearch ^[36] function in Matlab. The objective function was chosen to be the relative error, Eq (4.2). fminsearch is an optimization algorithm based on the Nelder-Mead method. This method depends on initial condition and derivative of objective function. After various numerical experiments with different initial conditions we notice that the optimal value for $\mu$ close to $3/2$. Thus, we concentrate on the sensitivity of five parameters $(D_1, D_2, a, n, \nu_i)$. The results are illustrated in Table 1. From Figure 2 it can be seen that when $\mu = 1.4907$, the model is able to fit the experimental data in tumor cell density and invasive radius better compare to the case with $\mu = 2$.

4.2.   Wave speed computation

Figure 3 shows numerical solutions of our model with optimal parameters for $5$, $10$, $15$, $20$, $25$ and $30$ days. Front of tumor moves linearly as time of simulation increases. The computation of wave speed is performed with optimized parameter values in Table 1. We find $x$ location such that the cell density is last over $4.2 \times 10^2$ cells/cm$^3$ (very small in clinical applications). Based on these $x$ values we fit to time $T$ linearly. Numerical values of slope of the line is $k = 0.0252$ cm/day. Equation 3.7 gives $k_{min} = 0.025$ cm/day from the optimal values, see Table 1.

5.   Discussions

In this study, we considered a density-dependent reaction-diffusion equation with a general growth function to describe GBM model governed by the Eq (2.1). The diffusion function is assumed to be differentiable (the condition (A1)), positive and decreasing functions (the condition (A2)) due to experimental study in ^[18]. A general growth function satisfies the conditions specified as (C1) and (C2) to capture both the logistic and Bernoulli growth functions. In the theoretical part, we carried out traveling wave analysis of the model with density dependent diffusion and general growth functions and found that the minimum speed of a traveling wave satisfies

This result agrees with the traveling wave analysis of Stepien et al. ^[17] which was carried out in the special case of the model (2.1) with the logistic growth function, i.e., $f(u^*) = 1-u^*$, and the following diffusion function

Thus, we extended the main result of Stepien et al. ^[17] since model (2.1) is more general than in ^[17].

In the numerical part, we considered the same diffusion function as in Stepien et al. ^[17], that is, the Eq (4.1), since in vitro experimental results studied by Stein et al. ^[18] suggests that diffusion is inversely proportional to the cell density, i.e., diffusion is larger in areas where the cell density is smaller, but diffusion is smaller where the cell density is larger. This could possibly be explained by cell–cell mutual interference in high density area and cell–cell adhesion in medium density regions ^[17,37]. The numerical findings suggest that the numerical values of the diffusion parameters ($D_1$ and $D_2$) are of the same order of magnitude and is much larger compare to the advection parameter in the model equation. In contrast, in the previous study ^[17], there was significant difference between these parameters $D_1$ and $D_2$ and $D_1 < \nu$. Model (2.1) is able to improve previous fittings to the realistic data sets, which were digitized in the one-dimensional setting. The main reasons to choose one dimensional model include its ability of capturing the proliferating tumor core and the invading migratory cells in a single equation and providing simplicity in analysis and faster convergence in numerical computations. Naturally, further study of two- or three-dimensional model with multiple cell populations would be of interest ^[9,16,21]. Another possible direction for the future is to consider different diffusion functions that satisfy the conditions (A1)–(A2) and different growth functions such as Gompertz and von Bertalanffy.

Acknowledgements

The authors would like to thank Professor Erica Rutter for providing Matlab codes of the previous study ^[17]. AK was supported in part by Nazarbayev University FDCRG No. 090118FD5353. YA would like to acknowledge support from the Nazarbayev University FDCRG No. 110119FD4502 the Ministry of Education and Science of the Republic of Kazakhstan Grant No. AP08052762. BS is partially supported by the Ministry of Education and Science of the Republic of Kazakhstan Grant No. АР08052345. YK is partially supported by NSF grants DMS-1615879, DEB-1930728 and an NIH grant 5R01GM131405-02.

Conflict of interest

No potential conflict of interest was reported by the authors.