BIFURCATION AND SPATIOTEMPORAL PATTERNS IN A BAZYKIN PREDATOR-PREY MODEL WITH SELF AND CROSS DIFFUSION AND BEDDINGTON-DEANGELIS RESPONSE

. In this article, the clasical Bazykin model is modifed with Bedding– ton–DeAngelis functional response, subject to self and cross-diﬀusion, in order to study the spatial dynamics of the model.We perform a detailed stability and Hopf bifurcation analysis of the spatial model system and determine the direc- tion of Hopf bifurcation and stability of the bifurcating periodic solutions. We present some numerical simulations of time evolution of patterns to show the important role played by self and cross-diﬀusion as well as other parameters leading to complex patterns in the plane.


1.
Introduction. Alan Turing [17], first proposed the reaction-diffusion theory for pattern formation in his seminal work on the chemical basis of morphogenesis. He suggested that, under certain conditions, chemicals can react and diffuse in such a way as to produce steady state heterogeneous spatial patterns of chemical or morphogen concentration [11]. With the pass of years the creation of such patterns, which we now call Turing patterns, has been playing significant roles in mathematical biology and other areas [11]. Particularly, the existence of Turing patterns in models of population dynamics, where individuals move not only through time but a space, has been extensively studied by many authors, taking into account the effect of self-as well as cross-diffusion, see for example [4,12,3,23,7]. The term self-diffusion implies the movement of individuals from a higher to a lower concentration region, influenced only by its own density, i.e. there is no response to the density of the other one. On the other hand, cross-diffusion implies that the per capita diffusion rate of each species be influenced by the other ones.
The self diffusion coefficient is a positive value, while cross-diffusion coefficient may be positive, negative or zero. A positive cross-diffusion coefficient denotes the movement of the species in the direction of lower concentration of another species while the negative cross-diffusion coefficient for one species tends to diffuse in the direction of higher concentration of another species [23].
Even though the Lotka-Volterra predator prey model has been extensively studied by many authors, it has been questioned by several biologists. Thus, the predator-prey model with the Beddington-DeAngelis functional response (first proposed by Beddington [2] and DeAngelis [5]) has been well studied [13]. The model proposed by Beddington and DeAngelis was the following: where f12x2 b2+x2+W12x1 is the functional response. It is similar to the well-known Holling type-II functional response but has an extra term in the denominator which models mutual interference between predators [13]. The ecological background and parameters can be seen in [5].
Zhang et al [23] considered a predator prey model with Beddington-DeAngelis functional response, self-and cross-diffusion. The model was the following: Where u(t, x), v(t, x) stand for the population densities of the prey and predator species at time t and space location x and ∆ is the Laplacian operator. The functions f 1 and f 2 were the following : Functions f 1 , f 2 are a particular case of the original ones used by DeAngelis in (1) under a change of variable.
Constants d 1 , d 2 > 0 are the diffusion rates, while α i,j are the cross diffusion rates. Functionsû =û(x, y, t),v =v(x, y, t) represent the population of prey and predator respectively at time t in the spatial coordinate (x, y) ∈ Ω; we consider a rectangular spatial domain Ω = [0, L x ] × [0, L y ]. ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 is the Laplacian operator,f andĝ are the kinetics defined aŝ System (4) is subject to zero flux boundary conditions : where n is the normal outward vector to ∂Ω. There are a lot of reasons for choosing zero flux conditions, the main one is that we are interested in self-organization of pattern; these conditions imply that there is not external input. If we put fixed boundary conditions then patterns formed for u and v could be a consequence of the fixed conditions [11]. Letû = ru e ,v = rpv e , T = rt (see [13]), then (4) is reduced to the following non-dimensional system: Where = bp Ae , α = r Ae , β = Brp Ae , γ = c r , δ = hp e , T = rt (using t instead T for notational convenience), δ 1 = d1 r , δ 2 = d2 r , α 11 = a11 e , α 21 = a21 e , α 12 = a12p e , α 22 = a22p e . a ij could be positive, negative or zero and other parameters are positive. The functions f, g are the following From now on and in the following chapters we will work with system (8).

