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

: Mathematical modeling for cancerous disease has attracted increasing attention from the researchers around the world. Being an e ﬀ ective tool, it helps to describe the processes that happen to the tumour as the diverse treatment scenarios. In this paper, a density-dependent reaction-di ﬀ usion equation is applied to the most aggressive type of brain cancer, Glioblastoma multiforme. The model contains the terms responsible for the growth, migration and proliferation of the malignant tumour. The traveling wave solution used is justiﬁed by stability analysis. Numerical simulation of the model is provided and the results are compared with the experimental data obtained from the reference papers.


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 aggressive 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 "goor-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 ≥ c 0 = √ 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)) → w(x) for t → ∞ and m(t) taken appropriately, and m(t) → 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 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) = φ(x + ct) is also a wave solution with the boundary conditions φ(−∞) = 0, φ(∞) = 1. If we consider the radius instead of x then the Fisher's equation in radially asymmetric disk case looks as follows (1. 2) The term 1 r ∂u ∂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, ρ 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, δ 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 convectivereaction 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 onedimensional 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, ρ is the intrinsic growth rate and ν i is cell migration rate. Finally, for the numerical simulations we choose the Bernoulli's growth function, i.e., g(u) = ρu 1 − u u M µ−1 , µ > 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 µ is close to 3/2 compared to the logistic equation, i.e., when µ = 2.

Model formulation
Motivated by the work of Stepien et al. [17], we consider the following density-dependent reactiondiffusion equation with a general reaction function. (2.1) Here u represents the invasive cells of the tumor at a position x and time t, u M is the carrying capacity, ρ > 0 is the intrinsic growth rate, ν 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:
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: There are many growth functions that satisfy these conditions. Let us define g is the logistic growth function and the model (2.1) is the same as in [17]. We can choose f (u * ) = 1 − (u * ) µ−1 , µ > 1, then g(u * ) = ρu M u * 1 − (u * ) µ−1 is the Bernoulli growth function. One can show that all of these choices satisfy above conditions on the function f .

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.
and define a new parameter p = ν i √ ρ 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 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.
One can easily see that det J 2 = 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 = 1 D(0) > 0 and trJ 1 = 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.
c ≥ c min = 2 D(0) + p, 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 ∈ (0, 1). Thus, the direction fields on the line L 1 is in the negative k direction. Next, let us consider the line where λ 1 is the more negative eigenvalue of J 1 , i.e., One can verify that λ 1 satisfies the following equation The normal vector to the line L 2 in the positive k direction is given by (−λ 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 ∈ (0, 1]. Thus, we have h D(h) ≥ 0 and λ 2 Further, since f (0) = 1 and f is a decreasing function we have f (h) − 1 ≤ 0 for h ∈ (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.  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, 1 2 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 ω−limit set of this orbit is also in L. Since h = k ≤ 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 ω−limit set cannot be in the interior of L and it must be on the boundary of L. So, one can see that the ω−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].

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, D(u) = D 1 − D 2 u n a n + u n , (4.1) where the constants are chosen so that D 1 ≥ D 2 , n ≥ 1, 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) ∂u(x,t) ∂x = 0 for x = 0 cm. Time step is 0.01 and ∆x = 1/1000. Equation (3.2) with f (u) = 1 − u µ−1 was solved numerically by using pdepe tool in Matlab. We note that ρ 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).

Parameter estimation
We have conducted numerical experiments on finding optimal parameters (D 1 , D 2 , a, n, ν i , µ) 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 µ close to 3/2. Thus, we concentrate on the sensitivity of five parameters (D 1 , D 2 , a, n, ν i ). The results are illustrated in Table 1. From Figure 2 it can be seen that when µ = 1.4907, the model is able to fit the experimental data in tumor cell density and invasive radius better compare to the case with µ = 2.   Table 1. We find x location such that the cell density is last over 4.2 × 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.

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 c min = 2 D(0)ρ + ν i . 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 D(u) = D 1 − D 2 u n a n + u n . 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 < ν. 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.