Stationary and non-stationary pattern formation over fragmented habitat

We use Routh-Hurwitz criteria to determine the local asymptotic stability of coexisting equilibrium point [22]. The Jacobian matrix of the system (4) evaluated at $E_{*}$ is given by

$\displaystyle J_{*}\,=\,\left[\begin{array}[]{cc}g(u_{*})+u_{*}g^{\prime}(u_{*})-h^{\prime}(u_{*})v_{*}&-h(u_{*})\\ sh^{\prime}(u_{*})v_{*}&-spv_{*}\\ \end{array}\right]\,\equiv\,\left[\begin{array}[]{cc}\alpha_{11}&\alpha_{12}\\ \alpha_{21}&\alpha_{22}\\ \end{array}\right].$

(9)

According to the Routh-Hurwitz criteria [23], $E_{*}$ is locally asymptotically stable if the following two conditions are satisfied:

$\displaystyle\textrm{tr}(J_{*})\,=\,\alpha_{11}+\alpha_{22}\,<\,0,\,\,\textrm{det}(J_{*})\,=\,\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21}\,>\,0.$

(10)

Components of $E_{*}$ are positive, parameters are positive and $h^{\prime}(u)>0$ for all $u>0$. Except $\alpha_{11}$, signs of other entries of $J_{*}$ are fixed and hence the sufficient condition for local asymptotic stability of $E_{*}$ is given by $\alpha_{11}<0$.

The coexisting steady-state $E_{*}$ loses stability through Hopf bifurcation whenever $\alpha_{11}>0$. The Hopf bifurcation condition for the coexistence equilibrium $E_{*}$ is given by

$\displaystyle\textrm{tr}(J_{*})\,=\,0,\,\,\textrm{det}(J_{*})\,>\,0,\,\,\frac{d}{ds}\left[\textrm{tr}(J_{*})\right]\,\neq\,0.$

(11)

It is interesting to mention that the components of $E_{*}$ can not be found explicitly, however, the model formulation allow us to obtain the Hopf-bifurcation threshold explicitly in terms of $s$ as:

$\displaystyle s_{H}=\frac{g(u_{*})+u_{*}g^{\prime}(u_{*})-h^{\prime}(u_{*})v_{*}}{pv_{*}}.$

(12)

Note that the Hopf bifurcation threshold is feasible when $\alpha_{11}>0$ and $\alpha_{11}\alpha_{22}>\alpha_{12}\alpha_{21}$. One can find a limit cycle in the vicinity of Hopf bifurcation threshold, and the stability of Hopf bifurcating limit cycle depends on the sign of the first Lyapunov number [22, 24]. To obtain the expression for the first Lyapunov number, we first expand the Taylor series expansion of the system (4) around $E_{*}$, $s=s_{H}$ as the following

	$\displaystyle\frac{du_{1}}{dt}$	$\displaystyle=$	$\displaystyle\beta_{10}u_{1}+\beta_{01}v_{1}+\beta_{20}u_{1}^{2}+\beta_{11}u_{1}v_{1}+\beta_{02}v_{2}^{2}+\beta_{30}u_{1}^{3}+\beta_{21}u_{1}^{2}v_{1}+\beta_{12}u_{1}v_{1}^{2}+\beta_{03}v_{2}^{3}+\cdots,$		(13a)
	$\displaystyle\frac{dv_{1}}{dt}$	$\displaystyle=$	$\displaystyle\gamma_{10}u_{1}+\gamma_{01}v_{1}+\gamma_{20}u_{1}^{2}+\gamma_{11}u_{1}v_{1}+\gamma_{02}v_{2}^{2}+\gamma_{30}u_{1}^{3}+\gamma_{21}u_{1}^{2}v_{1}+\gamma_{12}u_{1}v_{1}^{2}+\gamma_{03}v_{2}^{3}+\cdots.$		(13b)