2.
Basic properties of the model. The model has three equilibria points, the trivial equilibrium E 0 = (0, 0), E 1 = (1, 0), and an internal equilibrium, namely whenever − β + βu * = 0. u * is given by the positive root of equation It is clear that equation (10) has at least one real root, because it is a cubic equation, moreover the leading coefficient is positive. If > β, the constant term of (10) is negative, so if we consider condition > β we have that equation (10) possesses at least one positive root u * + . For notational convenience we set u * = u * + . However, there is the possibility that v * takes negative values, even when u * is non-negative, to avoid this situation we need also Clearly, 1 + αu * > 0 for u * ≥ 0 . If u * ≥ 1 then we have v * ≤ 0, so we ask for u * < 1. From equation (10) we have Then, if we guarantee, − γα − γ > 0 the positive root u * will be in interval [0, 1]. Taking into account that 0 < u * < 1 we have u * − 1 > −1, or equivalently: Therefore, we can enunciate the following theorem: Theorem 2.1. Let > β and γ < 1+α , there exists a positive internal equilibrium E * .
We denote w = (u, v) T , F (w) = (f, g) T , Φ(w) = (Φ 1 (w), Φ 2 (w)), where: Then system (8) can be rewritten as: The linearization of system (12) is given by ( [23]): F w (w * ) denotes de Jacobian matrix of F evaluated in (u * , v * ). Similarly for Φ w (w * ). In particular, where: From (13), the linearization can be written in terms of u and v as the following form: Any solution of system (15) can be expanded into a Fourier series [9] u(r, t) = ∞ n,m=0 where r = (x, y), and 0 < x < L x , 0 < y < L y , k n,m = (k n , k m ) and k n = nπ/L x , k m = mπ/L y are the corresponding wave numbers. By substituting u nm and v nm into (15), we obtain where k 2 = k 2 n + k 2 m . A general solution of (17) has the form C 1 e λ1t + C 2 e λ2t , where the eigenvalues λ 1 and λ 2 are given by the determinant of the following matrix From the definition of k 2 , for n = 0, 1, ... and m = 0, 1, ..., the set of all k 2 is numerable, ie, we have k 2 1 , k 2 2 , ..., k 2 i , ... for each i = 0, 1, 2, ..., and a pair of eigenvalues λ i 1 , λ i 2 for each i. λ i 1 , λ i 2 are solutions of the following equation According to Routh-Hurwitz, we get the following.  • If det(J L ) > 0 for all k 2 and there exist a k 2 > 0 such that tr(J L ) > 0 then the positive equilibrium E * is unstable. • If det(J L ) < 0 for a k 2 > 0 then the equilibrium is a saddle.
When at least one of the following conditions is violated, the system is unstable. For having Turing instability, suppose (u * , v * ) is stable in the absence of any spatial effects, i.e., the conditions tr(J) < 0 and det(J) > 0 hold. If d 11 + d 22 < 0 then there exists always a k 2 large enough such that 0 > J 11 + J 22 > k 2 (d 11 + d 22 ), so we have Turing instability.
If d 11 + d 22 > 0, the first condition tr(J L ) < 0 is not violated. Hence the only condition that could not be true is that det(J L ) > 0, which would gives rise to diffusion instability. Thus, the condition for diffusive instability is that there exist a k 2 > 0 such that: H(k 2 ) is a parabolic function and its sign depends on the cross diffusion parameters. • is open down and, no matter the value of other parameters, there exists always a large enough k 2 N > 0 such that H(k 2 N ) < 0 , therefore we have Turing instability.
is on the other quadrants, eventually the function will fall in fourth quadrant ( Turing instability ). For each k 2 we have a pair of eigenvalues, λ 1 , λ 2 , figure (1)(b) shows the graph of the real part of one of these eigenvalues, say λ 1 ,which is sufficient to ensure the existence of Turing instability. Again, dotted line is for α 12 = −1 and solid line is for α 12 = −2. We can see that λ 1 has positive real part and therefore we have instability.
is open up, and its minimum occurs at k 2 = k 2 m , where: According to this, a sufficient condition for Turing instability is that H(k 2 m ) is negative for a k 2 m > 0. Therefore or equivalently Figure (2) shows this case for same values as used in figure (1). det(Φ w (w * )) > 0 for α 12 ∈ (−0.04298, ∞) and we have that H(k 2 m ) = 0 for two values, −0.0361... = α 1 12 and 0.2237... = α 2 12 . If α 12 = α 1 12 then we have that (λ 1 ) touches the x axis (solid line) and its real part is negative for all k 2 . When α 2 12 > α 12 > α 1 12 then H is above the x axis and the eigenvalue has always negative real part (dotted line), but when α 12 < α 1 12 the parabola takes negative values and the real part of λ 1 is positive, so we have Turing instability (dashed line). According with previous results we have the following.
3. Hopf bifurcation. The Hopf bifurcation theorem gives a set of sufficient conditions to ensure that an autonomous differential equation with a parameter exhibits non-trivial time-periodic solutions for certain values of the parameter. Diffusive predator prey models with suitable functional response can generate temporal periodic patterns via Hopf bifurcation [15]. Our goal in this section is to obtain necessary conditions for the existence of a Hopf bifurcation for the reaction diffusion system with cross diffusion, (8). For the sake of completeness, we first consider the case of the local system where the diffusion parameters are zero.
3.1. Hopf bifurcation of the local system. The local system is given by the following system of ordinary differential equations: The equilibria points of the ODE's system are the same as the ones for the diffusive system, so we consider the positive endemic equilibrium E * = (u * , v * ) given in section 2. From the calculation in section 2 we know that the linearization of (27), is given by the Jacobian matrix: where the J ij 's are the ones in section 2. Then the characteristic equation is given by the quadratic polynomial: We are interested in the existence of a pair of complex conjugate roots of (29) with zero real part, and this happens if and only if det(J) > 0 and T r(J) = 0.
Then tr(J) = 0 iff s = 0. We analyze the existence of a Hopf bifurcation using s as the bifurcation parameter. To understand the detailed property of Hopf bifurcation, we need a further analysis of the normal form.
Assume det(J) > 0 and denote Λ = √ det J, then at s = 0 the Jacobian matrix (28) has two imaginary eigenvalues λ = ±Λi. The eigenvalues λ(s),λ(s) vary smoothly with s so we can affirm that near s = 0 the eigenvalues are still complex and its real part is given by: Therefore, by Hopf bifurcation theorem ([6] Theorem 3.4.2) system (27) has a small amplitude non constant periodic solution bifurcated from the interior equilibrium E * when s crosses through 0, ie, system (27) undergoes a Hopf bifurcation.
To understand the detailed behaviour of the system around s = 0 we need a further analysis of the normal form [15], readers can see chapter 3 of [6] for a deep study of bifurcation and normal forms.
We make a change of coordinates x = u − u * , y = v − v * to obtain an equivalent system to (27) with an equilibrium at (0, 0) in the x-y plane. Under this change of variables, the system obtained is: Where: It is not difficult to show that J is the Jacobian matrix of (32) at (0, 0), then by direct calculation: and s = 0 is equivalent to a 11 = −a 22 . Separating linear and non-linear terms of (32) we have: Where From the equality (34) we can compute the complex eigenvectors of matrix J and take its real and imaginary part as follows: Remark 1. Recall that we have assumed that > β, and in this case a 12 = 0.

