THE NORMAL VELOCITY OF THE POPULATION FRONT IN THE “PREDATOR–PREY” MODEL

. The propagation of one and two-dimensional waves of populations are numerically investigated in the framework of the “predator–prey” model with the Arditi - Ginzburg trophic function. The propagation of prey and predator population waves and the propagation of co-existing pop-ulations’ waves are considered. The simulations demonstrate that even in the case of an unstable quasi-equilibrium state of the system, which is established behind the front of a traveling wave, the propagation velocity of the joint population wave is a well-deﬁned function. The calculated average propagation velocity of a cellular non-stationary wave front is determined uniquely for a given set of problem parameters. The estimations of the wave propagation velocity are obtained for both the case of a plane and cellular wave fronts of populations. The structure and velocity of outward propagating circular cellular wave are investigated to clarify the local curvature and scaling eﬀects on the wave dynamics.


Introduction
The predator-prey models are applied to describe nonlinear dynamics of population of many living organisms, in particular, of plankton.Modeling of the phyto-and zooplankton population dynamics is important for estimations of the marine bio stocks since plankton is part of the basic chain of aquatic ecosystems [9,15,18,21,26].Classical models of plankton dynamics [5,7,8,11] describe time evolution of the concentration of phyto-, zooplankton and the nutrients concentrations in the surroundings.
Existing models differ in specification of the processes that determine the development of plankton, consisting of various types of microorganisms, the interactions among microorganisms, the distributions of nutrient concentrations and the influence of physical properties of the environment, such as solar radiation, convective transfer, and others [4].
In the simplest case, the dynamic behavior of the phytoplankton population can be described by the nonlinear Kolmogorov-Petrovsky-Piskunov-Fisher equation (Fisher-KPP equation) [10,17], which is encountered in heat and mass transfer problems, combustion theory, biology and ecology, in plasma physics and in the theory of phase transitions.This model includes only the growth rate, extinction and diffusion transfer of the plankton concentration.One of the most important properties of Fisher-KPP equation is the existence of a solution describing the front of a growing plankton population, which propagates at a constant speed along the normal to the surface of the colony.The population growth autowave propagating through the medium has a characteristic internal structure or "wave thickness" with a scale of l p = D p /r, where D p is the diffusion coefficient of plankton, and r is the characteristic time of the plankton growth.When the characteristic spatial size of the population is much larger than the internal structure of the wave, the population boundary can be considered as a surface (front) propagating relative to the environment with a given speed.In this case, the growth of the colony surface in a stationary or moving medium can be described by a simple kinematic equation for the front propagating with constant normal velocity in a given flow field, similar to the G-equation in combustion theory [14,29].
Since the plankton density inside the colony is usually constant (or can be appropriate determined), data on the temporal evolution of the population boundary allow predicting the growth of total plankton mass with an arbitrary initial shape and to describe plankton evolution in a given flow field of the environment.Note that this approach can be applied only in the case when the size of the population is much larger than the internal structure of the autowave of population.
Population autowaves also appear in other more complex models of plankton growth, for example, in the Nutrient, Phytoplankton, Zooplankton (NPZ) model, which takes into account the changes of the nutrients, the phytoplankton, and the zooplankton concentrations [22].The main goal of the present study is to analyze the dynamics of the autowaves describing by the partial case of the NPZ model, when the nutrients concentration is constant and large enough.The characteristic feature of the model is appearance, in some cases, of the diffusive instability or Turing instability of the plankton population [2,16,19,23,28].The diffusive instability manifests in formation of patching structures ahead of the propagating front and in appearance of the separate cells or "fingers" at the front of population.In this study, we show that at a certain choice of spatial and temporal scales, the propagation dynamics of the ecological communities can be described by using the front concept.The population's front is considered boundary separating areas with a distributed zoo-and phytoplankton populations and the areas where plankton is absent.
A similar approach was used earlier in works on the propagation of sporadic combustion waves consisting of individual flame cells.In a sporadic combustion wave, there is a collective spread of the foci due to the exchange of heat and the reagents among the individual foci.As it is shown in [13] the propagation speed of the sporadic combustion wave front is a well-defined value, despite the complex unsteady internal structure of the sporadic combustion wave.Note that combustion waves and the autowaves of biological population growth belong to the same class of the reaction-diffusion systems, differing only in the type of nonlinear interaction among the reagents of the system.
In this study, using the "predator-prey" model as an example, we introduce the concept of the normal propagation speed of the population front, which depends on the problem parameters such as the diffusion coefficients, the plankton growth rates, the mortality and the predation rates.Using numerical modeling, the features of the populations front propagation in a system with diffusion instability are investigated, the propagation velocities of the population's front with sporadic structure are determined, and the regularities are revealed that make it possible to estimate the propagation velocities for a given set of problem parameters.

