﻿Travelling Wave Solution of the Fisher-Kolmogorov Equation with Non-Linear Diffusion

Applied Mathematics
Vol. 4  No. 8A (2013) , Article ID: 35319 , 13 pages DOI:10.4236/am.2013.48A021

Travelling Wave Solution of the Fisher-Kolmogorov Equation with Non-Linear Diffusion

Department of Mathematics, Air University, Islamabad, Pakistan

Email: apshakeel2004@yahoo.com

Copyright © 2013 Muhammad Shakeel. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received May 11, 2013; revised June 11, 2013; accepted June 19, 2013

Keywords: Fisher-Kolmogorov Equation; Non-Linear Diffusion; Travelling Wave; Wave Speed; Pulled Front; Pushed Front

ABSTRACT

In this paper we study one-dimensional Fisher-Kolmogorov equation with density dependent non-linear diffusion. We choose the diffusion as a function of cell density such that it is high in highly cell populated areas and it is small in the regions of fewer cells. The Fisher equation with non-linear diffusion is known as modified Fisher equation. We study the travelling wave solution of modified Fisher equation and find the approximation of minimum wave speed analytically, by using the eigenvalues of the stationary states, and numerically by using COMSOL (a commercial finite element solver). The results reveal that the minimum wave speed depends on the parameter values involved in the model. We observe that when diffusion is moderately non-linear, the eigenvalue method correctly predicts the minimum wave speed in our numerical calculations, but when diffusion is strongly non-linear the eigenvalues method gives the wrong answer.

1. Introduction

For a long time diffusion has been used as a well known model for spatial spread in many biological systems. Okubo [1] and Murray [2] used diffusion in invasion and pattern formation and Skellam [3] studied diffusion in the field of ecology. For motile cell populations diffusion has been used in different situations. Sherrat and Murray [4] used diffusion to model wound healing and Chaplain and Stuart [5] used it to model the capillary growth network.

To study the spatial movement of cell populations, linear diffusion is one of the well known model. However, for closely packed cells such as epithelia, linear diffusion model is unsuitable to examine many biological systems including the movement of cell populations [6]. For closely packed cells, a single cell population can be modeled by a reaction diffusion equation [4,7], but for dissimilar cell populations, a diffusion term would mean that different cell populations can mix together completely and movement of one cells type can be influenced by the cells of other type. But actually the situation is totally opposite, different cell populations cannot move through one another: instead the cell will stop moving when it suddenly comes across another cell. This phenomenon is known as “contact inhibition of migration”. This process has been well documented in many types of cells [8].

In all biological systems the exchange of information at both interand intra-cellular level is almost continuous. In order to get sequential development and generation of the required pattern such communication is necessary e.g. for development and growth of an embryo. Propagating waveforms are one of the ways of conveying such biological information between the cells. Let us consider a simple one-dimensional diffusion equation

(1)

where is chemical (cell or nutrient) concentration and is diffusion coefficient. The time to exchange information in the form of changed concentration is, where is the length of domain. We can get this order by dimensional arguments of Equation (1). At the early stage of growth of an organism the diffusion coefficient can be very small: values of order 109 to 1011 cm2∙sec1. If diffusion is the main process to convey the biological information then to cover macroscopic distances of several millimeters requires a very long time. When the diffusion coefficient is O (109 to 1011 cm2∙sec1) and is order of 1 mm then the time required to convey the information is O (107 to 109 sec), which is very large in the early stages of growth of an organism. This means that simple diffusion is unlikely to be the main means of exchanging the information during embryogenesis. Karevia [9] and Tilman and Karevia [10] estimated the diffusion coefficient for insect dispersal in interacting population.

About seventy years ago Fisher [11] and Kolmogorov et al. [12] introduced a classical model to describe the propagation of an advantageous gene in a one-dimensional habitat. The equation describing the phenomenon is a one-dimensional non-linear reaction-diffusion equation,

(2)

where is chemical concentration, is the diffusion coefficient and the positive constant represents the growth rate of the chemical reaction. Since then a great deal of work has been carried out to extend their model to take into account the other biological, chemical and physical factors. The Equation (2) is also used in flame propagation [13], nuclear reactor theory [14], autocatalytic chemical reactions [15,16], logistic growth models [17] and neurophysiology [18].

