Dynamics and spatio-temporal patterns in a prey–predator system with aposematic prey

1.   Introduction

The intra- and inter-species mutualistic and antagonistic interactions have a significant impact on the persistence and extinction of constituent populations in the real-world ecological systems. One of the possible ways by which we can quantitatively understand the effect of cooperative and competitive interactions on the co-evolution of species is the intensive analysis of the underlying prey–predator model ^[1,2,3]. The evidence of combative interactions is widespread in unicellular to multicellular organisms ^[4,5]. The opportunistic human pathogen Pseudomonas aeruginosa uses a quorum-controlled mechanism, called policing effect ^[6], to actively mitigate the unrestricted emanation of cheating. In the course of strategic policing the defectors, the cooperative subpopulation of P. aeruginosa secrets hydrogen cyanide in sequence with the production of public goods, such as elastase, to combat the unimpeded exploitation of common goods by the toxin-sensitive defectors ^[6,7]. Multicellular organisms often use perceptible color signals to avoid the predation, commonly known as aposematism ^[8,9], which is conceptually similar to the policing effect in bacteria.

Aposematism or warning coloration is one of the strategic defensive mechanisms, which is often incurred by prey to reduce the predation-induced mortality rate ^[9]. Aposematic preys utilize their conspicuous color signals to express the deleterious effects of the consumption of the chemically-defended preys to the predators. However, from the present phenotypic state of defended prey, predators explicitly measure the trade-off between the nutritional benefit from the consumption and the cost of toxin ingestion, and regulate their rate of intake of toxic preys ^[10,11,12]. The undefended or defended species share the warning color pattern of the sympatric aposematic species which are described as Batesian and M$\ddot{\text{u}}$llerian mimicry, respectively ^[13]. It needs to be mentioned that the avoidance learning of predator on the basis of distinctive warning coloration of defensive preys fundamentally supports the evolution of aposematism and mimicry in an ecological system ^[10,11,12,9]. In a general resource-consumer model, the intrinsic characteristics of resource is to facilitate the abundance of consumer, whereas, consumer hinders the abundance of available resources; however, an inclusion of natural defense against consumer (e.g., prey–predator models with aposematic prey) additionally includes a negative feedback to its abundance.

The existence of aposematism is ubiquitous among animal taxa, including invertebrates, fishes, amphibians, snakes, and birds ^[14,15,16]. As the example, dendrobatid frog Oophaga pumilio ^[15], venomous coral snake M. fulvius ^[17] and mealworm Tenebrio molitor ^[10] often use the bright colors and/or patterns to represent their unpalatability or noxiousness to the potential predators. In La Selva Biological Station, Costa Rica, Saporito et al. have experimentally validated that the red colored plasticine models of O. pumilio receive lesser attacks from avian predators as compared to the brown frog models——indicates that the conspicuous bright color (i.e., aposematic signal) in real alkaloid-contained dendrobatid frog O. pumilio acts as a chemical-defense mechanism against predation ^[15]. The tricolor banded pattern of venomous coral snakes serves as an aposematic signal for the free-ranging avian predators as experimentally demonstrated by Edmund Brodie ^[17]. The experiments in Refs. ^[15] and ^[17] contain plasticine models to conceptualize the aposematic signal indicator in poison frog and coral snakes for the avian predators; however, Rowe et al. ^[10] have done the in-vitro experimentation with live mealworm larvae Tenebrio molitor to validate the impact of warning coloration on the avian predators——European starlings Sturnus vulgaris. In Ref. ^[10], the primary research question is: How is the evolution of advertising the unprofitability in the consumption of aposematic prey and the cognitive avoidance learning process of naive predators interrelated?

The sustenance of species primarily depends on the self-protection and food accumulation adjoined to mating ^[18]. In the present manuscript, our logical assumption is that any species primarily invests its total energy or time into two acts: (a) defense against predation and (b) search for food added with reproduction. Therefore, if any species invests more time in self-defense (e.g., time for aposematism), then the rate of reproduction naturally diminishes; however, it increases the prey's survival rate via reducing the rate of consumption of defended prey by predator. On the basis of this assumption, here, we construct a diffusive prey–predator model with defended (i.e., aposematic) prey, and we assume that the prey's reproduction rate and the rate of predation negatively linearly vary with the level of toxicity, which is proportional to the time invested for aposematism. Moreover, we consider the unlimited resource environment for prey which ensures the non-existence of intra-prey competition-induced death of prey. The primary finding of the current manuscript is that the aposematic time and food searching efficiency of prey can substantially regulate the temporal as well as spatio-temporal population dynamics of prey–predator systems.

