We present an analysis of a mathematical model describing the key features of the most frequent and aggressive type of primary brain tumor: glioblastoma. The model captures the salient physiopathological characteristics of this type of tumor: invasion of the normal brain tissue, cell proliferation and the formation of a necrotic core. Our study, based on phase space analysis, geometric perturbation theory, exact solutions and numerical simulations, proves the existence of bright solitary waves in the tumor coupled with kink and anti-kink fronts for the normal tissue and the necrotic core. Finally, we study the linear stability of the solutions to calculate the time of tumor recurrence.

#### Keywords

- Solitary waves
- bright solitons
- dark solitons
- mathematical oncology

#### MSC 2010

- 92B05
- 35C07
- 35C08
- 34D10
- 35B35

The application of tools from nonlinear dynamical systems to describe in a quantitative manner the underlying mechanisms encountered in a number of biologically- and/or medically-oriented scenarios has attracted an increasing interest in recent years [2, 11, 12, 13, 19, 32, 44, 45, 56]. As multidisciplinary collaboration between mathematical modelers with experimental biologists and/or clinicians is becoming a more common practice, it is possible to tackle complex problems and develop new armamentaria that could potentially display a great impact in public health, with cancer being a quite prominent example [38].

Here, we analyze in detail a mathematical model that was proposed in [39] to describe the salient physiopathological hallmarks of a very aggressive type of primary brain tumor: glioblastoma (GBM), the most frequent high-grade glioma in adults. At present, GBM carries a dismal prognosis; the median overall survival from diagnosis is 14.6 months with the current standard therapy [46], which includes neurosurgery followed by radiotherapy in combination with the alkylating agent temozolomide, and various salvage therapies (e.g., antiangiogenic) once tumor progression recurs [51]. Large-scale research endeavors are ongoing to provide a comprehensive understanding of all the genetic and epigenetic alterations that underlie GBM [48]. Despite all of these efforts, GBM remains incurable (the five-year overall survival rate is less than 10%), although it should be mentioned that, albeit still very small, a slow but steadily increase in the number of long-term surviving patients (more than five years) has been reported [47].

Various mathematical models have been proposed to describe specific aspects of GBM in vivo [3, 4, 5, 6, 7, 18, 23, 28, 40, 49, 50, 54]. Most of these models assume that the tumor evolves according to a simple reaction-diffusion equation (or a combination of): the Fisher-Kolmogorov (FK) equation [33], which epitomizes the simplest form of a nonlinear reaction-diffusion process

where _{∗} are positive parameters accounting, in the case of GBM, for the diffusion constant (cell migration), the proliferation rate and the maximum cell density (carrying capacity). This equation was first proposed by Fisher [17] in a completely different scenario; as a deterministic version of a stochastic model for the spatial spread of a favoured gene in a population. It is also the natural extension of the logistic growth population model when the population disperses via linear diffusion [10, 34]. The same equation appears in flame propagation [55], neurophysiology [53], and Brownian motion processes [8] among other situations. This equation and its variants have attracted the attention of many mathematicians, which have calculated different exact solutions under various conditions (see for instance [41, 42] and references therein).

In one-dimensional scenarios the FK equation has solitary wave solutions of kink-type [1, 29], which account for the progression of the tumor front edge. In dimensions higher than one, the FK equation does not display any kind of travelling wave solution [9], thus its analysis must resort to numerical methods. However, some essential features of brain tumor lesions in their transit to high-grade gliomas (i.e. from low-grade astrocytoma to GBM) are neglected by the FK model, such as the formation of a necrotic core responsible for the increase of the interstitial pressure and the subsequent intracranial deformation that may result in patient death. Furthermore, they do not incorporate the interaction of the tumor with the adjacent normal tissue (brain parenchyma).

Consequently, in order to understand glioma progression from low to high grade and thinking of connections with clinical practice there is a need for: (i) models accounting for the key features of the tumor microenviroment dynamics without involving too many details on any of the intervening specific processes, and (ii) models that are simple enough to allow for some quantitative understanding, e.g. using the tools of nonlinear dynamical systems and/or nonlinear waves well developed in other scenarios [45].

