Travelling wave solutions of the reaction-diffusion mathematical model of glioblastoma growth: An Abel equation based approach

We consider quasi-stationary (travelling wave type) solutions to a nonlinear reaction-diffusion equation with arbitrary, autonomous coefficients, describing the evolution of glioblastomas, aggressive primary brain tumors that are characterized by extensive infiltration into the brain and are highly resistant to treatment. The second order nonlinear equation describing the glioblastoma growth through travelling waves can be reduced to a first order Abel type equation. By using the integrability conditions for the Abel equation several classes of exact travelling wave solutions of the general reaction-diffusion equation that describes glioblastoma growth are obtained, corresponding to different forms of the product of the diffusion and reaction functions. The solutions are obtained by using the Chiellini lemma and the Lemke transformation, respectively, and the corresponding equations represent generalizations of the classical Fisher--Kolmogorov equation. The biological implications of two classes of solutions are also investigated by using both numerical and semi-analytical methods for realistic values of the biological parameters.

1. Introduction.Glioblastoma, the most common central nervous system (CNS) tumor, accounting for 50% of the 17,000 primary brain tumors diagnosed annually in the US [1], is also the most malignant form of brain cancer, having an extremely poor outcome [2].For the most aggressive grade of gliomas, known as glioblastoma multiforme (GBM), the life expectancies are from 6 to 12 months [3].Amongst patients treated with surgery and a radiation-containing regimen, median survival was 12.0 months in the period 2000-2003, and 14.2 months in 2005-2008, respectively.In the temozolomide era, median survival times varied from a high of 31.9 months, for patients in the age group 20-29 years, to a low of 5.6 months in patients age 80 years, and older [4].
One factor that makes glioblastoma extremely difficult to treat is its high invasiveness, enabling tumor cells to disperse from the main tumor mass into the surrounding normal brain.Glioblastoma is highly diffuse, and can invade a large portion of the cerebral cortex in a short period of time, making complete surgical excision impossible, so that dispersed glioma cells are out of reach of surgery, radiation, and chemotherapy, so that recurrence becomes inevitable [5].There are many factors determining the prognosis for patients with gliomas, like the histologic type, the grade of malignancy, the patient's age and the level of neurological functioning, respectively [6].The grade of malignancy for glioblastoma includes at least two factors, the net proliferation rate, and the invasiveness, respectively, which are estimated histologically (for the World Health Organization (WHO) classification and grading of brain tumors see [7]).However, a practical accurate definition of these factors is still missing [8].Unlike solid tumors, for which simple exponential or geometric expansion represents expansion of volume (equivalent to the number of cells in the tumor), gliomas consist of motile cells that can migrate as well as proliferate [8].Indeed, the invasiveness makes it almost impossible to define the growth rate as a classical volume-doubling time (DT) [8].DT is a parameter widely used for the measurement of the tumor growth rate, which can also be quantified by the specific growth rate (SGR), defined as follows.If the tumor volume V is measured at times t 1 and t 2 , then SGR can be obtained as SGR = ln (V 2 /V 1 ) / (t 2 − t 1 ) [9].DT is related to SGR by the relation DT = ln 2/SGR [9].
Cancer research has been a fertile ground for developing mathematical models that can describe the proliferation of malignant cells.The early models of cell proliferation were based on a simple exponential growth of solid (usually benign) tumors, so that n(t) = N 0 exp(λt), where n is the number of cells at time t, N 0 is the initial cell number, and λ is a constant [10].Another model used to describe tumor dynamics is based on the Gompertz curve, a mathematical model for a time series, where growth is slowest at the end of a time period, where N ∞ is the plateau cell number, which is reached at large values of the time.
The parameter b is related to the initial tumor growth rate [11].An alternative model describing the time variation of a mass m of any organism, including solid tumors, is given by [12] dm dt = am 3/4 − bm, where a, b = constant.However, such models do not take into account the spatial arrangement and distribution of the cells at a specific anatomical location, or the spatial spread of the cancerous cells.These spatial aspects are crucial in estimating tumor growth, since they determine the invasiveness of the tumor and the sharpness of the apparent border of the tumor [13].Explicitly taking into account the extent of infiltration of the tumor is necessary in different situations, like, for example, in estimating the benefit of surgical resection.A simple mathematical model for the proliferation and infiltration of the glioma cells was introduced in [13].The model was based in part on quantitative image analysis of histological sections of a human brain glioma and especially on crosssectional area/volume measurements of serial Computer Tomograph (CT) images while the patient was undergoing chemotherapy.In its general form, from a mathematical point of view the model represents a reaction-diffusion system, with the growth rate and the diffusion rate representing the key model parameters.An extension of the mathematical model based on proliferation and infiltration of neoplastic cells introduced in [13] that allows predictions to be made concerning the life expectancies following various extents of surgical resection of gliomas of all grades of malignancy was considered in [14].Numerical simulations using the model allow to estimate what would happen to patients if various extents of surgical resection, rather than chemotherapies, would have been used.It has been shown that the shell of the infiltrating tumor that remains after 'gross total removal' or even a maximal excision continues to grow and regenerates the tumor mass remarkably rapidly [14].The model initially introduced in [12,13] has been extended and used in different clinical situations for the mathematical study of glioma proliferation and invasion (for very informative reviews see [8] and [15]) and to include the effects of radiotherapy [16], where an extension to Swanson's reaction-diffusion model to include the effects of radiation therapy using the classic linear-quadratic radiobiological model was presented.Moreover, in [15] it was shown that the defining and essential characteristics of gliomas in terms of net rates of proliferation and invasion can be determined from serial Magnetic Resonance Images (MRIs) of individual patients.To gain some insight into glioblastoma invasion, experiments were conducted on the patterns of growth and dispersion of U87 glioblastoma tumor spheroids in a threedimensional collagen gel in [17].A continuum mathematical model of the dispersion behaviors was developed.The mathematical model quantitatively reproduces the experimental data.
Chemotherapy in a spatial model of tumor growth was considered in [18].The model, which is of reaction-diffusion type, takes into account the complex interactions between the tumor and surrounding stromal cells by including densities of endothelial cells and the extra-cellular matrix.When no treatment is applied the model reproduces the typical dynamics of early tumor growth.A mathematical model using a realistic three-dimensional brain geometry, and which and considers migrating and proliferating cells as separate classes was analyzed in [19].Several mechanisms for infiltrative migration were considered, and methods for simulating surgical resection, radiotherapy and chemotherapy were developed.It was shown that the model provides clinically realistic predictions of tumor growth and recurrence following therapeutic intervention.An important aspect of the biological modeling of glioblastoma is the mathematical handling of boundary conditions.An explicit and thorough numerical formulation of the adiabatic Neumann boundary conditions imposed by the skull on the diffusive growth of gliomas and in particular on glioblastoma multiforme was considered in [20].A detailed exposition of the numerical solution process for a homogeneous approximation of glioma invasion using the Crank-Nicolson technique in conjunction with the Conjugate Gradient system solver was also provided.
An individual-based stochastic model that analyses how the phenotypic switching between proliferative and migratory states of individual cells affects the macroscopic growth of the tumour was proposed in [21].The glioblastoma cells are either in a proliferative state, where they are stationary and divide, or in motile state, in which they are subject to random motion.This model may find some clinical applications for designing relevant cell screens for glioblastoma and cytometry-based patient prognostic.
The simple reaction-diffusion model of gliomas [13,14] predicts that the "edge" of the detectable glioma should behave as a travelling wave and follow Fisher's approximation that the velocity equals twice the square root of the product of the growth rate and of the diffusion coefficient.This predictions has been confirmed by using data derived from medical imaging techniques [22].The implication of this result is that the prognosis of any individual patient can be predicted if two sets of scans, separated in time so that a significant increase in growth can be measured, can be obtained before treatment is begun [22].
Reaction-diffusion equations usually do admit travelling wave solutions [23,24].It has been shown in [25] that wave solutions, which correspond to moving bands of concentration do exist in the system of equations describing the Belousov-Zhabotinskii reaction.Various relevant properties of the solutions have also been obtained.With parameter values obtained from experiment, numerical results have been given for the travelling wave solutions.A semi-inverse method, which renders exact static solutions of one-component, one-dimensional reaction-diffusion equations with variable diffusion coefficient D, requiring at most qualitative information on the spatial dependence of the latter was introduced in [26].Through a simple ansatz the reaction-diffusion equations can be mapped onto the stationary Schrödinger equations, having the form of the potential still at our disposal.A new substitution was used in [27] to reduce the problem of a quasi-stationary solution to a nonlinear reaction-diffusion equation with arbitrary, autonomous coefficients to either a linear ordinary differential equation, or the first order first kind Abel differential equation of the form where f i (x), i = 0, 1, 2, 3, are arbitrary real functions of x, defined on a real interval [28].The second kind Abel differential equation is defined as where g i (x), i = 0, 1 are real functions of x.The second kind Abel differential equation can be reduced to the first kind Abel equation by means of the substitution g 0 (x) + g 1 (x)y = 1/z [28].Exact periodical and solitary wave solutions were obtained by solving the corresponding Abel equations [27].In particular, it was shown that the problem of the kinetics of thin film growth (or of wire-like nanostructures) has bounded solutions in terms of elliptic functions.A direct method to obtain travelling-wave solutions of some nonlinear partial differential equations expressed in terms of solutions of the Abel differential equation of the first type with constant coefficients was proposed in [29].Exact solutions to the modified Benjamin, Bona, and Mahony (BBM) equation by viscosity were found in [30], by including the effect of a small dissipation on waves.Using Lyapunov functions and dynamical systems theory, it was proven that when viscosity is added to the BBM equation, in certain regions there still exist bounded travelling wave solutions in the form of solitary waves, periodic, and elliptic functions.
Travelling-wave solutions of different mathematical model describing the growth of tumors have been considered in several publications.Spreading cell fronts play an essential role in many physiological and biological processes.Classically, models of this process are based on the Fisher -Kolmogorov -Petrovsky -Piscounov equation (Fisher -Kolmogorov equation for short) [31,32,33]; however, such continuum representations are not always suitable as they do not explicitly represent behaviour at the level of individual cells.Additionally, many models examine only the large time asymptotic behaviour, where a travelling wave front with a constant speed has been established [34].
A travelling-wave analysis of a mathematical model describing the growth of a solid tumor in the presence of an immune system response was analyzed in [35].From a modelling perspective, attention was focused upon the attack of tumor cells by Tumor Infiltrating Cytotoxic Lymphocytes (TICLs), in a small multicellular tumor, without necrosis and at some stage prior to (tumor-induced) angiogenesis.The existence of travelling-wave solutions for the system was established.
A lattice-gas cellular automaton model of tumor cell proliferation, necrosis and tumor cell migration was introduced in [36], with the main aim of predicting the velocity of the traveling invasion front, which depends upon fluctuations that arise from the motion of the discrete cells at the front.An analytical estimate of the velocity was derived in the cut-off mean-field approximation via the discrete Lattice Boltzmann equation and its linearization.The front velocity scales with the square root of the product of probabilities for mitosis and the migration coefficient, while the width of the traveling front is found to be proportional to its velocity.In [37] it was shown, with the help of a simple growth model, that the short time required for the recurrence of a glioblastoma multiforme tumour after a gross total resection cannot be explained solely by a mutation-based theory.It was proposed that the transition to invasive tumour phenotypes can be explained on the basis of the microscopic Go or Grow mechanism (migration/proliferation dichotomy) and the oxygen shortage, i.e. hypoxia, in the environment of a growing tumour.This hypothesis was tested with the help of the lattice-gas cellular automaton model [36].Possible therapies that could help prevent the progression towards malignancy and invasiveness of benign tumours were also suggested.A mathematical model that incorporates the interplay among two tumor cell phenotypes, a necrotic core and the oxygen distribution, and which reveals the formation of a traveling wave of tumor cells was analyzed in [38].The model reproduces the observed histologic patterns of pseudopalisades.The simulations of the model equations also show that preventing the collapse of tumor microvessels leads to slower glioma invasion.
A cell-based model of glioblastoma growth, which is based on the assumption that the cancer cells switch phenotypes between a proliferative and motile state, was analyzed in [39].The dynamics of this model can be described by a system of partial differential equations, which exhibits travelling wave solutions whose wave speed depends crucially on the rates of phenotypic switching.Under certain conditions on the model parameters, a closed form expression of the wave speed can be obtained.By using singular perturbation methods an approximate expression of the wave front shape can be derived.
A simple free boundary model formed of a Hele-Shaw equation for the cell number density coupled to a diffusion equation for a nutrient was studied in [40].In this model a travelling wave solution does exist, with a healthy region separated from the progressing tumor by a sharp front (the free boundary), while the transition to the necrotic core is smoother.The pressure distribution vanishes at the boundary of the proliferative rim, with a vanishing derivative at the transition point to the necrotic core.
The Fisher-Kolmogorov equation belongs to a more general class of reactiondiffusion equations, given by [23,24]  It is the goal of the present paper to study exact traveling wave solutions in the one-dimensional reaction-diffusion models of glioblastoma tumor growth introduced in [13]- [16].More exactly, we will consider the possibility of the description of the tumor growth by a general reaction-diffusion system, containing two arbitrary tumor concentration dependent functions, called the diffusion function and the reaction function, respectively.By considering travelling wave solutions of the reactiondiffusion equation, we show that the second order non-linear differential equation describing wave propagation can be reduced to a first kind non-linear Abel equation.Note that the mathematical properties of the Abel equation and its applications have been intensively investigated in [41]- [48].
By using some standard integrability conditions of the Abel equation, we obtain several classes of exact solutions of the tumor growth equation.More exactly, the first class of solutions is obtained by using the integrability condition of Chiellini [49,50], which can be formulated as a differential condition relating the diffusion and the reaction functions.The use of this condition allows to obtain exact travelling wave solutions to some reaction-diffusion equations representing a generalization of the Fisher-Kolmogorov equation.The second class of solutions is obtained by transforming the Abel equation to a second order equation [50,51], which allows the construction of exact travelling wave solutions for several choices of the diffusion and reaction functions.The biological implications of two tumor growth models by travelling wave propagation, both representing generalizations of the Fisher-Kolmogorov model, are investigated in detail by using some realistic numerical values for the free parameters of the model.
The present paper is organized as follows.In Section 2 we briefly review the basic mathematical model of glioblastoma growth, and the reduction of the general diffusion-reaction equation to an Abel equation is presented.Exact solutions of the glioblastoma growth equation based on the Chiellini lemma are presented in Section 3. Further integrability cases of the tumor growth equation by diffusion are obtained in Section 4. A number of exact travelling wave solutions of the tumor growth equation, corresponding to different functional forms of the product of the diffusion and reaction functions are presented in Section 5.The biological implications of our models are briefly investigated in Section 6.We discuss and conclude our results in Section 7.
2. The mathematical model of tumor growth.The first model of the growth of an infiltrating glioma as a reaction-diffusion system was initially formulated as a conservation equation [52].From a phenomenological point of view, the model can be formulated in words as [13]- [16], [52]: "The rate of change of tumor cell population = the diffusion (motility) of tumor cells + the net proliferation of tumor cells." Mathematically, the growth of the gliomas can be described as a diffusion equation for the tumor cell density c ( r, t) at the position r and time t, where J is the cell flux, and ρ denotes the net proliferation rate.By assuming that the flux obeys the standard Fick law, J = D∇c, where D is the diffusion coefficient representing the active motility of glioma cells, Eq. ( 7) can be written as a reaction-diffusion equation, The model formulation is completed by boundary conditions, which impose no migration of cells beyond the brain boundary, and initial conditions c ( r, 0) = f ( r), where f ( r) defines the initial spatial distribution of malignant cells.As boundary conditions, it is required that there is no flux of cells to the outside of the brain, or into the ventricles, so at the boundaries of the two-dimensional domain n • ∇c = 0, where n is the unit vector normal to the boundary.However, more general models can also be formulated.In [13] and [14] glioblastoma modelling was done by using the Fisher -Kolmogorov equation [31,32,33], which describes randomly moving cells, and which simultaneously divide at rate ρ.
Cells throughout the tumor are assumed to proliferate at a constant rate ρ until they reach a limiting density c max .The Fisher equation exhibits travelling wave solutions, which from biological point of view describe a tumor invading the healthy tissue, with the velocity of the invading front given by v = 2 √ Dρ [39].The inclusion of the chemotherapy in the model leads to a modification of the free evolution of the cells according to the phenomenological prescription [13]: "The rate of change of tumor cell population = diffusion (motility) of tumor cells + net proliferation of tumor cells -loss of tumor cells due to chemotherapy." Mathematically, the simplest model of glioblastoma evolution in the presence of chemotherapy is formulated as where G(t) is a function describing the effects of chemotherapy.When chemotherapy is administered, G(t) is constant, and G(t) = 0 otherwise.The effect of the radiotherapy can be modelled as where R ( r, t) represents the effect of external beam radiation therapy at location r and time t [16].
In the following we propose to model the evolution of gliomas by means of 1+1 dimensional general reaction -diffusion system of the form where the dissipative (diffusion) function D(c) = 0 and the reaction term Q(c) = 0 both depend explicitly upon the cell density c only, and not on space x and time t variables.By introducing a phase variable where V f ≥ 0 is a constant wave velocity, Eq. ( 12) takes the form of a second order non-linear differential equation of the form where Eq. ( 14) must be considered together with the initial conditions c(0) = c 0 and (dc/dξ)| ξ=0 = c ′ 0 , respectively.By means of the transformations Eq. ( 14) can be transformed to the general form of the first order first kind Abel equation, given by Eq. ( 19) must be integrated with the initial condition v allows us to write Eq. ( 19) in the standard form of the Abel equation, which must be solved with the initial condition 3. The Chiellini integrability condition for general reaction-diffusion systems.An exact integrability condition for the Abel equation Eq. ( 21) was obtained by Chiellini [49] (see also [48] and [50]), and can be formulated as the Chiellini Lemma as follows: Chiellini Lemma.If the coefficients F (x) and L(x) of a first kind Abel type differential equation of the form where k = constant = 0, then the Abel Eq. ( 19) can be exactly integrated.
As applied to the Abel Eq. ( 21), the Chiellini lemma states that if the coefficients of the equation satisfy the condition the Abel equation is exactly integrable.In this case the Abel Eq. ( 21) takes the form With the help of the substitution Eq. ( 26) becomes which must be integrated with the initial condition Eq. ( 28) has the general solution where C −1 is an arbitrary constant of integration, and respectively.In order to obtain the explicit form of the function H(W, k) we have used the algebraic identities where a and b are arbitrary constants, with the particular condition a = b = 1, and 4 respectively.Therefore Eq. ( 26) has the general solution while for v we obtain Therefore, by using the Chiellini Lemma, we have obtained the following general solution for the general reaction-diffusion equation Eq. ( 12): Theorem 1.If the diffusion function D(c) and the reaction function Q(c) in a general one-dimensional reaction-diffusion equation satisfy the condition given by Eq. ( 25), then the equation admits exact travelling wave solutions, expressed in parametric form as where ξ 0 is an arbitrary constant of integration, and we have taken W as a parameter.
The constant k can be determined once the explicit functional dependence of D and Q on c is known as while the integration constant C is determined as For The integrability condition given by Eq. ( 25) fixes the constant k as k = ρD 0 /V 2 f > 1/4.Then from Eq. ( 37) we obtain the function ξ(W ) as giving where we have denoted ∆ = V 2 f − 4D 0 ρ.The general solution of Eq. ( 41) is given by where C n , n = +, − are arbitrary constants of integration.Note that the general solution of Eq. ( 41) as given by Eq. ( 44) can be obtained directly from Eq. ( 38), which is a linear second order homogeneous differential equation.