The expression for the first Lyapunov number is given by

	$\displaystyle\sigma$	$\displaystyle=$	$\displaystyle-\frac{3\pi}{2\beta_{01}\Delta^{3/2}}\left[\beta_{10}\gamma_{10}(\beta_{11}^{2}+\beta_{11}\gamma_{02}+\beta_{02}\gamma_{11})+\beta_{10}\beta_{01}(\gamma_{11}^{2}+\beta_{20}\gamma_{11}+\beta_{11}\gamma_{02})\right.$		(14)
			$\displaystyle+\gamma_{10}^{2}(\beta_{11}\beta_{02}+2\beta_{02}\gamma_{02})-2\beta_{10}\gamma_{10}(\gamma_{02}^{2}-\beta_{02}\beta_{20})-2\beta_{10}\beta_{01}(\beta_{20}^{2}-\gamma_{20}\gamma_{02})$
			$\displaystyle-\beta_{01}^{2}(2\beta_{20}\gamma_{20}+\gamma_{11}\gamma_{20})+(\beta_{01}\gamma_{10}-2\beta_{10}^{2})(\gamma_{11}\gamma_{02}-\beta_{11}\beta_{20})$
			$\displaystyle\left.-(\beta_{10}^{2}+\beta_{01}\gamma_{10})(3(\gamma_{10}\gamma_{03}-\beta_{01}\beta_{30})+2\beta_{10}(\beta_{21}+\gamma_{12})+(\gamma_{10}\beta_{12}-\beta_{01}\gamma_{21}))\right],$

where $\Delta=\beta_{10}\gamma_{01}-\beta_{01}\gamma_{10}$. The Hopf bifurcation is super-critial if $\sigma<0$ and is sub-critical for $\sigma>0$. We verify the stability of Hopf bifurcating limit cycle with the help of numerical example.

Here we derive the Turing instability condition for the system (1) around the homogeneous steady-state $E_{*}$, where $u(t,x,y)=u_{*}$, $v(t,x,y)=v_{*}$ is the homogeneous steady-state for (1). The model (1) is subject to the non-negative initial condition

$$u(0,x,y)\,=\,u_{0}(x,y)\,\geq\,0,\,\,v(0,x,y)\,=\,v_{0}(x,y)\,\geq\,0,\,\,0\,\leq\,x,y\,\leq\,L,$$

(15)

and zero-flux boundary condition along the boundary of the square domain.

After linearizing the system (1) around the homogeneous steady-state $E*$, we obtain the following system of linear PDEs in terms of the perturbation variables $u_{1}(t,x,y)$ and $v_{1}(t,x,y)$:

	$\displaystyle\dfrac{\partial u_{1}}{\partial t}$	$\displaystyle=d_{1}\bigg{(}\dfrac{\partial^{2}u_{1}}{\partial x^{2}}+\dfrac{\partial^{2}u_{1}}{\partial y^{2}}\bigg{)}+\alpha_{11}u_{1}+\alpha_{12}v_{1},$		(16)
	$\displaystyle\dfrac{\partial v_{1}}{\partial t}$	$\displaystyle=d_{2}\bigg{(}\dfrac{\partial^{2}v_{1}}{\partial x^{2}}+\dfrac{\partial^{2}v_{1}}{\partial y^{2}}\bigg{)}+\alpha_{21}u_{1}+\alpha_{22}v_{1}.$

We assume that the solution of the system of linear equations (16) satisfying the no-flux boundary conditions and the solution is in the form

$\displaystyle u_{1}(t,x,y)\,=\,\epsilon_{1}e^{\lambda t}\cos(k_{1}x)\cos(k_{2}y),\,\,v_{1}(t,x,y)\,=\,\epsilon_{2}e^{\lambda t}\cos(k_{1}x)\cos(k_{2}y),$

(17)

where $\epsilon_{1},\epsilon_{2}\ll 1$ are two arbitrary constants and $k=\sqrt{k_{1}^{2}+k_{2}^{2}}$ is wave number with $k_{1}=n_{1}\pi/L$, $k_{2}=n_{2}\pi/L$, $n_{1}$ and $n_{2}$ are two natural numbers. For a non-trivial solution (17) of the linear system of equations (16), $\lambda$ satisfies the following characteristics equation

$$|M(k)-\lambda I_{2}|\,=\,0,$$

(18)

where