One of the extensions of the Fisher and Kolmogorov model is to introduce a non-linear diffusion coefficient, which can be taken as a non-linear Fick’s law. The nonlinearity can arise in terms of space, time or density-dependent diffusion coefficient. In the Fisher and Kolmogorov model the reaction kinetics are coupled to diffusion which gives travelling waves of chemical concentration and it can affect biological change very much faster as compared to the processes governed by simple diffusion without the kinetic term. In many biological populations density-dependent dispersal has been observed e.g. Carl [19] observed that ground squirrels move from highly populated area to sparsely populated areas, Myers and Krebs [20] studied the population density cycles in small rodents. Several mathematical models have been developed to describe the density-dependent dispersal systems. Gurney and Nisbet [21,22] developed a first density-dependent diffusion model in ecological context by using a random walk approach and Montroll and West [23] developed a one-dimensional model for a single species by using the same approach. But there is not much work on the Fisher equation with non-linear diffusion. Sánchez and Maini [24] studied a travelling wave solution in a degenerate Fisher equation with non-linear diffusion.

Although Equation (2) is known as the Fisher-Kolmogorov equation, the discovery, investigation and analysis of travelling waves in chemical reactions was first reported by Luther [25]. He found that the wave speed is a simple consequence of the differential equations. He obtained the wave speed in terms of parameters associated with the reactions he was studying. The analytical form is the same as that found by Kolmogorov [12] and Fisher [11].

Seeking the travelling wave solution of Fisher equation has been a challenging and difficult task. Several authors used different methods to find the travelling wave solution of Fisher equation. Ablowitz and Zeppetella [26] obtained the first explicit form of travelling wave solution of Fisher equation using Painlevé analysis. The authors found that kink wave propagates from left to right with a speed. Wazwaz [27] used tanh method to find the travelling wave solution of non-linear partial differential equation. The author also extended this study to the equations which do not have tanh polynomials. He also demonstrated the efficiency of the method by applying it to variety of equations e.g. Fisher equation. Mansour [28] found a travelling wave solution of class of degenerate reaction diffusion equations in which non-linearity is assumed to be singular at zero. He employed two methods to find the wave profile and speed of the front. These method include travelling wave equation and initial boundary value problem with an adaptive travelling wave condition for partial differential equation. Tan et al. [29] used an analytical method namely Homotopy analysis method to solve the Fisher equation and they described a family of travelling wave solutions. Jiaqi et al. [30] found a travelling wave solution of non-linear reaction diffusion equation by using the homotopic method and theory of travelling wave transform. Feng et al. [31] found the travelling wave solution of reaction diffusion equation by two methods. Firstly they applied Divisor theorem for two variables in complex domain to find a quasi-polynomial first integral of an explicit form to an equivalent autonomous system. Then through this first integral, they reduced the reaction-diffusion equation to a first-order integrable ordinary differential equation, and found a class of traveling wave solutions. They compared their results with the existing results and found some errors in analytic results in literature. They clarified and introduced the refined results. Hou et al. [32] studied the travelling wave solution of reaction diffusion equation with double degenerate nonlinearities. They investigated existence, uniqueness and stability of wave solution.

In this paper we study the growth of cells in one-dimensional scaffold in the presence of constant nutrients. We model a system in which cell growth and diffusion takes place simultaneously but the cell diffusion depends on cell density. It increases with increase in cell density. The equation governing this system is a non-linear Fisher equation in which the diffusion coefficient depends on cell density. The non-linearity arises in the diffusion coefficient. We represent the diffusion coefficient as an exponential function of cell density. The form of nonlinear diffusion captures the feature that it is very small for small cell density and increases with the increase in cell density and it is maximum when cells stop growing. Results are presented for different values of the parameters. The Fisher equation has such wide applicability in itself but also it is the prototype equation which admits travelling wavefront solutions. We study the travelling wave solution of this equation and find the approximation for the minimum wave speed. We find that the minimum wave speed depends on the model parameters. For moderately non-linear systems the analytical method correctly predicts the wave speed in our numerical calculations. However for strongly non-linear diffusion the numerical wave speed becomes larger than the analytic minimum speed.

The paper is organized as follows in Sections 2 and 3 we present the dimensional and nondimensional model equations respectively. The parameter values and travelling wave solution are described in Sections 4 and 5, respectively. We present the numerical solution in Section 6 and calculate the numerical minimum wave speed in Section 7.

2. Mathematical Model