The rest of the manuscript is organized as follows: Section 2 presents the formulation of temporal and spatio-temporal models. Mathematical analysis of the temporal model is described in Section 3. Section 4 deals with the stability analysis of the spatio-temporal system. Extensive numerical simulations are presented in Section 5. Section 6 concludes the manuscript with a few remarks on the analytical and numerical results of the proposed model and future scope.

2.   Description of model

To analyze the impact of aposematism on a resource-consumer system, we consider a simple prey–predator model, where prey population shows anti-predation defense mechanism for their survival. Our two-fold modeling study consists of the analysis of temporal dynamics and the subsequent investigation of the spatio-temporal system. Suppose, $N$ is the density of prey individuals, and $P$ is the density of predator individuals. We divide the total time invested by aposematic prey into two parts: they regulate the degree of aposematism (i.e., the level of toxicity) to avoid predation and they search food for their growth-induced reproduction. Let $a$ be the fraction of time for aposematism in one unit time, and $b\; ( = 1-a)$ be the fraction of food and/or mate searching time in one unit time. Although the aposematic time $a$ is a behavioral strategy of prey, in the present study, our primary motivation is to analyze the impact of non-adaptive nature of $a$ on the prey–predator model first. We hope to find the influences of adaptation of aposematic time on population dynamics in the further studies.

To derive the recruitment rate of prey population, we assume that the resource input rate is constant, $r$. This assumption will simplify the dynamics and is also reasonable if the resource recovers very fast to remain in steady state ^[18]. Therefore, one mature prey admits the resource $\frac{r}{N}$. Suppose, $s_0$ is the searching efficiency of the resource for prey individuals without the aposematism. Since, prey invest only fraction $b$ time in one unit time, the searching efficiency is reduced from $s_0$ to $s_0b$, i.e., $s_0(1-a)$. Now, the uptake rate of resource for one individual prey is $w = s_0(1-a)\frac{r}{N}$. We consider the simple Beverton-Holt recruitment function $\frac{k_1 w}{1+k_2 w}$ to obtain per-capita birth rate of prey (which includes the saturation effect of converting nutrient into offsprings):

where $s = k_1s_0r$ and $k = k_2s_0r$. Now, with a consideration of Holling type $II$ predators' response, the proposed prey–predator model with aposematic prey is given by

where $N = N(t)$ and $P = P(t)$ denote the population densities of prey and predator at time $t$, respectively. $s$, $a$, $k$, $d$, $p$, $h$, $c$, $m_{1}$ and $m_{2}$ are positive constants. $d$ is the intrinsic death rate of prey, $p$ is the rate of predation, $h$ is the half-saturation constant, $c$ is the conversion efficiency for converting preys' biomass into predators' biomass, $m_1$ and $m_2$ are the natural death rate and intra-predator competition-induced death rate of predator, respectively. From the temporal dynamics (1), we can conclude that if prey invests its entire time for aposematism (i.e., $a = 1$), then its growth-induced reproduction rate diminishes to zero, and due to the excessive toxicity in prey, predator naturally avoids the consumption of aposematic prey which subsequently ensures the population extinction.

If we consider the diffusion in prey and predator populations in the above Model 1, then we have the following reaction-diffusion system

where $N \equiv N(x, y, t)$, $P \equiv P(x, y, t)$ and $(x, y)$ is the spatial coordinate. $\frac{\partial}{\partial t}$ denotes the partial derivative with respect to time $t$. $\nabla^2 \equiv \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ is the Laplacian Operator in two dimensional space. $D_N$ and $D_P$ are the self-diffusion coefficients of prey and predator, respectively, which describe the individual movement from higher to lower concentration. As we choose the uniform environment, the model parameters do not depend on space and time.

The spatio-temporal dynamics (2) is subject to the non-zero initial conditions and zero-flux (or Neumann) boundary conditions ^[19]:

Here, $\widehat{n}$ is the outward unit normal vector of the boundary $\partial \Omega$ which is assume to be smooth. We consider zero-flux boundary conditions (i.e., there is no external input) as in the present manuscript, we are primarily interested in the self-organization of patterns.

3.   Mathematical analysis of temporal model

In this section, first we will analyze our non-spatial Model 1.