"Predator-prey" model
To study the dynamics of populations, a reaction-diffusion "predator-prey" model is used.The two-component "predator-prey" model in the two-dimensional case has the form Here T is the time; X, Y are the spatial coordinates; P (X, Y, T ), Z(X, Y, T ) are the distribution densities of the prey and predator populations; D P , D Z are the corresponding diffusion coefficients and F , G are functions describing local kinetics of the system (2.1).This study focuses on models that admit autowave solutions that describe the propagation of population fronts even in the case of diffusion instability.One of the known models in ecology, where dissipative structures can form, is the Arditi-Ginzburg model [1] with the following functions: Here r is the growth rate of the prey; P 0 is the maximum achievable prey population density in the absence of a predator.In the formula for the G function, the conversion factor σ determines the efficiency of food assimilation, and µ is the mortality rate of a predator (zooplankton).g(P, Z) is the trophic function (functional response) depending on the ratio of prey to predator densities.Function g quantitatively characterizes the dependence of prey absorption by predators per unit of time.Parameter α is the coefficient of predation, and the parameter β is the time spent by the predator in searching for prey and processing food.There are many nonlinear functional responses represented, for example, in [24], where approximating functions consider the saturation of predators, prey searching time, and other factors.Present model allows the detailed description of the population dynamics in such systems as plankton communities (phytoplankton-zooplankton) in comparison with the classical twocomponent models with functional responses of Holling type or models of the Lotka-Volterra type [27].The functional dependence g(P, Z) has a significant drawback associated with the uncertainty of the functional response at zero values of population density.In the present model, the additional assumption is used to eliminate the uncertainty of the g function at the point P = 0, Z = 0: Note that in the more complex Beddington-DeAngelis model (BDA) [3,6], the function g(P, Z) is written in the form In the BDA model, the function g and its partial derivatives with respect to P and Z have no singularities and vanish at the point P = 0, Z = 0.For small values of the parameter h → 0, the BDA model asymptotically transforms into the Arditi-Ginzburg model.In this study, using the Arditi-Ginzburg model, we will assume that the derivatives and the function g(P, Z) itself at the point P = 0, Z = 0 are equal to zero.
In dimensionless variables p = P/P 0 , z = Z/(αβP 0 ), x = X r/D P , y = Y r/D P , t = rT , the system of equations (2.1) has the form Here a, d, m and s are the following dimensionless parameters: The system of equations (2.4) is supplemented by the boundary conditions for the absence of flows at the boundaries of the computational domain and the initial conditions.