In this Section we model cell growth in a bioreactor subject to uniform nutrient concentration. Consider the cells are seeded onto a porous scaffold with interconnected pores, which is placed in the bioreactor. The scaffold extends from, where denotes the spatial coordinate. In this model we represent the cell density by, initial cell density by, maximum carrying capacity of the system by and diffusion coefficient by. We assume that the environment is inhomogeneous i.e. the cell density and diffusion coefficient depends on spatial coordinate and time, i.e. and. We know that the cells require some nutrients such as oxygen and glucose etc. to live and perform specific functions. Suppose that the concentration of such nutrients remains uniform everywhere in the entire domain at all times. We want to model a system in which the change in cell density is due to cell proliferation and cells disperse in the entire domain due to diffusion. We assume that when the cell density is small diffusion is also very small and when the cell density reaches its maximum carrying capacity the cell proliferation stops and cells spread in the entire domain via diffusion. For that we consider a logistic growth model in which the cell population spreads via diffusion. That means we have a coupled system of reaction kinetics and diffusion. Our modelling of cell proliferation [33,34] has led us to the following Fisher equation with non-linear diffusion.

(3)

where is the maximum cell growth rate and is the non-linear diffusion. We observe that cells grow in numbers due to second term on the right hand side of Equation (3) and they disperse in the entire domain due to first term on the right hand side of Equation (3). To model cell proliferation we need to choose a form of so that, in regions of low cell density, cells grow by increasing their density with little or no spatial motion, whereas in regions of high cell density the newly formed cells diffuse rapidly towards regions of lower cell density. Thus to model this behaviour on a continuum level we choose the following form for the non-linear diffusion.

(4)

where parameter is the cell diffusion when cell density reaches its maximum carrying capacity, so we can call this the maximum value of cell diffusion and the parameter controls how fast cells are spreading in the domain. The parameters and have dimensions m2/sec and m3/cell respectively. The exponential function in Equation (4) ensures that cell diffusion will remain positive. The non-linear cell diffusion also has maximum value when. So we can say that non-linear cell diffusion is maximum either when cell density reaches its maximum carrying capacity or the value of parameter is zero. The non-linear diffusion becomes linear when. It is also clear from the Equation (4) that when there are no cells then in that case cell diffusion is, which is the minimum value of diffusion. So cell diffusion varies in the range,

(5)

Initially when the cell density is small, the diffusion coefficient is also small and increases with increasing cell density but it always remains positive.

The initial cell density is given by

(6)

3. Nondimensionalization

For convenience and to reduce the parameters in the equation we rescale all the variables to analyze the nondimensional form of the equation. We nondimensionalize all lengths by which we take as a typical length scale for the problem and all cell densities by the maximum carrying capacity.

(7)

(8)

We nondimensionalize time by the speed of the growth front (which is given by Equation (30)),

(9)

In dimensionless form Equation (3) can be written as

(10)

where and. The parameters and are dimensionless parameters in the model which are given by

(11)

The parameter is the ratio of cell growth rate to speed of growth front and parameter is the ratio of maximum value of non-linear diffusion to the speed of growth front. The parameter controls the diffusion and the parameter controls the growth term. The initial condition (6) in dimensionless form becomes

(12)

4. Parameter Values

The modified Fisher Equation (10) includes a number of parameters. Some parameters depend on the cell type and some parameters depend on scaffold geometry. Table 1 shows the values of the parameters used in the simulation. We assume that the length of the scaffold is 0.02 m. Some quantities such as cell growth rate depend on the cell type cultured in the bioreactor. The above proposed model is a generic model and can be applied to any cell type. To compare the model with the experimental data the cells used in the simulations are Murine immortalized rat cell. The maximum cell growth rate for cells is [35]. Since cell growth is a slow process we can choose that speed of growth front

is very small e.g.. The value of dimensionless parameter can be obtained by using the values of dimensional parameters, and. The values of parameters and are not available in the literature. To estimate the value of these parameters we use the expression for dimensionless wave speed i.e. (which we will calculate later in Section 5.2). We choose that the theoretical wave speed. In the expression for there are two unknowns and. In order to find the value of and we fix one of the parameters or in the expression for. If we fix the parameter then using the expression for we can find the value of parameter. We observe that the value of dimensionless parameter depends on the value of parameter. If the value of parameter is high then value of is also high. This means that cells will diffuse more quickly for high values of parameter. Table 1 shows the values of the dimensional parameters used in the model. We know that the initial cell density. We can use any form of.

5. Travelling Wave Solution

We are interested in travelling wave solutions because Equation (10) has two steady states and, which are, respectively, unstable and stable. This suggests that we should look for a travelling wave solution of Equation (10) for which; since negative cell density has no physical meaning. If a travelling wave solution exists it can be written in the form

(13)

where is the wave speed. Since Equation (10) is invariant if, may be negative or positive. To be specific we assume that. We have

(14)

Substituting the travelling wave solution (13) into Equation (10) and using relations from Equation (14) we get a second order ordinary differential equation,

Table 1. Model parameters used in this work.

(15)

