Waves in a Hyperbolic Predator–Prey System

: We address a hyperbolic predator–prey model, which we formulate with the use of the Cattaneo model for chemosensitive movement. We put a special focus on the case when the Cattaneo equation for the ﬂux of species takes the form of conservation law—that is, we assume a special relation between the diffusivity and sensitivity coefﬁcients. Regarding this relation, there are pieces arguing for its relevance to some real-life populations, e.g., the copepods ( Harpacticoida ), in the biological literature (see the reference list). Thanks to the mentioned conservatism, we get exact solutions describing the travelling shock waves in some limited cases. Next, we employ the numerical analysis for continuing these waves to a wider parametric domain. As a result, we discover smooth solitary waves, which turn out to be quite sustainable with small and moderate initial perturbations. Nevertheless, the perturbations cause shedding of the predators from the main core of the wave, which can be treated as a settling mechanism. Besides, the localized perturbations make waves, colliding with the main core and demonstrating peculiar quasi-soliton phenomena sometimes resembling the leapfrog playing. An interesting side result is the onset of the migration waves due to the explosion of overpopulated cores.


Introduction
The Patlak-Keller-Segel (PKS) law provides a simple macroscopic model for a perceptual motion (taxis) of the particles ensemble in response to the spatial gradient of a stimulus (or signal) field. The commonly recognized formulation of the PKS flux reads as: χp∇s, where the notations of p, s, and χ stand for the density of the medium that moves in response to the signal, the intensity of the signal itself and the sensitivity coefficient, correspondingly. A multitude of parabolic PKS systems resulting from the summation of the diffusive and PKS fluxes stand as the subject of intensive research for last five decades. A considerable number of reviews expose the advances in this area, e.g., [1][2][3].
However, usage of the parabolic framework is not a unique way of treating the PKS law. For example, Dolak and Hillen [4] proposed a different formulation known as the Cattaneo model for chemosensitive movement. In contrast with the parabolic models, this one takes into account the inertia of the response and becomes hyperbolic. At that, the flux has a contribution from the local time derivative of itself in addition to those mentioned above. From the articles [5,6], it follows that the common parabolic model, the Cattaneomodel and several other hyperbolic models represent the approximations of the kinetics equations under different hypotheses. There is a concise review by R. Eftimie [7] of this subject that covers deriving the models and the issues of exact solutions, stability and bifurcations, etc. It makes clear that the hyperbolic models are not less natural than their parabolic counterparts, e.g., while modelling the aggregation processes in the active media. Nevertheless, the former receive much less attention than the latter. The mentioned

The Governing Equations and Scaling
The governing equations read as: s t = G(p, s) + D s s xx . (3) In this system, the dependent variables, x and t, stand for the spatial and temporal coordinates, correspondingly, x ∈ R, t > 0. The dependent variables are p = p(x, t), q = q(x, t) and s = s(x, t). The first and the last one play the parts of the densities of the species, say, the predators and the prey, correspondingly. The first two equations constitute the Cattaneo model for the prey-sensitive movement of the predators, so that the remaining dependent variable, q, stands for the predators' flux. The prey spreads itself purely by diffusion, and the notation of D s stands for its diffusivity. In what follows, D s ≡ const by assumption. We also assume that the predators' diffusivity, µ(p, s), and the sensitivity, χ(p, s), are specified, and pχ(p, s) → 0, p → +0, ∀ s ≥ 0.
We assume that the reaction terms, G and F, are prescribed, but postpone further detailing them to Section 5. The mechanical analogy suggests treating the second term on the left hand side of the second equation as a contribution of the resistance of the environment to the predator's motion. So, we will be calling the correspondent coefficient, ν, as the resistivity.
Let the notations T, X, P, S, Q, D p , C, J p , J s stand for characteristic scales for the time, length, predator's density, prey's density, diffusivity, sensitivity, predator's and prey's sources densities correspondingly. The resistivity coefficient, ν, is naturally dimensionless.
Since the values of D p , D s and C depend on the concrete species, it is natural to consider them as anyhow given. In contrast, the values of X, T, S, P, J p and J s are free to choose. Given this, let us set and postpone defining the values of P and S. We also set µ(p,s) = D −1 p µ(Pp, Ss),χ(p,s) = C −1 χ(Pp, Ss).
In what follows, every variable employed is dimensionless by default. Upon the above scaling, the dimensionless form of the governing equations reads: For the forthcoming analysis, it is important to distinguish the case when the righthand side in the flux equation is integrable in the sense that: For such an integrability, it is necessary and sufficient to link the diffusivity to the sensitivity as follows:

Sedentary Prey and Hyperbolicity
Throughout this section, we consider the system (7)-(9) with δ = 0. Then it becomes a first order quasi-linear PDE system, which we put into the form: The eigenvalues of matrix A are ± µ(p, s), 0. They are real and distinct one from another as long as µ(p, s) > 0. It is always true by assumption. Hence the system of Equations (13) and (14) is strictly hyperbolic. The above triple of eigenvalues determines the triple of characteristic speeds-that is, every characteristic of system (13) allows a parametrization by the mapping t → (t, X(t)) that satisfies an equationẊ(t) = λ(p(X(t), t), s(X(t), t)), where λ ∈ {0, ± µ(p, s)}.
A question to ask about a first-order hyperbolic system is whether or not it allows diagonalizing by a pointwise transform ρ = R(z). Such an ansatz generally does not exist for a system that includes more than two equations. The diagonalizing is feasible, provided that the system matrix, A = A(z), satisfies an integrability criterion, which allows a straight algebraic formulation. Another formulation requires satisfying the Frobenius condition dω i ∧ ω i = 0 with every 1-form ω i = i · dz generated by the vectors of the dual eigenbasis { 1 , 2 , ...} of the matrix A (this is the eigenbasis of the transposed matrix). Then for every i there exists a factor α i = α i (z) such that α i ω i = dρ i , and ρ i = R i . The last criterion is handy to use as long as there is a handy dual eigenbasis, as in the case of matrix (14), for example. Then a routine calculation reduces the diagonalization criterion to the following condition: Under condition (12) the obtained criterion simplifies to: where c stands for an arbitrary function. Thus, assuming the diagonalization entails restrictions that are too artificial even in the simplified form. That is why we will not be considering this option in this study anymore. Let the condition (12) hold throughout all subsequent considerations. Then there exists a single-valued function ϕ = ϕ(p, s) such that: Hence, the system (13) and (14) consists of conservation laws, namely: This feature makes feasible the generalized solutions, e.g., the shock waves. Consider now a shock wave with discontinuities at some curve x = X(t). Then the velocity of the shock,Ẋ, has to satisfy the Rankine-Hugoniot conditions entailed by the conservations laws (18). They read as: Here the notation of a dependent variable put in the square brackets stands for the jump of this variable across the discontinuity-that is, the difference between the limit values evaluated for (x, t) → (X(t) + 0, t) and (x, t) → (X(t) − 0, t). At that, where the superscripts + or − stand for the unilateral limit values at the left and right shores of the discontinuity (when the observer looks forward alongside the time axis).
There are two possibilities. The first is the standing wave-that is, where the value of [s] remains undetermined. The second is the travelling wave-that is, where s * = s(X(t), t) is the value of the prey density directly on the shock (note that function s remains continuous across the shock), and p * ∈ (p − , p + ) is an unknown quantity that for each t satisfies the equation: The speed at which such a wave propagates can take two opposite values, namely: Generically, the shock wave speeds determined in (24) differ from the characteristic ones, which are discontinuous across the shocks. It depends on the inequalities between the velocity of a specific shock wave and the limit values of the characteristic velocities whether or not this shock propagates. These inequalities are known as Lax conditions. Checking them yields a conclusion that shock waves can propagate provided that: Moreover, exactly one of the two possible shock waves can propagate, and the propagation speed is equal to µ(p * , A peculiar degeneration occurs when a predator's diffusivity is independent of its density-that is, µ = µ(s). Given this, the expressions for the shock wave speeds read as either:Ẋ = ± µ(s * ) (26) orẊ = 0. It is worth recalling that s * stands for the trace of the prey's density, s, right on the shock. It is defined well due to the continuity of the function s across the shock. Hence, every possible wave's speed coincides with a characteristic one-that is, the shocks spread along the characteristics.
Matching the predator-independent diffusivity to the integrability condition, (12) and assumption (4) make the sensitivity, χ, predator-independent too. Given this, the integrability condition simplifies as follows: Modulo the scaling, relation (27) is equivalent to that deduced by Tyutyunov et al. from biological rationales in the aforementioned articles. Also, Tyutyunov et al. addressed some issues of stability and pattern formation for the corresponding PKS systems in the parabolic form [22]. Throughout the rest of the present article, we will be studying Cattaneo's systems that arise from relation (27).

The Travelling Shock Waves for the Predator-Independent Diffusivity
Let d/dt stand for the total derivative along a characteristic. From system (18) it follows that: where the notation of λ = ± µ(s) stands for the corresponding characteristic speed. If this characteristic supports a shock, then the same equations hold on the shores of the shock with the limit values of every quantity involved. Subtracting them and eliminating the variable q with the use of the Rankine-Hugoniot conditions (19) lead to the following equation on the shock where subscript * indicates the quantities, which depend on the values of s * only. Furthermore, from the continuity of the prey's density, s, it follows that the gap of ∇s across the shock is normal to it everywhere-that is, So, we are discussing the travelling shock waves. By definition, such a wave propagates at a constant speed c, and the corresponding solution depends on only one variable, Then the characteristics at which the shocks occur are the parallel lines determined by equations ξ = ξ * for the values of ξ * varying over some set, Σ, which we assume to be finite. This set gathers all the discontinuities of the solution in variable ξ. On the complement of the singular set, the system (18) reduces to an ODE system. Namely, The conditions for matching the solutions are defined on the adjacent intervals separated by a discontinuity following from Equations (29) and (30). The latter holds trivially for every function in only variable ξ, while the former turns to a system of functional equations as follows: The last equation in this set is equivalent to [q] = c[p]. The following constraint is for the wave speed, c, to obey: Alluding to the values of the unknown, s, taken on Σ is correct by the continuity that is consistent with the matching condition (33). The same is true regarding another unknown, r.
The first equation in system (32) degenerates when the independent variable, ξ, is approaching a discontinuity. Let us restrict ourselves within the solutions obeying the following condition: where the superscripts, + and − , are to distinguish the solutions settled on the right and left intervals adjacent to the discontinuity point, ξ * . This estimate entails one more condition for matching the solutions at the discontinuity. Namely, the following equations are to obey: where the notations of p ± * , s * and r * stand for the unilateral limits of variable p and the bilateral limits of variables s and r correspondingly. It is worth noting that subtracting the equalities (36) gives the first of matching condition (33).
At first, let us address the waves which allow a single shock only. Then there must be only one discontinuity on the ξ-axis, so let us place the origin there, and put ξ * = 0. Consider the equation: By conditions (36), a nonzero jump across the singularity in variable p due to a travelling wave presumes that Equation (37) has several distinct solutions (x i , y, z) with the same y, z, e.g., z = s * , y = r * , x i = p + * and x j = p − * for some i = j. Assume there exists a nonempty set filled with solutions (x, y, z) to Equation (37) such that x ≥ 0, z ≥ 0, and let the projection of this set onto the yz-plane alongside the x-axis cover some domain with multiplicity N ≥ 2. For every point (z, y) in this domain, let P i = P i (y, z), i = 1, 2, . . . , N be the coordinates of projection of its pre-image on the x-axis alongside the yz-plane. Then constructing the travelling shock waves can go the following way: set and then try to find the solutions p ± = p ± (ξ), r ± = r + (ξ), s ± = s ± (ξ) to Equation (32), that are defined and bounded for ±ξ > 0 and match the data listed above when ξ = 0. Note that the values of y, z are free to manipulate. One can try to obtain some pairs of conjugated waves by transposing i and j. Since we have been assuming the projection to be the N-leaf covering mapping, there are N(N − 1)/2 conjugated pairs of the datasets.
In the next subsection, we will be considering a case of a very simple implementation of the approach outlined above.

The Inertial Limit for the Sedentary Prey
At this point, we need more details regarding the reaction terms, F = F(p, s) and G = G(p, s). Henceforth, we will be assuming the following: The last two assumptions are equivalent to inequalities g p | p=0 < 0, g s | s=0 < 0 modulo the scaling of variables p, s since the characteristic scales of these variables, P and S, are indefinite still. For example, the Lotka-Volterra kinetics while being normalized in accord with (39) reads as Thus, parameter α plays the part of the carrying capacity given the scaling adopted here.
Consider the Cauchy problem Assume there exists a number a 0 > 0 such that for every a ∈ [0, a 0 ] problem (41) has a solution S = S(τ, a) defined on R, and such that the mapping a → S(τ, a) sends the segment [0, a 0 ] to itself for every τ ∈ R. Consider also the following system of functional equations: and assume there exists a positive solution s = s e > 0, p = p e > 0. The last assumption entails the existence of a strictly positive equilibrium for general system (7)-(9). By equilibrium, we mean a particular solution specified as follows: q = 0, p = const, s = const, so that both species distribute themselves homogeneously. A nonnegative equilibrium not being strictly positive can make sense too. For example, the Lotka-Volterra kinetics (40) allows an equilibrium with p = 0, s = α for every parameter setting, and there exists the strictly positive equilibrium s e = β/γ, p e = α − s e provided that β < αγ.
Throughout the rest of this section, we will be considering the inertial limit and sedentary prey. So, we put ν = δ = 0. We also assume all the listed assertions on the kinetics to be true. The Lotka-Volterra kinetics defined in (40) meet such an assumption for β < αγ; at that, the ODE involved in the Cauchy problem (41) reads as S τ = S(α − S).
Given the assumption made, there are at least two solutions, x = P i (y, z), i = 1, 2 to Equation (37) defined in some vicinity of every point (y, s e ). At that, P 1 ≡ 0, while function P 2 (in fact, in one variable, z) is that defined implicitly by equations µ(z) f (x, z) = κχ(z)g(x, z) = 0, P i (s e ) = p e . However, we do not need much manipulating with the values y, z here, since the choice is evident; namely, x 1 = 0, x 2 = p e , z = s e , and the values of y are arbitrary. Then there are two pairs of the travelling shock waves, one pair propagates at speed c = µ(s e ), and the other one propagates at the opposite speed. The formulae for both are the same, and they read as: p = 0, r = y, s = S(ξ, s e ), ξ ∈ (−∞, 0), p = p e , r = y, s = s e , ξ ∈ (0, ∞), (43) p = 0, r = y, s = S(ξ, s e ), ξ ∈ (0, ∞), p = p e , r = y, s = s e , ξ ∈ (−∞, 0), (44) where y is an arbitrary constant. Waves (43) and (44) are conjugated in the sense explained in the previous section but not mirroring one by the other one. Although both species in both waves take the same equilibrium values within the areas settled by predators, outside them, prey spreads itself differently. Indeed, the solutions to the Cauchy problem (41) generically behave differently for ±τ > 0 (except for the equilibria).
In Figure 1, the left and middle frames show profiles of waves (43) and (44) correspondingly. For both waves, the front of the predator's invasion moves towards the smaller concentration of prey. This feature is not as paradoxical as it can seem given the locality of the predator-prey interactions.
An overlay of waves (44) and (43) gives examples of finite predators' mass localized within a patch that moves uniformly. Such a wave propagates with two shocks at the speed c = µ(s e ) or at the opposite speed. The formulae read as: where the values of y and M > 0 are arbitrary constants. At that, the values of M stand as the mass of the patch. It is worth noting that the shapes of the patchy waves differ depending on the sign of the wave speed, c. The prey are dying out to the left (right) of the predators patch for negative (positive) c, so that the patch moves towards the smaller concentration of prey (see the right panel in Figure 1). In Section 7, we return to these waves to discuss their relevance to the population dynamics. . Each frame shows three instantaneous profiles for both the predator and the prey densities. In the left and middle frames, the blue (green) coloured lines are for the former (latter), while the solid, dashed and dotted lines picture the profiles taken at t = 0, t = 3/2 and t = 3. The right frame addresses the wave that transports a patch filled with the unit mass of predators. At that, the dashed (solid) lines mark out the densities of the predators (prey). The colours green, red and blue are for the shots taken at t = 0, t = 3/2 and t = 3 (so that the wave speed is negative). All three panels correspond to the Lotka-Volterra kinetics (40), where α = 1, β = 0.2, γ = 0.8. At that, the diffusivity function reads as µ = (1 + s) −1 , so that κ = 1, and the equality (27) determines the sensitivity, χ. These definitions give the wave speed, c = µ(s e ) ≈ 0.9. Finally, recall that the figure regards the diffusionless and inertial limit, and δ = ν = 0, hence.

Numerical Experiments
In this section, we present the numerical solutions to system (7)-(9) that we formulate with the use of the Lotka-Volterra kinetics (40) and predators' diffusivity µ = (1 + κs) −1 . The diffusivity determines the sensitivity, χ(s), by the simplified integrability condition (27). For the numerical implementation, we eliminate the flux, q, from this system with the use of ansatz by K. Hadeler (see references provided above). As a result, we arrive at the following equations: Substituting this second-order system for the original one (which is of order 1) requires ad hoc initial conditions. These are: where the notations of p 0 , q 0 , s 0 stand for the functions which determine the initial values of the dependent variables of the original system (7)- (9). The numerical implementation also requires restricting the solution within a finite spatial domain and formulating suitable boundary conditions. So, we set It turns out that the numerical solving of the initial-boundary problems (46)-(49) is quite feasible with the Maple built-in PDE solver.
Implementation of the numerical solution should cover a spatial interval wide enough to put the artificial boundaries far out of the domain in which the phenomena of interest occurs. We had controlled all the results below by widening the spatial area and found them reproducing themselves sustainably for values of L higher than 8. In particular, L = 10 for all the figures presented in this section. In a similar way, we checked the influence of the mesh refining and varying the level of smoothing used for preparing the initial data and found the results of this inspection quite satisfactory.
The first set of numerical experiments is for answering the question of whether the shock waves persist for the positive values of the resistivity, ν and the prey diffusivity, δ. To this end, we have been taking the profiles of the species' densities, s and p from the waves (43)-(45) and putting them as the initial profiles, s 0 and p 0 , which enter the boundary conditions (48). At that, we have been smoothing the shocks slightly. Further, we have been putting q 0 = c * p 0 , where c * = µ(s e ) is the shock wave speed. This choice is consistent with the definition of the shock waves. Appendix A provides the concrete details of setting initial conditions, the control parameters values, etc. The answer to the question formulated above is fully affirmative. The shocks become a bit smoother, but keep propagating at an almost constant speed that is nearly equal to the value of c * . Figure 2 shows the typical behavior of the slightly smoothed counterpart of the patchy shock wave. This figure tells us that the smoothed patch spreads as a kind of soliton, which is shaped rather sharply for the small resistivity and prey diffusivity. An increase in the resistivity produces scattering of the predators behind the rear front of the wave.
The next piece of computing addresses the interactions of the observed solitary waves with some perturbations applied initially. These are: (a) a small displacement of the species density profiles one relative to the other one with no deformation; (b) a small deformation of the species density profiles; (c) a small droplet of predators localized behind the main core; (d) a small droplet of predators localized ahead of the main core.
In cases (a) and (b), the smallness means that the magnitudes of mutual displacements (deformations) are approximately ten percent of the magnitude of the main patch. The smallness of the droplet means that its mass is ten times smaller than the mass of the main patch, while they both are localized in the intervals of nearly equal lengths. In Appendix A, there is the exact formulation of the initial conditions and the concrete values of the control parameters.
In cases (a) and (b), the effects of the initial perturbations manifest themselves mostly by the predators scattering, which goes almost the same way as shown in the bottom row of frames in Figure 2, with no any qualitative distinctions. So, we do not illustrate these cases. Case (c) is similar to the above, but shedding the predators is rather intensive. We illustrate this case in the top row of frames in Figure 3. Case (d) demonstrates a very peculiar interplay of solitary waves. Namely, the main patch and the droplet attract one to the other one until they clash. Then they play leapfrog: droplet climbs over the patch and rolls down to the other side of it. The patch drops some mass due to the scattering but keeps moving at almost the same speed. The droplet keeps moving too but in the opposite direction while getting bigger and sharper. We illustrate this case in the bottom row of frames in Figure 3. The results presented above confirm the stability of the travelling patchy waves. Figure 4 demonstrates the results of a more extreme crash-test. Appendix A provides the detailed description of the initial states and other settings used for this piece of computing. It is worth recalling that every travelling patchy wave has a counterpart that propagates at the opposite speed. Gluing this pair provides the initial state for computing the next five patterns pictured in Figure 4. Overall, they illustrate the collision of two patchy waves. The remarkable feature is that the collision gives rise to an expansion wave, at the fronts of which peaks are arising and sharpening in the course of its propagation. The mirror symmetry of the patterns is due to the mirror symmetry of the problem and initial data.  Further, we proceed with an asymmetrical initial configuration that arises from perturbing the travelling patchy wave in a way similar to that used in case (d) above. However, the perturbation is not small this time as its mass is equal to the patch mass. Appendix A provides a detailed description of the initial states and other settings used for this piece of computing. In Figure 5, the two rows of images display two ways of the evolution of this configuration for two different values of the resistivity, ν. The top row regards the smaller resistivity. The initial stage of evolution is qualitatively like that reported for case (d) above. New features arise after the leapfrog playing. Among them, the most notable is the spike that springs up suddenly at the foremost bound of the main patch. The other core becomes sharper too. The bottom row shows the changes due to increasing the resistivity to a considerable extent. It is easy to see a powerful smoothing, which is emergent from forcing out the waves by the equilibrium state. The rightmost frame indeed allows us to see that the values of the densities are close to the equilibrium ones near the origin. Here it is worth recalling that, for the kinetics (40), the equilibrium densities are p e = α − s e , s e = β/γ. So, p e ≈ 0.43, s e ≈ 0.57 for the parameters values adopted for computing the patterns presented herein.
The numerical experiments illustrated above are about the travelling shock waves presented in Section 5. The system we deal with, however, hides a multitude of interesting features, one of which is the formation of peaks after colliding the travelling patchy wave with another pattern, which is not necessarily a wave of the same type too. In Figure 6, we illustrate the occurence of a similar feature irrespective of the patchy waves. Namely, we consider an explosion wave due to a unit mass of the predators smoothly localized in a compact area at the initial time moment with the zero initial flux, q 0 . At that, the initial density of prey is everywhere equal to the carrying capacity, which, effectively, takes the value of the parameter α given the scaling applied here. At the same time, it is the value of carrying capacity that corresponds to the prey density at the equilibrium state with no predators. Then the initial core of the predators stands as a perturbation, which is not small though localized. The bottom (top) row of frames corresponds to the parameters values such that the equilibrium state with a positive predators' density is possible (impossible). The mirror symmetry of the patterns is due to the mirror symmetry of the problem and initial data, the detailed description of which is in Appendix A. Both rows demonstrate an explosion of the initial core accompanied by spreading the predators across a widening areal extent with peaks at the bounds. Outside these boundary layers, both densities tend to the equilibrium values. If the positive equilibrium density of the predators is not feasible then smoothing and decaying of the boundary peaks takes place by degrees. Otherwise, these peaks become sharper and higher.

Discussion
We have addressed a Cattaneo-type dynamics of the predator-prey system with the Lotka-Volterra kinetics term in one spatial dimension. By assumption, only the predators are capable of the perceptual motions, and the flux of them generally reads as pχ(p, s)s x − µ(p, s)p x , where µ(p, s) and χ(p, s) are some prescribed functions. To start with, we have considered the sedentary prey. In this approximation, the governing equations turn to form a strictly hyperbolic system, for which we have formulated the criterion for reducing to Riemann's invariants in terms of sensitivity, χ(p, s), and diffusivity, µ(p, s), explicitly. It has turned out, however, that the class of systems obeying this criterion looks rather artificial. At the same time, reducing the governing equations to the conservation laws happened to be less restrictive. In particular, such a reduction is possible provided that µ = µ(s), κχ(s) = −dµ/ds . Since there are biological rationales for this structural relation (see the reference above), we have accepted it.
For the systems of conservation laws, the shock waves are natural, and we have derived a system of conditions on the shocks. In the inertial limit-that is, for ν = 0-we have discovered a very simple exact solutions that describe the travelling shock waves. The wave pattern represents a semi-infinite or even finite patch of predators that propagates at a constant speed. There are no predators outside the patch while both species coexist in the equilibrium state inside, see Figure 1. The wave speed is equal to ± µ(s e ), where s e is the prey density at the equilibrium state.
The waves carrying semi-infinite patches describes the transitions between two equilibria states. The community goes either from the extinction of both species to the coexistence at the equilibrium or from the coexistence at the equilibrium to the extinction of the predators and restoring the prey up to the carrying capacity. In this sense they resemble the KPP-Fisher waves. Here, however, the extinction of predators means that the wave either has not brought them to the areal extent under consideration yet or has moved them out already. The interesting feature is that the front of the predators' invasion always moves towards the smaller concentration of prey. This is not as paradoxical as it can seem given the locality of the predator-prey interactions. Thus, two patterns of behavior are emergent from propagating these waves. The first is restoring the resource up to the carrying capacity after the migrating predators' withdrawal. The second is the transition from the prey dying out to the mutual equilibrium with the spreading predators. The travelling shock waves carrying a finite mass of the predators combine both features. Indeed, the profile of such a wave involves the transitions from the mutual extinction via the coexistence to the extinction of predators. Besides, the patch of predators moves towards the smaller concentration of prey again. So, the corresponding pattern of behavior looks like preventing the prey from dying out and even restoring the prey population due to migrating a compact core of the predators.
We have extended the travelling shock waves to the positive values of the prey diffusivity, δ, and the resistivity, ν, numerically. It turns out that they withstand such an extension, at least while the values of δ and ν remain small enough. The shocks smooth themselves, but the speed at which they propagate is nearly equal to the value of c = ± µ(s e ).
The shock waves carrying a finite mass of predators transform themselves into the smooth soliton-like waves, to which we have paid particular attention, see Figure 2. An interesting feature is shedding the predators behind the rear front of the wave due to an increase in the resistivity. As a result, the core of migration leaves behind itself a populated areal extent, which remains settled even when the migration core moves far ahead. Thus, propagating the patchy wave can play a part in settling the predators due to migration.
A relatively small predator's droplet that occurs instantly ahead of or behind the main core of the soliton-like wave enhances the above scattering but does not cause any other noticeable changes (see Figure 3). The main core never absorbs the droplet in the case of collision, but they both keep moving after a peculiar interaction resembling the leapfrog play. Quite a different pattern of behavior arises from colliding the cores, which have gathered the equal masses of the predators (see . We have examined colliding for several pairs of cores, which belong to the following classes: (i) two identical smoothed travelling patchy waves which propagate at opposite speeds; (ii) a compact group of predators that occurs suddenly ahead of the predators' patch of a travelling wave and get moving towards it; (iii) a finite mass of predators that suddenly have landed on a compact part of a greater areal extent where the resource density has attained the value of the carrying capacity. Merging the cores immediately leads to local overpopulation and a lack of prey, which, in turn, gives rise to an explosive migration. The explosion wave spreads the predators uniformly across a widening area while boundary layers of high density occur near the wave fronts. Deep in this area, the species coexist at equilibrium if feasible, or the predators become extinct, and the prey density goes back to the carrying capacity. Various numerical experiments have been reproducing this pattern of behavior sustainably though not literally, of course. Thus, we conclude that the explosive waves that we have been discussing deliver a mechanism for overcoming the local overpopulation and lack of resources. setting the initial conditions for the concrete computations. We also recall that, throughout all the computations, the predators' diffusivity reads as: At that, the relation (27) determines the predators sensitivity, χ.