3.2.
Hopf bifurcation of the reaction-diffusion model. In section 3.1 we obtained the conditions for the existence of temporal periodic solutions for system (8) when the diffusion coefficients are all zero (which gives us an ordinary differential equations system) and the results are given in theorem (3.1). Now we want to obtain a similar theorem for the case when the diffusion coefficients are not all zero.
Using the characteristic equation computed in (19) and the definition of s used in (30), we can rewrite the trace and determinant of J L as: For each k 2 in a numerable set {k 2 i , i = 0, 1, 2, ...}. For each i we have a pair of eigenvalues of matrix J L , λ section 2), and all these eigenvalues are eigenvalues for the original system with diffusion. Assume that all eigenvalues λ j are simple and denote its corresponding eigenfunction by φ j (r). We use the bifurcation parameter s (like in section 3.1) and recall the following necessary and sufficient condition to obtain a Hopf bifurcation at a value s = s H [20,24].
(H1) There exist i ∈ {0, 1, 2, ...} such that  (23), (24) we have that a sufficient condition for having det(J L )(0, k 2 ) = 0, ∀k 2 = 0 is det(Φ w (w * )) > 0 and H(k 2 m ) > 0, which is equivalent to 2 Under these conditions we have a single pair of complex eigenvalues with zero real part given by λ = ±iΛ and near the imaginary axis (see equation (31) and (47)): Therefore we have the following theorem: then system (8) undergoes a Hopf bifurcation at s = 0, the periodic solutions are spatially homogeneous and coincide with the orbits of the corresponding local system.
To determine the stability of bifurcated periodic solutions we need to know the behaviour of the system (8) in its center manifold at the bifurcation point ( as with ODE's systems ). Many authors have studied the center manifold of a reactiondiffusion system with 1D and 2D diffusion, as well as its normal form and the existence hopf bifurcation, see for example [8,15,18,20,24]. For this model we follow the methodology used for [8,15,18,20].
Under the change of coordinates used in the previous section, X = u − u * , Y = v − v * , and setting u = X, v = Y for notational convenience, the system (8) can be rewritten as: Separating linear and nonlinear part: Where D is a matrix, D = (d ij ), J is the jacobian matrix taken as (34) and F 3 is given by Functions F 1 2 and F 2 2 are given by (36), (37). Define the real-valued Sobolev space And X C = X + iX. Let , be the complex valued inner product in L 2 (Ω) × L 2 (Ω) 1 defined as with α = (α 1 , α 2 ) T and β = (β 1 , β 2 ) T in X C . Define the adjoint operator of L as With domain D L * = X C = D L . It is not difficult to check that For any eigenvalue λ j (which depends of k 2 ) of L we can find an eigenfunction of the form φ j = (q 1 , q 2 ) T cos kr where q 1 and q 2 satisfy the condition: Particularly, when k 2 = 0 (at s H = 0) and λ = ±Λi, we have a single pair of complex eigenvectors (q 1 , q 2 ), (q 1 ,q 2 ) , where q 1 , q 2 satisfy (57).

Numerical simulations.
4.1. Local system. We give some numerical simulations to illustrate the results obtained in theorem (3.1). Set the parameters as follows: α = 10, = 5, γ = 1/11, δ = 0.3 and let β be arbitrary. Then the system is: We have that γ < /(1 + α) = 5/11, therefore, if 5 > β, by Theorem (2.1) there exists a positive internal equilibrium By Theorem (3.1), the Hopf bifurcation occurs around the interior equilibrium E * at s = 0, but E * depends on β, and s depends on E * and β, however we do not have an explicit expression for E * ( specifically for u * ). To find a possible Hopf bifurcation value , define the following functions: (69) With v * given above. For each β, implicit function g = 0 provide us with points where u * is an equilibrium, while function s = 0 gives, for each u * , the value where Hopf bifurcation occurs, therefore, the intersection of both curves provide us with the value where u * is an equilibrium and there is a Hopf bifurcation. From figure (3) we can see that both curves intersect each other between β = 1 and β = 2.
Let β = 1.1, the equilibrium is given by u * = 0.16290828513399125233, v * = 0.53951361805065861059. 4.2. Reaction diffusion system. We perform numerical simulations of model (8) in two-dimensional space with parameters satisfying Turing conditions, to obtain the so called Turing patterns. We classify the patterns obtained in spots (isolated zones with high population densities), holes (isolated zones with low population density), stripes (like zebra's skin) or mixed, where we have the spot-stripe mixture or the hole-stripe mixture [19]; however this classification is based on human perception and there is not an exact definition of each class of pattern [14]. To solve differential equations by computer, one has to discretize the space and time of the problem. The continuous problem defined by the reaction-diffusion system in two-dimensional space is solved in a discrete domain (x i , y j ) and t k . The spacing between the lattice points of space is defined by the lattice constant ∆h, so x i = i∆h, y j = j∆h. The time evolution is also discrete, the time goes in steps of ∆t, so t k = k∆t. Define the value of functions u and v at each point of grid as u(x i , y j , tk) = u k ij and v(x i , y j , t k ) = v k ij . In the discrete system the laplacian describing diffusion is calculated using finite differences, i.e., the derivatives are approximated by differences over ∆h. First, we use the following properties: Where ∆(u 2 ) = 2u∆u + 2(∇u) 2 , ∆(v 2 ) = 2v∆v + 2(∇v) 2 , ∆(uv) = (∆u)v + 2(∇u)(∇v) + u∇v.
To approximate the laplacian ∆, we use the standard five point approximation with zero flux boundary conditions, while for gradient ∇ we use central difference approximation: 2∆h .
The time evolution can be solved by using the finite difference approach: Let x = y = (0, ∆h, 2∆h, ..., N ∆h) and t = (0, ∆t, 2∆t, ..., t max ) with N = 200. Fix also the step-size ∆ h = 0.25 and ∆ t = 0.01. We use as initial conditions the value of equilibrium (u * , v * ) with a perturbation of order 10 −3 and run the program for multiple time steps t max and different values of α 12 . We plot the results obtained for u and v with a color bar graphic where the x and y labels are the position (x, y) in the space, and the color is the density of population (the value of u, v). Each figure represents the population at a time step t max , so if we want the population density for 0 < t < a, we need a set of figures for each time step t i between 0 and a. Some of the patterns were obtained for large times t max , we run the simulations for each value of α 12 until a time size t max after which, we did not observe a significant change in the figure. Moreover, when det(Φ w (w * )) > 0, the second condition of theorem (2.4) for this case is valid for α 12 ∈ (−0.028292, 0.11354), so, we use values for α 12 in the interval above. To see the initial condition, figure (5) shows the results for α 12 = 0.0001 with t max = 0 (a and b), where we are located in a perturbation of the equilibrium. When we go forward in time, for t max = 50, we observe that some figures begin to appear (c and d). These figures become more visible when t max = 200 (e and f). The patterns obtained in this example are called spots, zones with high concentration of population with a visible boundary in a low concentration background. Once we have the first two conditions of theorem (2.4), we can see that det(Φ w (w * )) > 0, iff α 12 ∈ (−0.013849, ∞), condition (26) is valid for α 12 ∈ (−0.013849, 0.024327). We take α 12 = 0.0001 and run the program for t max = 400, the pattern obtained here is showed in figure (6) and it is called holes. They are similar to spots, but now, there are isolated zones with low population density in a background of high concentration, like holes of population. Example 3. Now, we take α 12 = 8 and fix other parameters as in previous example. The equilibria point for these values is (0.12929, 0.47952), with det(J) = 0.053714, tr(J) = −0.00741, det(Φ w (w * )) = 0.47952 α 12 + 0.01 and it is positive iff α 12 ∈ (−0.020854, ∞). Moreover, condition (26) for Turing instability is equivalent to α 12 ∈ (−0.020854, 0.22219). If we take α 12 = 0.0001, then all conditions for Turing instability are satisfied and some patterns appear. Figure (7) shows the pattern obtained for this example, with time t max = 1000, this pattern is formed by stripes (like zebra's skin) and we can observe also some spots in it, so we could call it a spot-stripe mixture pattern. From theory, we have that there exists Turing patterns for all α 12 that satisfies conditions of theorem (2.4), but numerically we can not obtain results for all values. For instance, in example 2, the patterns should appear for α 12 ∈ (−0.013849, 0.024327), but when we run program for α 12 = −0.01 or less, the We have similar situation for α 12 > 0.01 and for example 1 and 3, that is the reason why some authors prefer to use the implicit methods more than the explicit ones [10].

5.
Conclusions and remarks. In this paper we study the spatial dynamics of the Bazykin predator prey model, modified with Beddington-DeAngelis functional response, subject to self as well as cross-diffusion and zero flux boundary conditions. We also study the conditions for Turing instability and pattern formation. The most important result obtained, in theory and numerical simulations, is that pattern formation can be induced by cross diffusion parameters (theorem (2.4)). These results are very important in real life, and there exists a lot of examples for the application of them (see for instance [1,11]). A structure with spots in prey and predator means there exist isolated zones with high concentration of population, while the hole structure shows the existence of zones with low or none population in the territory; therefore this model can be used to describe the population dynamics for a predator and prey population. However, in this paper we assume that parameters are all constant and are not affected by natural elements of the zone, which is not true in real life, because there exists a lot of environmental factors that could affect them, for example temperature, dampness, migration, etc. Considering these situations the parameters are not constant and present some fluctuation due to this natural noise [19], it could be a topic for future work. The patterns obtained with the Turing bifurcation in section 2 are of spatial type, i.e. for a fixed value of t we observe periodic solutions of u in the variables (x, y) (it can be observed in the numerical simulations). We obtained also the conditions for existence of temporal patterns in section 3, where Hopf bifurcation theorem is used to find periodical solutions in time. Now, an interesting topic is, when Turing and Hopf bifurcation occur simultaneously? Is that possible at all? What does it mean biologically? For example in [15] authors make an analysis of the existence of Turing-Hopf bifurcation for a 1-D reaction-diffusion system and put some numerical examples where patterns are observed in time and space.
A useful tool for study of pattern formation is the amplitude equation, which can be used to determine the amplitude of the patterns that are formed when Turing or Hopf bifurcation occur. This equation is obtained via the weakly nonlinear analysis, a technique based on the multiples scales method, where the bifurcation parameter (Turing or Hopf, whatever the case) is perturbed to analyse the behaviour of solutions where the parameter is near to a bifurcation threshold. Reader can see the work of [16,21,22], for a more detailed analysis of these methods.
Numerically, we used a finite difference scheme for our cross diffusion system and the graphs obtained are very informative, but it could be possible to use another numerical method, for example finite element method for space scheme or an implicit method for time scheme [10] to obtain results for larger values of t max and maybe with less computation time. Recall that we run the simulations for negative values of the cross diffusion parameter α 12 , and not only positive ones.
We believe that our results in the Bazykin's model could help mathematicians, ecologists and interested people to understand the spatial dynamics of predator and prey populations and motivate them for more complex models and investigations.