In this paper we present in detail one of those models allowing both for a description of features seen in real GBM and for a theoretical analysis. Interestingly, this model leads to the existence of multicomponent solitary waves encompassing bright solitary waves in the glioma compartment coupled with kink and antikink waves for the normal tissue and necrotic core respectively. We will study the model’s phenomenology using a combination of numerical simulation methods, dynamical systems theory and qualitative analysis. In the course of the study we wish to bring out a number of distinguishing features with respect to the FK model that are also relevant when comparing their predictions of how GBMs evolve.

Our plan is as follows, first in Sec. 2 we present the model in detail. Next in Sec. 3 we analyze the existence of traveling wave solutions of the system of equations using qualitative analysis and geometric perturbation theory. Then in Sec. 4 we calculate bright soliton solutions and determine numerically the minimum speed of these solitary waves solutions. The linear stability analysis of small initial data is discussed in Sec. 5 together with its application to the computation of the relapse time after surgery. Finally Sec. 6 summarizes our conclusions.

The model to be discussed in detail in this paper for describing the progression of GBM was proposed in [39]. To account for the spatial structure occupied by the malignant glioma cells and their environment, let Ω be a domain of ℝ^{n}^{+}, which are nonnegative, twice differentiable in the interior of their domain and satisfy the initial-value problem

∀^{+}

for some nonnegative continuous functions _{0}, _{0}, _{0} :

The biological meaning of all the above quantities is the following: _{∗}(

Eqs. (2) add a number of quantities with respect to the Fisher-Kolmogorov Eq. (1). The standard FK equation is recovered by setting

In Eqs. (2), tumor cell spreading is taken into consideration by using a standard Fickian diffusion term in (2a). This is the simplest transport mechanism and the one employed in most of the continuous mathematical models on cell motility. More realistic and complicated diffusion terms in gliomas should probably be governed by fractional (anomalous) diffusion [14, 16, 35] or other more elaborate approaches [15, 27, 52] to account for the high infiltration observed in this type of tumors [21, 36] and the fact that cells do not behave like purely random walkers and may actually remain immobile for a significant amount of time before compelled to migrate to a more favorable place in terms of resources. Other possibilities are to build in (at least) two cell phenotypes observed in malignant gliomas [26, 43]: one of invasive type migrating through the brain parenchyma and the other, more proliferative, dependent on angiogenesis and migrating due to the effect caused by the destruction of the anomalous tumor vasculature. A number of mathematical models have been proposed to incorporate such a tumor heterogeneity [24, 30, 31], however, here we will maintain a single tumor cell phenotype that encompasses many of the features observed in GBM progression.

The nonlinear term in (2a) corresponds to proliferation and is mediated by a competition for space that is occupied by the two cell types and the necrotic core, which comprises apoptotic cells and cell debris. There is a _{∗} that may depend on space without loss of generality but, here, we will assume _{∗} to be a constant. We also explicitly add a tumor cell death term to include the fact that tumor cells, although generally lacking the apoptotic (programmed cell death) mechanisms [22], may succumb by means of a number of mechanisms that include the interaction with the immune system (e.g. microglia) in the normal tissue, hypoxia and acidosis generated by the anomalous metabolism of the glioma cells in the high density tumor areas, and deficiency of nutrients and physical support in the necrotic core. This tumor cell death term assumes that, on average, the characteristic tumor cell life time is 1/

As to the normal cell dynamics, via Eq. (2b), due to their differentiation state, the proliferation of normal brain tissue is almost negligible within the characteristic time scales of GBM progression, thus we will assume that the normal tissue will not be able to regenerate. This is why in Eq. (2b) we have represented the cell loss due to the interaction with the tumor by means of an arbitrary form 𝓕(

To simplify the analysis and to focus on the main dynamical features, we will restrict ourselves to a one spatial dimension and assume that _{∗}(

It is convenient to introduce the new functions

and the rescaled variables

Notice that

with 0 ≤

Travelling waves with constant speed are solutions of Eq. (5) depending on

Using Eq. (6b) to get

Let us first consider the case

In the limit

Defining 𝓥_{ξ} = 𝓦, Eq. (9) can be written as the autonomous system

The equilibria of Eqs. (10) are of the form (𝓥, 𝓦) = (𝓥_{∗}, 0), with 𝓥_{∗} being an arbitrary real number. We are interested in values of 𝓥_{∗} ≥ 0. Thus, linearizing Eq. (10) around (𝓥_{∗}, 0), we obtain that the eigenvalues of the Jacobian matrix are given by

and its corresponding eigenvectors are (1, 0) and (1, 𝓥_{∗} + _{2} ≠ 0). These points are nonhyperbolic points. If 𝓥_{∗} > 1 – _{∗},0) possesses a local unstable manifold and a local center manifold. Otherwise, (𝓥_{∗}, 0) has a local stable manifold and a local center manifold. Thus, to get a heteroclinic orbit joining two points, say (𝓥_{–}, 0) and (𝓥_{+}, 0), with 𝓥_{–} > 𝓥_{+}, it is a necessary condition that 𝓥_{+} < 1 – _{–} > 1 – _{–} < 1). These two conditions can be summarized in the physical constraint imposing that the normal-necrotic compartment density is always strictly larger than the normal compartment density; the difference originating in the fraction of tumor cell loss.

Bright solitons, when existing, correspond to positive

Its solution is

for any arbitrary _{±} with 𝓥_{∗} = 𝓥_{–} > 𝓥_{+} (cf. Eq. (6b)) and 𝓦(± ∞) = 0, thus

to have two positive roots 𝓥_{–} > 𝓥_{+} > 0. Thus, for the initial condition 𝓦(𝓥_{+}) = 0, Eq. (13) becomes

It is straightforward to prove that, given the initial condition 𝓦(𝓥_{+}) = 0, with 0 < 𝓥_{+} < 1 – _{–}, with 𝓥_{+} < 1 – _{–} < 1, such that there exists a heteroclinic orbit joining both values, 𝓥_{–} and 𝓥_{+}; the function 𝓦(𝓥) has a single minimum at

which is larger than 𝓥_{+} for 𝓥_{+} < 1 – _{+}) = 0 and the function 𝓦(𝓥) is decreasing in the interval (𝓥_{+}, 𝓥_{m}), it is clear that