3.1.   Boundedness of positive solutions

All non-negative solutions of the Model 1 that start in $\mathbb{R}^2_+$ are uniformly bounded. It is easy to check that the solution $N(t)$ of 1 tends to $0$ (hence also $P(t)$ tends to $0$) if $s\leq dk$. In fact $\frac{dN}{dt} < 0$ if $s\leq dk$ and $N > 0$.

Theorem 3.1. The system 1 is positively invariant and uniformly ultimately bounded in $\mathbb{R}^2_+$, with the following properties

Proof. For any $N > 0$ and $P > 0$, we have $\frac{dN}{dt}\Big|_{N = 0} = 0$ and $\frac{dP}{dt}\Big|_{P = 0} = 0$, which implies that $N = 0$ and $P = 0$ are invariant manifolds, respectively. Due to the continuity of the system, we can conclude that the system 1 is positively invariant in $\mathbb{R}^2_+$.

Now, choose any point $(N, P)\in\mathbb{R}^2_+$, such that $N > \frac{(1-a)(s-dk)}{d}$, then due to the positive invariant property of 1, we have

In addition, we have $\frac{dN}{dt}\Big|_{N = \frac{(1-a)(s-dk)}{d}} = 0$. Therefore, we can conclude that

Now define, $X(t) = N(t) + \frac{1}{c}P(t)$, then we have

Case Ⅰ: $d\geq m_1$

For sufficiently large $t\geq0$, we have

By applying the theory of differential inequality, we have

Case Ⅱ: $d < m_1$

For sufficiently large $t\geq0$, we have

Hence,

Consequently, we have

3.2.   Equilibrium analysis

The Model 1 has three non-negative equilibria:

$1.$ The extinction equilibrium $E_0 = (0, 0)$, which always exists.

$2.$ The axial equilibrium, where only prey population survive is given by $E_n = \Big(\frac{1}{d}(1-a)(s-dk), 0\Big)$, and which exists if $s > dk$.

$3.$ The interior equilibrium is given by

and which exists if

therefore, we must have $m_1 < cp(1-a)$. Here, $N_*$ is the real positive root of cubic polynomial

where

Theorem 3.2. There exists a unique positive equilibrium if $s > dk$ and $m_1 < \frac{cp(1-a)^2(s-dk)}{(1-a)(s-dk)+hd}$.

Proof. Since, $s > dk$, then the coefficient $C_4 < 0$. Now, by Descartes' rule of signs the polynomial $\Phi(z)$, can have more than one real positive solution if and only if $C_2 < 0$, but $C_3 > 0$. Therefore, for the existence of unique positive real solution, it is sufficient to show that, under the condition $s > dk$, $C_3$ must be negative, when $C_2$ is negative. Now, $C_3$ can be written as

Therefore, there exists at most one positive real root $N_*$ of $\Phi(z) = 0$.

Now, since $\frac{cp(1-a)^2(s-dk)}{(1-a)(s-dk)+hd} < cp(1-a)$, and $m_1 < \frac{cp(1-a)^2(s-dk)}{(1-a)(s-dk)+hd}$, we have $m_1 < cp(1-a)$. We will prove that this unique solution $N_*$ satisfies

It can be easily verified that $cp^2(1-a)^2\overline{N}(k(1-a)+\overline{N}) - m_1(1-a)p(h+\overline{N})(k(1-a)+\overline{N}) = 0$.

Therefore,

This means $N_* > \overline{N}$.

3.3.   Stability analysis

$1.$ The extinction equilibrium $E_0$ is locally asymptotically stable if $s < dk$, as the eigenvalues associated with the equilibrium $E_0$ are given by $\lambda_1 = \frac{1}{k}(s-dk)$ $< 0 \mbox{ if } s < dk$, and $\lambda_2 = -m_1$ $(<0)$.

In fact, $E_0$ is globally asymptotically stable under this condition. In fact, we have

Consider the system, $\frac{dq}{dt} = \Big(\frac{s}{k}-d\Big) q$; with $q(0)\geq N(0)$. Since, $s < dk$, we define $\Big(\frac{s}{k}-d\Big) = -\varepsilon \; (<0)$. Therefore,

which gives $q(t) \rightarrow 0$ as $t \rightarrow \infty$.

By comparison theory, we have $N(t) \leq q(t)$, for $t\geq0$. Hence, $N(t) \rightarrow 0$ as $t \rightarrow \infty$, i.e., prey population goes to extinction.

