2007 International Conference in Honor of Claude Lobry

Complexity in a prey predator model

Baba I. Camara* — Moulay A. Aziz-Alaoui

Laboratoire de Mathématiques Appliquées du Havre,
Université du Havre, 25 rue Philippe Lebon
76600 Le Havre, FRANCE
camarab@univ-lehavre.fr, aziz.alaoui@univ-lehavre.fr

* Corresponding author

**ABSTRACT.** In this paper we consider a predator-prey model given by a reaction-diffusion system. It incorporates the Holling-type-II and a modified Leslie-Gower functional response. We focus on qualitative analysis, bifurcation mechanisms and patterns formation.

**RÉSUMÉ.** Nous considérons un modèle proie-prédateur exprimé sous forme de système de réaction-diffusion. En absence de diffusion, le système étudié est de type Holling-type-II et la réponse fonctionnelle une forme modifiée du terme de Leslie-Gower. Dans cet article, nous nous intéressons à l'analyse qualitative des solutions, l'étude des bifurcations et la formation de motifs spatio-temporels.

**KEYWORDS :** Predator-prey model, Stability, Bifurcations, Spatiotemporal chaos, Self-organization.

**MOTS-CLÉS :** Modèle prédateur-proie, Stabilité, bifurcations, chaos spatio-temporel, auto-organisation
# 1. Introduction

The dynamical relationships between species and their complex properties are at the heart of many important ecological and biological processes. Predator-prey dynamics are well-studied in the process of control and conservation of some ecosystems. We assume that only basic qualitative features of the system are known, namely the invasion of a prey population by predators. The local dynamics has been studied in [4, 8]. This model incorporates the Holling-type-II and a modified Leslie-Gower functional responses. Without diffusion it reads as,

$$ egin{cases} \frac{dH}{dT} = \left(a_1 - b_1 H - \frac{c_1 P}{H + k_1}\right) H \\ \frac{dP}{dT} = \left(a_2 - \frac{c_2 P}{H + k_2}\right) P \end{cases} \qquad (1) $$

with,

$$ H(0) \geq 0, P(0) \geq 0. $$

H and P represent the population densities at time T. $r_1$, $a_1$, $b_1$, $k_1$, $r_2$, $a_2$, and $k_2$ are model parameters assuming only positive values. $a_1$ is the growth rate of preys H. $a_2$ describes the growth rate of predators P. $b_1$ measures the strength of competition among individuals of species H. $c_1$ is the maximum value of the per capita reduction of H due to P. $c_2$ has a similar meaning to $c_1$. $k_1$ measures the extent to which environment provides protection to prey H. $k_2$ has a similar meaning to $k_1$ relatively to the predator P.

The historical origin and applicability of this model is discussed in details in [4, 8, 15, 16, 6]. The corresponding PDE version has been first done and partially studied in [9].

This paper is organized as follows. In Section 2, we prove the global existence of solutions, we also study the stability of the positive steady states. We investigate complex pattern formation and spatiotemporal Chaos emergence.

# 2. Global existence of solutions

The mathematical model we consider here consists of reaction-diffusion equations which express conservation of predator and prey densities. It has the following form,

$$ egin{cases} \frac{\partial H}{\partial T} = D_1 \Delta H + \left(a_1 - b_1 H - \frac{c_1 P}{H + k_1}\right) H \\ \frac{\partial P}{\partial T} = D_2 \Delta P + \left(a_2 - \frac{c_2 P}{H + k_2}\right) P \end{cases} \qquad (2) $$

$H = H(T, X)$ and $P = P(T, X)$ are the densities of preys and predators, respectively. $\Delta$ is the laplacian operator. $D_1$ and $D_2$ are the diffusion coefficients of prey and predator

## A R I M A
respectively.

To investigate problem (2), we introduce the following scaling tranformations,

$$t = a_1 T, \quad x = X\left(\frac{a_1}{D_1}\right)^{\frac{1}{2}}, \quad y = Y\left(\frac{a_1}{D_1}\right)^{\frac{1}{2}}, \quad u(t) = \frac{b_1}{a_1}H(T), \quad v(t) = \frac{c_2b_1}{a_1a_2}P(T)$$