Travelling wave solutions for
As a second example of application of Theorem 1 we consider that the diffusion and the reaction functions are given by The travelling wave equation for the reaction-diffusion system with these forms of D and Q is given by Eq. ( 46) represents the travelling wave form for a generalized Fisher equation given by Since the functions D(c) and Q(c) satisfy the integrability condition given by Eq. ( 25), with D 0 ρ = kV 2 f , the general solution of Eq. ( 46) can be obtained in an exact parametric form.

3.2.1.
The case k = 1/4.We consider first the case k = 1/4.Then the general solution of Eq. ( 46) is given in parametric form as where C is an arbitrary constant of integration.

3.2.2.
The case k = 1/4.For k = 1/4, we obtain where the constant k is determined by the model parameters V f , D 0 and ρ.

3.3.
Travelling wave solutions of the first generalized Fisher-Kolmogorov equation.Exact travelling wave solutions of more general equations of the form where D 0 , α and ρ are arbitrary constants, are given, in a parametric form, by the following equations, and respectively, with k = ρD 0 /V 2 f .In the following we will call Eq. ( 54) the first generalized Fisher-Kolmogorov equation.

4.
Further exact travelling wave solutions of the general reaction-diffusion equation.In the following we introduce first a new variable θ(c) through the Lemke transformation, defined as Therefore Eq. ( 21) becomes By taking into account the differential identity Eq. ( 58) takes the form [50, 51] Eq. ( 60) must be integrated with the initial conditions c (θ 0 ) = c 0 and Therefore we have obtained the following Theorem 2. If a solution c = c(θ) of Eq. ( 60) is known, then the general reaction-diffusion equation Eq. ( 14) admits travelling wave solutions, given in parametric form as where ξ 0 is an arbitrary constant of integration.
If the solution of Eq. ( 60) is obtained in a parametric form, with parameter τ , then the general travelling wave solution of the reaction -diffusion equation Eq. ( 14) is obtained as 5. Travelling wave solutions of the tumor growth equation.Depending on the functional form of the product of the diffusion and reaction functions D(c) and Q(c), Eq. ( 60) can be integrated exactly in a number of cases, thus leading, with the help of Theorem 2, to exact travelling wave solutions of the general diffusionreaction Eq. ( 14), which we present in the following.

