Effective Particle Methods for Fisher-Kolmogorov Equations: Theory and Applications to Brain Tumor Dynamics

Effective Particle Methods for Fisher-Kolmogorov Equations: Theory and Applications to Brain Tumor Dynamics

Juan Belmonte-Beitia Departamento de Matemáticas, E. T. S. I. Industriales and Instituto de Matemática Aplicada a la Ciencia y la Ingeniería (IMACI), Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain (juan.belmonte@uclm.es). Gabriel F. Calvo Departamento de Matemáticas, E. T. S. I. Caminos, Canales y Puertos and Instituto de Matemática Aplicada a la Ciencia y la Ingeniería (IMACI), Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain (gabriel.fernandez@uclm.es). Víctor M. Pérez-García Departamento de Matemáticas, E. T. S. I. Industriales and Instituto de Matemática Aplicada a la Ciencia y la Ingeniería (IMACI), Universidad de Castilla-La Mancha, 13071 Ciudad Real, Spain (victor.perezgarcia@uclm.es).
July 15, 2019

Extended systems governed by partial differential equations can, under suitable conditions, be approximated by means of sets of ordinary differential equations for global quantities capturing the essential features of the systems dynamics. Here we obtain a small number of effective equations describing the dynamics of single-front and localized solutions of Fisher-Kolmogorov type equations. These solutions are parametrized by means of a minimal set of time-dependent quantities for which ordinary differential equations ruling their dynamics are found. A comparison of the finite dimensional equations and the dynamics of the full partial differential equation is made showing a very good quantitative agreement with the dynamics of the partial differential equation. We also discuss some implications of our findings for the understanding of the growth progression of certain types of primary brain tumors and discuss possible extensions of our results to related equations arising in different modelling scenarios.

Fisher-Kolmogorov equations, brain tumors, effective particle methods

1 Introduction

Many partial differential equations of relevance in applied sciences have robust localized solutions displaying particle-like behavior. A wealth of distinct behaviors are encountered and in general they are referred to as coherent structures and/or solitary waves[1, 2, 3, 4].

In a limited number of prominent cases the equations are known to be integrable. When that happens, there is an infinite number of conserved quantities and the solution of the initial value problem can be constructed using different mathematical methods such as the inverse scattering transform. In integrable systems initial data can be rigorously decomposed into solitons plus radiation (linear modes) and a complete analysis of the asymptotic dynamics can be made.

While there are several physically relevant systems ruled by integrable partial differential equations, there is a vast majority of problems that are nonintegrable. Remarkably, some of them consist of small perturbations to integrable problems. In those cases one can still construct a rigorous theory for the dynamics of solitons that allows for an analytical description of the dynamics [4, 5]. However, in many other instances the perturbations are not “small” and/or the basic underlying problem is not integrable but coherent structures still persist and constitute a basic elemente in the dynamics.

In many of those problems a variational formulation can be written and then a very popular method is the, so-called, effective Lagrangian method, collective coordinate method or effective particle method. This method assumes the profile of the solution to be given by a specific ansatz depending on a small number of time dependent parameters. The specific choice of the ansatz depends on the equation under study and in many cases is suggested by physical considerations. The names “effective particle” and “collective coordinates” come from the fact that, in the framework of this approach, one simplifies the dynamics of an extended field with spatio-temporal dependencies, i.e. having “infinite” degrees of freedom” to a finite (small) number of time-dependent quantities (coordinates). The name comes from analogy with classical mechanics that provides a simple description (coordinates) of a typically extended object.

While the ansatz does not provide an exact solution of the PDE it is tipically used through the variational formulation as a test function and equations are obtained for the evolution of the solitary wave parameters. This approach works remarkably well for many relevant problems having Hamiltonian structure and provides a way to describe the infinite-dimensional dynamics in a simple form when the dynamics is dominated by coherent structures. This method has been exploited extensively in applied sciences in a large number of works for a broad variety of equations having a variational formulation (see e.g. Refs. [1, 6, 7, 8, 9, 10, 11, 12] and references therein). When coherent structures are robust it furnishes a simple description of the dynamics. The main weaknesses of the effective particle method are that: (i) it requires some experience to select appropriate ansatzes that capture adequately the dynamics, and (ii) the reduction to finite dimensions is provided without a measure of the error of the approximation that is estimated a posteriori on the basis of numerical simulations of the parent PDE. While the method has been used in hundreds of papers dealing with the dynamics of nonlinear waves in non-integrable systems, to our knowledge there are no papers obtaining a priori error bounds for such types of approximation.

Of particular interest are those equations that cannot be derived from a variational principle as it happens e.g. in dissipative systems. In that context there has been a great interest on different types of finite-dimensional descriptions of the dynamics of systems ruled by PDEs in a variety of contexts (see e.g. [13, 14, 15, 16] and references therein). However a simple procedure such as the one provided by effective particle methods that allows applied scientists to reduce the dynamics of a partial differential equation with solitary waves to a set of finite dimensional simple equations for the solitary wave parameters is not available yet. In this paper we present a very simple methodology that allows to obtain those types of approximations for the Fisher-Kolmogorov (FK) and related reaction-diffusion equations. The simplest version of the FK equation is


and describes the evolution of a population density measured in units of a maximal population on a given spatial domain. This equation is the simplest reaction diffusion model incorporating two effects: dispersion with a dispersal rate and proliferation or population growth with rate . In Eq. (1) the population growth is of the so-called logistic type although other terms with similar qualitative form have been used in the literature. Eq. (1) is written here in dimensional form in order to connect better with applications, although one can rescale the spatial and temporal variables to get rid of the coefficients and .

The FK equation and its extensions are a family of ubiquous reaction-diffusion models arising in population dynamics problems [17, 18, 19], most prominently in cancer modelling [20, 21, 22], in the description of propagating crystallization/polymerization fronts [23], chemical kinetics [24], geochemistry [25] and many other fields (see e.g. [26, 27] and references therein). These equations do not admit a Lagrangian density depending on the field [2] and thus the variational formulation for the effective particle parameters cannot be written in the usual way.