$$a = \frac{a_2c_1}{a_1c_2}, \quad b = \frac{a_2}{a_1}, \quad e_1 = \frac{b_1k_1}{a_1}, \quad e_2 = \frac{b_1k_2}{r_1}, \quad \delta = \frac{D_2}{D_1}$$

We obtain the following equations, for the local model,

$$\left\{ \begin{aligned} \frac{dE_1}{dt} &= E_1 \left( 1 - E_1 - \frac{aE_2}{E_1 + e_1} \right) = f(E_1, E_2) \\ \frac{dE_2}{dt} &= bE_2 \left( 1 - \frac{E_2}{E_1 + e_2} \right) = g(E_1, E_2) \end{aligned} \right. \quad (3)$$

and the following system for the spatio-temporal equations :

$$\left\{ \begin{aligned} \frac{\partial u(t, x)}{\partial t} &= \Delta u + u \left(1 - u - \frac{av}{u + e_1}\right) = \Delta u + f(u, v), \quad x \in \Omega, \ t > 0 \\ \frac{\partial v(t, x)}{\partial t} &= \delta \Delta v + bv \left(1 - \frac{v}{u + e_2}\right) = \delta \Delta v + g(u, v), \quad x \in \Omega, \ t > 0 \end{aligned} \right. \quad (4)$$

We consider the Neumann boundary conditions given by,

$$\frac{\partial u}{\partial \nu} = \frac{\partial v}{\partial \nu} = 0$$ $$u(0, x) = u_0(x), \quad v(0, x) = v_0(x), \quad x \in \Omega \subset \mathbb{R}^n, \ n = 1, 2, 3.$$

Here, $\Omega$ is a bounded domain, the initial data $u_0$ and $v_0$, are non-negative functions.

$\frac{\partial u}{\partial \nu}$ and $\frac{\partial v}{\partial \nu}$ are respectively the normal derivatives of $u$ and $v$ on $\partial \Omega$.

## 2.1. Global existence

By standard existence theory, e.g. see [1], [3], and [2], it is not difficult to establish the local existence of the unique solution $(u(., t), v(., t))$ of (4) for $0 = t < T_{max}$, where $T_{max}$ is determined by $u_0(x)$ and $v_0(x)$. Now we establish the global existence by proving that for any finite time $T$, $||u(., t)||_{L^\infty}$, $||v(., t)||_{L^\infty}$ are bounded for $0 \le t < T$.

**Theorem 1** For any smooth nonnegative functions $u_0(x)$ and $v_0(x)$, such that,

$$\left\{ \begin{array}{l} u_0(x) \le 1 \\ \max_{\Omega} (u_0(x) + v_0(x)) \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b}, \end{array} \right. \quad (5)$$

system (4) has a unique smooth global solution for $t > 0$.

*Proof.* First, it is easily seen that $u(x, t) \ge 0$ and $v(x, t) \ge 0$ since $(0, 0)$ is a sub-solution of each equation of (4). We have

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial u(t, x)}{\partial t} \le \Delta u + u(1-u) \\ \displaystyle \frac{\partial u}{\partial \nu} = 0, \quad t > 0 \\ u(x, 0) = u_0(x) \le u_{01} \equiv \max_{\Omega} u_0(x). \end{array} \right. \quad (6)$$
$\frac{\partial u}{\partial v}$ is the normal derivative of $u$ on $\partial\Omega$.

By the comparison principle, we have $u(x, t) \le u_1(t) \le 1$,

where, $u_1(t) = \frac{u_{01}}{u_{01} + (1 - u_{01})e^{-t}}$ is the solution of the initial value problem:

$$ \left\{ \begin{aligned} \frac{du_1}{dt} &= u_1(1 - u_1) \\ u_1(0) &= u_{01} \le 1 \end{aligned} \right. \qquad (7) $$

From the second equation of system (4) we have

$$ \frac{\partial v}{\partial t} = \delta \Delta v + bv \left( 1 - \frac{v}{u + e_2} \right), $$

By the comparison principle, we deduce that

$$ \frac{\partial v}{\partial t} \le \frac{dE_2}{dt} $$

where $E_2$ is the solution of the second equation of system (3) satisfying $E_2(0) = \max_{\Omega} v_0(x)$. Thus by comparison principle and from [4, 8],

$$ \frac{\partial v}{\partial t} \le \frac{dE_2}{dt} + \frac{dE_1}{dt}. $$

Let us denote by $\sigma = E_2 + E_1$,

$$ \frac{\partial v}{\partial t} \le \frac{d\sigma}{dt} $$

From [4, 8] we have,

$$ \frac{d\sigma}{dt} \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b} - \sigma $$

Since $\sigma(0) \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b}$ and by Gronwall lemma we obtain,