where

(16)

where primes denote the differentiation with respect to. A typical wave front solution is a solution such that approaches one steady state as and approaches the other steady state as. Therefore the boundary conditions for travelling wave solution are

(17)

5.1. Phase Plane Analysis

From Equation (15) we observe that we have to determine the values of the wave speed such that a nonnegative solution of Equation (15) exists. We use phase plane analysis to characterize solutions of Equation (15) in the phase plane where,

(18a)

(18b)

The system of ordinary differential Equations (18) is nearly singular at, since for high values of parameter. To remove the near singularity we introduce a new parameter [24], in such a way that

(19)

Except at, where is not defined,.

Thus we have

(20)

and we obtain

(21)

Substituting and into the system of Equations (18) we get

(22a)

(22b)

where dot denotes the derivative with respect to.

The phase trajectories of (22) are solutions of

(23)

The fixed points are the points where, these are steady states. So in this case fixed points are and. The local behaviour of the trajectories of system (22) can be obtained by analyzing the linear approximation of system (22) around each fixed point.

5.2. Stability of Fixed Points

The system of Equations (22) has two fixed points, but which of these fixed points are stable? The local stability of a fixed point is determined by linearization of the dynamics at the intersection. So with the linear approximations the system of Equation (22) becomes

(24a)

(24b)

which can be written in the matrix form as

(25)

where

(26)

Let and be eigenvalues of then we have

(27)

where

Eigenvalues for the fixed points (0, 0) and (1, 0) are

(28a)

(28b)

We know that and. Using values of and in Equation (28) we get

(29a)

(29b)

It is clear from the Equation (29a) that fixed point is a stable node if

, with the case when giving a degenerate node. The fixed point (0, 0) is a stable spiral if or;

i.e. oscillates in the vicinity of origin. When

then it is not physically realistic because cannot be negative. The fixed point is a saddle point. The solution of our modified Fisher equation evolve to a travelling wave if the fixed point is a stable node and minimum wave speed of wave front is. The condition that is a stable node is a necessary condition for travelling wave propagation but not sufficient. Sometimes it gives wrong answer as we shall in see Section 7. If the propagation speed of the front is determined by the leading edge of the population distribution the non-linear fronts whose asymptotic propagation speed equals is called a “pulled front”. In this case the eigenvalues give the right answer because the wave speed is determined by what happens at the front edge where. As above is determined by the eigenvalues thus we call speed of pulled front. On the other hand if the speed of front is determined by the whole wave-front but not the behaviour of the leading edge the front is called the “Pushed front” [36,37]. In this case simply computing the eigenvalues gives the wrong answer for the minimum wave speed and the minimum wave speed is.

Thus the wave speed for a pulled front is given by

. The wave speed depends on the parameters, and, is directly proportional to the square root of the product of and. In terms of the original dimensional Equation (3) the wave speed is

(30)

We observe from the system of equations (28) that for a general function the wave speed is determined by. So in general the wave speed of the Pulled front is.

Figure 1 shows the phase plane sketch of the trajectories of Equation (22) when. We see that when

Figure 1. Phase plane trajectories of Equation (22). Here parameter values are χ = 132.1739, δ = 0.0139 γ = 2 and v = 1.5 > vc. For explanation of parameter values see Section 4.

the fixed point is a stable node because all the trajectories from to have the same limiting direction towards and the fixed point point is a saddle point because there are two incoming trajectories and two out going trajectories and all the other trajectories in the neighborhood of the critical point bypass. Similarly Figure 2 shows the phase plane sketch of the trajectories of Equation (22) when. We observe that when the fixed point is a stable spiral because all the trajectories from to spiral around the point.

Figures 3 and 4 show the phase plane sketch of trajectories of Equation (22) for various wave speeds and respectively. We observe from the Figure 3 that when all the trajectories in the phase plane from to remain entirely in the quadrant where and, with for all wave speeds. Similarly from the Figure 4 we see that for all wave speeds the phase trajectories from to spiral around the fixed point. In this case oscillates in the vicinity of the origin giving negative which is unphysical.

Figure 2. Phase plane trajectories of Equation (22). Here parameter values are χ = 132.1739, δ = 0.0139, γ = 2 and v = 0.5 < vc. For explanation of parameter values see Section 4.

Figure 3. Phase plane trajectories of Equation (22) for different values of v ≥ vc. The other parameter values are same as in Figure 1. Each curve (starting from bottom) represents the trajectory for various values of speed v ≥ vc, i.e. v = 1, v = 1.5, v = 2, v = 2.5 and v = 3 respectively.

