On Waves and Distributions in Population Dynamics

We discuss dynamics of spatially distributed interacting populations. The following two cases are discussed: (I.) migration waves for the case or negligible random fluctuations of the populations densities and (II.) probability distributions of the spatially averaged populations densities for the case of significant random fluctuations of these densities. For each of the cases we obtain in general a system of differential equations that is treated analytically (when possible) or numerically. The obtained results are discussed from the point of view of population dynamics. Keywords-interacting populations; nonlinear waves; method of simplest equation; density fluctuations; probability density functions for populations densities I. I NTRODUCTION In many cases the dynamics of interacting populations is studied on the basis of models consisting of equations that contain only time dependence of the population densities [1]-[4]. These models are very useful for understanding the complex dynamics of the interacting populations but they neglect two important aspects of this dynamics: (I.) the possible influence of spatial characteristics of the environment; and (II.) the possible fluctuations of the population densities caused by different factors. Below we shall investigate two kinds of population dynamics models that account for each of these effects. First of all we shall discuss the dynamics of spatially distributed populations and this will be a continuation of our previous work [5], [6]. Then we shall show that by appropriate averaging the spatial model can be reduced to model in which the population densities depend only on the time. This model is valid not only for the cases of small values of the densities of the populations and because of this we shall discuss the influence of random fluctuations of the population densities for arbitrary values of the densities. The result of the action of these fluctuations is that instead of equations for the trajectories of the populations in the phase space of the population densities we shall write and solve equations for the probability density functions of the densities of the interacting populations. II. M ODEL EQUATIONS A. Spatially distributed populations Let us consider an two-dimensional area S where N competing populations are present. The density of each population isρi(x, y, t) = ∆Ni ∆S , where∆Ni is the number of the individuals of the i-th population that are present in the small area ∆S at the moment. Now let a movement of population members through the borders of Citation: N. Vitanov, Z. Dimitrova, On Waves and Distributions in Population Dynamics, Biomath 1 (2012), 1209253, http://dx.doi.org/10.11145/j.biomath.2012.09.253 Page 1 of 7 N. Vitanov et al., On Waves and Distributions in Population Dynamics the area∆S be possible and let ~ji(x, y, t) be the current of this movement. Then(~ji · ~n)δl is the net number of members fromi-th population, crossing a small border line δl with normal vectorn. Let the density changes (other then these caused by the border crossings) be summarized by the functionCi(ρ1, ρ2, . . . , ρN , x, y, t). Then the change of the density of members of the ith population in the studied area is described by the equation ∂ρi ∂t + div~ji = Ci. (1) Below we shall discuss the case where ~ji has the form of linear multicomponent diffusion. In this case


I. INTRODUCTION
In many cases the dynamics of interacting populations is studied on the basis of models consisting of equations that contain only time dependence of the population densities [1]- [4].These models are very useful for understanding the complex dynamics of the interacting populations but they neglect two important aspects of this dynamics: (I.) the possible influence of spatial characteristics of the environment; and (II.) the possible fluctuations of the population densities caused by different factors.Below we shall investigate two kinds of population dynamics models that account for each of these effects.First of all we shall discuss the dynamics of spatially distributed populations and this will be a continuation of our previous work [5], [6].Then we shall show that by appropriate averaging the spatial model can be reduced to model in which the population densities depend only on the time.This model is valid not only for the cases of small values of the densities of the populations and because of this we shall discuss the influence of random fluctuations of the population densities for arbitrary values of the densities.The result of the action of these fluctuations is that instead of equations for the trajectories of the populations in the phase space of the population densities we shall write and solve equations for the probability density functions of the densities of the interacting populations.