When preys goes to extinction, then we have $\frac{1}{P}\frac{dP}{dt} < -m_1 < 0$, and by similar way we can show that $P(t) \rightarrow 0$ as $t \rightarrow \infty$, i.e., predator population also goes to extinction.

$2.$ The axial equilibrium $E_n$, exists if $s > dk$, and locally asymptotically stable if $m_1 > \frac{cp(s-dk)(1-a)^2}{(s-dk)(1-a)+hd}$. Since the eigenvalues at the equilibrium $E_n$ are given by $\lambda_1 = -\frac{d}{s}(s-dk)$ $< 0$ as $s > dk$, and $\lambda_2 = \frac{cp(s-dk)(1-a)^2-m_1(s-dk)(1-a)-m_1hd}{(s-dk)(1-a)+hd} < 0$, if $cp(s-dk)(1-a)^2-m_1(s-dk)(1-a)-m_1hd < 0$, i.e., if $m_1 > \frac{cp(s-dk)(1-a)^2}{(s-dk)(1-a)+hd}$.

$3.$ The Jacobian matrix at the interior equilibrium $E_{*}$ is given by

where,

The characteristic equation at the interior equilibrium is given by

where $A = -\mathrm{tr}(J) = -(a_{11}+a_{22})$, and $B = \det(J) = a_{11}a_{22} - a_{12}a_{21}$. Thus, by Routh-Hurwitz criterion, the interior equilibrium $E_{*}$ is locally asymptotically stable iff $A > 0$ and $B > 0$.

The existence and stability conditions of equilibria for the Model 1 are summarized in Table 2.

Next, we investigate the possibility of Hopf-bifurcation at the interior equilibrium $E_*$ by considering the parameter $s$, the searching efficiency of prey for resource, as the bifurcation parameter.

The interior equilibrium $E_*$ loses its stability through Hopf-bifurcation when the eigenvalues are complex conjugate with zero real parts. We consider the parameter $s$ as the bifurcation parameter. Let, $\lambda(s) = \lambda_r(s) + i\lambda_i(s)$ be the eigenvalues of the characteristic equation 3. After substituting the value of $\lambda$ in equation 3, and separating real and imaginary parts, we get

At the Hopf-bifurcation point, we have $\lambda_r(s) = 0$. We set at $s = s_H$, $\lambda_r(s_H) = 0$, and put $\lambda_r = 0$ in 4. Therefore,

Therefore, from the above equations, we have $A(s_H) = 0$, and $\lambda_i(s_H) = \sqrt{B(s_H)} > 0$, i.e., $\det{J}\Big|_{s = s_H} = B(s_H) > 0$. Thus, at the Hopf-bifurcation point, we have

Differentiating equations 4, w.r.t. $s$, and put $\lambda_r(s) = 0$, we have

By solving the above system of equations, we have

provided $[2\lambda_i^2\frac{dA}{ds}+A\frac{dB}{ds}]\Big|_{s = s_H} \neq 0$. Hence, we can write the following theorem:

Theorem 3.3. If $\det{J}\Big|_{s = s_H} > 0$ and $\frac{d\Re(\lambda(s))}{ds}\Big|_{s = s_H} \neq 0$ hold, then the interior equilibrium $E_{*}$ of Model (1) is locally asymptotically stable when $s < s_{H}$, and the Model (1) undergoes Hopf-bifurcation at $E_{*}$ when $s = s_H$.

4.   Stability analysis of the spatial model

In this section, we will investigate the Turing instability of the Model 2, where the spatially-homogeneous steady state $E_*$ of the system is stable without diffusion, but unstable in the presence of diffusion. First, we consider the linearized form of the system 2 about the positive equilibrium $E_* = (N_*, P_*)$ as follows:

where $N = N_*+u$ and $P = P_*+v$, and $(u, v)$ are small perturbation of $(N, P)$ about the interior equilibrium point $E_* = (N_*, P_*)$. Now we can consider the solution of the system 5, as

Here $\xi$ is the growth rate of perturbation in time $t$, and $\kappa_x$, $\kappa_y$ are the wave numbers of the solutions. The Jacobian matrix of the linearized system can be written as:

In spatial model, the growth rate of perturbation $\xi$ depends on the sum of square of the wave numbers $(\kappa_x^2+\kappa_y^2)$. Therefore, both the wave numbers will affect the eigenvalues of $\overline{J}$. For simplicity, we can consider $\xi$ as rotational symmetric function on the $\kappa_x \kappa_y$ - plane and substitute $\kappa^2 = \kappa_x^2+\kappa_y^2$ and obtain the results for the two-dimensional case from the one-dimensional formulation. The corresponding characteristic equation is given by

where

Thus using the Routh-Hurwitz criterion, we have the following theorem for the stability of spatially-homogeneous positive equilibrium of the Model 2.

Theorem 4.1. The positive equilibrium point $E_*$ of the system 2 is locally asymptotically stable iff $\overline{A} > 0$ and $\overline{B} > 0$.

From the expressions of $A$ and $\overline{A}$, it is clear that if $A > 0$, then $\overline{A} > 0$. Thus the only possibility for the occurrence of diffusion-induced instability is the case when $B > 0$, but $\overline{B} < 0$. Thus the condition, for diffusive instability is given by

Therefore, the diffusion can induce the loss of stability with respect to the perturbation of certain wave numbers. Here, $H$ is a quadratic function of $\kappa^2$ and the graph of $H(\kappa^2)$ is a parabola. Let the minimum of $H(\kappa^2)$ is reached at $\kappa^2 = \kappa_c^2$, where $\kappa_c^2$ is given by

As $a_{11}+a_{22} < 0$, and $\kappa_c$ is real then we must have $a_{11}a_{22} < 0$. The sufficient condition for instability is that $H(\kappa_c^2) < 0$. Thus, with the critical value of $\kappa = \kappa_c$, the condition for diffusive instability can be written as

It is easy to see that the change of sign in $H(k^2)$ occurs when $k^2$ enters and leave the interval $(\kappa_{-}^2, \kappa_{+}^2)$ where

In particular, we have $H(\kappa^2) < 0$ for $\kappa_{-}^2 < \kappa^2 < \kappa_{+}^2$. Also, the diffusive instability cannot occur unless the diffusivity ratio is sufficiently away from unity. In fact, we have $a_{11}+a_{22} < 0$, and therefore $a_{11} < -a_{22}$ (where $a_{22} < 0$). Then, from the condition 8, we have

The above conditions are the necessary conditions for the Turing instability, which is applicable to any two-species activator-inhibitor system. In particular, we must have $D_P\neq D_N$.

From the inequality 8, the explicit condition for the occurrence of Turing instability becomes:

5.   Numerical simulations

In this section, we perform some numerical simulations for both temporal and spatio-temporal dynamics. For our numerical simulations, the fixed parameter values are given in Table 1.

5.1.   Temporal dynamics

First, we study the temporal dynamics of the Model 1. We draw the bifurcation diagrams of both prey and predator with respect to the bifurcation parameter $a$ (the fraction of time for aposematism in one unit time). From the Figure 2, we can observe that when $a$ is small $(a\approx0.08)$ system shows the oscillatory coexistence between prey and predator. As we increase the parameter $a$, system becomes asymptotically stable, and the abundance of both prey and predator increases up to a certain threshold value of $a$. Predator and prey abundance reaches maximum at $a\approx0.45$ and $a\approx0.75$, respectively and then decreases to zero. The predator becomes extinct at very high value of $a$ $(a\approx0.9)$, when preys are too much defensive such that predator may not capture sufficient prey to survive, and the predator is specialist. At $a = 1$, the extinction of population occurs which is also clear from the dynamics (1). Therefore, from the bifurcation diagrams, we can say that aposematism time can increase the stability and abundance of both the populations up to a certain range, after that the abundance decreases.

Similarly, we examine the effect of the parameter $s$ (the searching efficiency of prey) for the Model 1. From the bifurcation Figure 3, we can observe that as we increase the parameter $s$, the system behavior changes from stable to limit cycle (at $s\approx0.48$) through Hopf-bifurcation and again changes from limit cycle to stable behavior (at $s\approx1.26$). Therefore, for some intermediate values of $s$, the system shows oscillatory behavior, otherwise system become stable for a long range of the parameter values of $s$. Also, as we increase the parameter $s$, the abundance of both prey and predator increases.

5.2.   Spatial dynamics

5.2.1.   Existence of spatial heterogeneity