Therefore, as 𝓦(1) > 0, the existence of a value 𝓥_{–} in the interval (𝓥_{m}, 1), such that 𝓦(𝓥_{–}) = 0, follows from a direct application of Bolzano’s Intermediate Value Theorem.

Moreover, it is straightforward to verify that no orbits can cross the triangle

from the outside. Therefore, _{–} such that 1 – _{–} < 1, since these points have an unstable manifold inside _{+}, which satisfy 𝓥_{+} < 1 – _{+} ∈ _{+} < 1 – _{+} is a _{–} ∈ _{–} < 1, such that there is always an orbit connecting 𝓥_{–} to 𝓥_{+} inside

Thus, we have shown the following result: there exists a heteroclinic orbit that connects 𝓥_{–} and 𝓥_{+} and, therefore, 𝓥 > 0, (i.e. to remain entirely in the fourth quadrant).

In general, these fast moving solitons are small ones as Fig. 1(c) shows. Let us notice that

that allows for a solution for 𝓥 in the form of quadratures, once initial data is specified,

Finally, in the limit _{+} + ^{–(1–β–𝓥∗–𝓥+)(ξ–ξ0)}. In a similar way, for _{–} + ^{–(1–β–𝓥∗–𝓥–)(ξ–ξ0)}.

Our previous analysis proves the existence of the heteroclinic orbit in the limit

which has the critical points (𝓥_{∗}, 0, 0), for 𝓥_{∗} ∈ ℝ.

On the other hand, writing ^{2}, system (21) becomes

which is the dual fast system associated to (21) (see, for instance, [20]). For

If in the system (21)

which is a two-dimensional submanifold of ℝ^{3}. We choose any compact subset of such manifold and designate this compact subset as _{0}. Geometric perturbation theory uses both the above systems: (22) provides us with an invariant manifold _{ε}_{0} and we study the flow of (21) restricted to this manifold.