Figure 4. Phase plane trajectories of Equation (5.1) for different values of v < vc. The other parameter values are same as in Figure 1. Each curve (starting from bottom) represents the trajectory for various values of speed v < vc, i.e. v = 0.5, v = 0.6, v = 0.7 and v = 0.8 respectively.

5.3. Selection of Initial Condition

A very important question at this stage is what kind of initial condition will evolve into the travelling wave solution and if the travelling wave solution exists what is its wave speed? Fisher [11] found that Equation (10) has an infinite number of travelling wave solutions for which for all wave speeds. Kolmogorov [12] proved that Equation (10) has a travelling wavefront solution and the wave speed is, if has compact support. A function is said to have a compact support if

(31)

where

where and is continuous in. If the initial condition is other than (36) then solution depends on the behaviour of.

If then Equation (15) reduces to

(32)

A travelling wave solution of Equation (32) in explicit form for was found by Ablowitz and Zeppetella [26] for the special wave speed,

(33)

But if is not a constant then it is not possible to find the exact solution of Equation (15). Solution of such non-linear problem can be approximated by perturbation theory or numerical investigation.

6. Numerical Solution

Numerically we solve the modified Fisher Equation (10) by using the commercial finite element solver COMSOL. We subdivide the domain into a suitable number of mesh elements (intervals) of length. The end points of each interval are called node points and the elements do not have to have the same length. But in this case the length of each element is same. To obtain meaningful results care is required in the choice of a suitable number of mesh elements, finite element approximation and model parameters. Convergence is achieved by successively refining the mesh elements. The refined mesh contains 30,721 mesh vertices and 30,720 mesh elements. The dependent variable is approximated by a quadratic shape function and solved for 61,441 degrees of freedom. We assume that at time the cell density is and after the time the cell density is. We start with initial cell density and after each time we replace by and solve the Equation (10) again for updated cell density. The time from to is subdivided as, where is the time step size from to and is the time when we update the cell density. The cell density at each mesh point is obtained for different times. To estimate the wave speed numerically we look for the point after each time where the cell density is half of its maximum value i.e.. When the cell density is half of its maximum value at time, then and at time,. So we can estimate the total distance traveled by the wave in the time interval. Hence

(34)

Results are plotted for different values of the dimensionless parameters and in Sections 6.1 to 6.4. In the next Section we consider various cases in which we fix the value of dimensionless parameter χ and vary the values of dimensionless parameter and parameter in such a way that the theoretical wave speed

remains 1. Let us assume thatwhere N0 and r are constants and is the Heaviside step function.

6.1. Case I: χ = 1.3217

In this case we fix the value of dimensionless parameter and find the values of parameter and such that theoretical speed of growth front is 1. Table 2 shows the values of dimensionless parameter for corresponding values of parameter.

Figures 5 and 6 show the numerical solution of the modified Fisher Equation (10) for the parameter values given in the Table 2. We observe from Figures 5 and 6 that the solution does not evolve into a travelling wave solution before it reaches the domain edges. The cell density drops down because the rate of diffusion is larger than the growth rate, which means that the dimensionless parameter is larger than the dimensionless parameter. Diffusion dominates in this simulation so we need to look at the case where diffusion is smaller and cell growth is larger, to find a travelling wave. Numerical results are plotted for and in Figures 5 and 6. In a larger domain and with a longer time of solution the system will eventually evolve into a travelling wave but we are interested in a finite domain.

6.2. Case II: χ = 13.2173

We consider the case when the dimensionless parameter

Table 2. Values of and for.

Figure 5. Numerical results of profile of cell density at different times when γ = 1, δ = 0.5141, χ = 1.3217. Initial cell density is, where, and r2 = 0.1. The time step size Δt = 0.001 and cell update time tnew = 0.01. The figure shows the cell distribution after each time tnew and final time is t = 0.3.

Figure 6. Numerical results of profile of cell density N at different times when γ = 2, δ = 1.3976, χ = 1.3217. Ninit, Δt and tnew are same as in Figure 5.

and find the value of parameter and such that theoretical speed of growth front is 1. Table 3 shows the values of dimensionless parameter for the corresponding values of parameter.

Numerical results of modified Fisher Equation (10) are plotted in Figures 7 and 8 for parameter values given in Table 3. It is clear from the Figures 7 and 8 that in this case the solution evolves into travelling wave fronts. So in this case the growth term is somewhat larger than the diffusion term. But the diffusion term is not too small because the cells also diffuse very quickly. This feature is evident from the Figures 7 and 8 because at the edges the front are smooth.

6.3. Case III: χ = 132.1739