First, we prove the possibility of the existence of Turing instability in our Model 2. The diffusion-driven instability of an equilibrium means if it is asymptotically stable without diffusion, but unstable in the presence of diffusion. To confirm the diffusion-driven instability around the interior equilibrium, here we check the analytic conditions numerically. First, we check the diffusion-induced instability for the parameter $a$. In the Figure 4(a), we plot the polynomial $H(\kappa^2)$ 7 with respect to $\kappa^2$ for different values of $a$ $(a = 0.05, \; 0.15, \; 0.25, \; 0.30, \; \mbox{and }\;0.35)$. From the Figure 4(a), we can observe that as the parameter $a$ increases, the range of $\kappa^2$, for which the polynomial $H(\kappa^2)$ remains negative decreases, and finally $H(\kappa^2)$ becomes completely positive for $a > 0.30$. Thus, by increasing the aposematism time $a$, the possibility of the occurrence of Turing instability decreases. As the largest real part of $\xi$ 6 also provides information about the existence of Turing instability. Largest real part of $\xi$ is positive, implies the possibilities of the Turing instability. We can observe that as $a$ increases the maximum of the real part of $\xi$ decreases and become completely negative for $a > 0.30$.

Similarly, in the Figure 5, we draw the polynomial $H(\kappa^2)$ and the largest real part of $\xi$ with respect to $\kappa^2$ and check the possibility of the occurrence of Turing patterns for different values of $s$. Here, we can observe that Turing instability can occur only at the intermediate values of $s$. For the fixed set of parameter values as in Table 1, with $a = 0.25$, if $s$ is between $0.4$ and $1.1$ then only pattern can be observed. The polynomial $H(\kappa^2)$ remains positive and real part of $\xi$ remains negative if $s < 0.4$ or $s > 1.1$.

5.2.2.   Effect of aposematic time and searching efficiency

Next, we investigate how the stability of the system changes by changing both aposematism time $(a)$ and searching efficiency of prey $(s)$ in $(s, \; a)$ parameter plane (Figure 6). Here, black region is for the extinction of both prey and predator populations. In this region the extinction equilibrium $E_0$ is globally asymptotically stable. In the magenta region, only prey population survive and the extinction of predator occurs. In this region the axial equilibrium $E_n$ is locally asymptotically stable. The white region stands for the stability region of the coexistence equilibrium $E_*$ in both spatial and non-spatial cases (both with and without diffusion). The yellow region represents the region of Turing instability, where the system is stable without diffusion but unstable and shows different spatial pattern in the presence of diffusion. The red region represnets the Turing-Hopf region, i.e., system shows oscillatory coexistence without diffusion but there is a possibility of inhomogeneous stationary patterns due to the interaction of Turing instability. From the Figure 6, it is clear that both the parameters $s$ and $a$ are very important for determining different dynamics of the system. For example, for this set of fixed parameter values as in the Table 1, when $s$ is very small, both the prey and predator extinct, and for large value of $a$, predator population extinct. For some intermediate values of $s$ and low value of $a$, different spatial patterns (both Turing-Hopf and Turing patterns) occur. For large value of $a$ $(a > \approx0.30)$ there is no possibility of Turing patterns and the system is stable at the coexistence equilibrium in a large parameter domain.

5.2.3.   Different types of pattern formations

Now, after finding the possibilities of getting Turing instability and range of parameter values for the parameters $s$ and $a$, we investigate the different types of Turing patterns by varying $s$ and $a$ for the Model 2. For obtaining the spatio-temporal dynamics in two dimensional spatial domain, the system of equations in 2 is solved numerically by using a finite difference method. The forward difference Euler scheme is used for the reaction part and the standard five point explicit finite difference scheme is used for two dimensional diffusion terms. We use positive initial conditions and homogeneous Neumann boundary conditions at the boundary of the domain. Numerical simulations for all the stationary patterns are carried out over a square domain $\Omega = [0,200]\times[0,200]$ with the spatial grid sizes $\Delta x = \Delta y = 1$ and the temporal grid size $\Delta t = 1/30$. For the simulations, the initial conditions $(N_0, P_0)$ are always a small random perturbation around the interior equilibrium $E_*$. $(N_0, P_0) = (N_*+\delta \zeta_{i, j}, P_*+\delta \eta_{i, j})$, where $\delta = 0.0001$, $\zeta_{i, j}$ and $\eta_{i, j}$ represent the usual Gaussian white noise. After the initial period of perturbation spreads, the system goes into either a time-dependent state or an essentially steady state solution, and we stopped our simulations after sufficeint time to assume that the patterns attained the stationary state and they do not change further with time. In all the simulations of pattern formation, we run our simulations up to $t = 50, 00,000$ time points, so that the patterns do not change further with time. While doing the simulations, different types of spatial dynamics have been observed, and almost one-to-one correspondence is observed between the stationary patterns of prey and predator. Thus, we have only shown the spatial distribution of prey.