$$M(k)=\begin{pmatrix}\alpha_{11}-d_{1}k^{2}&\alpha_{12}\\ \alpha_{21}&\alpha_{22}-d_{2}k^{2}\end{pmatrix},$$

(19)

$\alpha_{ij}$’s $(i,j=1,2)$ are given in (9). After expanding the characteristic equation (18), we obtain

$$\lambda^{2}-((\alpha_{11}+\alpha_{22})-(d_{1}+d_{2})k^{2})\lambda+h(k^{2})=0,$$

(20)

where

$$h(k^{2})=d_{1}d_{2}k^{4}-(d_{2}\alpha_{11}+d_{1}\alpha_{22})k^{2}+\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21}.$$

(21)

The homogeneous steady-state $E_{*}$ is stable under small amplitude heterogeneous perturbation if both the eigenvalues of the characteristic equation (20) has a negative real part. Turing instability sets in due to the destabilization of homogeneous steady-state $E_{*}$, when one root of the characteristic equation (20) passes through zero at some critical wave number. At the same time, the other root of the characteristic equation remains negative. The critical wave number $k_{c}$ can obtain by solving the equation $\frac{dh(k^{2})}{dk^{2}}=0$ for $k^{2}$, which is given by

$$k^{2}_{c}=\dfrac{1}{2d_{1}d_{2}}(d_{2}\alpha_{11}+d_{1}\alpha_{22}).$$

(22)

As the critical wave number $k_{c}$ is a real number, the dimensionless diffusion coefficients $d_{1},d_{2}>0$, feasible existence of $k_{c}$ demands the satisfaction of the implicit parametric restriction $d_{2}\alpha_{11}+d_{1}\alpha_{22}>0$. Turing instability occurs when a stable homogeneous steady-state (which is stable under small amplitude homogeneous perturbation) becomes unstable under small amplitude heterogeneous perturbation. We first assume that the homogeneous steady-state $E_{*}$ is locally asymptotically stable, i.e., the conditions (10) are satisfied. Hence, for positive $d_{1}$ and $d_{2}$, $\alpha_{11}+\alpha_{22}<0$ and $d_{2}\alpha_{11}+d_{1}\alpha_{22}>0$ satisfied simultaneously when $\alpha_{11}\alpha_{22}<0$. The equation of Turing bifurcation boundary can be obtained by eliminating $k^{2}$ from $h(k^{2})=0$ with the help of $dh(k^{2})/dk^{2}=0$ as follows,

$$d_{2}\alpha_{11}+d_{1}\alpha_{22}\,=\,2\sqrt{d_{1}d_{2}}\sqrt{\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21}}.$$

(23)

The Turing bifurcation curve is feasible under the implicit parametric restrictions $\alpha_{11}\alpha_{22}<0$ and $\alpha_{11}\alpha_{22}>\alpha_{12}\alpha_{21}$. These two conditions can be written in terms of $M(k)$ as follows

$$\textrm{tr}(M(0))\,<\,0,\,\,\textrm{det}(M(0))\,>\,0.$$

(24)

The wavenumber $k$ is given by $k=\sqrt{\dfrac{n_{1}^{2}\pi^{2}}{L^{2}}+\dfrac{n_{2}^{2}\pi^{2}}{L^{2}}}$, where $L$ is the length of the square domain, $n_{1}$ and $n_{2}$ are two natural numbers. We use the notation $\omega=\pi/L$. After substituting $k^{2}=(n_{1}^{2}+n_{2}^{2})\omega^{2}$ in $h(k^{2})=0$ and then solving for $d_{2}$, we find the explicit expression for Turing bifurcation boundary corresponding to the admissible unstable eigenmodes $(n_{1},n_{2})$ as

$\displaystyle d_{2T}(n_{1},n_{2})\,=\,\dfrac{(n_{1}^{2}+n_{2}^{2})\omega^{2}d_{1}\alpha_{22}-\alpha_{11}\alpha_{22}+\alpha_{12}\alpha_{21}}{d_{1}(n_{1}^{2}+n_{2}^{2})\omega^{2}((n_{1}^{2}+n_{2}^{2})\omega^{2}-a_{11})}.$

(25)