5.1.
Travelling wave solutions for D(c)Q(c) = constant.In the case the arbitrary function subsequently Eq. ( 60) has the general solution where C 1 and C 2 are arbitrary constants of integration.Hence for w we obtain Eq. (68) identically satisfies Eq. ( 21), and therefore we can take the arbitrary integration constant C 1 as zero, C 1 = 0. Using the transformations v (θ) = D (θ) w (θ), dc/dξ = 1/v (θ), and with the help of Eq. (68), we obtain the following expression for ξ, giving, together with Eq. (67), the general solution of the general reactiondiffusion Eq. ( 12) with the coefficients D(c) and Q(c) satisfying the condition (66) as, where ξ 1 is an arbitrary constant of integration, and we have taken θ as a parameter.Eqs. ( 67) and (69) give the general solution of the general reaction-diffusion equation Eq. ( 14) with diffusion and reaction functions satisfying the condition (66).In order to determine the integration constants C 1 and C 2 , we use the initial conditions that give and giving and respectively.Eqs.(73) determines the value of θ 0 as a function of the initial conditions (c 0 , c ′ 0 ) and the free parameters {α, V f , D (c 0 )}.

5.2.
Travelling wave solutions with a linear dependence of D(c)Q(c).Next, in order to obtain another exact general solution of Eq. ( 60), we assume that the diffusion function D (c) and the reaction function Q (c) satisfy the condition where β is an arbitrary constant.By inserting Eq. (74) into Eq.( 60), the latter equation becomes Eq. ( 75) can be integrated to yield where C 3 and C 4 are arbitrary constants of integration, and The function v (θ) is given by Thus we obtain where ξ 2 is an arbitrary constant of integration, and we have taken θ as a parameter.Thus the general solution of the reaction -diffusion equation Eq. ( 14) with diffusion and reaction functions satisfying the condition given by Eq. ( 74) is given, in a parametric form, by Eqs.(76) and Eq.(79), respectively.The numerical values of the integration constants C 3 and C 4 are determined from the initial conditions as and respectively.
If the arbitrary constant α vanishes, then we obtain the condition D(c)Q(c) = βc given by Eq. ( 25), thus we regain the solution obtained in the previous Section by using the Chiellini lemma.
where K is a constant, Eq. ( 60) takes the form of an Emden-Fowler equation [53], Eq. ( 83) has the exact solution, given in parametric form, where C 5 and C 6 are arbitrary constants of integration, and erf(z) = (2/ √ π) z 0 e −t 2 dt is the error function, representing the integral of the Gaussian distribution.For the parametric dependence of ξ we obtain Eqs. ( 84) and (86) give the general solution of the reaction -diffusion equation with diffusion and reaction functions satisfying the condition (82).