First, we find different types of patterns formation in Turing domain (parameter values are chosen from the yellow region of the Figure 6). As we vary the aposematic time investment of prey in one unit of time $(a)$, and keeping all the other parameters fixed as mentioned in Table 1, we get different spatial patterns (Figure 7). When the aposematic time investment is small, the system shows spot patterns for a long range of values of $a$. Here, in Figure 7a, we have drawn the spot patterns for $a = 0.22$, where the abundance of prey is higher in isolated red zones. In view of population dynamics of the proposed model, this 'spot' pattern signifies that when prey aposematism time is small, then prey population is outcompeted by the predator population, and fixes to a very low density in majority of the spatial region due to its weak defense level. It results in the formation of patches of high prey density surrounded by areas of low prey densities ^[20].

Now, as we increase the value of aposematism time, at $a = 0.25$ we get a mixture of spot and stripe patterns (Figure 7b). Further increasing the value of $a$, at $a = 0.28$, there is a regular peaks and troughs of prey density, i.e., the stripe pattern emerges (Figure 7c). At $a = 0.29$, we observe a mixture of stripes and holes (Figure 7d). Finally, for relatively large value of $a$, in the Turing domain, at $a = 0.295$, the system dynamics changes to holes (Figure 7e), i.e., the prey population is at low density in the isolated zone and exist at high density in the remaining region. Thus by increasing the value of aposematism time $a$, the patterns changes from spots to holes through spot-stripe mixture $\rightarrow$ stripes $\rightarrow$ stripe-hole mixture.

Next, in the Turing domain, we look at the different patterns formed by varying the prey searching efficiency parameter $(s)$ and keeping all the other parameters fixed as mentioned in the Table 1. At low values of s in the Turing domain, spots patterns is observed (at $s = 1.1$, Figure 8a). Thus, as prey searching efficiency is small, prey population with very low density is observed in the majority of the spatial region. As we increase the searching efficiency $s$, at $s = 1.21$ a mixture of spots and strips pattern is observed (Figure 8b). Further increase of $s$ results stripes pattern (at $s = 1.4$, Figure 8c). At $s = 1.45$, a mixture of stripes and holes pattern can be observed (Figure 8d). Finally, at large value of $s$ in the Turing domain, we observe the holes pattern only (at $s = 1.5$, Figure 8e). Thus with high searching efficiency, higher density of prey has been observed in the majority of the spatial domain. Therefore, in this case also, as we increase the value of searching efficiency the pattern sequence 'spots $\rightarrow$ spot-stripe mixture $\rightarrow$ stripes $\rightarrow$ stripe-hole mixture $\rightarrow$ holes' is observed. In both the cases, as we increase the aposematism time or searching efficiency of prey, the system changes from predator predominance to prey predominance.

5.2.4.   Turing-Hopf patterns

Finally, we find the types of patterns observed in the Turing-Hopf domain (parameters are chosen from the red region of the Figure 6). In the Turing-Hopf domain, only spots patterns have been observed. In Figure 9, we have shown the spots patterns for both prey and predator. Thus, for our model, when temporal dynamics shows limit cycle oscillations, then the spatio-temporal dynamics shows predator predominance in the system. From Figures 7, 8, and 9, we can observe that the density variation is higher in Turing-Hopf domain (non-Turing patterns) as compared to corresponding Turing patterns.

6.   Discussion

Self-defense and growth-induced reproduction are two basic intrinsic characteristics of any species in an ecosystem ^[18]. Specifically, in a prey–predator model, the prey population most often invests its total time in two actions——(ⅰ) self-protection against predation (such as camouflage, warning coloration or aposematism) and (ⅱ) search food for growth-induced mating; However, predators are primarily involved in searching of preys for their diet, which is adjoined to mating of mature predators. The aposematic time and time invested in searching of food for prey are complementary of each other. An increase in time for aposematism naturally decreases the reproduction rate of prey, whereas, predators often diminish the rate of intake of highly defended prey ^[11]. With the consideration of these logical assumptions, in the present manuscript, we formulate a diffusive prey–predator model with aposematic prey, which follows Beverton-Holt recruitment function ^[18] for its per-capita growth rate. The proposed model (1) contemplates unlimited-resource environment for prey which justifies the non-existence of intra-prey competition; however, the Model (1) considers intra-predator competition-induced death of predators. To include the saturation effect in the processing of food, we use Holling type-Ⅱ functional response for predation and corresponding per-capita growth of predator.