A. Spatially distributed populations
Let us consider an two-dimensional area S where N competing populations are present.The density of each population is ρ i (x, y, t) = ∆Ni ∆S , where ∆N i is the number of the individuals of the i-th population that are present in the small area ∆S at the moment t.Now let a movement of population members through the borders of the area ∆S be possible and let j i (x, y, t) be the current of this movement.Then ( j i • n)δl is the net number of members from i-th population, crossing a small border line δl with normal vector n.Let the density changes (other then these caused by the border crossings) be summarized by the function C i (ρ 1 , ρ 2 , . . ., ρ N , x, y, t).Then the change of the density of members of the ith population in the studied area is described by the equation Below we shall discuss the case where j i has the form of linear multicomponent diffusion.In this case where D ik is the diffusion coefficient.Eq. ( 2) accounts for the possibility that the spatial motion of the population members happens not only as consequence of gradient of density of the own population but also as consequence of gradients of the densities of the other populations.We shall not specify the kind of the function C i as we shall consider the general case of relative small populations densities that allow us to write C i as Taylor series expansion around the zero values of all population densities as follows where the constant coefficients α (i) n1,n2,...,nN are as follows We shall discuss the one-dimensional case and in addition we shall assume that the diffusion coefficients D ik are constants.Then the substitution of Eqs. ( 2) and ( 3) in (1) leads to the following system of nonlinear PDEs for the studied N interacting populations: For the case of 1 population the system (5) is reduced to the equation

B. Spatially Averaged Equations
Below we shall apply spatial averaging to the system of equations ( 5) similar to the averaging used in the optimum theory of turbulence [7] - [8].In the general two-dimensional case let a quantity q(x, y, t) be defined in an large two-dimensional plane area S with acreage | S |.Then by definition the spatial average of q is Then q(x, y, t) can be separated in spatial averaged part q and the rest quantity Q(x, y, t): q(x, y, t) = q(t) + Q(x, y, t).
Let | S | be large enough so that the plane average of any product of the rest quantities vanish: In addition we shall assume that S dx dy∇ 2 Q has finite and small value such that The application of the averaging to Eq. ( 1) in presence of the assumptions given by Eqs. ( 2), ( 3) and (4) (note that in this case we have two spatial dimensions) leads to the system of ODEs as follows (i = 1, 2, . . ., N ): We note that Eq. ( 9) follows directly from Eq. ( 5) in the spatially homogeneous case.For the case of one population the equation becomes We note that the equations of the kind of Eqs.( 9) and ( 10) are often used as model equations in population dynamics not only for small values of population densities but also for large values of these densities, i.e., for large ρ i .

