Cosymmetry approach and mathematical modeling of species coexistence in a heterogeneous habitat

We explore an approach based on the theory of cosymmetry to model interaction of predators and prey in a two-dimensional habitat. The model under consideration is formulated as a system of nonlinear parabolic equations with spatial heterogeneity of resources and species. Firstly, we analytically determine system parameters, for which the problem has a nontrivial cosymmetry. To this end, we formulate cosymmetry relations. Next, we employ numerical computations to reveal that under said cosymmetry relations, a one-parameter family of steady states is formed, which may be characterized by different proportions of predators and prey. The numerical analysis is based on the finite difference method (FDM) and staggered grids. It allows to follow the transformation of spatial patterns with time. Eventually, the destruction of the continuous family of equilibria due to mistuned parameters is analyzed. To this end, we derive the so-called cosymmetric selective equation. Investigation of the selective equation gives an insight into scenarios of local competition and coexistence of species, together with their connection to the cosymmetry relations. When the cosymmetry relation is only slightly violated, an effect we call 'memory on the lost family' may be observed. Indeed, in this case, a slow evolution takes place in the vicinity of the lost states of equilibrium.


1.
Introduction. The preservation of biodiversity requires a scientific understanding of the processes driving the changes of spatial distributions of species, e.g. predators and prey, both inhabiting a common biotope. In this context, the application of systems of parabolic partial differential equations describing the dynamics of populations is widely used [6,22]. Most approaches to modeling the spatial migration of species are based on homogeneous diffusion, [1,12,20,24,28]. However, from the viewpoint of biology, this simple concept of diffusion, i.e. free roaming of individuals, is deficient. It is necessary to relate the direction of migration to the living activity of the whole community, cf. [2,11]. Studies of migration under inhomogeneous living conditions, including availability of resources, growth conditions and mobility, were undertaken in recent years e.g. by [4,5,17,18,19,21].
It is an important observation that changes in the habitat and the migration of species may lead to a shift and even the loss of the ecological balance. Hence, models of the coexistence of species must take into account states of equilibrium, states of very slow evolution, and the sensitivity to data. Such characteristics are principally found in cosymmetric systems, [25,26], which exhibit non-unique stationary solutions. In the case of perturbed cosymmetry, a very slow approach to equilibrium states is typical [27].
Recently, cosymmetric models of competing populations were analyzed in [3,7,9,10]. In [9,10] the case of a system with homogeneous diffusion and a specific nonlinearity was investigated, while in [3] a nonlinear diffusion approach was considered together with a logistic growth function. Thus, the consideration of predator-prey systems with a cosymmetry relation, along with such featuring a slight violation of cosymmetry, seems very interesting.
The aim of the present paper is to apply the cosymmetry theory to an analyzis of scenarios of predator-prey interaction in a two-dimensional habitat.