Travelling waves solutions with a power law dependence of D(c)Q(c).
If D(c) and Q(c) satisfy the relation then the basic equation determining the travelling wave solutions of the corresponding general reaction-diffusion equation is and it has the exact parametric solution given by [53] θ and C 7 and C 8 are arbitrary constants of integration.The parametric dependence of ξ is obtained as 90) and (92) give the general solution of the general reaction-diffusion equation with diffusion and reaction functions satisfying the condition (87).
6. Biological applications.In the present Section we briefly point out some possibilities of biological applications of the obtained results to model the growth of glioblastomas.In order to obtain some numerical results we need to chose some values of the parameters D and ρ that characterize the tumor dynamics.These parameters can be computed observationally from as few as two pre-treatment MRI observations [16], and current data from 29 tumors show a range of 6 -324 mm 2 /year for D and 1 -32 /year for ρ.However, in the following, we describe the virtual tumor with parameter values as follows: D = 1.43 cm 2 /year, and ρ = 16.25/year, which serve as representative means of published ranges [16].The quantity c represent the number of tumor cells within the volume V .For the sake of simplicity we fix the initial concentration of tumoral cells as c(0) = c 0 = 1000 cells/cm 3 [54], and we take (dc/dξ) | ξ=0 ≈ 7 × 10 8 cells/cm 4 .To estimate c max , we assume that the volume of a typical cell is 1200 µm 3 , as it is for EMT6/Ro tumor cells [55].Assuming that half the volume of the spheroid is made up of tumor cells, the maximum density is c max = 4.2 × 10 8 cells/cm 3 [17].This tumor model is only applicable for tumors having a volume greater than 1 mm 3 [8,17].
In the following we will consider two tumor growth toy models, described by the travelling wave solutions of the generalized Fisher-Kolmogorov equations Eq. ( 54), and by Eq. ( 88), respectively.6.1.Tumor growth in the first generalized Fisher-Kolmogorov equation model.We will investigate first the properties of the travelling wave model for tumor growth in the first generalized Fisher-Kolmogorov equation, whose solutions are given by Eqs. ( 55) and (56), respectively.In the present Section for D 0 we adopt the numerical value D 0 = D = 1.43 cm 2 /year.