$$ \sigma \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b} $$

Thus,

$$ v \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b}. $$

The proof is complete. $\square$

**Theorem 2** The domain given by

(i) $A = [0, 1] \times \left[ 0, \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b} \right]$ is a positively invariant region for the global solutions of system (4).

(ii) The solutions of problem (4), which initial conditions are in $\mathbb{R}^+ \times \mathbb{R}^+$ converge towards A.
*Proof.* For any initial condition $(u_0(x), v_0(x))$ of system (4) we have by comparison principle,

$$u \le E_1, v \le E_2, \text{ with } E_1(0) = \max_{\Omega} u_0(x) \text{ and } E_2(0) = \max_{\Omega} v_0(x)$$

and from ([4, 8]) we have

$$\overline{\lim}_{t \to +\infty} E_1(t) \le 1,$$

$$\overline{\lim}_{t \to +\infty} (E_1(t) + E_2(t)) \le \frac{5}{4} + \frac{(1+b)^2(1+e_2)}{4b}.$$

This completes the proof □

## 2.2. Stability of steady states

The steady states $(u(x), v(x))$ of system (4) satisfy,

$$\left\{
\begin{array}{l}
\Delta u + u \left(1 - u - \frac{av}{u + e_1}\right) = 0, \quad x \in \Omega \\
\delta \Delta v + bv \left(1 - \frac{v}{u + e_2}\right) = 0, \quad x \in \Omega \\
\frac{\partial u}{\partial \nu} = \frac{\partial v}{\partial \nu} = 0.
\end{array}
\right.
\qquad (8)$$

From analysis above and from [4, 9, 8], it is easily seen that system (8) has the following nonnegative solutions:

$$S_0 = (0, 0)$$

$$S_1 = (1, 0)$$

$$S_2 = (0, e_2)$$

(9)

$S_3 = (u^*, v^*)$, where $(u^*, v^*)$ is the positive solution of the system

$$\left\{
\begin{array}{l}
1 - u^* - \frac{av^*}{u^* + e_1} = 0 \\
1 - \frac{v^*}{u^* + e_2} = 0.
\end{array}
\right.$$

$S_4 = (u(x), v(x))$, where $u(x)$ and $v(x)$ are two positive functions.

In this section, we are going to investigate the linear stability of the above equilibrium solutions $S_i$ of system (8) whose existence has been proved in last section. It is well-known (see [11]) that the stability question for $S_i$ is answered by considering the corresponding eigenvalue problem for the linearized operator around $S_i$. Let us substitute

$(u(x, t), v(x, t)) = S_i + W(x, t) = S_i + (w_1(x, t), w_2(x, t))$ into system (4) and then pick up all the terms which are linear in $W$:

$$\frac{\partial W}{\partial t} = D \Delta W + L(S_i)W, \qquad (10)$$

where

$$D = \mathrm{diag}(1, \delta),$$
$$L(u, v) = \left( \begin{array}{cc} 1 - 2u - \frac{ae_1 v}{(u+e_1)^2} & -\frac{au}{u+e_1} \\ \frac{bv^2}{(u+e_2)^2} & b - \frac{2bv}{u+e_2} \end{array} \right). \qquad (11)$$

**Proposition 1** $S_0 = (0,0)$ is unstable.

*Proof.* From equation (10), the linearized system of equation (4) around $S_0$ is

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial w_1}{\partial t} = \Delta w_1 + w_1, \quad x \in \Omega \\[1em] \displaystyle \frac{\partial w_2}{\partial t} = \delta \Delta w_2 + bw_2, \quad x \in \Omega \\[1em] \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0. \end{array} \right. \qquad (12)$$

Now we study the following eigenvalue problem :

$$\left\{ \begin{array}{l} \Delta w_1 + w_1 = \eta w_1, \quad x \in \Omega \\ \delta \Delta w_2 + bw_2 = \eta w_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0. \end{array} \right. \qquad (13)$$

To prove Proposition 1, we need to prove that the largest eigenvalue of system (13) is positive. Let $\eta$ be an eigenvalue of system (13) with eigenfunction $(w_1, w_2)$. If $w_1 \neq 0$, then $\eta$ is an eigenvalue of $\Delta + 1$ with homogeneous Neumann boundary condition. Therefore, $\eta$ must be real. Similarly, if $w_2 \neq 0$, $\eta$ is also real. Hence all eigenvalues of system (13) are real. Let $\eta_1$ be the largest eigenvalue of system (13). The principal eigenvalue $\lambda_1$ of

$$\left\{ \begin{array}{l} \Delta w_1 + w_1 = \lambda w_1, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (14)$$

is positive and the associated eigenfunction $\tilde{w}_1 > 0$. We claim that $\lambda_1$ is also an eigenvalue of system (13). In fact, we take $w_2 \equiv 0$ then $(w_1, w_2) = (\tilde{w}_1, 0)$ satisfies system (13) with $\eta = \lambda_1$. So $\lambda_1 > 0$ is an eigenvalue of system (13). Therefore we must have $\eta_1 \ge \lambda_1 > 0$. Hence $S_0$ is unstable. $\square$

**Proposition 2** $S_1 = (1,0)$ is unstable.

*Proof.*

From (10), the linearized system of (4) around $S_1$ is

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial w_1}{\partial t} = \Delta w_1 - w_1 - \frac{a}{1+e_1} w_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_2}{\partial t} = \delta \Delta w_2 + bw_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (15)$$

Now we study the following eigenvalue problem:

$$\left\{ \begin{array}{l} \Delta w_1 - w_1 - \frac{a}{1+e_1} w_2 = \eta w_1, \quad x \in \Omega \\ \delta \Delta w_2 + bw_2 = \eta w_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (16)$$
We need to prove that the largest eigenvalue of system (16) is positive. First, same as in Proposition 1, all eigenvalues of system (16) are real. Let $\eta_1$ be the largest eigenvalue of system (16). Since $b > 0$, the principal eigenvalue $\lambda_1$ of

$$ \left\{ \begin{array}{l} \delta\Delta w_2 + bw_2 = \lambda w_2, \quad x \in \Omega \\ \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (17) $$

is positive and the associated eigenfunction $\tilde{w}_2 > 0$. Let us prove that $\lambda_1$ is also an eigenvalue of system (16). In fact, we take $\tilde{w}_1$ to be the solution of linear problem

$$ \left\{ \begin{array}{l} \Delta w_1 - (1 + \lambda_1)w_1 = \frac{a}{1 + e_1} \tilde{w}_2, \quad x \in \Omega \\ \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (18) $$

then $(w_1, w_2) = (\tilde{w}_1, \tilde{w}_2)$ is a solution of problem (16) with $\eta = \lambda_1$. So $\lambda_1 > 0$ is an eigenvalue of problem (16). Therefore we must have $\eta_1 \ge \lambda_1 > 0$. Hence $S_1$ is unstable.

<figure><img src="image_1.png" /></figure>

**Proposition 3**

If $e_1 > ae_2$ then $S_2 = (0, e_2)$ is unstable. If $e_1 < ae_2$ then $S_2 = (0, e_2)$ is stable.

*Proof.*

From equation (10), the linearized system of equation (4) around $S_2$ is

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial w_1}{\partial t} = \Delta w_1 + \left(1 - \frac{ae_2}{e_1}\right)w_1, \quad x \in \Omega \\ \displaystyle \frac{\partial w_2}{\partial t} = \delta \Delta w_2 + bw_1 - bw_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (19) $$

Now we study the following eigenvalue problem :

$$ \left\{ \begin{array}{l} \Delta w_1 + \left(1 - \frac{ae_2}{e_1}\right)w_1 = \eta w_1, \quad x \in \Omega \\ \delta \Delta w_2 + bw_1 - bw_2 = \eta w_2, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (20) $$

We need to prove that the largest eigenvalue of system (20) is positive if $e_1 > ae_2$. First, same as before, all eigenvalues of problem (20) are real. Let $\eta_1$ be the largest eigenvalue of problem (20). Since $e_1 > ae_2$, the principal eigenvalue $\lambda_1$ of

$$ \left\{ \begin{array}{l} \Delta w_1 + \left(1 - \frac{ae_2}{e_1}\right)w_1 = \lambda w_1, \quad x \in \Omega \\ \displaystyle \frac{\partial w_1}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (21) $$

is positive and the associated eigenfunction $\tilde{w}_1 > 0$. Let us prove that $\lambda_1$ is also an eigenvalue of problem (20). In fact, we take $\tilde{w}_2 > 0$ to be the solution of linear problem

$$ \left\{ \begin{array}{l} \delta \Delta w_2 - (b + \lambda_1)w_2 = -b \tilde{w}_1, \quad x \in \Omega \\ \displaystyle \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (22) $$
then $(w_1, w_2) = (\tilde{w}_1, \tilde{w}_2)$ satisfies problem (20) with $\eta = \lambda_1$. So $\lambda_1 > 0$ is an eigenvalue of problem (20). Therefore we must have $\eta_1 \ge \lambda_1 > 0$. Hence $S_2$ is unstable.

Let $(\tilde{w}_1, \tilde{w}_2)$ be the principal eigenfunction of problem (20) corresponding to the largest eigenvalue $\eta_1$. If $\tilde{w}_1 \ne 0$, then $\eta_1$ is also an eigenvalue of problem (21). Then we must have $\eta_1 < 0$ if $e_1 < ae_2$ because, in this case, the largest eigenvalue of problem (21) is $\lambda_1 = 1 - \frac{ae_2}{e_1} < 0$.

If $\tilde{w}_1 \equiv 0$, then we have $\tilde{w}_2 \ne 0$. Therefore $\eta_1$ is an eigenvalue of

$$ \left\{ \begin{array}{l} \delta \Delta w_2 - b w_2 = \lambda w_2, \quad x \in \Omega \\ \frac{\partial w_2}{\partial \nu} \Big|_{\partial \Omega} = 0 \end{array} \right. \qquad (23) $$

Obviously the largest eigenvalue of problem (23) is $-b < 0$. Therefore we also have $\eta_1 < 0$. Thus we know that if $\lambda_1 = 1 - \frac{ae_2}{e_1} < 0$, then $S_2 = (0, e_2)$ is stable. $\square$

**Proposition 4** Assume that $a \ge \frac{1}{2}$ and $0 < e_1 < \bar{e}_1$ with

$$ \bar{e}_1 = -(a+1) + \sqrt{(a+1)^2 + 2a(1+2a) - 1}. $$

Then $(u^*, v^*)$ is stable.

*Proof.* Let $\phi_j$ denote the j-th eigenfunction of $-\Delta$ on $\Omega$ with homogenous Neumann boundary condition. That is,

$$ \left\{ \begin{array}{l} -\Delta\phi_j = \lambda_j\phi_j, \quad \text{in } \Omega \\ \frac{\partial\phi_j}{\partial\nu} = 0, \quad \text{on } \partial\Omega, \end{array} \right. \qquad (24) $$

for scalar $\lambda_j$ satisfying $0 = \lambda_0 < \lambda_1 < \lambda_2...$

From (10), the linearized system of (4) around $(u^*, v^*)$ is

$$ \frac{\partial W}{\partial t} = D \Delta W + \Sigma W, \qquad (25) $$

where $D = \text{diag}(1, \delta)$, and

$$ \Sigma = \begin{pmatrix} A & B \\ Q & R \end{pmatrix} = \begin{pmatrix} 1 - 2u^* - \frac{e_1(1-u^*)}{u^*+e_1} & -\frac{au^*}{u^*+e_1} \\ b & -b \end{pmatrix}. \qquad (26) $$

We expand the solution $W$ of (25) via

$$ W = \sum_{j=0}^{\infty} z_j(t) \phi_j(x), \qquad (27) $$

where each $z_j(t) \in \mathbb{R}^2$. Substituting (27) into (25) and equating the coefficients of each $\phi_j$, we have

$$ \frac{dz_j}{dt} = C_j z_j, $$

where $C_j$ is the matrix

$$ C_j = \Sigma - \lambda_j D. $$
Now the solution $(u^*, v^*)$ is stable if and only if each $z_j(t)$ decays to zero. This is equivalent to the condition that each $C_j$ has two eigenvalues with negative real parts. The eigenvalues $\eta_{1,2}$ of $C_j$ are determined by

$$ \eta^2 - \eta[A + R - \lambda_j(1 + \delta)] + \lambda_j^2\delta - \lambda_j(R + \delta A) + AR - BQ = 0. $$

Therefore the fact that each $C_j$ has two eigenvalues with negative real parts is guaranteed by

$$ A + R - \lambda_j(1 + \delta) < 0, \qquad (28) $$

and

$$ \lambda_j^2\delta - \lambda_j(R + \delta A) + AR - BQ > 0. \qquad (29) $$

We have $R = -b < 0$. Therefore, observing that $\lambda_j \ge 0$ and $B < 0$, (28) and (29) hold if

$$ A \le 0, \text{ and } Q > 0. $$

Therefore to have $A \le 0$, we need

$$ 2u^* - 1 + e_1 \ge 0 $$

Let us recall that, $u^* = \frac{1}{2}(1 - a - e_1 + \Delta^{\frac{1}{2}})$, where $\Delta = (a + e_1 - 1)^2 - 4(ae_2 - e_1)$.

Consequently, we need that $1 - a - e_1 + \Delta^{\frac{1}{2}} - 1 + e_1 \ge 0$

Thus,

$$ (a + e_1 - 1)^2 - 4(ae_2 - e_1) - a^2 \ge a $$

$$ e_1^2 + (2a + 2)e_1 - 2a(1 + 2e_2) + 1 \ge 0. $$

This is true from the assumptions of the proposition, hence $A \le 0$. Obviously $Q > 0$, and this completes the proof. □

# 3. Complex pattern formation and spatiotemporal chaos

## 3.1. Local bifurcation in a one dimensional space

In this section, we present our numerical results in one dimensional space. We suppose that the two species diffuse on a line, $\Omega \subset \mathbb{R}$.

At boundaries we use the zero-flux condition. Let us consider the two following initial conditions :

$$ \begin{aligned} u(0, x) &= u_0 \quad \text{for } L_{1u} < x < L_{2u}, \quad \text{otherwise } u(0, x) = 0 \\ v(0, x) &= v_0 \quad \text{for } L_{1v} < x < L_{2v}, \quad \text{otherwise } v(0, x) = 0 \end{aligned} \qquad (30) $$

The initial domain, where the prey moves, is larger than that of the predator for making, during the simulation, the impact of the boundaries as small as possible. Thus, we assume,

$$ 0 < L_{1u} \le L_{1v} < L_{2v} \le L_{2u} < L $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** System (4) density of spatial distribution, the parameters and initial data are fixed as given in (31) and (32) and $b = 0.256$</figcaption></figure>

We choose the parameters so that the populations do not disappear when the growth or the decrease degree of prey, (that is (birth quantity of prey)/(death quantity of prey)) and the predator vary. In the following of this section, the parameters are fixed as follows,

$$L = 100, L_{1u} = 40, L_{1v} = 48, L_{2v} = 56, L_{2u} = 60, u_0 = 1.0, v_0 = 0.1 \quad (31)$$

$$e_1 = 0.08, e_2 = 0.01, a = 3.0, \delta = 1. \quad (32)$$

Figure 1 is an example of species spatial distribution, observed for $b = 0.256$, at (a) $t = 250$, (b) $t = 750$ and (c) $t = 1200$. With these fixed parameters, and these initial distributions, there are two patches at the beginning and the end of the field which are formed at the first moments of simulation. Between these patches the densities remain constant with respect to the time parameter.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** A cascade of bifurcations leading to the onset of chaotic oscillations in the phase plane (U,V) for different values of b : (a) $b = 0.197$, (b) $b = 0.203$, (c) $b = 0.23$, (d) $b = 0.26$. The other parameters are given in (31) and (32)</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Bifurcation diagram when the parameter b varies.</figcaption></figure>

To explore the properties of the population dynamics as a whole we estimate the species size of prey and predator by,

$$U(t) = \int_{0}^{L} u(t, x)dx \text{ and } V(t) = \int_{0}^{L} v(t, x)dx \quad (33)$$

The aim is to study the properties of the oscillations of the dynamics of the populations when one varies the control parameter, the choice of this parameter is then important. We will choose b as control parameter since it determines the ratio of two factors which are the birth rates of the prey and the predator. Therefore, while b varies between 0.195 and 0.26, other parameters will be fixed as in equation (32) and initial conditions are given as in (31). For each value of b system (4) is solved with the initial conditions given in equation (32). We leave a rather large transitory time, so that the total quantities of the species U and V are in the attractor domain. We start with b = 0.26 and make b decreasing. For b = 0.26, the system presents an attractor focus in (U, V), see figure 2(d). The same phase plane is obtained as long as b is higher than 0.255. We have a first bifurcation when this ratio is equal to 0.255. When b belongs to [0.208, 0.255] the system exhibits periodic attractors, see figure 2(c). A second bifurcation leads to the dynamics of the species in quasi-periodic attractors, for b between 0.208 and 0.199, see figure 2(b). Finally for b between 0.199 and 0.195, it becomes chaotic, see figure 2(a). These results are summarized by the bifurcation diagram given in figure 3.

## 3.2. Complex pattern formation in two dimensional space

Biological systems are in general far-from-equilibrium systems. Self-organization can operate via symmetry-breaking instabilities. Structures such as target patterns and spiral waves are observed in a wide set of chemical and biological system, see [12, 13]. For model (4), we observe spiral waves, see figure 4. Another kind of spatial pattern is spatial chaos, in which concentration waves change aperiodically through space, as in turbulence. These structures and spatiotemporal phenomena are emergent and arise from very simple causes, in the sense that they are macroscopic patterns emerging from simple coupling interaction and local diffusion rules. The last has huge implications for our understanding of complexity. In other words, they are not programmed in the equations structure, see
[5, 10, 14, 13].

In this last subsection we study the problem in the limited field $\Omega = [0, 900] \times [0, 900]$ of $\mathbb{R}^2$. We are interested in the emerging structures when the homogeneous equilibrium $(u^*, v^*)$ of system (4) is unstable and the species diffuse in the same way. Let us remark

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Spatial spiral waves type distribution of species at $t = 1200$ for respectively the prey and the predator. Initial conditions are given in (35). $b = 0.042$, $e_2 = 0.2$, the other Parameters are given in (36)</figcaption></figure>

that, if $(u^*, v^*)$ is unstable for equation (3) it becomes also unstable for (4) as soon as the coefficient of diffusion is equal to one. In the case of Turing instability, developed in ([13, 17]), all the eigenvalues of $(f(u, v), g(u, v))$ Jacobian matrix at $(u^*, v^*)$ have negative real parts. A necessary condition to observe a manifestation of instability is that predator must diffuse faster than prey. In this part, we suppose that the prey and the predator diffuse in the same way ($\delta = 1$). The global emerging structures are given only by the local interactions of the functional response $(f(u, v), g(u, v))$. We start from an initial condition rather close to $(u^*, v^*)$, having a low disparity of the space distribution. These initial conditions have been chosen as,

$$u(0, x, y) = u^* - 2 \cdot 10^{-7} (x - 0.1y - 231)(x - 0.1y - 632) \quad (34)$$

$$v(0, x, y) = v^* - 3 \cdot 10^{-5} (x - 450) - 1.2 \cdot 10^{-4} (y - 150) \quad (35)$$

we now set the parameters so that the positive homogeneous solution is unstable,

$$e_1 = 0.3, e_2 = 0.1, b = 0.02, a = 1.1, \delta = 1 \quad (36)$$

We observe the following time evolution of spatial distributions. The left figures are the evolution of the prey spatial distribution and the right are the predator's.

<figure><img src="image_6.png" /><figcaption>Spatial distribution of species at $t = 400$</figcaption></figure>
<figure><img src="image_7.png" /></figure>

$t = 800$

<figure><img src="image_8.png" /></figure>

$t = 3000$

<figure><img src="image_9.png" /></figure>

$t = 10000$

**Acknowledgments** I would like to say my deepest thanks to the referees of this article.

## 4. References

[1] AMANN H., "Dynamics theory of quasilinear parabolic equations-I. Abstract evolution equations", *Nonlinear Anal*, vol. 12, (1988), 895-919.

[2] AMANN H., "Dynamics theory of quasilinear parabolic equations-II. Reaction-diffusion system", *Diff. Int. Eq.* vol. 3, (1990), 13-75.

[3] AMANN H., "Dynamics theory of quasilinear parabolic equations-III. Global Existence", *Math. Z.*, vol. 202, (1989), 219-250.

[4] AZIZ-ALAOUI M.A., DAHER OKIYE M., "Boundedness and global stability for a predator-prey model with modified Leslie-Gower and Holling type II themes", *Applied Math. Lets.* vol. 16, (2003), 1069-1075.

[5] BASCOMPTE J., SOLÉ R. V., "Rethinking and complexity : modelling spatiotemporal dynamics in ecology", *Trends in Ecology and Evolution*, vol. 10, (1995), 361-366.

[6] CAMARA B. I., AZIZ-ALAOUI M.A., "Dynamics of predator-prey model with diffusion", *To appear in Dynamics of Continuous, Discrete and Impulsive System, series A*, (2008)

[8] DAHER OKIYE M., AZIZ-ALAOUI M.A., "On the dynamics of a predator-prey model with the Holling-Tanner functional response", *MIRIAM Editions, Editor V. Capasso, Proc. ESMTB conf.* vol. (2002), 270-278.
[9] DAHER OKIYE M., "Study and asymptotic analysis of some nonLinear dynamical systems : Application to predator-prey problems", *PHD Thesis, Le Havre University, France*, (2004).

[10] KAWATA M., TOQUENAGA Y., "From artificial individuals to global patterns", *Trends Ecol. Evol.* vol. 9, (1994), 417-421.

[11] KIELHÖFER H., "Stability and semilinear evolution equations in Hilbert space", *Arch. Rational Mech. Anal.*, vol. 57, (1974), 150-165.

[12] MAINI P.K., PAINTER K.J., CHAU H.N.P., "Spiral pattern formation in chemical and biological systems", *Journal of Chim. Soc.*, vol. 93, num. 20, (1997) 3601-3610.

[13] MURRAY J. D., "Mathematical Biology: II. Spatial Models and Biomedical Applications", Springer-Verlag Berlin Heidelberg., 2003

[14] MURRAY J. D., "Mathematical Biology", Springer Verlag, 1993

[15] NINDJIN A.F., AZIZ-ALAOUI M.A., CADIVEL M., "Analysis of a Predator-Prey Model with Modified Leslie-Gower and Holling-Type II Schemes with Time Delay", *NonLinear Analysis Real World Applications*, vol. 7, num. 5, (2006) 1104-1118.

[16] NINDJIN A. K., AZIZ-ALAOUI M.A., "Persistence and global stability in a delayed Leslie-Gower type three species food chain", *Journal of Mathematical Analysis and Applications*, vol. 340, num. 1, (2008), 340-357.

[17] TURING A., "The Chemical Basis of Morphogenesis", Phil. Trans. R. Soc. London B, vol. 237, (1952), 37-72.