Stationary solutions and their linear stability
Stationary, space-homogeneous states of the system (2.4) are described by solutions of the following algebraic equations: The system of algebraic equation (3.1) has three solutions.The first solution (I) describes the absence of population p 1s = 0, z 1s = 0.The second solution (II), p 2s = p * , z 2s = z * , where describes the case of coexistence of both types of species.From (3.2) it follows that the inequality s > m is a necessary condition for the existence of predator in the system: the mortality of a predator should be lower than its birth rate.The third stationary state (III) corresponds to the case when the prey density is maximum p 3s = 1, and there are no predators z 3s = 0.
To study the linear stability of the stationary solutions of (3.1), the perturbed variables are written as follows: All perturbed quantities contain the factor φ = exp(Ωt + ik x x + ik y y) , where Ω is the growth rate of perturbations with a perturbation wavelength k = k 2 x + k 2 y , i is the imaginary unit.Substituting expressions (3.3) into equations (2.4) and linearizing them with respect to perturbations, one can obtain a homogeneous system of equations for p and z.The condition for the existence of a solution can be written in the form of a dispersion equation for the increment Ω Here g = zp/(z + p) , g p = ∂g/∂p and g z = ∂g/∂z are partial derivatives of g determined at the stationary point p = p s , z = z s .At the stationary point p s = 0, z s = 0, where g p = 0 and g z = 0, the solutions of dispersion equation At the stationary pointp s = 1, z s = 0, where g p = 0 and g z = 1, the solutions of dispersion equation (3.4) read Ω 1 = −1 − k 2 and Ω 2 = s − m − dk 2 .The coefficient s related with generation of predator has to be larger than coefficient of mortality m.In this case the Ω 2 is positive in some range of the wave numbers near the point k = 0. Thus, the state (III) is unstable.
At the stationary point p s = p * , z s = z * the derivatives are the following g , where c = m/s.The dispersion equation, in this case, can be written in the form of quadratic equation Here The two roots of equation (3.5) read .
The instability occurs if the real part of increment Ω 1,2 is positive Re(Ω 1,2 ) > 0. It can be realized if A is negative A < 0 or in the case of negative B (B < 0) and any sign of A. The conditions of instability A < 0 can be written as It follows from this condition that the greatest value of the growth rate of disturbances is achieved in the absence of spatial disturbances (k = 0).Disturbances with a wavelength greater than the critical value 2π/k * , which is determined by the condition A(k * ) = 0, are unstable.The conditions of instability B < 0 reads The minimal value of the expression on the right-hand side of the inequality is achieved at k = k cr , where This type of instability is possible only if B(k cr ) < 0. Conditions (3.6) and (3.7) determine the set of parameters a, d, m, a under which the formation of spatial structures in the distribution of concentrations can be observed.
The four typical dependencies of real (Ω 1 (k)) and imaginary parts (Ω 1 (k)) of growth rates on the wave number k evaluated for different parameters a, d, m, s are shown at Figure 1.
Figure 2 shows the regions of different instability modes of the stationary statep s = p * , z s = z * , calculated for fixed values s = 1 and m = 0.6 in the plane (a, d).
In the case a) the concentrations of species remain uniform in the space and equal to stationary values p s = p * , z s = z * .The one-dimensional instability (case b)) can lead to extinction of one or both types of populations.The development of oscillating instability of spatial perturbations in the cases c) and d) can manifest itself as appearance of spots or other nonstationary patterns in the distributions of the both types of species [2].In these nonstationary states, the averaged concentrations in space and time are constant values.Therefore such state can be considered as quasi-stationary states, because the mean population concentrations are constant values and both types of species exist infinitely long time.If the populations in stationary or quasi-stationary states are restricted by a domain located in the surrounding medium, where population distribution is in unstable states of first (I) type, the domain of populations can expand and occupy territory with unstable populations.In the next sections, we consider the propagation of fronts of such populations in the surroundings.

Numerical method
The nonstationary problem of the wave propagation is numerically investigated within the framework of equations (2.4).The system of equations (2.4) is solved by using an explicit finite-difference forward time centered space numerical scheme (FTCS) [25].One dimensional simulations are conducted in the domain 0 < x < L, where L = 1000 with a uniform grid.
Calculations of two-dimensional waves are carried out for a rectangular computational domain 0 < x < L, 0 < y < H. Integration over time is carried out with a step τ = 10 −4 , which ensure the convergence of solutions.The accuracy of the algorithm is checked by comparison of the simulation results on the sequence of the grid refinement.Simulations performed on the grids of mesh sizes δ = 1, 0.1 and 0.01 have demonstrated qualitatively identical behavior in the wide range of problem parameters.Herewith, the average speeds of the wave differ by less than 0.1% for the grids with meshes 0.1 and 0.01.We choose the grid with mesh δ= 0.1 so that at least 10 nodes identify the internal structure of the wave.
The zero boundary conditions for concentrations fluxes are used at the boundary of the domain.At the initial moment, the distributions of both types of species are set constant in two regions of the computational domain, in such a way that the equilibrium or quasi-equilibrium state is located in the left region, and the unstable state on the right.The detailed description of the initial conditions one can find in the next sections.During evolution, a wave of population be formed, which propagated from left to right.
Let J(t i ) is the total number of nodes of the computational grid, in which the concentration of organisms, at least once for the entire calculation time t i , take on a value greater than some small threshold value ε.In the two-dimensional case, the area of such a colony at time t i is equal to J(t i ) δ 2 .If the colony spreads in a rectangular channel, then the ratio of the colony area to the lateral dimension of the channel H will be equal to the average position of the wave front of the population XF (t i ) = J(t i )δ 2 /H.In one-dimensional calculations, the product of the number of nodes and the spatial step of the computational grid δ specifies the coordinate of the front X F (t i ) = J(t i )δ at the time t = t i .
The mean front propagation speed versus time is calculated by the formula where ∆t is the time interval.The choice of the time interval has to be such that the difference between the averages positions of the front XF (t i + ∆t) − XF (t i ) would be greater than the characteristic size of the internal structure of the wave.Noting that in the case of a cellular front, the interval XF (t i + ∆t) − XF (t i ) has to be at least greater than several characteristic cell sizes.
In the case of a diverging circular wave, expression (4.1) is used to calculate the average radial velocity, in which, instead of XF (t i ) the average colony radius RF (t i ) is substituted.
The average colony radius RF (t i ) is determined by the expressions