In this paper we get a small set of ordinary differential equations mimicking, not only the asymptotic dynamics of fronts arising in the Fisher-Kolmogorov (FK) equation but also describing their transient evolution towards the asymptotic regime. The plan of the paper is as follows: First, in section 2, we present the theoretical approach and discuss its application to the find the evolution of kink-like initial data classical FK equation. A comparison of the results of the effective particle method with the numerical solution of the FK equation is made. Next, in section 3, we apply the method to get effective equations for the dynamics of initially localized solutions to the FK equation. We identify different dynamical regimes for three relevant quantities associated to the spatio-temporal evolution of such localized profiles of the FK equation which are not apparent from the usual numerical solution. Secs. 4 and 5 are devoted to several applications of the method relevant for the understanding of the growth dynamics of certain types of brain tumors. Next, in Sec. 6 we present an example of applications to models beyond the FK equation by adding a spatial dependence to the diffusion coefficient. Finally, in section 7, we discuss the implications of our results and summarize the conclusions.

2 Effective-particle method description of front-type solutions

2.1 Derivation of effective equations for the front parameters

Nonnegative single-front-type travelling wave solutions of Eq. (1) satisfying the boundary conditions


obey the ODE


and have been studied in detail [17]. It is well known that such fronts can be constructed whenever and a celebrated result by Kolmogorov and coworkers [28] states that compact support initial data decay asymptotically into this type of waves with . However, less is known on the transient dynamics until the asymptotic regime is reached (see e.g. [29] and references therein).

Following the basic idea behind effective particle methods, our aim is to find an approximation for the dynamics of Eq. (3) by means of a simple finite-dimensional expression of the form


The “effective particle” describing the front is parametrized by three quantities depending only on time, namely, the wave amplitude , the front position and the width .

A key point of the method is the choice of a suitable profile function in Eq. (4) approximating the spatial profile of the solution. In our case, we will take the profile in (4) to be inspired by the Ablowitz solution [30]


The solution given by Eq. (5) is the only simple explicit solution known for the Fisher-Kolmogorov equation (in its adimensional version, i.e. with ), but corresponds to the specific speed , slightly larger than the minimal speed solution to the FK equation. Thus, a natural choice for our front profile is


that has the expected asymptotic exponential decay for large values of . In what follows we will try to obtain equations for the dynamics of the “effective particle” defined by the parameters .

To proceed with the method, let us define the integral quantities:


which are related to the -norm (number of particles), center of mass and width of the gradient of the density , respectively.

Then, introducing (6) in integral (7a), it follows that . The evolution of can be obtained by a direct formal calculation


where we have used Eq. (1) and the fact that .

We may calculate in a similar way to get


On the other hand, differentiating with respect to time and using Eq. (1), we find


where we have taken into account that . Thus, differentiating (9) and combining it with (10), it follows that


To get an equation for the time evolution of the width we can first use Eq. (7c) to obtain


The time evolution of the width can be derived once more from the FK equation


Therefore, using expression (12) for together with Eq. (13), the time evolution of the width easily follows.

Summarizing, we get the following set of differential equations for the evolution of the front parameters as described by Eq. (4)


Thus, our method provides a simple set of differential equations governing the evolution of a few (but very noteworthy) quantities describing the propagation of the front. As it will be described in Sec. 2.2, explicit solutions for Eqs. (14) can be found what means that the dynamics of the reduced system can be easily found in closed form. Since Eq. (1) is time invariant and consequently Eqs. (14) autonomous we will choose without loss of generality in what follows without loss of generality. Since “a priori” error estimates are not available for our approach, we will later verify the quality of the approximation by resorting to numerical simulations to compare the results obtained from the PDE (1) with the reduced ODE model of Eqs. (14).

2.2 Analytical solutions

Although Eqs. (14) are a set of coupled nonlinear evolution equations their exact solutions can be found in closed form. First, Eq. (14a) is a logistic equation whose solution is given by


where . The expression for given by Eq. (15) can be inserted into Eq. (14c), which is linear in , to obtain the explicit solution for which reads (with )




is the Gaussian hypergeometric function with , , and the Euler gamma function. Notice that for .

Finally, we can also obtain the analytical expression for the velocity by combining (14b) and (16) to get


Thus Eqs. (15-18) provide the full dynamics of the front for all times. We can easily find the asymptotic behavior of these solutions to be


Also, it is worth mentioning that despite the crudeness of the approximation of the effective particle method, i.e. that the front maintains its basic shape during the evolution, the velocity in the asymptotic limit is very close to the exact one , with the percentual relative error being of the order of 2%, that is the difference between the real asymptotic speed and the one of the Ablowitz solution. However, the main strength of the method is that it provides qualitative information on the transient evolution of the front parameters.

2.3 Comparison with the FK equation

Due to the approximate nature of Eqs. (14) and the lack of a priori error bounds, it is necessary to validate the predictions of the effective particle description through a direct comparison with the numerical solution to Eq. (1). To do so, we have solved numerically the Eq. (1) with initial data given by (6) and a particular set of initial parameters . From the solution and using Eqs. (7), we obtain the soliton parameters numerically in terms of integral quantities, i.e.


These values are to be compared with the solutions of Eqs. (14), i.e. with Eqs. (15-18).

Figure 1 displays typical results of the comparative evolution of fronts according to our reduced model (14) and from the Eq. (1) through (20). The agreement between the reduced ODE set and the full PDE is excellent. To exclude the possibility of this result being the consequence of a fortunate choice of the initial data, we have explored different sets of initial conditions in the range , , and and model parameters in the range and . In all cases a very good quantitative agreement among the three sets of curves is observed for all times, the percent relative errors being always smaller than . The small (and expected) discrepancy of the asymptotic values for the speed and the width is always present in our calculations and is a result of our ansatz choice.

Figure 1: [Color Online]. Comparison of the evolution of front solutions of the Fisher-Kolmogorov equation described by Eq. (1) with the analytical solutions of the effective particle method given by Eqs. (15-18). The FK equation is solved numerically using a standard second order in space and time finite difference method with zero derivative boundary conditions and constants and . The initial data is given by Eq. (6) with , , and . In all subplots (a)-(c) the solid lines correspond to the parameters (20) extracted from the numerical solution of the FK equation. The dashed lines correspond to the analytical solutions of the ODEs. The subplots show the: (a) amplitude (dashed) versus (solid), (b) velocity of the front (dashed) versus (solid), (c) width of the solution (dashed) versus (solid).