For a given set of parameter values, suppose the homogeneous coexisting steady-state $E_{*}$ is stable under temporal perturbation (i.e. the condition (24) holds), then $E_{*}$ is stable under heterogeneous perturbation for $d_{2}<d_{2T}$ and is unstable for $d_{2}>d_{2T}$. Depending upon the chosen parameter values involved with the reaction kinetics, one can find more than one set of values for $(n_{1},n_{2})$ such that $d_{2T}(n_{1},n_{2})$ is feasible. Consequently, for a fixed set of parameter values, we find multiple feasible thresholds $d_{2T}(n_{1},n_{2})$, which indicates the existence of multiple heterogeneous stationary solutions. However, it is not easy to determine the number of heterogeneous stationary solutions. We obtain the pattern through numerical simulation is determined by the feasible value of $d_{2}>d_{2T}(n_{1},n_{2})$ and the admissible value of $(n_{1},n_{2})$ corresponding to the maximum real part of $\lambda(k)$ along with the size of the domain and the choice of the initial condition. We will explain this idea in detail with the help of a numerical example in the following section.

The critical Turing bifurcation threshold, in terms of $d_{2}$, is given by a positive solution of the equation (23) when other parameter values are given, and let us denote this solution as $d_{2T}$. Now if we choose a value of $d_{2}>d_{2T}$, keeping other parameters fixed, then we can find a range of values of $k$, say $(k_{1},k_{2})$, such that $h(k^{2})<0$ for $k_{1}<k<k_{2}$, where $k_{1}$ and $k_{2}$ are given by

	$\displaystyle k_{1}^{2}$	$\displaystyle=$	$\displaystyle\frac{d_{2}\alpha_{11}+d_{1}\alpha_{22}-\sqrt{(d_{2}\alpha_{11}+d_{1}\alpha_{22})^{2}-4d_{1}d_{2}(\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21})}}{2d_{1}d_{2}},$		(26a)
	$\displaystyle k_{2}^{2}$	$\displaystyle=$	$\displaystyle\frac{d_{2}\alpha_{11}+d_{1}\alpha_{22}+\sqrt{(d_{2}\alpha_{11}+d_{1}\alpha_{22})^{2}-4d_{1}d_{2}(\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21})}}{2d_{1}d_{2}}.$		(26b)

For $d_{2}>d_{2T}$, one can verify that the condition $d_{2}\alpha_{11}+d_{1}\alpha_{22}>2\sqrt{d_{1}d_{2}}(\alpha_{11}\alpha_{22}-\alpha_{12}\alpha_{21})$ holds, and the positivity of $d_{2}\alpha_{11}+d_{1}\alpha_{22}$ ensures the feasibility of $k_{1}$ and $k_{2}$. Depending upon the choice of parameter values and the size of the domain, we can find more than one combination of $(n_{1},n_{2})$ for which the following result holds

$\displaystyle k_{1}<\frac{\pi}{L}\sqrt{n_{1}^{2}+n_{2}^{2}}<k_{2}.$

(27)

For a given set of parameter values, we can find a different admissible set of values for $(n_{1},n_{2})$ such that the above inequality is satisfied. But, the resulting stationary pattern is the combination of $(n_{1},n_{2})$ which corresponds to the maximum value of the real part of $\lambda(k)$ which is a root of the characteristic equation (20).

Here we consider the pattern formation by the model (1) over a square domain $[0,L]\times[0,L]$, subject to zero-flux boundary condition and non-negative initial condition

$\displaystyle u(0,x,y)=u_{*}+\epsilon\xi(x,y),\,\,v(0,x,y)=v_{*}+\epsilon\eta(x,y),$

(28)