One-dimensional population wave
Three cases of the traveling waves formation are considered below.

Propagation of prey population front in the absence of predators
In the case when there are no predators in the environment (z = 0), system (2.4) is simplified to one equation for the prey concentration p, which has the form of the KPP-Fisher equation: The KPP-Fisher equation [17] has self-similar solutions p(x − vt) describing the propagation of the population front at a constant velocity v over a medium where prey population is absent.The only solution with minimal possible velocity v = 2 is stable among others with larger velocities.Far behind the wave front, the prey concentration is equal to the maximum possible value p = 1.Within the framework of the KPP-Fisher model, this state with the maximum concentration of prey is stable.In the case of a more complex model (2.4), assuming the possibility of the existence of predators eating preys, the state with a constant prey concentration p = 1 becomes unstable.In this case, model (2.4) has solutions that describe the propagation of a wave of predators over a prey colony with a concentration of p = 1.

Propagation of predator population wave along the prey field
Suppose that at the initial moment the population's distribution is given by step functions of the form Here w ≥ 0 is the width of the area separating the boundary between the prey and predator populations.Distributions (5.2) correspond to the case when in a certain area both types of populations coexist in equilibrium (p * , z * ) (3.2), and outside this area there is only prey population with the constant concentration p = 1.
The numerical solution of system (2.4) with the initial conditions (5.2) showed that after a certain transitional state, a predator wave is formed, which propagates at a constant speed over the area occupied by prey colony.The numerical calculations of the propagation velocity of the predator population front by formula (4.1) is in good agreement with the theoretical formula This formula follows from [17] for the propagation velocity of the autowave front in the model, which is described by the equation If predators propagate over an area with prey concentration p = p 0 , then the function W (z) in the equation for the predators concentration z in system (2.4) has the form From here, using formula (5.5), one can obtain an estimate for the speed (5.3).In the case when the concentration distribution p = p * , z = z * behind the traveling wave is a stable state, the wave propagates at a constant velocity described by formula (5.5).Typical spatial distributions of preys and predators concentrations are shown in Figure 3.
In the case when the concentration distribution p = p * , z = z * behind the traveling wave is unstable, two characteristic modes of wave propagation are realized.In the first case, when the state is unstable and the growth rate of Ω has maximum at k = 0 (see Fig. 1b), the predators population wave is localized in space and it has a bell-like shape.The concentration distributions in the wave for this case are shown in Figure 4.
In the second case when the final state is unstable with respect to disturbances with a nonzero wavelength, as, for example, in Figure 1c and d, the concentration distribution in the wave is shown in Figure 5. Behind the wave front, the population concentrations experience damped fluctuations, which turn into regular periodic fluctuations in the concentrations in the far field.
Note that the average speed of the front in all considered cases tends to a constant value determined by formula (5.3) and the difference between the theoretical and the calculated speeds is less than 0.1% in all cases.