3 Effective particle methods for localized solutions

3.1 Motivation

While front solutions of the FK equation have relevance in many practical scenarios, there are cases where the solutions are initially localized. A typical example are models related to the propagation of tumors, that start from the onset as localized low amplitude cell densities and extend through the healthy tissue as localized solutions.

To derive finite-dimensional simple models able to tackle these questions we may extend the effective particle method to obtain approximate localized solutions of the Fisher-Kolmogorov equation. The procedure, however, is less straightforward. The first aspect to be addressed is the choice of a proper ansatz.

3.2 The choice of the ansatz

In contrast with the profile given by Eq. (6), we now look for a nonnegative two-front wave satisfying the boundary conditions


Our simple finite-dimensional approximation will be of the form


with representing a single front. As before, the two counter-propagating fronts are parameterized by three quantities; the amplitude , the right-front position and the front widths . We will further assume that the resulting profile is spatially symmetric although this restriction can be lifted in systems without spatial symmetries. To this end, we resort to an extension of our previous ansatz


The above ansatz (23) satisfies the following properties: first, if , then ; next , for ; also , for all . Finally the amplitude at is given by , for all .

3.3 Evolution equations for the parameters of the effective particle

We now proceed to define the integral quantities:


These integral quantities are different from the ones () used to characterize the front solutions due to the fact that now we are dealing with localized solutions. The parameter represents the total “mass” and in population dynamics applications represents the normalized number of individuals. gives the variance of the density distribution having the biological meaning of spatial width or spatial extension occupied by the population. Finally, the parameter provides the right-front size, giving an estimate of the size of the infiltration region in applications.

Upon substitution of Eq. (23) in Eqs. (24a)-(24c), and after integration, we get

Figure 2: [Color Online]. Comparison of the evolution of front solutions of the Fisher-Kolmogorov equation described by Eq. (1) (solid curves) with that provided by the ansatz (23) (dasher curves) with parameters given by Eqs. (42)-(44) for various times: from the innermost to the outermost profiles. The diffusion coefficient is and the growth rate . The initial data is given by Eq. (23) with (a) , , and ; (b) , , and .

The evolution of , and is obtained via the FK equation (1) as in the case of single-fronts. For , we find


where we have used the fact that . Let us now consider , for which we get


where we have made use of Eqs. (25a), (25b) and (26).

Finally, for , we obtain


Now, in order to get a system of ordinary differential equations for the dynamically relevant quantities , and we can differentiate Eqs. (25a)-(25c) and use Eqs. (26)-(28). Long calculations lead to a final closed set of ordinary differential equations that are written in Appendix A.

Figure 3: [Color Online]. Comparison of the evolution of front solutions of the Fisher-Kolmogorov equation described by Eq. (1) with the solutions of the effective particle method given by Eqs. (42)-(44). The FK equation is solved numerically using a standard second order in time finite difference method with zero derivative boundary conditions and constants and . The initial data is given by Eq. (23) with , , and . The corresponding ODEs are solved taking the later as initial values for . In all subplots (a)-(c) the solid lines correspond to the parameters , and , extracted from the numerical solution of the FK equation and Eqs. (25a)-(25c). The dashed lines correspond to the results of the ODEs. The subplots shown are: (a) amplitude (dashed) versus (solid), (b) velocity of the front (dashed) versus (solid), (c) width of the solution (dashed) versus (solid).
Figure 4: [Color Online]. Comparison of the evolution of front solutions of the Fisher-Kolmogorov equation described by Eq. (1) with the solutions of the effective particle method given by Eqs. (42)-(44). The FK equation is solved numerically using a standard second order in time finite difference method with zero derivative boundary conditions and constants and . The initial data is given by a Gaussian profile . In all subplots (a)-(c) the solid lines correspond to the parameters , and , extracted from the numerical solution of the FK equation and Eqs. (25a)-(25c). The dashed lines correspond to the results of the ODEs. The figures shown are: (a) amplitude (dashed) versus (solid), (b) velocity of the front (dashed) versus (solid), (c) width of the solution (dashed) versus (solid).

3.4 Comparison of the effective particle dynamics with the FK equation

Let , and as in Sec. denote the initial values of the right-front position, amplitude and width, respectively. To test the validity of Eqs. (42-44) as approximations to the full PDE dynamics for localized initial data, we have run extensive series of simulations for different parameter values. As in the case of the front solutions discussed in Sec. 2 we have found that the effective particle model approximates with a very good accuracy the dynamics of Eqs. (1).

As an example in Fig. 2 we compare the profiles of the numerical solution to the FK equation (1) with the dynamics provided by the ansatz (23) together with Eqs. (42-44) for various times and different initial conditions. Figure 2(a) corresponds to the case where , whereas in Fig. 2(b) . Notice that in Fig. 2(b) the right and left fronts advance at a faster pace when compared to Fig. 2(a) despite the fact that the initial profile is much smaller. In both cases we get an excellent agreement since the profile chosen closely resembles the one arising spontaneously from the partial differential equation (1).

To get a more direct comparison between the parameter evolution computed from the ODEs (42-44) and the PDE (1) we have compared also the evolution of the parameters in many simulations. Figures 3 and 4 provide two typical examples. In Fig. 3 the initial condition is given by the ansatz (23), whereas in Fig. 4 the initial condition is a Gaussian profile that does not have initially the expected exponential decay of for large values of . Despite this deviation of the initial data there is generally a very good agreement of the approximate effective particle equations with the simulations of the PDEs.

3.5 Asymptotic regime

Despite Eqs. (42-44) are ordinary differential equations allowing to get the evolution of the parameters in a more direct way than the PDE, the fact that they have a complicated structure originated in the interactions between the two fronts used to approximate the solution, makes its qualitative analysis complicated. However, there are several limits in which it is possible to elucidate the dynamics of in terms of simpler expressions.

Let us consider the case when the population density is much more extended than the infiltrative zone, i.e. , a situation that is always verified for long enough times. In that case, the corresponding system of ordinary differential equations reduces to