where $\xi(x,y)$ and $\eta(x,y)$ are two Gaussian white noises $\delta$-correlated in space [16, 17] and $\epsilon$ measures the amplitude of small heterogeneous perturbation from the homogeneous steady-state. As mentioned above, henceforth, we mention only the values of $s$, $a$, $d_{1}$ and $d_{2}$ to obtain various stationary and dynamic patterns. We present three different stationary patterns in Fig. 2. For stationary patterns we fix $s=3$, $d_{1}=0.15$, $d_{2}=10$, we find hot spot pattern for $a=1.65$, labyrinthine pattern for $a=2$ and cold spot pattern for $a=2.2$. We choose the length of the square domain as $L=200$ for the numerical solution of stationary patterns. It is important to mention here that the obtaining stationary patterns presented in Fig. 2 are at a different time as the transient pattern persists for different duration of time. Three different stationary patterns are obtained at $t=500$ (hot spot), $t=1800$ (labyrinthine) and $t=2500$ (cold spot) respectively. For these choices of parameter values, homogeneous coexistence steady-state is stable, and the magnitude of diffusion coefficients satisfy the Turing instability condition. Hence, the resulting patterns are Turing patterns and are stationary. However, the duration of the transient pattern depends on the strength of temporal interactions and the magnitude of the diffusion parameters.

The dispersion relations corresponding to three stationary Turing patterns presented in Fig. 2, are presented in Fig. 3. For $a=1.65$, the largest real part of the eigenvalues of the characteristic equation (20) is positive for $k\in(k_{1},k_{2})=(0.49,1.21)$ and is maximum at $k=0.773$. There are plenty of choices for $n_{1}$ and $n_{2}$ such that the inequality $0.49<\frac{\pi}{L}\sqrt{n_{1}^{2}+n_{2}^{2}}<1.21$ holds, for $L=200$. In particular, the value $k=0.773$, corresponding to the maximum real part of the eigenvalue, can be achieved from the expression $\frac{\pi}{L}\sqrt{n_{1}^{2}+n_{2}^{2}}$ for a couple of choices of $n_{1}$, $n_{2}$, namely $(n_{1},n_{2})=(35,35)$, $(n_{1},n_{2})=(34,36)$, $(n_{1},n_{2})=(36,34)$ and so on. The choice of $(n_{1},n_{2})$ for which $\frac{\pi}{200}\sqrt{n_{1}^{2}+n_{2}^{2}}$ is close to $k=0.773$ is $(n_{1},n_{2})=(35,35)$. The resulting pattern presented in Fig. 2(2(a)) shows equal number of peaks within the domain $[200\times 200]$ as the number of peaks we find for the function $\cos\frac{n_{1}\pi x}{200}\cos\frac{n_{2}\pi x}{200}$ when plotted over a same domain size. Similar argument holds for the patterns presented in other panels of Fig. 2.

We present some critical curves for the Turing instability of the spatially homogeneous steady state in Fig. 4. We plot these curves in thin black. For each $(n_{1},n_{2})$, we find these curves by solving the equation (21) for $d_{2}$ with varying $a$. Specifically, two critical curves are marked (thick magenta curves) with the mode numbers, while the rest are not (thin black curves). For $a=1.65$, the Turing bifurcation threshold lies close to the critical curve mode $(19,55)$. It is the only exciting mode that determines the spatial pattern of the emerging wave dynamics. For $a=1.65$, $d_{2}=4.986$, the largest real part of the eigenvalues of the characteristic equation (20) is positive for $k\in(k_{1},k_{2})=(0.911,0.921)$ and is maximum at $k=0.916$. Therefore, for $n_{1}=19$, $n_{2}=55$ and $L=200$, the inequality $0.911<\frac{\pi}{L}\sqrt{n_{1}^{2}+n_{2}^{2}}<0.921$ holds as $k=0.916$ corresponds to the maximum real part of the eigenvalue, and it is close to the value of the expression $\frac{\pi}{L}\sqrt{n_{1}^{2}+n_{2}^{2}}=0.914$ while $(n_{1},n_{2})=(19,55)$. Hence, $(19,55)$ is a reasonable critical mode and the corresponding mode curve is close to the Turing bifurcation threshold [25]. Similarly, the mode curve for $a=2.25$ and close to Turing bifurcation curve is determined by the mode numbers $(n_{1},n_{2})=(9,47)$.