Propagation of coexisting prey and predator populations to an area where both populations are absent
This formulation of the problem corresponds to the case when in a certain area, there is an equilibrium or quasi-equilibrium state between a predator and a prey, and outside this area, both populations are absent.The characteristic feature the populations' fronts propagation in this case is the formation of two population waves -the wave of predators and the wave of prey, moving with different speeds.At the initial moment, the If the width of the region w is very large, then two, practically not interacting, waves of the preys and the predators can be formed.In this case, the preys wave propagates in a space free from preys with a dimensionless speed 2. The predators wave follows the preys wave and propagates through the preys colony with a speed specified by expression (5.3).If the speed of the predators wave is less than the speed of the preys wave then the distance between the waves is constantly increasing and the predators wave lags behind the preys wave.
Figure 6 shows the concentrations profiles in diverging waves at different points in time.With a sufficiently large distance between the fronts, the concentration distributions and the wave structures will be the same as for the single waves considered before.
Calculations have shown that if the speed of the prey wave is greater than the speed of the predator wave, then one wave overtakes the other, and in the future, both waves propagate with the same speed.It occurs when the following condition is fulfilled The speed of the two wave's joint propagation is always less than the speed of prey wave propagation in the absence of predators.
Figure 7 shows the distributions of populations when the two waves of biological organisms propagate together.The results of calculations demonstrating the effect of the coefficients a, d on the speed of two waves joint propagation, at fixed values of the parameters s = 1, m = 0.6, are shown in Figure 8.It is found that at a ≤ 2.1, the front propagation speed does not depend on time, and at a > 2.1, the speed of the wave of joint propagation of populations undergoes fluctuations around a certain average value v, the amplitude of which increases with an increase in the parameter a.
When a ≥ 2.5, the extinction of populations occurs over the entire computational domain.Simulations show that at large diffusion coefficients d > 10, the front propagation speed weakly changes.The hollow squares correspond to the average values of the velocities calculated at a > 2.1 when the wave propagates with oscillations.Figure 9 shows the dependencies of the prey wave front coordinate X r and the speed of its propagation on time in the case of joint propagation of population waves.The dotted line indicates the average speed.The time interval in calculations of instant velocity by formula (4.1) is ∆t = 2.
The time dependencies of instant velocity calculated for different time interval ∆t are shown in Figure 10.Increase in ∆t leads to smoothing of the oscillations and the averaged velocity is irrespective to ∆t.

Propagation of planar two-dimensional wave
At the initial moment, the concentration distributions are given as follows where ψ(y) = θ sin(πny/H) is the small perturbation of the boundaries along y axis; δ is the mean distance between populations and boundary and θ is the amplitude of initial perturbations.
If the wave of the prey population moves faster than the wave of predators, then the waves diverge, as it occurs in the one-dimensional case shown in Figure 8.
During evolution, the initial disturbances of the wave fronts are smoothed out and the waves become flat.The propagation velocity of the prey front tends to v = 2 and the propagation velocity of the predator front tends to that determined by the formula (5.3).At d(s − m) = 1, both waves propagate at the same speed at some distance from each other.
In the case d(s − m) > 1, regardless of the values of the parameters δ and a, the wave of predator populations moves faster than the wave of the prey population and, at some instant of time, both waves begin to propagate together.The waveform in this case can be flat or cellular.
The transition of the perturbed initial distributions into stable plane wave is shown in Figure 11.During evolution, the initial perturbations at the wave front decay and the wave becomes flat.At the same time, cells at the wave front can form even when the state behind the wave front is stable.This instability is not associated with the instability of the state behind the wave front, and the formation of cells occurs due to the difference in the diffusion coefficients of the predator and prey.Figure 12 shows an example of the prey population cellular front formation at successive moments of time when the stationary state behind the front is stable.If the concentration distribution (p * , z * ) in the region behind the wave is spatially unstable, the wave front is always cellular.Figure 13 shows the distributions of concentrations during the joint propagation of cellular waves of a predator and a prey.
In the case Figure 13a, the concentrations distribution behind the front are in the form of stationary spatial structure that can be called as "spotty" state.The labyrinth-like spatial structures are formed behind the travelling wave that is shown in Figure 13b.In the case shown in Figure 13c, the nonstationary patterns are formed with appearing and disappearing spots of concentrations.In the all considered examples, the cells at the wave front are in continuous motion; the smaller ones increase in size to a certain critical value, and then fall apart again into smaller ones.
The chaotic motion of the cells occurs in such a way that the average size of the cells is approximately constant.One can assume that the average size of the cells approximately corresponds to the characteristic  scales of disturbances at which the maximum growth rate, obtained from the linear analysis of stability, is reached.Calculations have shown that the number of local maxima of the front is on average preserved during wave propagation and their number linearly depends on the transverse size of the computational domain H. On the other hand, the constancy of the average cell size at the wave front makes it possible to explain the existence of the average propagation velocity of the cellular wave, which is irrespective to the size of the computational domain and the type of initial perturbations.
The time dependence of the joint wave propagation velocity corresponding to the case shown in Figure 13c is shown in Figure 14.The deceleration of the wave propagation speed at the stage 0 < t < 200 is associated with the transition from the initial distribution to the developed cellular front, which propagates with a constant average speed.
The results of calculations are presented in Figure 15, allowing to compare the values of the average velocities of the fronts in the one-dimensional and two-dimensional cases.An increase in the propagation velocity in the two-dimensional case occurs due to an increase in the total surface of the cellular wave front.The cells have small amplitudes at the a values close to unit, therefore the difference in the average velocities of 1D and 2D waves is small.The difference in the velocities of one-dimensional and two-dimensional waves increases with an increase of the diffusion coefficients d.