III. TRAVELING WAVES
Let us discuss the simplest case of one population described by Eq. ( 6).First we introduce the travelingwave coordinate ξ = x−vt where v is the velocity of the wave.In addition we shall assume that the polynomial non-linearity in Eq.( 6) is up to order L. We rescale the coefficients in Eq.( 6) as follows: Then Eq.( 6) becomes: Below we shall obtain exact solution of Eq. ( 12) by application of the modified method of simplest equation for obtaining exact solutions of nonlinear PDEs.
A. Application of the Modified Method of Simplest Equation to Eq. ( 12) The method of simplest equation [9]- [11] is based on the fact that after application of appropriate ansatz some NPDEs can be reduced to ODEs of the kind and for some equations of the kind ( 13) particular solutions can be obtained which are finite series constructed by solution φ(ξ) of more simple equation referred to as simplest equation.The simplest equation can be the equation of Bernoulli, equation of Riccati, etc.The substitution of Eq. ( 14) in Eq. ( 13) leads to the polynomial equation where the coefficients σ i , i = 0, 1, . . ., r depend on the parameters of the equation and on the parameters of the solutions.Equating all these coefficients to 0, i.e., by setting one obtains a system of nonlinear algebraic equations.Each solution of this system leads to a solution of kind (14) of Eq. ( 13).
In order to ensure non-trivial solution by the above method we have to ensure that σ r contains at least two terms.To do this we have to balance the highest powers of φ that are obtained from the different terms of the solved equation of kind (13).As a result of this we obtain an additional equation between some of the parameters of the equation and the solution.This equation is called balance equation [12]- [15].
Let us now apply the methodology to Eq. (12).Let us search a solution as finite series where φ(ξ) is a solution of the Riccati equation i.e., The substitution of Eq.( 17) in Eq.( 12) and the balance of the largest powers of φ that arise from the different terms of Eq.( 12) lead to the balance equation Thus we have the possibilities: P = L = 2 or P = 1; L = 3. Below we discuss the first possibility, namely In such a way we shall obtain exact traveling-wave solution of the equation which will be of the kind The parameters of the solutions are determined by the following system of kind ( 16) Now we have 2 possibilities: α † 0 = 0 and α † 0 = 0 (which is more close to the classical population dynamics models that usually do not possess terms independent on the population density).
1) Case α † 0 = 0: For this case the solution of the system (23) is as follows: 2) Case α † 0 = 0: For this case the solution of the system (23) is as follows: The obtained solutions describe kink waves that can be considered as traveling waves of change of the value of the population density of the studied population.Appropriate values of the boundary conditions at ξ = ±∞ can ensure that ρ(ξ) is non-negative everywhere.We note that the parameters of the solved Eq. (21) (after the rescallings ( 11)) are D † and α † 0,1,2 .The first relationship from (24) connects these 4 parameters.Then Eq.( 22) is exact solution of Eq. ( 21) only when the mentioned above 4 parameters are bounded by the corresponding relationship.This means that only 3 of these 4 parameters are free.

IV. STATISTICAL DISTRIBUTIONS AND EXIT TIME
Eqs. ( 9) and ( 10) are typical equations for description of dynamics of dynamical systems.The general case of such equations is For such kind of systems there exists a theory that allows us to characterize some system properties in the case when the system is under the action of random perturbations.Pontryagin, Andronov and Vitt [16] developed such theory for random impulses that occur after every interval of time τ and each impulse causes the phase point of the dynamical system described by Eqs. ( 26) to jump through a distance a along a random direction.Let us first consider the case of single population and one spatial dimension.For the case when a tends to 0 together with τ in such a way that the ratio a/τ tends to finite limit b it is possible to obtain an equation for the probability density function f (x, t) as follows: For the general case of N populations the equation for the probability density function becomes where b ij are again coefficients that characterize the random impulses.
Another kind of problem that can be solved by this approach is to calculate the mathematical expectation of the exit time.Let us again first discuss the case of one population and one spatial dimension.We have a phase point that is inside the interval [ 1 , 2 ] ( 1 < 2 ) and the system is under the influence of the same random perturbations as described above.The exit time is the time for which the phase point that was inside the above interval at t = 0 will leave this interval through 1 or through 2 .If we denote as F (x) the mathematical expectation for the exit time then F (x) is a solution of the equation [16] b 2 with boundary conditions F ( 1 ) = F ( 2 ) = 0.For the case of many populations the zero boundary conditions are on the entire border of the multidimensional phase space area that has to be exited and the equation for the probability density function of the exit time is Here the fixed points are 3.The two maxima of the p.d.f.distribution are centered around the two stable fixed points and the minimum is centered on the unstable fixed point ρ = 0.Here the fixed points are 3 again but one of the two stable fixed points is more preferred which can be seen from the larger peak of p.d.f.function at that point.
Let us now apply this theory to Eqs. ( 9) and (10).For the case of single population let us be interested in the case when after a long time the probability density function f becomes stationary and depends only on the spatial coordinate x.In this case X(x) = L n1=0 α n1 ρ ni .The integration of Eq.( 27) leads to where C 1 is a constant of integration.This constant is equal to 0 if following condition is fulfilled With C 1 = 0 we can continue the integration of Eq.( 7) and the result is where the constant of integration C is determined by the normalization condition Eq. (32) in combination with Eq. (33) mean that L n1=0 α n1 = 0.In addition f (ρ) must tend to 0 when ρ → ±∞.The dominant term at large values of ρ is α L ρ L .Then α L must be negative (to ensure f → 0 at large positive values of ρ and L must be odd (to ensure f → 0 at large negative values of ρ).Several examples of f (ρ) are shown in Fig. 1.
Let us now calculate the exit time expectation on the basis of Eq.(29).We calculate the distribution for exit from the initial position ρ to a position q < ρ.One integration of Eq.(29) leads to the equation dF The solution is searched in presence of two requirements: 1) F (ρ = q) = 0 , 2) F (ρ, q) increases in the slowest possible manner as ρ → ∞ (i.e.dF dρ from Eq. ( 35) is as small as possible).The second requirement leads to C 1 = 0, and the integration of Eq.( 35) leads to the result Fig. 2 shows the dependence of the exit time expectation on the population density and coefficients of the model equation for the case L = 3.The theory can be easily applied for the case of system of many interacting populations but even in the simplest one-dimensional case the integral from Eq.( 36) must be calculated numerically.