Notice that Eqs. (29a) and (29b) are exactly the same as the first two equations in (14). Combining Eqs. (29b) and (29c) we arrive at the third equation in (14). This is consistent with the intuitively expected fact that if becomes very broad then the left and right fronts no longer interact and behave as independent single-propagating particles whose time evolution is described by Eqs. (14). It is clear also that in the limit , the parameters , and tend to the asymptotic values given by Eqs. (19). Therefore, if the dynamics of the localized profile is the same as the simpler single front wave.

4 Applications to brain tumor dynamics (I): Transition to malignancy and time of birth of low grade gliomas

4.1 Motivation

Low grade glioma (LGG) is a term used to describe World Health Organization grade II primary brain tumors of astrocytic and/or oligodendroglial origin [31]. These tumors are highly infiltrative and generally incurable but have a median survival time of years because of low proliferation [32, 33, 34]. While most patients remain clinically asymptomatic besides seizures, the tumor transformation to aggressive high grade glioma is eventually seen in most patients.

Management of LGG has historically been controversial because these patients are typically young, with few, if any, neurological disorders. Historically, a wait and see approach was often favored in most cases of LGG, due to the lack of symptoms in these mostly young and otherwise healthy adults. The support for this practice came from several retrospective studies showing no difference in outcome (survival, quality of life) if therapy was deferred [35, 36]. Other investigations have suggested a prolonged survival through surgery [37]. In absence of a randomized controlled trial, recently published studies may provide the most convincing evidence in support of an early surgery strategy [38] and waiting for the use of other therapeutical options such as radiotherapy and chemotherapy. However, the decision on the individual treatment strategy is based on a number of factors including patient preference, age, performance status, and location of tumor [33, 34].

The FK equation arises as a basic model of the dynamics of low grade brain tumors, accounting in a simple way for the two main features of tumor cells: infiltration and proliferation (this type of tumors do not metastasize to other organs). In this context the density describes a wave of invasive tumor cells as it has been studied in many papers [39, 40, 41, 42]. One of the main characteristics of those tumors is their low cellular density. However, when the cellular density becomes too high, hipoxia arises triggering the hypoxic response [43, 44] with a cascade of metabolic alterations in the cell, including genetic instability, and has a potential effect on the acceleration of the progression which may have an impact in what is called the transition to malignancy.

4.2 Estimates for the time of transition to malignancy

Assuming that one of the events determining the transition to malignancy is the high local tumor density, we can get an estimate of the time taken by the tumor cell density to progress from an initial maximal density to a critical threshold density beyond which hypoxia and tissue damage trigger the cell changes leading to the malignant transformation.

The exact analytical solutions provided by the effective particle method allow us to provide an estimate for the time required for that process. Finding this time from the FK generally requires numerical computation. In contrast, it is straightforward to get it from Eq. (15).


It is interesting to note that according to Eq. (30), the time of transition to malignancy does not depend on the diffusion coefficient . Eq. (30) should be taken as an estimate or order of magnitude upper bound since tumor density profiles are not expected to be as smooth as our ansatz functions (e.g. Eq. (6)) and the transition to malignancy may depend on the highest initial cell density present throughout the tumor.

To use (30) in clinical scenarios it is necessary to estimate the parameters , and . Firstly, would correspond to the detection amplitude and has been discussed in several papers (see e.g. [20] and references therein), is less well known but can be estimated to be in the range between 0.5 and 0.8. However more work with real data is necessary to confirm this choice. Finally the estimate of for individual patients is very difficult in the case of low-grade gliomas because of the slow and irregular growth of these tumors. Recent work has suggested that this number can be obtained from the response of the tumor to radiotherapy [44] and it seems reasonable that those estimates may correlate also with the histology results for proliferation markers when available (e.g. MIB/Ki-67 labelling index).

4.3 Estimates for the time of birth of the tumor

Running the equations (42-44) backwards in time it is also possible to propose an estimate for the time of birth of the tumor, i.e. the time for which the tumor maximum density corresponds to a single cell (assuming that the tumor starts from one mutated cell). In our case, since astrocytic cells have a typical size of 10 m and thus the maximal linear density is about 100 cells/mm, we would get . The equation


gives the tumor time of birth. As with Eq. (30), Eq. (31) depends inversely on and thus, assessing this parameter from the available patient’s data may result in a value dependent on each patient. Getting estimates for is very relevant clinically in order to correlate the prediction with the patient’s clinical history and to extract information on possible causes for the origin of this type of tumors. Up to now, they are thought to be sporadic, but there is a strong interest among clinicians to investigate possible causes for these tumors, for which Eq. (31) may be helpful. In fact, this problem has been considered in the framework of simulations of the full FK equation in [42]. Our approach complements that work providing an equation that allows to circumvent the direct simulation of the PDE.

5 Applications to brain tumor dynamics (II): Simple description of the response to radiotherapy

5.1 Radiotherapy of low grade gliomas

The effectiveness of radiation therapy against low grade gliomas was proven in the clinical trial of Ref. [45]. However, the timing of radiotherapy after biopsy or debulking is debated. It is now well known that immediate radiotherapy after surgery increases the time of response to the therapy known as progression-free survival, but does not seem to improve the overall survival time. At the same time, radiotherapy may contribute to the observed neurological deficit of this patients as a result of the damage to the normal brain [46]. This is why radiotherapy is usually deferred in time until an increase in the tumor grade is diagnosed or, else, offered to selected patients with a combination of low risk factors such as age, sub-total resection, and diffuse astrocytoma pathology [47].

It is interesting that slight modifications of the radiotherapy protocol have been found to have little or no effect on overall survival [48]. Mathematical modelling has the potential to select patients that may benefit from early radiotherapy. Also, it may help in developing specific optimal fractionation schemes for selected patient subgroups. However, despite its enormous potential, mathematical modelling has had a very limited use with strong focus on some aspects of radiation therapy (RT) for high-grade gliomas [49, 50, 51, 52, 53, 54]. Up to now, no ideas coming from mathematical modelling have been found useful for clinical application in any of these contexts.

There is thus a need for models accounting for the fundamental features of low-grade glioma dynamics and their response to radiation therapy without involving excessive details on the -often unknown- specific processes. The increasing availability of systematic and quantitative measurements of LGG growth rates provides key information for the development and validation of such models [55, 56].