Propagation of the circular diverging two-dimensional wave
Radial wave propagation may be of interest for the experimental study of population wave propagation, since this configuration makes it possible to exclude the influence of the external boundaries.In addition, creating in experiments the initial distribution of a colony of organisms in a small circular area seems to be an easier task than creating a uniform distribution of organisms in a rectangular area, which is necessary for studying plane waves of a population.An expanding circular wave with a cellular structure of the front is a convenient object for studying the dynamic behavior of the cells and studying the effect of the cellular front structure on the wave propagation speed.In the case, when the radial increase in the wave front leads to an increase in the size of the formed cells, followed by the growth of secondary cells on their surface, a self-acceleration of the radial wave is possible due to an increase of its surface.
This phenomenon occurs, for example, with the development of instability on the surface of a radially diverging flame front [12].If the size of the cells remains constant on average, then the speed of propagation of the radial wave will asymptotically approach a constant value, which is equal to the propagation speed of the cellular front of the plane wave.This section presents the results of a numerical study of a radially diverging wave of prey and predator joint propagation under conditions when the cellularity of the wave front can form.
The numerical algorithm is implemented in the OpenFoam environment by using the resources of the computing cluster of the IAM FEB RAS with 400 Intel Xeon Gold 6230R (2.1 GHz) processors.
Figure 16 shows the prey concentration distributions at various points in time.
The state (p * ,z * ) behind the wave front corresponds to the instability case shown in Figure 1d.Calculations show that during evolution, no significant change in the size of the cells is observed, and the larger cells are divided into smaller ones at the initial transition stage.In this case, the total surface of the flame front increases linearly with the increasing radius, like in the case of a uniformly expanding circular front.Figure 17 shows the dependence of the average front radius on time, which is calculated with a time step ∆t = 2.The wave radius grows linearly with time after a transitional stage in the initial moment of the wave propagation.
Figure 18 shows the dependence of the front speed on the radius, which is calculated by formula (4.2) with a time step ∆t = 4.
As it follows from Figure 18, the speed tends to a constant value at large radiuses.The value of the average radial velocity at large radiuses is close in magnitude to the average propagation velocity of the plane wave with cellular front calculated for the same values of the parameters.At small values of the colony radius, the velocity depends on the local curvature of the wave front.This dependency appears by taking into account the finite thickness of the internal wave structure and it can be useful in evaluating of a colony growth having any  where v is the mean propagation velocity of the flat wave with cellular front and χ is the nondimensional length characterizing the scale of inner wave structure.Note that similar approach is applied in the estimations of the flame front curvature effects in combustion theory [20].The dependency (7.2) is plotted by dashed line in Figure 18 for χ = 11.24.The value of χ that was roughly estimated by matching dependency (7.2) to the calculated data.

Conclusion
The structure and propagation velocity of one and two-dimensional population's waves are numerically investigated in the framework of the "predator-prey" model with the Arditi-Ginzburg trophic function.The propagation of population's wave of living organisms occurs due to the growth of an already formed colony in an area where living organisms are absent, and there are conditions for their reproduction.Inside the colony, the distribution of the living organism's concentrations corresponds to a stationary or quasi-stationary state of the system.In the case of propagation of only a prey population in the absence of predators or distribution of only predators over an area with a constant concentration of prey, the equations describing wave propagation are similar to the Kolmogorov-Petrovsky-Piskunov-Fisher equation.In these cases, the wave propagation speed is the same as in these classical models.The dynamics of wave propagation becomes more complex when the populations of the prey and predator behind the wave front are in an unstable, quasi-equilibrium state.In this case, changes in the concentration of prey and predator in the wave are interrelated, and the two-dimensional wave sometimes assumes cellular shape with unstable chaotically moving cells.The region of the problem parameters corresponding a quasi-equilibrium state behind the wave front can be determined by the linear stability analysis of the homogeneous in space concentrations distribution.Calculations have shown that the spatial instability of the wave front, at which cells are formed, is observed in a wider range of parameters, compared with the range of parameters that determine the instability of a state homogeneous in space.The propagation of a joint wave of a predator and a prey can be non-monotonic, with a speed periodically varying in time.It is shown that the average propagation speed of a combined wave is a well definite value, even in the case of the formation of cellular wave front.The simulation results showed that the radial growth rate of the diverging circular population at large colony radii is equal to the propagation speed of plane populations' wave.This allows applying the concept of normal propagation velocity to estimate the mass growth of organisms in a large colony with an arbitrary initial shape.In conclusion, we note that the data on the normal velocity and the structure of the population's wave can be conveniently used to verify and to refine the model parameters by comparing the theoretical values with experimental data.