Let us introduce the following definition, which is useful for our purposess:

The manifold _{0} is said to be normally hyperbolic if the linearisation of (22) at each point in _{0}, restricted to _{0}, has exactly dim _{0} eigenvalues on the imaginary axis.

Now, we use the Fenichel’s invariant manifold theory [20, 25] to prove that for _{ε}^{3} which is within distance _{0} and which is invariant for the flow (22):

_{ε}_{0}_{0}. ^{r}

Thus, from the previous definition, such a perturbed invariant manifold _{ε}_{0} is normally hyperbolic.

The linearisation of the fast system (22), restricted to _{0} has the matrix

and the eigenvalues of this are _{0} is normally hyperbolic and, by Fenichel’s theory, for sufficiently small _{ε}

We would like to determine the dynamics on _{ε}_{ε}

where ^{ε}^{r}^{0}(𝓥, 𝓦) = 0, and the equations on _{ε}

Now, substituting (25) into (21) yields that ^{ε}

We Taylor expand ^{ε}

where

Equating powers of ^{i}

Thus

and the equations on _{ε}

These equations determine the dynamics on the manifold _{ε}

When _{–},0) and (𝓥_{+},0), 𝓥_{–} > 𝓥_{+} ≥ 0, of (10) was proved in the previous subsection. Note that when _{0},𝓦_{0}) be the solution of (33) for _{–},0) and (𝓥_{+},0).

To prove the existence of such a connection, we write

where 𝓥_{1}, 𝓦_{1} ∈ ^{∞}(ℝ). We assert that (34) is a heteroclinic connection of system (33). To see it, we substitute (34) in (33). To order ^{2} the coupled system of differential equations governing 𝓥_{1} and 𝓦_{1} is

and we wish to prove that this system has a solution satisfying

In matrix notation, we can write Eqs. (35) as

Let us denote

By using the Fredholm alternative, we state that problem (37) has a solution if, and only if,

where 〈⋅, ⋅〉 is the Euclidean inner product on ℝ^{2} and for all solutions _{1}, _{2}) of the adjoint problem

with boundary conditions
^{2}(ℝ). Now, our goal is to know the number of solutions of the adjoint problem. Although the matrix in (40) is non-constant, (𝓥_{0}, 𝓦_{0}) is the front of (10) about which we know a great deal. Letting _{0} → 𝓥_{+} and 𝓦_{0} → 0 and the matrix is then a constant matrix with eigenvalues

Thus, as ^{2} for the adjoint problem is therefore the zero solution _{ε}

_{–}_{+}.

Since the relationship between 𝓥(

Instead of considering kink solutions, via Eq. (7), we can work in the variable

From Eqs. (5) and defining

Combining Eqs. (41) into a single (and more general) equation for the tumor density we get

where subscripts denote partial differentiation. Notice that Eq. (42) admits (trivial) exact solutions of the form

where _{0}(

where _{1} and _{2} are arbitrary constants. (44) remains localized if

We will look for bright solitary solutions of Eq. (42) in the form _{η}

Just as we did with the kink-type solutions of Eq. (7) representing the normal-necrotic compartment, within the fast limit

with constants _{0} and _{0}, representing the amplitude and a shift, respectively, which, without loss of generality, can be set _{0} = 0. Figure 2 compares the exact numerical profile from Eq. (45) to the explicit form Eq. (46) and their corresponding homoclinic orbits.

We now wish to determine the minimum speed _{0} above which the solutions of Eq. (45) are positive for all _{ξ}

By probing the parameter space for 0 <

This autonomous differential equation possesses an exact explicit quadrature. Assuming that _{0} and

Positive solutions of Eq. (45) require that _{0}, _{+} < 0 (_{–} > 0). These asymptotic values can be found from Eq. (49) by means of the ancillary function

The nonlinear equation _{0}) = 0 has two real roots _{+} < 0 and _{–} > 0 if _{0}, whereas for 0 < _{0} the only root of _{–} > 0. At _{0}, _{+} satisfies

that is, _{+} becomes a double root of

By imposing the first two conditions of Eq. (51) to Eq. (50) we arrive at the two simultaneous transcendetal equations

Eqs. (52) give rise to two solutions for _{+}, _{+} = 0 and _{+} = −
_{+} = 0 leads to _{0} = 0). Upon substitution of _{+} = −
_{0} in terms of _{0}