V. CONCLUSION
In this paper we have discussed two aspects of population dynamics.First we have presented a model of the space-time dynamics of the system of interacting population in 2 spatial dimensions.For the most simple case of one spatial dimension and for one population we have obtained exact traveling-wave solution of the model nonlinear PDE by means of the recently developed modified method of simplest equation for obtaining exact and approximate solutions of nonlinear PDEs.The obtained exact solution describes the spreading of changes of the population density in the space.The generalization of this theory to the case of many populations is straightforward and describes the spreading of coupled waves of changes of densities of the studied populations.This research will be reported elsewhere.The second discussed aspect of the population dynamics was connected to the influence of the random fluctuations on the population densities.The presence of fluctuations leads to description in terms of probability density functions for the population densities.The discussed general theory is illustrated again for the simplest possible case of one population in two aspects: calculation of probability density functions and calculation of the expected extinction time.The minima and maxima of the obtained probability density functions are exactly at the fixed points of the corresponding non-perturbed model system of differential equations.The expected extinction time strongly depends on the coefficients of the model equations.Because of the lack of space we do not discuss the case of more than one population here.This case will be reported elsewhere.

Fig. 1 .
Fig. 1.Several profiles of f (ρ) from Eq.(33) (in the figures ρ is denoted as ρ). Figure (a): b = 2; α0 = α2 = α3 = 0; α1 = −0.3.In this case the model equation has a single fixed point ρ = 0 and at this fixed point the maximum of the probability density function is centered.Figure (b): α0 = α2 = 0; α1 = 0.3; α3 = −0.05;b = 3.Here the fixed points are 3.The two maxima of the p.d.f.distribution are centered around the two stable fixed points and the minimum is centered on the unstable fixed point ρ = 0. Figure (c): α0 = 0.05; α1 = 0.3; α2 = 0.01; α3 = −0.05;b = 3.Here the fixed points are 3 again but one of the two stable fixed points is more preferred which can be seen from the larger peak of p.d.f.function at that point.
Fig. 1.Several profiles of f (ρ) from Eq.(33) (in the figures ρ is denoted as ρ). Figure (a): b = 2; α0 = α2 = α3 = 0; α1 = −0.3.In this case the model equation has a single fixed point ρ = 0 and at this fixed point the maximum of the probability density function is centered.Figure (b): α0 = α2 = 0; α1 = 0.3; α3 = −0.05;b = 3.Here the fixed points are 3.The two maxima of the p.d.f.distribution are centered around the two stable fixed points and the minimum is centered on the unstable fixed point ρ = 0. Figure (c): α0 = 0.05; α1 = 0.3; α2 = 0.01; α3 = −0.05;b = 3.Here the fixed points are 3 again but one of the two stable fixed points is more preferred which can be seen from the larger peak of p.d.f.function at that point.