The stationary patterns presented in Fig. 2 correspond to the parameter values from the pure Turing domain as well as from Turing-Hopf domain. From the dispersion relations, we can identify that $a=1.65$ corresponds to the parameter value in the pure Turing domain as the real part of $\lambda(k)$ is negative at $k=0$. On the other hand, for $a=2$ and $a=2.5$, the real parts of $\lambda(k)$ is positive when $k=0$. For parameter values within the Turing-Hopf domain, the real part of the eigenvalue for $k=0$ is less than the maximum real part of the eigenvalue for some $k\neq 0$, and hence we find stationary patterns. A consolidated pattern diagram is presented in Fig. 5 for a range of values of the parameters $a$ and $d_{2}$. Stationary heterogeneous patterns exist for parameter values above the Turing bifurcation curve (magenta color curve). Below the Turing bifurcation curve, homogeneous steady-state is stable, and the level of stable homogeneous steady-state varies with the change of the value of $a$ from 1.5 to 2.4.

Now we consider a non-Turing pattern for $d_{1}=d_{2}=1$ which is in fact a non-stationary pattern. For $a=1.0$, $s=2.0$, and other parameter values as mentioned earlier, we find spatio-temporal chaos as shown in Fig. 6(a). We obtain this spatio-temporal chaotic pattern through the numerical simulation over a square domain of size $L=200$ by taking the initial condition as mentioned in (28) and the zero-flux boundary condition. The chaotic pattern is shown for $t=1000$, the nature of the chaotic pattern is interacting spiral chaos. We have verified the chaotic nature of the resulting pattern by plotting the chaotic oscillation of time evolution of spatial averages of prey and predator population densities $\langle u\rangle$ and $\langle v\rangle$ in Fig. 6(b). This spatio-temporal chaotic pattern also exhibits sensitivity to the initial condition, but we omit that result for the sake of brevity. A detailed discussion about the verification of sensitivity of spatial pattern to initial condition for the spatio-temporal model is available in [18].

In this section, we consider the effect of the shape of the domain on spatio-temporal pattern formation. To understand the effect of the complex shape of the domain on the pattern formation, we consider the U-shaped domain as shown in Fig. 7. Here, the connection between two large habitats is like a channel where the individuals can move from one domain to another. The U-shaped domain is characterized by the length ($L_{2}$) and width ($L_{x_{1}}$, $L_{x_{3}}$) of two habitats and the length ($L_{x_{2}}$) and width ($L_{y}$) of the connecting channel. First, we consider the pattern formation over the U-shaped domain when the size of the two patches are the same, and the connecting channel is reasonably wide. Next, we consider how the length and width of the connecting channel and the difference in the size of the two patches affect the various aspects of spatio-temporal pattern formation.

For the numerical simulations, we consider the following initial condition:

$\displaystyle u(0,x,y)=\left\{\begin{array}[]{ll}u_{*}+\epsilon\xi(x,y)&(x,y)\in\mathcal{D}_{1},\\ u_{*}&\textrm{elsewhere}\\ \end{array}\right.\mbox{and}\leavevmode\nobreak\ v(0,x,y)=\left\{\begin{array}[]{ll}v_{*}+\epsilon\eta(x,y)&(x,y)\in\mathcal{D}_{1},\\ v_{*}&\textrm{elsewhere}\\ \end{array}\right.,$

(33)

where $\mathcal{D}_{1}=[0,L_{x_{1}}]\times[0,L_{2}]$ and $\xi(x,y)$, $\eta(x,y)$ are spatially uncorrelated Gaussian white noise terms. The boundary condition is assumed to be no-flux along the boundaries of the entire domain. For our convenience of further discussion, we denote the domain consisting of the connecting corridor and the patch on the right-hand side as $\mathcal{D}_{2}$. It can be defined as $\mathcal{D}_{2}=\{[L_{x_{1}},L_{x_{1}}+L_{x_{2}}]\times[0,L_{y}]\}\cup\{[L_{x_{1}}+L_{x_{2}},L_{1}]\times[0,L_{2}]\}$.

Here, we first consider both the patches identical, and the parameter values are uniform over the entire U-shaped domain. Figure 8 depicts the pattern formation over the entire domain for different values of $a$, $s$, $d_{1}$ and $d_{2}$ keeping other parameters fixed as mentioned in the previous sub-section. We see that a similar type of stationary or dynamic pattern covers the entire U-shaped domain for fixed parameter values, as if the domain shape does not affect the resulting pattern. This observation is true if we perform the simulation over a reasonably long time and the width of the patches and the connecting channel is reasonably large. The characterization of the phrase ‘reasonably large’ is explained below in detail for the specific choices of $L_{x_{j}}$, $j=1,2,3$ and $L_{y}$. The domain size in Fig. 8 is $L_{2}=200$, $L_{x_{1}}=L_{x_{3}}=80$, $L_{x_{2}}=40$ and $L_{y}=40$. Depending on the choice of the parameters, the stationary and dynamic pattern engulf the entire U-shaped domain for the case of a reasonable large connecting channel, like the shape of the domain does not affect the pattern formation.

Now we explain how the size of the patches and connecting channel affect the spatio-temporal pattern formation. To explain this idea, we start with the parameter values $s=3$, $a=1.65$, $d_{1}=0.15$ and $d_{2}=10$. Specific values pertaining to the domain size are $L_{2}=200$, $L_{x_{1}}=L_{x_{3}}=80$, $L_{x_{2}}=40$ and $L_{y}=2$ and we use the initial condition (33). We present the pattern formation at different time steps in Fig. 9. First, we observe a hot-spot pattern in the domain $\mathcal{D}_{1}$ and an almost uniform distribution of the populations in the rest of the domain. Then a periodic travelling wave appears in the connecting channel and spread the population in the domain $\mathcal{D}_{2}$ where $\mathcal{D}_{2}=[L_{x_{1}}+L_{x_{2}},L_{1}]\times[0,L_{2}]$. The periodic travelling wave-like pattern breaks down once it hit the boundary, and then a stationary hot spot covers the entire U-shaped domain $\mathcal{D}_{2}$ [see Fig. 9(c)]. It is interesting to note that the formation of stationary pattern in $\mathcal{D}_{2}$ is regulated by the width of the connecting channel (i.e. $L_{y}$).

The initial wave propagation profile through the connecting channel to the domain $\mathcal{D}_{2}$ solely depend upon the width of the channel and the choice of the initial condition in $\mathcal{D}_{2}$. With the same choice of parameter values as mentioned in the above paragraph, we choose the initial condition in the domain $\mathcal{D}_{1}$ is the stationary hot-spot pattern as shown in earlier figures [cf. Fig. 9(c)]. We define two different initial conditions as follows:

$\displaystyle u(0,x,y)=\left\{\begin{array}[]{ll}\bar{u}(x,y)&(x,y)\in\mathcal{D}_{1},\\ u_{*}&(x,y)\in\mathcal{D}_{2}\\ \end{array}\right.,v(0,x,y)=\left\{\begin{array}[]{ll}\bar{v}(x,y)&(x,y)\in\mathcal{D}_{1},\\ v_{*}&(x,y)\in\mathcal{D}_{2}\\ \end{array}\right.,$

(38)

and

$\displaystyle u(0,x,y)=\left\{\begin{array}[]{ll}\bar{u}(x,y)&(x,y)\in\mathcal{D}_{1},\\ 0&(x,y)\in\mathcal{D}_{2}\\ \end{array}\right.,v(0,x,y)=\left\{\begin{array}[]{ll}\bar{v}(x,y)&(x,y)\in\mathcal{D}_{1},\\ 0&(x,y)\in\mathcal{D}_{2}\\ \end{array}\right.,$

(43)

where $\bar{u}(x,y)$ and $\bar{v}(x,y)$ is some stationary or dynamic pattern obtained through prior simulation. To be specific, we consider the stationary pattern for the prey in domain $\mathcal{D}_{1}$ as shown in Fig. 9(c) is $\bar{u}(x,y)$ and the corresponding pattern for the predator is chosen as $\bar{v}(x,y)$. We choose the U-shaped domain as $L_{2}=200$, $L_{x_{1}}=L_{x_{3}}=80$, $L_{x_{2}}=40$ and $L_{y}=25$. The simulation results at $t=500$ for two different initial conditions (38) and (43) are shown in Fig. 10 in the left and right panel, respectively. For a large time, a non-homogeneous stationary pattern occurs for both the initial conditions. But, the absence of the prey and predator species on the right side patch and the connecting channel (i.e. in $\mathcal{D}_{2}$) favours the rapid settlement of the stationary pattern in $\mathcal{D}_{2}$. We find significantly excess number of hot-spots in $\mathcal{D}_{2}$ corresponding to the initial condition (43) while compared to the pattern obtained through the initial condition (38). This is not only true for the two-dimensional case, but it is also true for the one-dimensional case. We choose a similar type of initial condition for the one-dimensional case. Here, the population distribution in the left half of the domain is the non-homogeneous stationary solution obtained through prior numerical simulation, and the right half is either a homogeneous steady-state (corresponding to the choice of parameter values) or zero population steady-state for both populations. We see that the zero population steady-state on the right half domain takes less time to settle the non-homogeneous stationary pattern over the entire domain than the homogeneous steady-state [see Fig. 11].

The width of the connecting channel and the size of the patches are responsible for the type of stationary pattern and the time required to reach the stationary distribution. The resulting stationary pattern is related to the domain size, as the domain size allows the heterogeneous distribution to settle down to the corresponding most unstable eigenmode. To illustrate this idea, we consider the parameter values which corresponds to the labyrinthine like pattern as shown in Fig. 8(b). In order to have the right-hand patch with smaller width compared to the patch on the left, we consider the domain with $L_{2}=200$, $L_{x_{1}}=80$, $L_{x_{2}}=90$, $L_{x_{3}}=30$ and $L_{y}=5$. Here, we choose the initial condition as (33) for the numerical simulation. We find a labyrinthine pattern in the left patch $\mathcal{D}_{1}$. However, a stripe pattern appears in the right patch [see Fig. 12]. We obtain this pattern at $t=1800$; this stationary pattern remains unaltered if we run the numerical simulations over a longer time.

In order to see the formation of spatio-temporal chaotic pattern over the U-shaped domain, we choose $a=1$, $s=2$, $d_{1}=d_{2}=1$. We take the initial condition (33) and domain size as $L_{1}=400$, $L_{2}=200$, $L_{x_{1}}=L_{x_{3}}=180$, $L_{x_{2}}=40$ and $L_{y}=2$. We plot the chaotic patterns at two different time steps in Fig. 13. As the width of the channel is narrow, the spatio-temporal chaotic pattern occurs on the right-hand patch after a considerable elapsed time interval. Initially, we observe a modulated periodic travelling wave-like pattern propagates, and then the chaotic spatial distribution appears later. Without any loss of generality, we can consider the modulated periodic travelling wave-like pattern in the right-hand patch as a transient pattern. If we look at the emerging pattern, it is challenging to realize the periodic nature of the transient pattern. However, we can interpret the transient pattern that appears at the right-hand patch from initial to some advancement of time as periodic with varying amplitude. For this, we plot the time evolution of the spatial average of prey densities in the domains $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$.

It is interesting to report that the time evolution of spatial average of the prey population density in $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ exhibit significantly different behaviour during some initial time interval. The time evolution of spatial averages of prey density in $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ are plotted in blue and magenta color, respectively. The numerical simulation is performed with the initial condition (33) and over a significantly large time interval $[0,1500]$. From Fig. 14, we can observe that the spatial average of prey density in $\mathcal{D}_{1}$ exhibits a short transient oscillation before entering the chaotic regime. However, the spatial average of prey density in $\mathcal{D}_{2}$ shows large-amplitude oscillation of varying magnitude for a considerable time before entering the chaotic regime. The high amplitude periodic oscillation of the spatial average sets in within a short time compared to the duration over which it declines. The width of the chaotic oscillation of the average prey density in $\mathcal{D}_{2}$ matches with the average prey density in $\mathcal{D}_{1}$ after a considerable time.

Finally, we like to remark that the duration of transient oscillation of average prey density in $\mathcal{D}_{2}$ depends on the width of the connecting channel. The time evolution of spatial averages of prey density in the domains $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ are presented in Fig. 14 for two different values of $L_{y}$. We have chosen $L_{y}=25$ and $L_{y}=50$, keeping other measures related to the domain size unaltered. With an increase in the magnitude of $L_{y}$, the duration of the periodic oscillation with varying amplitudes decreases.