Figure 1 .
Figure 1.The dependencies of real (Ω 1 (k)) (solid lines) and imaginary (Ω 1 (k)) (dashed lines) parts of growth rates on the wave number k obtained from linear stability analylis of stationary state p s = p * , z s = z * and evaluated for s = 1, m = 0.6 and different parameters a and d.(a) Stable state, d = 1, a = 1.(b) Instability with oscillations of long wave perturbations 0 < k < k 1 , d = 2, a = 2. (c) Instability with respect of spatial perturbations without oscillations in the diapason of wave numbers k 2 < k < k 3 , d = 10, a = 1.9.(d) Instability without oscillations in the diapason of wave numbers k 5 < k < k 6 accompanied by oscillating instability of long wave perturbations 0 < k < k 4 , d = 5, a = 2.1.

Figure 2 .
Figure 2. The regions of different instability modes of the stationary state p s = p * , z s = z * , calculated for the fixed values s = 1 and m = 0.6 in the plane (a, d).The red dotted line corresponds to the stability boundary calculated by the equation B = 0, provided k 2 cr > 0. The vertical line is the border of neutral stability k = 0 calculated by the equation A = 0.

Figure 3 .
Figure 3. Distribution of the predators (red curve) and the preys (blue curve) concentrations in the traveling wave, calculated for values corresponding to a steady state p = p * , z = z * far behind the wave front.s = 2, m = 0.6, a = 1, d = 2.The dotted lines show the initial distributions.

Figure 4 .
Figure 4. Distribution of the predators (red curve) and the preys (blue curve) concentrations in the traveling wave evaluated for s = 1, m = 0.6, a = 2.2, d = 1.The dotted lines show the initial distributions.

Figure 5 .
Figure 5. Distribution of the predators (red curve) and the preys (blue curve) concentrations in a traveling wave, calculated for the values s = 1, m = 0.6, a = 1.9, d = 10.

Figure 6 .
Figure 6.Distribution of the predator (red curve) and the prey (blue curve) concentrations in diverging traveling waves evaluated for s = 1, m = 0.6, a = 1, d = 1 at different times.The dotted lines are the initial distributions.

Figure 7 .
Figure 7. Distribution of the predators (red curve) and the preys (blue curve) concentrations in joint propagation of two waves evaluated for a = 1.5, d = 5, s = 1, m = 0.6.The dotted lines are the initial distributions.

Figure 8 .
Figure 8.The dependence of the front propagation speed v on the parameter a evaluated for d = 3 (green curve), d = 5 (blue curve) and d = 10 (red curve).Hollow squares correspond to the average values of the velocities calculated at a > 2.1.

Figure 9 .
Figure 9.The dependencies of the prey wave front coordinate X F (left figure) and the speed of its propagation v (right figure) depending on time in the case of joint propagation of population waves evaluated at d = 3, a = 2.3, m = 0.6, s = 1, ∆t = 2.The dotted line indicates the average speed v = 0.424.

Figure 14 .
Figure 14.The propagation velocity of the joint populations wave evaluated for a = 2.1, d = 5, s = 1, m = 0.6, ∆t = 10 The concentration distributions corresponding to this case is shown at Figure 13c.

Figure 15 .
Figure 15.The average velocities of the joint populations wave as function of the a parameter for d = 3 and d = 10 in one-dimensional (1D) and two-dimensional (2D) cases, calculated for s = 1, m = 0.6.Hollow squares denote the average velocity of the oscillating front in the one-dimensional case.

Figure 16 .
Figure 16.The dynamics of the prey population distribution propagating from the initial circle with the radius r 0 = 50, evaluated for parameters a = 2.1, d = 5, s = 1, m = 0.6.