The dependence of _{0} on _{0} is depicted in Fig. 3.

It is worth comparing the values of the minimum velocity that arise in the present model with those of the FK equation. It is well known [33] that the only solution of the FK equation evolving from positive compactly supported initial data that remains bounded and propagates as an antikink solitary wave is the one with the constant minimal speed
_{0} > 0, can actually be lower than the limit set by the FK equation. The consequence of this is that for the same values of

In this section we explore the possibility of calculating the time of tumor recurrence after the macroscopic tumor volume (detectable in MRI) has been surgically removed. To study this scenario from a mathematical point of view we examine the linear stability of the solutions

for the system (5), where _{*} is a positive constant. In fact, we are interested in the stability of the tumoral density

Now consider a small perturbation on (

with 0 <

Otherwise,

The fact that

Substituting Eqs. (55a) into Eq. (5) and keeping only the first-order terms in

We now look for solutions to the linear Eqs. (58a) by setting

which, upon substituting into (58a) and (58b), gives

From Eq. (60b) we obtain that

Therefore, one has the following condition for the instability of the solution

Thus, if

with _{s}_{1} and _{2} are the respective solutions to Eqs. (60a) and (60b). Notice that Eq. (60a), together with boundary conditions (63), is a regular Sturm-Liouville problem and has two infinite families of solutions for

and for

For the first case, the only solution of Eq. (60a) satisfying the boundary conditions (63) is the trivial solution. Thus, we focus on the second case: _{*} –1 <

Hence, the most unstable eigenvalue is achieved for _{1}, such that _{n1+1} > 0 and _{n1} < 0. For this value of _{1}, _{s}

Therefore, the eigenfunctions are given by

and, from (61), we obtain the function _{2}

Equivalently, another way to get the same result is by the following way: If we multiplying both sides of Eq. (58a) by exp(−(1−_{*})

This equation has the same form as the heat equation. The boundary conditions for such equation are

where

is the profile of the initial condition, where _{0} > 0.

The solution of Eq. (70) is given by

where _{n}_{n}

Now, we wish to calculate the time of tumor recurrence. We assume that the amplitude of the perturbation is appreciable for a value _{0}. It is clear that the most unstable modes occur when _{n}_{n}

where the symbol ⌊⋅⌋ denotes the integer part of a number. In the biological literature, there is a vast range of values for the diffusion and proliferation coefficients. To carry out the estimations, we resort to the following value for the proliferation ^{−1}, which is in the range [0.01–0.5] day^{−1}, taken from [28, 54] and ^{2}/day (which is in the range [0.0004–0.1] mm^{2}/day) [39].

Finally, we take ^{−1}, _{0} = 10 mm, _{*} = 0.3 in units of _{*}. Therefore, the recurrence time is

In this paper we have analyzed a simple model of glioma progression incorporating the normal tissue, tumor cells and the necrotic core. In comparison with the Fisher-Kolmogorov equation, widely used to model GBM, the model studied here only adds a single extra effective parameter accounting for the finite tumor cell lifetime in the tumor microenvironment due to the effect of vascular degeneration, competition for space and resources, hypoxia, acidosis and interaction with the immune system. Remarkably, with only this simple addition the model displays many of the signatures of GBM embodied as bright solitary wave solutions behaving as attractors of the dynamics of the tumor rim. We have studied the formation and propagation speed of this tumor front, we have proven numerically (and analytically in certain limits) that it takes the form of a solitary wave, indeed a vector soliton composed of a bright soliton in the malignant cell population coupled to a kink in the normal cell population and an antikink in a necrotic core compartment. We have also computed what is the relapse time of the tumor after surgical extirpation of the tumor even when it is possible to achieve total resection of the visible part of the tumor. It is worth mentioning that the outcomes of this model, such as the relation between the width of the bright solitons with the biological constants (diffusion, proliferation and death reate constants), have been used to predict novel imaging biomarkers for GBM patients [37]. In particular, the average size of the tumor rim, as observed in contrast enhanced T1-weighted magnetic resonance images, is a statistically significant metric for predicting the overall survival of GBM patients.

^{18}F-FMISO-PET, J. R. Soc. Interface, 12, 20141174. ^{18}F-FMISO-PET

Regarding new wave distributions of the non-linear integro-partial Ito differential and fifth-order integrable equations Nonlinear Mathematical Modelling of Bone Damage and Remodelling Behaviour in Human Femur Value Creation of Real Estate Company Spin-off Property Service Company Listing Entrepreneur's Passion and Entrepreneurial Opportunity Identification: A Moderated Mediation Effect Model Applications of the extended rational sine-cosine and sinh-cosh techniques to some nonlinear complex models arising in mathematical physics Study on the Classification of Forestry Infrastructure from the Perspective of Supply Based on the Classical Quartering Method A Modified Iterative Method for Solving Nonlinear Functional Equation New Principles of Non-Linear Integral Inequalities on Time Scales Has the belt and road initiative boosted the resident consumption in cities along the domestic route? – evidence from credit card consumption Analysis of the agglomeration of Chinese manufacturing industries and its effect on economic growth in different regions after entering the new normal Study on the social impact Assessment of Primary Land Development: Empirical Analysis of Public Opinion Survey on New Town Development in Pinggu District of Beijing Possible Relations between Brightest Central Galaxies and Their Host Galaxies Clusters and Groups Attitude control for the rigid spacecraft with the improved extended state observer An empirical investigation of physical literacy-based adolescent health promotion MHD 3-dimensional nanofluid flow induced by a power-law stretching sheet with thermal radiation, heat and mass fluxes The research of power allocation algorithm with lower computational complexity for non-orthogonal multiple access Research on the normalisation method of logging curves: taking XJ Oilfield as an example A Method of Directly Defining the inverse Mapping for a HIV infection of CD4+ T-cells model On the interaction of species capable of explosive growth Research on Evaluation of Intercultural Competence of Civil Aviation College Students Based on Language Operator Combustion stability control of gasoline compression ignition (GCI) under low-load conditions: A review Research on the Psychological Distribution Delay of Artificial Neural Network Based on the Analysis of Differential Equation by Inequality Expansion and Contraction Method The Comprehensive Diagnostic Method Combining Rough Sets and Evidence Theory Study on Establishment and Improvement Strategy of Aviation Equipment Design of software-defined network experimental teaching scheme based on virtualised Environment Research on Financial Risk Early Warning of Listed Companies Based on Stochastic Effect Mode System dynamics model of output of ball mill The Model of Sugar Metabolism and Exercise Energy Expenditure Based on Fractional Linear Regression Equation Constructing Artistic Surface Modeling Design Based on Nonlinear Over-limit Interpolation Equation Optimal allocation of microgrid using a differential multi-agent multi-objective evolution algorithm About one method of calculation in the arbitrary curvilinear basis of the Laplace operator and curl from the vector function Numerical Simulation Analysis Mathematics of Fluid Mechanics for Semiconductor Circuit Breaker Cartesian space robot manipulator clamping movement in ROS simulation and experiment Effects of internal/external EGR and combustion phase on gasoline compression ignition at low-load condition Research of urban waterfront space planning and design based on children-friendly idea Characteristics of Mathematical Statistics Model of Student Emotion in College Physical Education Human Body Movement Coupling Model in Physical Education Class in the Educational Mathematical Equation of Reasonable Exercise Course Sensitivity Analysis of the Waterproof Performance of Elastic Rubber Gasket in Shield Tunnel Impact of Web Page House Listing Cues on Internet Rental Research on management and control strategy of E-bikes based on attribute reduction method A study of aerial courtyard of super high-rise building based on optimisation of space structure Exact solutions of (2 + 1)-Ablowitz-Kaup-Newell-Segur equation