6.1.1.
The case k = 1/4.We will begin our analysis of the travelling wave solutions in the first generalized Fisher-Kolmogorov model by considering the case k = 1/4.The velocity of the travelling front is given in this model by Then the general travelling wave solution of the first generalized Fisher-Kolmogorov equation is given in parametric form as and respectively.By performing a series expansion of the integrand in Eq. ( 95) we obtain where K = Cc max , giving The cell concentration c(W ) can be obtained approximately as Eqs. ( 98) and (99) give an approximate parametric representation of the travelling wave solution of the generalized Fisher-Kolmogorov equation Eq. ( 54) for k = 1/4.In the first order of approximation, In the limit W → ∞ we have lim The variation of the cell number density c as a function of ξ is represented in Fig. 1.As one can see from the figure, the cell number density increases linearly with ξ, and reaches a constant value at a finite ξ.
where k = ρD 0 /V 2 f > 1/4.In the second order of approximation the integrand in Eq. (102) becomes giving In the same order of approximation for the cell density we obtain In the first order of approximation the cell density is given by which shows that for small ξ the cell number density is a linearly increasing function of ξ.
In the limit of large W the cell number density tends to a constant value, given by lim respectively.
The variation of the cell number density c with respect to ξ is represented in Fig. 2. The cell number linearly increases as a function of ξ, and reaches a maximum constant value, whose numerical value is dependent on the value of V f .6.1.3.Biological implications.One of the important predictions of the standard Fisher-Kolmogorov equation Eq. ( 9) is that the velocity v F K of the detectable tumor margin approximately satisfies the relation v F K ≈ 2 √ ρD [15].This result was obtained from the observation that a cell population with a dynamics determined by diffusion and growth alone expands at a constant velocity v F K for large times, thus expanding linearly for any given value of the (constant) diffusion coefficient, and ρ, respectively.An experimental study of the growth of low grade gliomas was performed in [56], and it indeed confirmed that low-grade gliomas do grow both slowly and linearly.As for the tumor growth velocity, in the first 27 patients studied in [56], the average velocity of the diameter is about 4 mm/year.On the other hand a much higher radial tumor velocity, of the order of 12 mm/year, was reported in a single rare patient with a glioblastoma that was followed with repeated MRIs for a year without intervening treatment [15].
In the glioblastoma growth models described by the first generalized Fisher-Kolmogorov equation, the assumption of the concentration depending diffusion coefficient, and of a generalized reaction function do allow a wide range of velocities for the tumor front.Interestingly, these velocities are independent of the constant α in Eq. ( 93), and they depend only on the constant value D 0 of the diffusion function The initial values of the cell concentration and its derivative are c(0) = 1000 cells/cm 3 and (dc/dξ)| ξ=0 = 7 × 10 8 cells/cm 4 .D(c) for c = 0, D 0 = D(0), on the proportionality coefficient ρ in the reaction function, and on the arbitrary constant k.For k = 1/4, and for the adopted values of D 0 and ρ, the velocity V f of the tumor front is of the order of V f ≈ 10 cm/year, rather different to the value of 12 mm/year adopted in [15].In this case we obtain for the velocity of the tumor growth front the same relation as for the approximate expression of the velocity in the standard Fisher-Kolmogorov equation, V f | k=1/4 = v F K .On the other hand, a very wide range of tumor front velocities V f = ρD 0 /k can be obtained for k = 1/4, without any need of modifying the numerical values of the basic model parameters D 0 and ρ.This shows that by adopting tumor growth models with concentration dependent diffusion and reaction functions, a wide variety of experimental/clinical observations could be modelled, and the differences in the tumor growth observations could be attributed to the variations in the diffusion and reaction properties of the tumors.
All the considered travelling wave solutions of the first generalized Fisher -Kolmogorov equation show a linear expansion of the tumor size, who reaches suddenly the plateau phase, where the cell density becomes a constant.The tumor size evolves linearly in both the small cell density and high cell density phases, and therefore it is a general property of the present models.This behavior is consistent with the experimental results presented in [56].As shown in Fig. 1, for k = 1/4, there is strong dependence of the cell concentration on the parameter α, determining the diffusion and the reaction laws.On the other hand, for a fixed α and for k = 1/4, the modifications in V f , due to its k-dependence, do affect significantly the maximum value of the cells in the plateau phase, as well the moment when this maximum is reached.6.2.Tumor growth in the second generalized Fisher-Kolmogorov model.
We consider now the model in which the diffusion and reaction functions satisfy the condition given by Eq. (87).Moreover, we assume D(c) = D 0 = constant, which fixes the reaction function as where A m = A(m + 3) 2 /(2m + 1).The corresponding reaction-diffusion equation describing tumor growth is given by We call Eq. ( 111) the second generalized Fisher-Kolmogorov equation.Its travelling wave solution is given, in a parametric form, by Eqs. ( 90) and (92).In Eq. ( 90) we take, without any loss of generality, C 7 = 1.The constant b is given by b The dependence of the cell number density on the parameter τ is obtained as while the dependence of ξ on the parameter τ is given by where the integration constant ξ 0 can be taken as zero without any loss of generality.The integral can be computed exactly to give By taking into account that for τ = τ 0 we have c (τ 0 ) = c 0 , from Eq. ( 113) we obtain the value of the integration constant C 8 as Therefore for c(τ ) we obtain .
By estimating at τ = τ 0 , and by using Eq. ( 22), it follows that the initial value τ 0 of the parameter τ is obtained as a solution of the algebraic equation In the first order of approximation we have With the use of the above expression for τ we obtain the number cell density in the first order approximation as where The variation of the ratio of the cell number density divided by c max is represented, in a logarithmic scale, and for different values of m and A m in Figs. 3 and  4, respectively.7. Discussions and final remarks.In the present paper we have investigated the travelling wave solutions of a general class of diffusion-reaction systems, with the diffusion and reaction functions given as functions of the concentration c of the diffusing component.In this case the second order differential equation, describing the travelling wave solution for the system can be reduced to a first order first kind Abel ordinary non-linear differential equation.The integrability conditions of the Abel equation allow to obtain several classes of exact travelling wave solutions of the general reaction-diffusion system, with the reaction and diffusion functions satisfying some compatibility conditions.We have considered the solutions that can be found by using the Chiellini integrability condition, and the Lemke transformation, respectively.From the large class of exactly integrable reaction-diffusion equations we did concentrate on two equations, both representing generalizations of the Fisher-Kolmogorov equation with more general diffusion and reaction functions.More exactly, with the choices The diffusion coefficient D 0 = 1.43 cm 2 /year, V f = 9 cm/year, while the initial value of the parameter τ is τ 0 = 10.The tumor concentration wave is travelling in the direction of negative x. and respectively, the travelling wave solutions of the corresponding one dimensional reaction-diffusion equations can be obtained in an exact parametric form.Both equations represent generalizations of the standard diffusion and Fisher-Kolmogorov equations, to which they reduce in some limiting case.For the first generalized Fisher-Kolmogorov equation, for α = 0 we reobtain the standard diffusion equation, while for m = 2 the equation of the generalized tumor model growth reduces to the Fisher-Kolmogorov equation.From a very general physical point of view the behavior of microscopic particles is called ideal if they form very dilute solutions.In this limit, their positions, orientations, and movements are considered to be non-correlated [57].The macroscopically measured quantities usually are averages over an ensemble of molecules.In real solutions, however, even weak inter-particle interactions will cause a concentration dependence of the observed properties [57].The study of such systems with weak inter-particle interactions is important in many processes.For examples, weak interactions determine the tendency for a suspension of particles to remain in solution, to aggregate, to overcome phase separation or to form crystals [57].Concentration-dependent diffusion coefficients are also use to model mass transfer phenomena in polymer membranes [58], as well as water sorption in a polymer matrix composite [59].Several functional forms for the concentration dependence of the diffusion function have been proposed.A very simple such form is a linear dependence of the "long time" diffusion function D(c), given by [57] where D 0 = constant is the "short time" diffusion coefficient, and k D is a constant that can be expressed quite generally in terms of the pair correlation function of the particles, and the equilibrium segment distribution about their centre of mass.A more general, both temperature T and concentration dependent diffusion coefficient was considered in [59], with where the activation energy E a (c) is given by with R the ideal gas constant, T b and T a two distinct temperatures, and D b , α, β and γ constants.
We have also investigated the possibility of modelling of the glioblastoma tumor growth by using the obtained exact solutions of the general reaction-diffusion equation, corresponding to the specific choices of the diffusion and reaction functions.In order to obtain the numerical results we have used realistic values of the biological and physical parameters determining the tumor growth.Similarly the to Gompertz curve, for the first model the growth slows done at the end of a time period, with the density of the cells reaching a constant value.A very different behavior characterizes the second model, which for a monotonically increasing cell number density requires V f < 0, that is, ξ = x + V f t.Moreover, for arbitrary m the cell number density increases very rapidly with ξ.Therefore, this model is not relevant for the study of the glioblastoma growth, and its biological dynamic.For V f > 0 we obtain a traveling wave solution of the second generalized Fisher-Kolmogorov equation traveling in the direction of negative x.
In the present paper we have considered the biological implications of the travelling wave solutions of the generalized tumor growth equations, which corresponds to some given forms of the diffusion and reaction functions D(c) and Q(c).In the standard model of glioblastoma growth [8], two major biological phenomena underlying the growth of gliomas at the cellular scale are taken into account: proliferation and diffusion.The simplest choice for the proliferation term is a constant growth rate ρ, leading to an exponentially growing total number of glioma cells.For the invasive properties of gliomas, cell migration is assumed to be a random walk, corresponding to a passive (Fickian) diffusion characterized by a single coefficient D.Besides the logistic form Q(c) = ρc (1 − c/c max ), other forms for the reaction function Q(c) have been considered in the literature, like, for example, Q(c) = ρc ln (c m /c), corresponding to the Gompertz law [60,61].Improved mathematical models by including anisotropic extension of gliomas have been also investigated, by adopting a cell diffusion tensor derived from the water diffusion tensor, as given by MRI diffusion tensor imaging [61].
In our analysis of the tumor growth model we have extended the allowed functional forms of D(c) and Q(c), and we have also considered the possibility of the existence of a functional dependence between the diffusion of the malignant cells and logistic growth.In the first considered case, corresponding to the choices given by Eq. (125) for D(c) and Q(c), the product of these functions satisfy the relation D(c)Q(c) = D 0 ρc, that is, the product of the diffusion and reaction functions is proportional to the cell concentration.This implies a cell concentration dependence of the diffusion function of the form D(c) ∝ c/Q(c), which imposes a tight relation between D(c) and Q(c).From a biological point of view one can assume that at the beginning of the tumor evolution its growth is dominated by the logistic reaction function Q(c) = 0, so that lim c→0 D(c) = lim c→0 c/Q(c) = 0, implying that at the early growth stages of the tumor one can neglect diffusion.On the other hand with the increase of the malign cell concentration, and with the decrease of Q(c), the diffusion becomes the main physical process governing the spread of the tumor.
The clinically distinct feature of glioblastoma lies within its infiltrative potential, rendering complete tumor resection nearly impossible [2].The glioblastoma cells have a great migratory and invasive potential of their surroundings, which can be related mainly to their diffusive properties.In the varying diffusion and reaction cell growth model considered in the present papers this invasive potential can be enhanced, as compared to the standard Fisher-Kolmogorov model.As one can see from Fig. 1, for k = 1/4, the increase of the coefficient α in the diffusion and reaction functions leads to a significant increase in the size of the tumor, from around ξ = 0.3 cm for α = 1/2 to ξ = 0.4 cm for α = 3/2.Therefore the invasive properties of the glioblastoma cells are strongly dependent on the functional forms of D and Q.On the other hand, in this case the maximum cell number is not affected by the variations of α.A strong dependence of the invasiveness of glioblastoma cells on the numerical values of k = 1/4, which fixes the front velocity V f , can be observed in Fig. 2. In this case, for fixed diffusion and reaction functions, the change in k determines a large variation of both the tumor size and the maximum cell number.Therefore in this model the invasiveness of the glioblastoma cells with fixed D and Q is determined by the coefficient k.For the case of the growth models described by the second generalized Fisher-Kolmogorov equation, presented in Fig. 4, since the diffusion coefficient is a constant, the tumor evolution and invasiveness is mostly determined by the reaction function Q(c), whose dependence on the concentration is described by the parameter m, and the constant A m .There is a very significant difference in invasiveness for the different models, with combinations of the free parameters resulting in tumoral extensions as large as ξ = 1 cm.
In the present paper we have ignored the effects of the radiotherapy and chemotherapy on the glioblastoma growth and evolution, described by Eqs.(10) and (11), respectively.In these cases the presence of the explicitly time and space dependent terms G(t) and R ( r, t) do not allow, in general, the reduction of the corresponding reaction-diffusion equations to an Abel type equation, unless some very particular functional forms of G and R are assumed.
Most of the mathematical studies of the reaction-diffusion equations, including those modelling the growth of glioblastomas, have been done by using either qualitative methods, or asymptotic evaluations.The usual way to study nonlinear reaction-diffusion equations similar to Eq. ( 14) consists of the analysis of the behavior of the system on a phase plane (dc/dξ, c), useful for qualitative analysis, however insufficient for finding any exact solution or testing a numerical one.Finding exact analytical solutions of the tumor growth models can considerably simplify the process of comparison of the theoretical predictions with the experimental data, without the need of using complicated numerical procedures.