Consider the case when the value of dimensionless parameter is very high i.e. and we find the values of parameters and such that theoretical speed of growth front is 1. Table 4 shows the values of dimensionless parameter for corresponding values of parameter.

In the modified Fisher Equation (10), growth and diffusion are taking place simultaneously. But if the diffusion is very small compared to growth, then first cells will grow quickly and when they reach maximum carrying capacity growth stops and they spread in the domain via diffusion. In the present case the growth term is very big as compared to the diffusion term.

Figures 9 and 10 show the numerical solution of the modified Fisher’s Equation (10) for and values of parameters and given in Table 4. The wave front takes some time to settle down to a travelling wave, which moves at a constant speed. When the number of cells reaches its maximum limit the proliferation

Table 3. Values of and for.

Figure 7. Numerical results of profile of cell density N at different times when, ,. Ninit, Δt and tnew are same as in Figure 5. In this case the final time is t = 0.6.

Figure 8. Numerical results of profile of cell density N at different times when, ,. Ninit, Δt and tnew are same as in Figure 5. In this case the final time is t = 0.6.

Table 4. Values of and for.

stops and then the cells spread via diffusion in the entire domain. In this case diffusion term is much smaller than the growth term, and due to this reason the shape of the front is very sharp. It is evident from the Figure 9 and 10 that when the front settles down then it moves with

Figure 9. Numerical results of profile of cell density N at different times when, ,. Ninit, Δt and tnew are same as in Figure 5. In this case the final time is t = 0.6.

. Numerical results of profile of cell density N at different times when, ,. Ninit, Δt and tnew are same as in Figure 5. In this case the final time is t = 0.6.

constant speed and shape.

6.4. Case IV: χ = 0

If then it represents the case when there is no growth of cells. The number of cells will not increase because of the absence of growth term. In that case for every value of the solution will not evolve into a travelling wave. Figures 11 and 12 show the cell density in the pure diffusion case i.e. Fisher equation without the growth term for the same times and parameter values used in Figure 7 except in this case. The behaviour of the solution is different from the Fisher equation

. Numerical results of profile of cell density N at different time without growth term when γ = 1, δ = 0.05141, χ = 0. Ninit, Δt and tnew are same as in Figure 5.

. Numerical results of profile of cell density at different time without growth term when γ = 2, δ = 0.139760, χ = 0. Ninit, Δt and tnew are same as in Figure 5.

with the growth term. Clearly the solution does not grow due to the absence of growth term and the behaviour of the solution is not wave-like.

7. Numerical Minimum Wave Speed

From the phase plane analysis it is clear that a travelling wave front solution exists for a range of wave speeds. We choose the values of parameters and such that the theoretically wave speed. If diffusion is linear in the modified Fisher Equation (10) then travelling wave move with the minimum wave speed [2]. In this model corresponds to linear diffusion. But when the value of then diffusion is no longer linear. Here we see that the numerical wave speed agrees with the theoretical value when. However at the numerical speed begins to diverge from the theoretical value and become increasingly large as increases.

Figures 13 and 14 show the phase plane sketch of the trajectories of Equation (10) for and respectively. We observe from the Table 5 that when and the numerical value of the minimum wave speed. We see that from that when then all trajectories from to remain entirely in the region where and for all wave speed. Similarly from we observe that when wave speed then for all the trajectories from to, becomes negative, which is unphysical. is still a stable node. We observe that method of finding by looking at the eigenvalues gives the wrong answer. It is beyond the scope of this work to find the analytical formula for the minimum wave speed, because the solution depends on the whole trajectory, but we have a very good agreement between numerical results and phase plane analysis.

. Phase plane trajectories of Equation (24) for different values of. The other parameter values are χ = 132.1739, γ = 3, and δ = 0.037990. Each curve (starting from bottom) represents the trajectory for various values of velocity i.e. v = 1.15, v = 1.2, v = 1.3 and v = 1.5.

. Phase plane trajectories of Equation (24) for different values of. The other parameter values areχ = 132.1739, γ = 3, and δ = 0.037990. Each curve (starting from bottom) represents the trajectory for various values of velocity i.e. v = 1, v = 1.02, v = 1.05 and v = 1.09.

Table 5. Numerical minimum wave speed.

shows the shape of growth front at time for fixed value of but different values of and. It is evident from that the speed of the growth front is not the same for all values of parameter. The speed of the growth front depends on the value of and increases as value of increases.

8. Conclusion

In this paper we studied the modified Fisher-Kolmogorov equation with non-linear diffusion, in detail. The diffusion coefficient in this case is non-linear, depending on the cell density. The form of non-linear diffusion is such that it produces similar behaviour to cell proliferation.