2.
Model. The population density on the two-dimensional domain Ω is described by the functions u i (x, y, t), i = 1, 2, . . . , m for the prey populations and i = m + 1, m + 2, . . . , n for the predators, (x, y) ∈ Ω, t ≥ 0 denotes time. The balance of species is expressed by the migrational fluxes q i and the growth functions f i = u i ψ i , cf. [22]: The migration fluxes are given by vectors q i = (q ix , q iy ), where q ix and q iy act in the respective directions of the spatial coordinate system, a · b is scalar product in R 2 . The fluxes of the prey populations are expressed via gradients of their densities, and the function describing their common resource -which is also called carrying capacity p(x, y), cf. [3,6,22]. The last term in (2) describes the combined impact of the inhomogeneity of all species, both prey and predators, on the directed migration of species. This corresponds to phenomena of taxis [13,15].
When it comes to the predator population, the inhomogeneity of the resourse p(x) has not to be taken into account: α i = 0, i = m + 1, . . . , n.
In expressions (2), (3), the diagonal matrices k i contain diffusion rates. The remaining contributions are defined using square diagonal matrices α i and β ij . Notice that β ij appears in (2) as well as in (3), but with i from a different range. The spatial heterogeneity of the resource influences the migration fluxes of the prey species via the gradient ∇p. Obviously, in (1), p itself is present in the growth functions f i (i = 1, 2, . . . , m). Hence, the model takes into account the influence of the distribution of each of the populations on the remaining ones. For instance, the transportation of the quantity u i caused by the inhomogeneity of the distribution u j (·, ·, t) is expressed as β ij u i (·, ·, t)∇u j (·, ·, t). The sign of the coefficient β ij is responsible for the mode of the reaction on the distribution of the species u j (x, y, t), zero meaning neutrality. For now, only one single resource is taken into account, the coefficient matrices are all constant, and the influence is local and instantaneous, i.e., there are neither delays nor non-local effects included.
The natural growth of a prey population is determined by a logistic type of law with coefficients µ i > 0, while the loss due to the activity of predators is assumed to be given by the following terms with the coefficients l ij > 0: The rate of change of a predator density in a point (x, y) of the domain satisfies (5), where according to [23], the phenomenological constants are positive, µ ij > 0, l i > 0: Further, we focus on a rectangular domain Ω = [0, a] × [0, b], a, b > 0. On its boundary ∂Ω, we choose Neumann conditions. In this case, the absence of fluxes on the edges is expressed by Finally, the system (1)-(6) is completed by initial conditions on the densities of all populations: The Fisher-Kolmogorov-Petrovsky-Piskunov equation [8,16] follows from (1)-(5) in the case of m=n=1, α 1 =0 and p=const. For m=1, n=2, α 1 =β 21 =0, in the limit of large p, i.e. abundance of resources, and homogeneous initial distributions (u 0 i = const), the system (1)-(7) implies the Lotka-Volterra model [23]. Further, a model of two coexisting populations, competing for a common resource as studied in [3], follows from (1)-(5) with m = n = 2.
3. Cosymmetry. For suitable choices of parameters, the considered model (1)- (6) belongs to the class of cosymmetric dynamical systems, cf. [25], which means that a continuous family of stationary states (equilibria) appears. An important difference between symmetry and cosymmetry lays in the stability spectrum of equilibria belonging to the mentioned family. While for a system featuring a classical symmetry, the stability spectrum is identical for all members of the family, cosymmetric systems have a variable spectrum. In particular, a continuous family of stationary states may contain stable and unstable stationary solutions simultaneously. For the problem (1)-(6), the existence of a cosymmetry is established by in the following statement.
is a cosymmetry of the system (1)-(2), (6), under the following conditions: The coefficients γ ij form an antisymmetric matrix, containing two nonzero diagonal blocks. These correspond, respectively, to prey (block of dimension m × m) and predators (skew matrix of dimension (n − m) × (n − m)). In the special cases of m = 1 (only one pray) or n − m = 1 (only one predator), such blocks will be just scalars and thus are equal to zero.
Proof. Due to the definition of cosymmetry, the vector L(u) should be orthogonal to the right-hand side F (u) defined by the system (1)-(6) for any arbitrary set of functions u i (·, ·, t), i = 1, 2, . . . n. This abstract formulation breaks down to the calculation of two integrals, one collecting the prey contributions, the other for the prey.
Using integration by parts, we obtain for I: where ds denoted the length measure along the boundary. Taking in account boundary conditions (6) and after substituting (9), we can write Using (10) and the antisymmetry of γ ij the integral I 1 can be rewritten in the form: Since the diffusion coefficients k i and the functions ψ ij are restricted by the relation (11), it follows that I 1 = 0. Analogically using relations (10) we find the equalities I 2 = I 3 = 0. Due to (10) we rewrite I 4 as Thus, we derive I 4 = 0 and finally I = 0.
In an analogous way the conditions (10), (11) imposed on the diffusion coefficients k i and functions φ i , ψ i lead to the nullification of the integral J. This way, the lemma is proved, i.e., the operator L is indeed a cosymmetry of the studied system in the Hilbert space L 2 (Ω) n .
Based on the specific expressions of functions ψ i given in (4), we derive conditions on parameters, for which cosymmetry (8)-(9) holds Analogically, considering φ i in (3), we obtain conditions on parameters k i , α i , β ir (r = 1, . . . , n) By its definition, the choice of a cosymmetry is not unique, cf. [25]. In particular, any product of the vector function L(u) by a nonzero real factor again gives a cosymmetry. It should be noted that there exists a cosymmetry on a subspace. Indeed, in the case of absence of predators (u i = 0, i = m+1, . . . n) the cosymmetry for the system (1)- (7) is given as (8) with ζ i = 0 and only conditions (16) must hold. This case was considered in [3]. The predator-prey system in one space dimension, with a slightly different form of the growth term, was analyzed in [7].
In any case, a necessary condition for the appearance of a family of equilibria is the presence of at least one stationary state u, on which the cosymmetry L(u) does not vanish, i.e., it holds F (u) = 0 and L(u) = 0. For the choice of parameters (16)-(21) the system (1)-(7) has nontrivial solutions, which form a one-dimensional continuous manifold M ⊂ H.
The numerical results in Sec. 4 affirm that a discretized version of the considered system has a family of stationary states, which can be constructed by a curvefollowing strategy of predictor-corrector type.
When we model real processes, constraints like (16)-(21) will rarely be satisfied exactly. Naturally, some parameters are based on experimental data and determined with substantial error. Moreover, in ecology situations are known, where some species can prevent their extinction by modifying their mobility (diffusion rates k i ), or applying mimicry (diminishing coefficients l ij ) or chemorepulsion (by varying coefficients β ij ), [22].
However, in real-world problems even in the absence of cosymmetry, a study of conditions for the existence of a cosymmetry turns out to be useful for analyzing sets of stationary distributions of species. For the study of nearly cosymmetric situations, the notions of the cosymmetric defect and the selective equation were introduced [27]. where L is a cosymmetry of F and Q(u, ν) is a disturbation term with Q(u, 0) = 0, breaking the cosymmetry for ν = 0.
To study solutions of perturbated system we consider selective equation [27]: S(θ, ν) = 0. By definition, it is automatically fulfilled for ν = 0. In [27] a theorem was proven, which states that a nondegenerate solution of the selective equation initiates a branch of solutions parametrized by the perturbation ν.