3. 1 .
Travelling wave solutions for D(c) = constant and Q(c) ∝ c.As a first example of the application of Theorem 1 we consider the case when the diffusion and the reaction functions are given by D(c) = D 0 = constant and Q(c) = ρc, ρ = constant.Therefore the equation describing the travelling wave propagation in the system takes the form

5. 3 .
Travelling wave solutions with D(c)Q(c) ∝ 1/c.If the diffusion and the reaction functions D(c) and Q(c) satisfy the relation

6. 1 . 2 .
The case k = 1/4.For k > 1/4, the general solution of the first generalized Fisher-Kolmogorov equation is given by

Figure 3 .
Figure 3. Variation of the ratio of the cell number density and the maximum cell number c max , c/c max , as a function of ξ (in a logarithmic scale) for the second generalized Fisher-Kolmogorov model for m = 1.2 and A m = 10 7.5 (solid curve), m = 1.3 and A m = 10 6.94 (dotted curve), m = 1.4 and A m = 10 6.4 (short dashed curve) and m = 1.5 and A m = 10 5.87 (dashed curve), respectively.The diffusion coefficient D 0 = 1.43 cm 2 /year, V f = −9 cm/year, while the initial value of the parameter τ is τ 0 = 10.The tumor concentration wave is travelling in the direction of positive x.

Figure 4 .
Figure 4. Variation of the ratio of the cell number density and the maximum cell number c max , c/c max , as a function of ξ (in a logarithmic scale) for the second generalized Fisher-Kolmogorov model for m = 1.2 and A m = 10 7.5 (solid curve), m = 1.3 and A m = 10 6.94 (dotted curve), m = 1.4 and A m = 10 6.4 (short dashed curve) and m = 1.5 and A m = 10 5.87 (dashed curve), respectively.The diffusion coefficient D 0 = 1.43 cm 2 /year, V f = 9 cm/year, while the initial value of the parameter τ is τ 0 = 10.The tumor concentration wave is travelling in the direction of negative x.