. Shape of growth front at time t = 0.3 for fixed and value of for corresponding value of are given in Table 3. Ninit, Δt and tnew are same as in Figure 5.

The diffusion is high where the cell density is high and it is low in the regions of low cell density but it remains positive. This leads us to choose the non-linear diffusion to be an exponential function of cell density. We found a travelling wave solution of the Fisher equation and found the theoretical minimum wave speed of growth front by using an eigenvalue analysis of stationary points. The results reveal that for moderately non-linear diffusion the wave speed found by eigenvalue method agrees with the results of the wave speed found numerically. The waves of this kind are called “Pulled front” [38]. However for highly non-linear diffusion, numerical minimum wave speed of growth front is greater than the minimum speed of growth front found by using the eigenvalue analysis. The waves of this kind are called “Pushed front” in which the wave speed is determined by non-linear effects. We tested the system for various values of parameters and found that the system exhibits the travelling wave solution when the growth term is dominant over the diffusion term. We also found that the minimum wave speed of growth front depends on the values of parameter and it increases as value of increases. The Fisher equation is widely used in modelling the physical phenomenon where diffusion and growth are taking place simultaneously. This equation is very useful in modelling the cell growth in in vitro tissue engineering, because the form of diffusion used in this model produces the similar behaviour as cell proliferation.

REFERENCES

1. A. Okubo, “Diffusion and Ecological Problems: Mathematical Models,” Springer-Verlag, Berlin, 1980.
2. J. D. Murray, “Mathematical Biology,” Biomathematics Texts, Vol. 19, Springer Verlag, New York, 1989.
3. J. G. Skellam, “Random Dispersal in Theoretical Populations,” Bulletin of Mathematical Biology, Vol. 53, No. 1, 1991, pp. 135-165.
4. J. A. Sherratt and J. D. Murray, “Models of Epidermal Wound Healing,” Proceedings: Biological Sciences, Vol. 241, No. 1300, 1990, pp. 29-36.
5. M. A. J. Chaplain and A. M. Stuart, “A Model Mechanism for the Chemotactic Response of Endothelial Cells to Tumour Angiogenesis Factor,” Mathematical Medicine and Biology, Vol. 10, No. 3, 1993, pp. 149-168.
6. J. A. Sherratt, “Wavefront Propagation in a Competition Equation with a New Motility Term Modelling Contact Inhibition between Cell Populations,” Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, Vol. 456, No. 2002, 2000, pp. 2365-2386.
7. M. A. J. Chaplain and A. M. Stuart, “A Mathematical Model for the Diffusion of Tumour Angiogenesis Factor into the Surrounding Host Tissue,” Mathematical Medicine and Biology, Vol. 8, No. 3, 1991, pp. 191-220.
8. M. Abercrombie, “Contact Inhibition in Tissue Culture,” In Vitro Cellular & Developmental Biology-Plant, Vol. 6, No. 2, 1970, pp. 128-142.
9. P. M. Kareiva, “Local Movement in Herbivorous Insects: Applying a Passive Diffusion Model to Mark-Recapture Field Experiments,” Oecologia, Vol. 57, No. 3, 1983, pp. 322-327.
10. D. Tilman and P. M. Kareiva, “Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions,” Princeton University Press, Princeton, 1997.
11. R. A. Fisher, “The Wave of Advance of Advantageous Genes,” Annals of Eugenics, Vol. 7, No. 4, 1937, pp. 353- 369.
12. A. Kolmogorov, I. Petrovsy and N. Piskounov, “Study of the Diffusion Equation with Growth of the Quantity of Matter and Its Applications to a Biological Problem,” Moscow University Mathematics Bulletin, 1989, p. 105.
13. D. A. Frank-Kamenetskii, “Diffusion and Heat Exchange in Chemical Kinetics,” Princeton University Press, Princeton, 1955.
14. J. Canosa, “Diffusion in Non-Linear Multiplicative Media,” Journal of Mathematical Physics, Vol. 10, No. 10, 1969, pp. 1862-1868.
15. H. Cohen, “Non-Linear Diffusion Problems,” In: A. H. Taub, Ed., Studies in Applied Mathematics, The Mathematical Association of America, 1971, pp. 27-64.
16. P. C. Fife and J. B. McLeod, “The Approach of Solutions of Non-Linear Diffusion Equations to Travelling Front Solutions,” Archive for Rational Mechanics and Analysis, Vol. 65, No. 4, 1977, pp. 335-361.
17. J. D. Murray, “Lectures on Non-Linear Differential-Equation Models in Biology,” Oxford University Press, Oxford, 1977.
18. H. C. Tuckwell, “Introduction to Theoretical Neurobiology: Non-Linear and Stochastic Theories,” Cambridge University Press, Cambridge, 1988.
19. E. A. Carl, “Population Control in Arctic Ground Squirrels,” Ecology, Vol. 52, No. 3, 1971, pp. 395-413.
20. J. H. Myers and C. J. Krebs, “Population Cycles in Rodents,” Scientific American, Vol. 230, No. 6, 1974, p. 38.
21. W. S. Gurney and R. M. Nisbet, “The Regulation of Inhomogeneous Populations,” Journal of Theoretical Biology, Vol. 52, No. 2, 1975, pp. 441-457.
22. W. S. Gurney and R. M. Nisbet, “A Note on Non-Linear Population Transport,” Journal of Theoretical Biology, Vol. 56, No. 1, 1976, pp. 249-251.
23. E. W. Montroll and B. J. West, “On an Enriched Collection of Stochastic Processes,” Fluctuation Phenomena, Vol. 66, 1979, p. 61. doi:10.1016/B978-0-444-85248-9.50005-4
24. F. Sánchez-Garduño and P. K. Maini, “Existence and Uniqueness of a Sharp Traveling Wave in Degenerate NonLinear Diffusion Fisher-KPP Equations,” Journal of Mathematical Biology, Vol. 33, No. 2, 1994, pp. 163-192.
25. R. Luther, “Propagation of Chemical Reactions in Space,” Zeitschrift für Elektrochemie und Angewandte Physikalische Chemie, Vol. 12, No. 32, 1906, p. 596. doi:10.1002/bbpc.19060123208
26. M. J. Ablowitz and A. Zeppetella, “Explicit Solutions of Fisher’s Equation for a Special Wave Speed,” Bulletin of Mathematical Biology, Vol. 41, No. 6, 1979, pp. 835-840.
27. A. M. Wazwaz, “The tanh Method for Traveling Wave Solutions of Nonlinear Equations,” Applied Mathematics and Computation, Vol. 154, No. 3, 2004, pp. 713-723.
28. M. B. A. Mansour, “Accurate Computation of Traveling Wave Solutions of Some Nonlinear Diffusion Equations,” Wave Motion, Vol. 44, No. 3, 2007, pp. 222-230.
29. Y. Tan, H. Xu and S. J. Liao, “Explicit Series Solution of Travelling Waves with a Front of Fisher Equation,” Chaos, Solitons & Fractals, Vol. 31, No. 2, 2007, pp. 462-472.
30. J. Q. Mo, W. J. Zhang and M. He, “Asymptotic Method of Travelling Wave Solutions for a Class of Nonlinear Reaction Diffusion Equation,” Acta Mathematica Scientia, Vol. 27, No. 4, 2007, pp. 777-780.
31. Z. Feng, S. Zheng and D. Y. Gao, “Traveling Wave Solutions to a Reaction Diffusion Equation,” Zeitschrift für Angewandte Mathematik und Physik, Vol. 60, No. 4, 2009, pp. 756-773.
32. X. Hou, Y. Li and K. R. Meyer, “Traveling Wave Solutions for a Reaction Diffusion Equation with Double Degenerate Nonlinearities,” Discrete and Continuous Dynamical Systems, Vol. 26, No. 1, 2010, pp. 265-290.
33. M. Shakeel, P. C. Matthews, S. L. Waters and R. Graham, “Continuum Modeling of Cell Growth and Nutrient Transport in a Perfusion Bioreactor,” 2011, pp. 90-117.
34. M. Shakeel, “Continuum Modelling of Cell Growth and Nutrient Transport in a Perfusion Bioreactor,” Ph.D. Thesis, The University of Nottingham, Nottingham, 2011.
35. F. Coletti, S. Macchietto and N. Elvassore, “Mathematical Modeling of Three Dimensional Cell Cultures in Perfusion Bioreactors,” Industrial & Engineering Chemistry Research, Vol. 45, No. 24, 2006, pp. 8158-8169.
36. F. Rothe, “Convergence to Pushed Fronts,” Rocky Mountain Journal of Mathematics, Vol. 11, No. 4, 1981, pp. 617-634.
37. W. Van Saarloos, “Front Propagation into Unstable States,” Physics Reports, Vol. 386, No. 2-6, 2003, pp. 29-222.
38. A. N. Stokes, “On Two Types of Moving Front in Quasilinear Diffusion, Mathematical Biosciences, Vol. 31, No. 3, 1976, pp. 307-315.