4.
One predator and two prey species on a rectangular domain. In this section, numerical results for the system with n = 3 are presented. The corresponding finite-difference scheme is described in Appendix and based on the approach [3]. The ordinary system of differential equations (28)-(31) may be written aṡ where the vector Y = (U 1 , U 2 , . . . , U n ) T comprises the values of the unknown variables in the nodes: U i = (u i,11 , u i,12 , . . . , u i,1ny , u i,21 , u i,22 , . . . , u i,nxny ), i = 1, . . . , n.
A finite-dimensional analog of the vector field L(Y ), which is obtained from (8), (9) by an appropriate discretization, turns out to be a cosymmetry of the system (23) with respect to the Euclidean scalar product in the space of discretized solutions.
Any nonzero stationary solution Y * , satisfying Φ(Y * ) = 0, does not nullify the cosymmetry, i.e. L(Y * ) = 0. Consequently, Y * belongs to the one-dimensional family of equilibria we are looking for.
We consider the Jacobian ∇Φ(Y * ) for the study of stability. Obviously, zero is an eigenvalue, it corresponds to the neutral direction tangential to the family of equilibria. If all remaining spectral values lay in the left complex half-plane, the state Y * is a stable equilibrium. This corresponds to asymptotic stability of the solution in the manifold transverse to theconsidered family of equilibria. The families of stationary distributions described here are calculated by the algorithm laid out in the paper [9].
As computational example, we consider now two populations of prey (m = 2) and a single predator on the rectangular domain Ω = [0, a] × [0, b], where a = 2, b = 1. The matrices k 1 , k 2 , k 3 , α 1 , α 2 , β 31 , β 32 are diagonal ones: Further, suppose the migration in the direction y is negligeble. Then the matrices β 12 and β 12 are given as In the considered special case, it is more convenient to change notation. The prey densities are denoted by u(x, y, t) ≡ u 1 (x, y, t) and v(x, y, t) ≡ u 2 (x, y, t), and for the predator we write w(x, y, t) ≡ u 3 (x, y, t).
Numerical experiments are performed for Neumann conditions on the boundary and different values of migration parameters α 1 , β 12 , β 21 and the growth and loss coefficients µ k , µ ij , l k , l ij . To demonstrate a concrete cosymmetric case and its destruction, we have frozen the following parameters: Further, we choose a resource distribution p, which has two local maxima at x ≈ a/6, y = b and x ≈ 5a/6, y = b: p(x, y) = p 0 + sin (πy/2b) e x sin (3πx/a) /2 , p 0 = 0.1 , One may consider the regions around those points as favorable zones for the prey, hence they will eventually also be attractive to the predator.
Finally, the initial conditions are given in the form of prey distributions centered around the points (a/6, b/4) and (a/2, b/4), when the predator is located in vicinity of the point (5a/6, b/4).
The possibility of cosymmetry in the problem of coexistence of two closely related populations without a predator (n=m=2) was demonstrated in [3]. In the present paper, the formation of stable heterogeneous spatial patterns in the population densities of prey and predators is under consideration.
As expected, in the cosymmetric case a family of steady states is formed, while destruction of cosymmetry leads to a few isolated stationary solutions. In Table 1 we present the integral mean values of final densities.
{U , V , W } = 1 (n x + 1)(n y + 1) nx r=0 ny s=0 {u rs , v rs , w rs }. and the three elements of the stability spectrum of the stationary solutions that are closest to the imaginary axis. In the case of cosymmetry of the system (two first rows of Table 1), it turns out that by choosing initial densities with distinguished distributions within the area, one can smoothly control the resulting stationary states. These have a (near) zero eigenvalue in their spectrum of stability (numerically, it is about λ = 10 −6 ). The other eigenvalues vary, since the problem is cosymmetric, not symmetric.
The resulting densities are presented in Fig. 1. The left column of plots corresponds to a density of the population v closer to the favorable spots on the habitat. These plots represent the first row of Table 1 and one can see a predominance of the specie v in that case.  In the case of an initial distribution u 0 (x, y) arranged closer to the left of the two favorable regions, the whole rectangular area is filled with both species without a dominance in either direction, see Fig. 1, right column. The mean density of the predator exceeds the value reached in the previous experiment (left column Fig. 1). This is due to the fact that the growth coefficient of the predator with respect to the population u(x, y, t) is bigger than that with respect to v(x, y, t): µ 31 > µ 32 .
Furthermore, the directed migration in the case of nonvanishing coefficients α 1 , β 12 and β 21 is analyzed.
The left column of Fig. 2 corresponds to the case, when for the growth of the prey only the heterogeneity of the living conditions is taken into account. This matches the data in the third row of Table 1.
A comparison of the left columns of Figs. 1 and 2 shows that under a nonzero migration coefficient α 1 the inhomogeneity of the resource p(x, y) leads to a concentration of the population u in the favorable zone. This is accompanied by a drop of the density v in the favorable zone, and a compensating surge away from the loci of higher values u.
Taking into account the influence of heterogeneity of resources as well as the impact of the species in the neigborhood, we obtain the results presented in Fig. 2, right column (see also row 4 of Table 1). The cross-dependence between the populations leads to more expressed local extrema of the distributions. In the vicinity of a local peak of the species u, the species v is indented, and vice versa. Such solutions model the occurrence of spatial separation of species and the creation of ecological niches in nature.
On Fig. 3 projections of the families of stationary states onto the plane of their mean deviations σ(u), respectively σ(v), are depicted for two growth coefficients for the predator: µ 31 and µ 32 . Here σ(u) is determined on the chosen two-dimensional grid according to: In the absence of a predator, the family of stationary distributions for the prey populations is given by the straight line P Q in Fig. 3. For sufficiently large values of µ 31 and µ 32 a cosymmetric family appears, where both prey populations coexist in the presence of the predator (see curve 1 in Fig. 3). When the growth coefficient µ 31 , respectively (µ 32 ), is reduced, in the case of dominance of the u-population (resp. v), the predator may suffer extinction as indicated by curves 2 (resp. 3) on Fig. 3.