In the process of analyzing the temporal and spatio-temporal version of the suggested model, we primarily concern about the effect of aposematic time $(a)$ and searching efficiency $(s)$ of prey on the dynamic outcomes. The Model (1), which is positively invariant and uniformly ultimately bounded (see Theorem 3.1), shows Hopf-bifurcation with respect to the bifurcation parameter $s$ (see Theorem 3.3 and Figure 3). Unlike the searching efficiency $s$, the proportion in which prey and predator coexist (see Theorem 3.2 and Figure 1) declines after a certain value of aposematic time $a$; moreover, at the highest level of aposematism $(a = 1)$, the prey population rarely gets time for mating, and on the other side, predators avoid the consumption of extremely defended prey which jointly results in the extinction of entire population as shown in Figure 2. Like temporal dynamics, the impact of aposematic time and searching efficiency on the formation of spatial patterns (as emerged from the reaction-diffusion model (2)) is notable. The steady-state spatial abundance of prey or predator population, in the form of regular Turing and non-Turing patterns, varies with respect to both $a$ and $s$ (see yellow and red regions in Figure 6). The individual effect of $a$ and $s$ on the organization of Turing patterns is conceptually similar. For the low value of aposematic time (i.e., less defended prey), the rate of predation increases with a less dominant increment in the reproduction rate of prey——results in 'spot' pattern, which confirms the presence of high density in discrete patches with low density in the surroundings as shown in Figure 7 (a). Further increase in the value of $a$ shows the transition from 'spots' to 'mixture of spots and stripes' to 'stripes' to 'mixture of stripes and holes' to 'holes' (low density in the center of patches and high density in the surroundings) pattern, which is generated due to declination in the predation rate of substantially defended prey (see Figure 7). Similarly, low food searching efficiency of prey results in 'spot' patterns, whereas, the abundance of prey increases ('spots' to 'stripes' to 'holes') with an increase in $s$ as shown in Figure 8. As, in the proposed reaction-diffusion model (2), the predator dynamics is coupled with the dynamics of prey, we get the same Turing patterns for both prey and predator. The combined effect of $a$ and $s$ in the formation of non-Turing patterns is shown in Figure 9, where, we get only 'spot' pattern with the dominance of predators.

The work presented in the current manuscript can give a clue towards understanding the impact of both prey defense and preys' food searching efficiency on the extinction, coexistence and spatial distribution of the constituent populations of a real ecosystem. Here, we consider that the reproduction rate of prey and the rate of predation vary linearly with the time for aposematism $(a)$; however, one can analyze the effect of nonlinearity in aposematic time on the dynamic outcomes of the proposed model. Our logical assumption is that the predators avoid the aposematic preys on the basis of its level of toxicity (proportional to the amount of time invested for aposematism), which can be perceived through the degree of warning coloration, although, we have not considered the dynamics of $a$. The consideration of adaptive dynamics for the toxicity level of aposematic prey can incline the dynamic outcomes more towards reality. In the existing literature ^[11], it has been shown that the predator's decision to include the aposematic preys in their diet is influenced by the presence of undefended (i.e., non-toxic) preys, therefore, one can further analyze the 3D dynamics with the subpopulations of undefended prey, defended prey and predator ^[21].

In the realm of evolutionary biology, the evolution of aposematism and its paradoxical effects have long been the topics of interest ^[22,23,24]. The advertisement of noxiousness and unpalatability of aposematic prey through conspicuous warning coloration often initially suffers from immense predation by naive predators——questioned those mechanisms, which are responsible for the conversion of cryptic prey into overtly defended prey ^[22].

Acknowledgment

SKS's research is supported by the JSPS post-doctoral fellowship. JB's research is supported by the Teaching Assistantship from Ministry of Human Resource Development, Government of India. YT's research is supported by the JSPS through the "Grant-in-aid 26400211." We would like to thank Ritwik Kumar Layek of IIT Kharagpur and Helen Byrne of Mathematical Institute, University of Oxford for the fruitful discussions.

Conflict of interest

The authors declare that they have no conflict of interest.