5.2 A simple mathematical model of tumor response to therapy

We will assume that the radiation doses (Gy) are given instantaneously at times for integer , what implies assuming that the total irradiation time is very short in comparison with the tumor response times. This is typically the case in external beam radiotherapy where treatment times are in the order of minutes while tumor response times are in the range of weeks or months. We will also assume that the radiation is spatially uniform through the tumor what is also a good approximation to clinical practice whenever possible. Following the standard practice in radiotherapy, we will take the damaged fraction of tumor cells as given by the classical linear-quadratic (LQ) model [57], i.e. the fraction of cells that are not lethally damaged by a dose , to be given by


where and are respectively the linear and quadratic coefficients for tumor cell damage of the LQ model. The precise values of the parameters remain unknown despite many studies, because what is clinically more relevant is the ratio which for glioma cells is around 10. The full treatment consists of a total dose split in a series of -typically equal- doses delivered at times . For instance, for high grade gliomas the standard protocol consists of doses of Gy for a total of 60 Gy, administered once per day (except for the weekends) during 6 weeks. For LGGs typical doses range from 45 to 54 Gy with Gy during 5 or 6 weeks of treatment.

Taking into account all this information, and depending on the cell death mechanism, different models can be written. In Ref. [44] a model has been proposed including delayed response to radiation. However the simplest possible model can be based on the assumption that damaged cells die instantaneously after the therapy (or in times very short compared with the characteristic natural history of the tumor. This leads to the simplest possible model where the evolution of the tumor cell density is given by Eq. (1) or its ODE effective particle reduction described by (14)) for each interval between doses . The initial data for each subinterval will be take then as given by


where is the density after the dose is given.

Figure 5: Tumor amplitude evolution under radiation therapy for several choices of the parameter values. In both cases the solid line correspond to the amplitude computed using the FK equation (1) using the equation and the dashed line (overlapping with the solid one) corresponds to the approximation of the effective particle method given by Eqs. (36) and (37). Shownare: (a) day, mm/day, and initial tumor cell density given by with measured in mm. Radiotherapy follows the standard scheme (6 weeks with 1.8 Gy doses from monday to friday) and starts at time days. The differences between both curves are minimal with . (b) day, mm/day, and initial tumor cell density given by with measured in mm. Radiotherapy follows a non-standard scheme with radiation sessions the first week of every two months for a total of 6 weeks with 1.8 Gy doses from monday to friday and starts at time day. The differences between both curves are minimal with .

This allows us to write a simple model for the response of the tumor amplitude to radiation given by Eq. (32) together with Eq. (15). For instance, starting from an initial amplitude at time , we get that before and after the first irradiation at time the tumor amplitudes will be given by

Thus, defining and taking by definition we arrive to the following recursive formula for the tumor amplitudes after each irradiation cycle


Thus, the amplitude of the tumor after a sequence of radiation sessions administered at times with survival fractions , in the framework of our simple model, is given for all times by the equations


where the constants , are obtained recursively using Eq. (36).

5.3 Validation and implications

Eqs. (36) and (37) provide a very simple model of the response of a low grade glioma to radiation therapy. However these simple formulae have been obtained in the framework of the simple effective particle method. To test if the approximation provided by this approach is acceptable when describing the solutions of Eq. (1) together with an instantaneous response to radiation given by Eq. (32) we have numerically simulated these equations, computed the amplitude using . The evolution of this “exact” amplitude has been compared with the approximation provided by Eqs. (36) and (37) obtained in the framework of the effective particle approximation to the dynamics.

We have compared the evolution in different scenarios. The first one corresponds to the standard radiation fractionation of a total of 54 Gy in 30 fractions of 1.8 Gy over a time range of 6 weeks (5 sessions per week from monday to friday). Radiation is started 1 week after the initial simulation time () and the evolution is followed for six years. The comparison between the approximation based on the effective particle method and the full Fisher-Kolmogorov equation is displayed in Fig. 5(a). A second example is provided in Fig. 5(b) where a faster growing tumor ( day versus day in the previous case) is irradiated following a non-standard radiotherapy scheme consisting of five consecutive sessions of 1.8 Gy every two months with the intention to control the tumor rather than focusing on minimizing the tumor mass. The agreement between both approaches is excellent.

We have explored several parameter regimes finding a very good accuracy in the approximation provided by the effective particle method meaning that one can follow the evolution of the tumor amplitude using the simple equations Eqs. (36) and (37). This approach allows us to reduce the problem of solving a partial differential equation to computing a finite number of iterations of a simple nonlinear mapping. This fact has very interesting implications since one may then pose optimization problems for finding optimal fractionation schemes in a much simpler context. One example is optimizing radiation for delaying the transition to malignancy, a problem that can be reformulated completely in terms of amplitudes assuming that there is a critical amplitude beyond which the transition to malignancy is triggered. Thus the problem becomes a discrete optimization one with a very simple non-differential model that is equivalent to finding the extrema of a function of several variables with several restrictions. While the complete analysis of this problem is beyond the scope of this paper and will be considered elsewhere, we want to highlight the potential of the effective particle method to obtain simple equations amenable to a complete analysis.

6 Extension to related Fisher-Kolmogorov equations

Figure 6: [Color Online]. Comparison of the evolution of front solutions of the Fisher-Kolmogorov equation with space-dependent diffusion described by Eq. (38) with the analytical solutions of the effective particle method given by Eqs. (14). The FK equation is solved numerically using a standard second order in time finite difference method with zero derivative boundary conditions with as (40), where , , and . The initial data is given by Eq. (6) with , , and . In all subplots (a)-(c) the solid lines correspond to the parameters (20) extracted from the numerical solution of the FK equation. The dashed lines correspond to the analytical solutions of the ODEs. The subplots show the: (a) amplitude (dashed) versus (solid), (b) velocity of the front (dashed) versus (solid), (c) width of the solution (dashed) versus (solid).