5.
Coexistence of species on one-dimensional area. For the study of a disturbance of the cosymmetry conditions, we investigate the case of a resource function that depends only on x, being constant along the y-direction of the domain: Due to the Neumann conditions on the boundary segments y = 0, y = b, we arrive at a spatially one-dimensional problem, where β 12 , β 21 and α 1 , α 2 are scalar quantities.
Firstly, we consider the system of two prey without predator. In Fig. 4 graphs of the selective function (22) are presented for various values of the migration coefficients. Here the continuous parameter θ points to the stationary state on the family. The figure illustrates the possibility of coexisting prey populations. Thus, loss of cosymmetry leads to the destruction of the family and the appearance of isolated sulutions.
For coefficients of the same sign, a stationary distribution of coexisting populations is obtained (see curves [1][2][3]. Depending on the ratio β 12 /β 21 we find a shift of the point, where the selective function attains the zero level. However, for different signs of the migration coefficients, i.e. β 12 β 21 < 0 (curve 4), there are only solutions with one of the prey species gone extinct.  For decreasing parameter β 12 = β 21 we observe the convergence to the family. The corresponding transient process is rather slow, so that the approach to the final state takes more time.
Perturbation of cosymmetry may lead to coexistence of populations as well. Fig. 6 displays regions on the parameter plane (α 1 , α 2 ) with different scenarios. There are three courses of events, at the end of which either both prey species are present (III), or one of them is present and the other extinct (I and II). The regions corresponding to the different possible scenarios intersect in a common point on the straight line k 1 α 2 = k 2 α 1 , which presents the cosymmetric case. An influence of the predator on the location of this common point can be inferred (the bottom of Fig. 6).  Figure 6. Map of migration parameters, corresponding to coexistence (III) or survival of only one of the species, u (I) or v (II), without predator (top) and with predator (bottom). The thick line indicates the existence of a continuous family of stationary states Conclusion. A model of migration of competing prey and predator species in an inhomogeneous habitat is under investigation. The case of migration fluxes depending on the spatial inhomogeneity of the densities of populations and resources is considered. The system is formulated in terms of nonlinear parabolic partial differential equations, which under certain additional conditions belong to the class of cosymmetric dynamical problems [25]. On the basis of the Method of Lines, combined with staggered grids, a numerical study of scenarios of extinction as well as of coexistence of species is performed. The influence of the predator on the formation of distributions of the prey is discussed. Ranges of parameters are found for which nonunique stationary solutions can be observed. In the case of cosymmetry, the asymptotic equilibrium state depends continuously on the initial distribution of populations. Effects of predominance of one of the species in a certain subdomain are studied. Special parameter values were distinguished, which lead to a division of the biotope and the formation of ecological niches.