The effective particle method can be extended to other types of reaction-diffusion equations having fronts as asymptotic attractors of the dynamics. For instance, for the cases where the diffusion coefficient either density dependent , space and/or time dependent and/or the growth rates also have extra dependencies , as well as for other situations where the system is multicomponent (e.g. in cancer modelling, there are several phenotypes). These examples are of interest in a broad range of biological processes [18] and particularly in cancer modelling problems [61, 60, 58, 59, 62]. Here, as an example of how this methodology is robust and can be extended beyond the most basic FK equation, we consider the equation with space-dependent diffusion


Thus, we are considering Fickian diffusion but with a diffusion coefficient dependent on the space coordinates [17]. Then, we can find the set of effective particle equations for the amplitude , the center of mass and the width of the travelling wave (6), for Eq. (38) following the same procedure of Section 2. The result is


To test the limits of the method we will use a piecewise constant discontinuous diffusion coefficient given by


Effective particle methods have tipically difficultities in dealing with discontinuous or very fast varying coefficients, and other effects that influence strongly the shape of the solution thus providing asymmetric deformations of the initial ansatz. However as we will see in what follows, in this case the method is able to follow qualitatively the details of the evolution of the deformed front solution.

The explicit choice for the diffusion coefficient given by Eq. (40) has also some interest in applications. Specifically, in brain tumor modelling it arises in situations in which the tumor invades the gray matter from the white matter. In that case Eq. (38) has been proposed as a toy model in which the diffusion coefficient is a piecewise constant function, corresponding to different tumor cell motilities in the white and grey matter [41, 60, 17]. Analogous mathematical problems arise in other application scenarios (see e.g. [27] and references therein). The explicit form of given by Eq. (40) allows us to compute explicitly the integral in Eq. (39c) to get


Eq. (41) together with Eqs (39a) and (39b) is again a closed system of ODEs ruling the dynamics of the front in the framework of the effective particle method. We have run extensive simulations to compare the dynamics of Eqs. (39) with the numerical solution calculated for Eq. (38) and the ansatz (6). Despite the potential problems that might be expected coming from the discontinuity of the diffusion coefficient the agreement is very in all of the cases studied. A typical example is shown in Fig. 6 where it is seen how the amplitude dynamics is fully captured by the effective particle method and the width and speed dynamics have only quantitative transient differences with the asymptotic behavior been again correctly described by the approximation method.

7 Conclusions

In this paper we have presented an extension of effective particle methods to deal with a non-Hamiltonian problem of relevance in applied science: the Fisher-Kolmogorov (FK) equation. The method provides a very simple picture in terms of ordinary differential equations of the behavior of both a single-front and a localized travelling wave. It yields direct information on three relevant parameters: the amplitude, the front position and the width of the wave, which turn out to be parameters more easily accessible to experimental measurement in application scenarios than the entire profile , yet furnishing sufficient insight on the characteristic dynamics of the dynamics of the partial differential equation in certain regimes.

In addition to presenting the method and quantifying its accuracy, we have also discussed, through the specific application to problems arising in the description of brain tumors, how it can be used to get very simple estimates useful for applied scientists. Specifically we have developed explicit formulae to estimate the times of transition to malignancy and of birth of a low grade glioma. We have also provided a way to transform the problem of optimizing radiation delivery on the PDE to a finite-dimensional problem involving only a discrete map.

The method presented in this paper has very broad implications and potential uses. As an example we have shown its appropriateness to deal with a problem with spatially dependent diffusion with discontinuous diffusion coefficient. However, it can be extended to many other reaction diffusion equations in order to get a simple qualitative understanding of the dynamics of coherent structures. Secondly extending it to higher dimensions, may allow to get simple models in situations were theoretical results are much more scarce and numerical simulations more difficult. In that case approximate front profiles may be used as tentative test functions for the method [63]. Finally, the set of ODEs provided by the effective particle method also allow simplifying optimal control problems, such as those involved in finding the optimal combinations of different therapies, that are much more difficult to cope within the framework of partial differential equations.

We hope that this paper would further stimulate the application of the method to get useful information for the many applications of the Fisher-Kolmogorov and related equations.


This work has been supported by grants MTM2009-13832 and MTM2012-31073 (Ministerio de Economía y Competitividad, Spain). We would like to acknowledge Alicia Martínez (Universidad de Castilla-La Mancha, Spain) and Philip Maini (Oxford University, UK) for discussions.

Appendix A Full form of the effective particle equations for localized initial data

For completeness, we detail the full expressions of the system of ordinary differential equations for , and which follow by differentiating Eqs. (25a)-(25c) and equating them to Eqs. (26)-(28).

The first equation corresponds to the total number


The second equation corresponds to the variance


Finally, the third equation corresponding to the right-front size is



  • [1] T. Dauxois, M. Peyrard, Physics of Solitons, Cambridge University Press, 2006.
  • [2] A. Scott, Nonlinear Science: Emergence and Dynamics of Coherent Structures, Oxford University Press, 2003.
  • [3] P. G. Drazin, R. S. Johnson, Solitons: An introduction, Cambridge University Press, 1989.
  • [4] J. Yang, Nonlinear Waves in Integrable and Non-integrable Systems, SIAM, 2010.
  • [5] Y. S. Kivshar, B. A. Malomed, Dynamics of solitons in nearly integrable systems, Rev. Mod. Phys. 61, 763-915 (1989).
  • [6] B. A. Malomed, Variational methods in nonlinear fiber optics and related fields, Prog. Opt. 43, 7-193 (2002).
  • [7] N. R. Quintero, E. Zamorano-Sillero, Lagrangian formalism in perturbed nonlinear Klein-Gordon equations, Physica D 197, 63 (2004).
  • [8] V. M. Pérez-García, Self-similar solutions and collective coordinate methods for nonlinear Schrödinger equations, Physica D 191, 211-218 (2004).
  • [9] R. Carretero-González, D. J. Frantzeskakis and P. G. Kevrekidis, Nonlinear waves in BoseÐEinstein condensates: physical relevance and mathematical techniques, Nonlinearity 21, R139 (2008).
  • [10] V. M. Pérez-García, P. Torres, G. D. Montesinos, The method of moments for Nonlinear ScheÁrödinger equations: Theory and Applications, SIAM J. Appl. Math. 67, 990 (2007).
  • [11] S. Cuenda. A. Sánchez, Length scale competition in nonlinear Klein-Gordon models: A collective coordinate approach, Chaos 15, 023502 (2005).
  • [12] N. R. Quintero, F. G. Mertens, A. R. Bishop, Generalized traveling-wave method, variational approach, and modified conserved quantities for the perturbed nonlinear Schršdinger equation, Phys. Rev. E 82, 016606 (2010).
  • [13] C. Foias, B. Nicolaenko, R. Temman (Eds.), The connection between infinite dimensional and finite dimensional dynamical systems, American Mathematical Society (1987).
  • [14] T. Bohr, M. H. Jensen, G. Paladin, A. Vulpiani, Dynamical systems approach to turbulence, Cambridge University Press, 1998.
  • [15] R. Joly, G. Raugel, A striking correspondence between the dynamics generated by the vector fields and by the scalar parabolic equations, Confluentes Mathematici 3, 471-493 (2011).
  • [16] A. N. Carvalho, J. W. Cholewa, G. Lozada-Cruz, M. R. T. Primo, Reduction of infinite dimensional systems to finite dimensions: Compact convergence approach, Cuadernos de Matemática 11, 281-330 (2010).
  • [17] J. Murray, Mathematical biology, Third Edition, Springer (2007).
  • [18] N. Shigesada, K. Kawasaki, Biological Invasions: Theory and Practice, Oxford University Press (1997).
  • [19] V. Volpert, S. Petrovskii, Reaction-difusion waves in biology, Physics of Life Reviews 6, 267Ð310 (2009).
  • [20] K. R. Swanson, R. Rostomily, E. C. Alvord, Jr, Predicting Survival of Patients with Glioblastoma by Combining a Mathematical Model and Pre-operative MR imaging Characteristics: A Proof of Principle, Br. J. Cancer, 98, 113-119 (2008).
  • [21] C. Wang, J. K. Rockhill, M. Mrugala, D.L. Peacock, A. Lai, K. Jusenius, J. M. Wardlaw, T. Cloughesy, A. M. Spence, R. Rockne, E. C. Alvord Jr., K. R. Swanson, Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model, Cancer Res. 69, 9133-40 (2009).
  • [22] V. M. Pérez-García, G. F. Calvo, J. Belmonte-Beitia, D. Diego, L. A. Pérez-Romasanta, Bright solitons in malignant gliomas, Phys. Rev. E 84, 021921 (2011).
  • [23] J.F. Douglas, K. Efimenko, D. A. Fischer, F. R. Phelan, J. Genzer, Propagating waves of self-assembly in organosilane monolayers, Proc. Nat. Acad. Sci. 104, 10324 (2007).
  • [24] I. Epstein and J. A. Pojman, An Introduction to Nonlinear Chemical Dynamics, Oxford University Press, New York (1998).
  • [25] P. Grindrod, The Theory and Applications of Reaction-Diffusion Equations, Oxford University Press, New York (1996).
  • [26] J. Xin, Front Propagation in Heterogeneous Media, SIAM Rev., 42, 161-230 (2000).
  • [27] C. W. Curtis, D. M. Bortz, Propagation of fronts in the Fisher-Kolmogorov equation with spatially varying diffusion, Phys. Rev. E 86, 066108 (2012).
  • [28] A. Kolmogoroff, I. Petrovsky, and N. Piscounoff, Etude de l’equation de la diffusion avec croissance de la quantite de matiere et son application a un probleme biologique. Moscow University, Bull. Math. 1, 1-25 (1937).
  • [29] J. A. Sherrat, On the transition from initial data to travelling waves in the Fisher-KPP equation, Dynamics and Stability of Systems 13, 167-174 (1998).
  • [30] M. J. Ablowitz, A. Zeppetella, Explicit solutions of Fisher’s equation for a special wave speed, Bull. Math. Biol. 41, 835 (1979).
  • [31] Louis, D. N., Ohgaki, H., Wiestler, O. D., Cavenee, W. K., Burger, P. C., Jouvet, A., Scheithauer, B. W. & Kleihues P, World health organization classification of tumours of the central nervous system, 4th ed., Renouf Publishing Co. Ltd., Geneva. pp. 33-46 (2007)
  • [32] Pignatti, F., Van den Bent, M., Curran, D., Debruyne, C., Sylvester, R., Therasse, P., Afra, D., Cornu, P., Bolla, M., Vecht, C. & Karim, A.B., Prognostic factors for survival in adult patients with cerebral low- grade glioma, J. Clin. Oncol., 20, 2076-84 (2002).
  • [33] Ruiz, J., & Lesser, G. J., Low-Grade Gliomas, Curr. Treat. Opt. Oncol., 10, 231-242 (2009).
  • [34] Pouratian, N., & Schiff, D. Management of low-grade glioma, Curr. Neurol. Neurosci. Rep. 10, 224-231 (2010)
  • [35] Olson, J.D., Riedel E., & DeAngelis, L. M. (2000) Long-term outcome of low-grade oligodendroglioma and mixed glioma. Neurology 54, 1442-1448.
  • [36] Grier, J. T. & Batchelor, T. (2006) Low-Grade Gliomas in Adults. The Oncologist, 11, 681-693.
  • [37] Smith, J.S., Chang, E.F., Lamborn, K.R., Chang, S.M., Prados, M.D., Cha, S., Tihan, T., Vandenberg, S., McDermott, M.W. & Berger, M.S. Role of extent of resection in the long-term outcome of low-grade hemispheric gliomas., J. Clin. Oncol. 26, 1338-1345 (2008).
  • [38] Jakola, A.S., Myrmel, K.S., Kloster, R., Torp, S.H., Lindal, S., Unsgard, G., Solheim, O., Comparison of a strategy favoring early surgical resection vs a strategy favoring watchful waiting in low-grade gliomas, JAMA 308, 1881-1888 (2012).
  • [39] J. D. Murray, Mathematical Biology I and II, 3rd ed, Springer (2003).
  • [40] E. Mandonnet, J-Y Delattre, M-L Tanguy, K. R. Swanson, A. F. Carpentier, H. Duffau, P. Cornu, R. Van Effenterre, E. C. Alvord Jr., and L. Capelle, Continuous growth of mean tumor diameter in a subset of WHO grade II gliomas, Annals of Neurology, 53, 524-528 (2003).
  • [41] S. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K. R. Swanson, M. Pelegrini-Issac, R. Guillevin, H. Benali, Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging, Magnetic Resonance in Medicine, 54, 616-624 (2005).
  • [42] C. Gerin, J. Pallud, B. Grammaticos, E. Mandonnet, C. Deroulers, P. Varlet, L. Capelle, L. Taillandier, L. Bauchet, H. Duffau, M. Badoual, Improving the time-machine: estimating date of birth of grade II gliomas. Cell. Prolif. 45, 76-90 (2012).
  • [43] K. R. Swanson, R. C. Rockne, J. Claridge, M. A. Chaplain, E. C. Alvord Jr, A. R. Anderson, Quantifying the role of angiogenesis in malignant progression of gliomas: in silico modeling integrates imaging and histology., Cancer Research 71, 7366-7375 (2011).
  • [44] V. M. Pérez-García, M. Bogdanska, A. Martínez-González, J. Belmonte, P. Schucht, L. A. Pérez-Romasanta, Delay effects in the response of low grade gliomas to radiotherapy: A mathematical model and its therapeutical implications, Mathematical Medicine and Biology (submitted) (2013).
  • [45] Garcia, D.M., Fulling, K.H. & Marks, J.E., The value of radiation therapy in addition to surgery for astrocytomas of the adult cerebrum, Cancer, 55, 919-927 (1985).
  • [46] Van den Bent, M.J., Afra, D., de Witte, O., Ben Hassel, M., Schraub, S., Hoang-Xuan, K., Malmström, P.O., Collette, L., Piérart, M., Mirimanoff, R. & Karim, A.B., Long-term efficacy of early versus delayed radiotherapy for low-grade astrocytoma and oligodendroglioma in adults: the EORTC 22845 randomised trial., Lancet 366, 985-990 (2005).
  • [47] Higuchi Y., Iwadate Y. & Yamaura A., Treatment of low-grade oligodendroglial tumors without radiotherapy., Neurology 63, 2384-2386 (2004).
  • [48] A. B. M. Karim, B. Maat, et al., A randomized trial on dose-response in radiation therapy of low-grade cerebral glioma: EORTC study 22844, Int. J. Rad. Oncol. Biol. Phys. 36, 549-556 (1996).
  • [49] G Powathil, M Kohandel, S Sivaloganathan, A Oza and M Milosevic, Mathematical modeling of brain tumors: effects of radiotherapy and chemotherapy, Phys. Med. Biol. 52, 3291-3306 (2007)
  • [50] Rockne, R., Hendrickson, K., Lai, A., Cloughesy, T., Alvord, E.C. Jr & Swanson, K.R. Predicting the efficacy of radiotherapy in individual glioblastoma patients in vivo: a mathematical modeling approach, Phys. Med. Biol. 55, 3271-3285 (2010).
  • [51] Bondiau, P.Y., Frenay, M. & Ayache, N., Biocomputing: numerical simulation of glioblastoma growth using diffusion tensor imaging, Phys. Med. Biol. 53, 879-893 (2008).
  • [52] Konukoglu, E., Clatz, O., Bondiau, P.Y., Delingette, H. & Ayache, N., Extrapolating glioma invasion margin in brain magnetic resonance images: suggesting new irradiation margins, Med. Image Anal. 14, 111-125 (2010).
  • [53] Barazzuol, L., Burnet, N. G., Jena, R., Jones, B., Jefferies, S. J. & Kirkby, N. F., A mathematical model of brain tumour response to radiotherapy and chemotherapy considering radiobiological aspects, Journal of Theoretical Biology, 262, 553-565 (2010).
  • [54] Stamatakos, G. S., Antipas, V. P. & Uzunoglu, N. K. (2006) Simulating chemotherapeutic schemes in the individualized treatment context: The paradigm of glioblastoma multiforme treated by temozolomide in vivo. Comp. Biol. Med. 36, 1216-1234.
  • [55] J. Pallud, L. Taillandier, L. Capelle, D. Fontaine, M. Peyre, F. Ducray, H. Duffau, & E. Mandonnet, Quantitative Morphological MRI Follow-up of Low-grade Glioma: A Plead for Systematic Measurement of Growth Rates, Neurosurgery 71, 729-740 (2012).
  • [56] J. Pallud, J. F. Llitjos, F. Dhermain, P. Varlet, E. Dezamis, B. Devaux, R. Souillard-Scemama, N. Sanai, M. Koziak, P. Page, M. Schlienger, C. Daumas-Duport, J. F. Meder, C. Oppenheim, & F. X. Roux, Dynamic imaging response following radiation therapy predicts long-term outcomes for diffuse low-grade gliomas, Neuro-Oncology, 14, 496-505 (2012).
  • [57] Van der Kogel, A., & Joiner, M., Basic clinical radiobiology, Oxford University Press, 2009
  • [58] K. R. Swanson, E. C. alvord, J. D. Murray, Dynamics of a model for brain tumors reveals a small window for therapeutic intervention, Disc. Cont. Dyn. Sys. B 4, 289-295 (2004).
  • [59] S. Jbabdi, E. Mandonnet, H. Duffau, L. Capelle, K. Swanson, M. Pelegrini- Issac, R. Guillevin, H. Benali, Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging, Magnetic Reson. in Med. 54, 616Ð624 (2005).
  • [60] E. Konukoglu, O. Clatz, P.Y. Bondiau, H. Delingette, N. Ayache, Extrapolating glioma invasion margin in brain magnetic resonance images: suggesting new irradiation margins, Med. Image Anal. 14, 111-125 (2010).
  • [61] F. Sánchez-Garduño, P. K. Maini, Travelling wave phenomena in some degenerate reaction-diffusion equations, J. Diff. Eqs. 117, 281-319 (1995).
  • [62] J. Belmonte-Beitia, T. E. Wooley, J. G. Scott, P. K. Maini, E. A. Gaffney, Modelling biological invasions: Individual to population scales at interfaces, Journal of Theoretical Biology (submitted)
  • [63] P. K. Braznik, J. J Tyson, On travelling wave solutions of Fisher’s equation in two spatial dimensions, SIAM J. Appl. Math. 60, 371-391 (1999).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters