Dynamics and pattern formation of a diffusive predator−prey model in the subdiffusive regime in presence of toxicity

Motivated by the fact that the restrictive conditions for a Turing instability are relaxed in subdiffusive regime, we investigate the effects of subdiffusion in the predator−prey model with toxins under the homogeneous Neumann boundary condition. First, the stability analysis of the corresponding ordinary differential equation is carried out. From this analysis, it follows that stability is closely related to the coefficient of toxicity. In addition, the temporal fractional derivative does not systematically widen the range of parameters to maintain a point in the stability domain. Furthermore, we derive the condition which links the Turing instability to the coefficient of toxicity in the subdiffusive regime. System parameters are varied in order to test our mathematical predictions while comparing them to ecological literature. It turns out that the memory effects, linked to the transport process can, depending on the parameters, either stabilize an ecosystem or make a completely different configuration. D. C. Bitang A Ziem University of Yaounde 1, 812 Yaounde, Cameroon E-mail: danielcassidi7@gmail.com C. L. Gninzanlong University of Yaounde 1, 812 Yaounde, Cameroon C. B. Tabi Botswana International University of Science and Technology, Private Bag 16, Palapaye, Botswana T. C. Kofane University of Yaounde 1, 812 Yaounde, Cameroon


Introduction
During evolution, venomous animals have produces a panoply of peptide toxins, formidable weapons for defense against predators or for the capture of prey. The targets of these toxins are proteins involved in nervous conduction as well as nervous transmission, mainly ion channels. These venoms are an invaluable source of very specific pharmacologicals agents for characterization of ion channel and receptor subtypes membranes. By their efficiency and selectivity of action, they offer significant therapeutic potential as evidenced by toxins whose effectiveness in neurological disorders is now being tested in therapeutic trials. Recent studies point out the fact that nerve conduction and transmission are overwhelmingly dominated by subdiffusion [1][2][3]. Subdiffusion has acquired relevance in the past decades since it has been experimentally detected in several systems such as porous media [4], glasses [5], transport through cell membranes [1,2], and other biological systems [3]. Studies of subdiffusion in reaction-diffusion systems, using time fractional derivatives, have shown that the conditions for a Turing instability are not different from the normal conditions for special cases of subdiffusion [6,7]. However, Hernandez et al. [8] came out with a new way of obtaining Turing patterns with subdiffusion. Since the well−known Turing reaction−diffusion model was discovered [9], enormous efforts have been made to investigate nonlinear self−organization phenomena in nature. Recently, interest was given to the link between toxic substances and Turing bifurcation, Hopf bifurcation and Steady − State bifurcation [10]. This lead to the conclusion that the toxic coefficient is important for the occurence of the complex spatial dynamics. Additionaly, they considered a typical reaction−diffusion system, which implies that their combined influence on the evolution of the density of a species is a simple sum, because reaction and diffusion processes in this way are separable [11]. When it comes to subdiffusive entities undergoing reactions, reaction and diffusion processes are no longer separable because of the strong memory effects in the transport mechanism [11], which means that the non-Markovian nature of subdiffusion results in a nontrivial combination of reactions and spatial dispersal. In this situation, can toxicity still be considered as an important parameter responsible for complex spatial dynamics?
The answer to this question requires to investigate how strong memory effects in the transport process influence Turing pattern induced by toxic substances, which is the purpose of this paper. The scheme of this paper is as follows. In Sec. 2, we start by examining the influence of toxic substances on the dynamics properties of the fractional−order ODE system. Sec. 3 is devoted to the analysis of a model describing an ecosystem where the prey produces substance that is toxic for the predator, and give general conditions for toxicity coefficient allowing Turing instability to appear. We verify the predictions of the analysis via numerical calculations in Sec. 4. In Sec. 5, we summarize our results and draw some conclusions.
2 Effect of toxic substances on the dynamics properties of the fractional−order ODE system In this work, the following parabolic activator-inhibitor model with subdiffusion consisting of two variables is investigated: The reaction terms in Eq. (1) were first proposed by Zhang and Zhao to describe an ecosystem where the prey produces substances that are toxic to the predators [10]. The term βuv 2 denotes the effect of toxic substances and β represents the efficiency of toxicity. u and v denote the sizes of the prey population and predator population, respectively. r, K, m, a, s, h, and β are all positive constants, and η ≤ 1. Note that in this form, the evolution of the uniform state is controlled by subdiffusion, even if the spatial term is not present. In here the time operator is the Caputo fractional derivative of order η [12]. We first examine the stability and Hopf bifurcation of the ODE version of system (1). Omitting the diffu-sion terms, the following ODE is obtained from (1) An obvious steady state of system (2) is E 1 = (K, 0), the predator−free equilibrium. Assuming that the other equilibrium E * = (u * , v * ) of system (2) satisfies Eliminating v yields The number and the signs of the roots can easily be predicted from (4), according to Descartes' rule of signs. In doing so, system (4) admits one or three positive roots u = u * if one of the following inequalities holds a < K, hs < Kaβ, Km + ahr < Khr, a < K, hs > Kaβ, Km + ahr < Khr, Once we know the number of equilibria, a word should be said about their stability. The stability of the steady states is determined by the eigenvalues of the Jacobian matrix and the associated characteistic equation is where T = a 11 + a 22 , ∆ = a 11 a 22 − a 12 a 21 , with If the trace T of the Jacobian matrix is negative, with ∆, its determinant, being positive, the steady state is stable. While looking at the sign of T and ∆, respectively, we should not forget to distinguish the case of a single positive root from the one of three positive roots. Then, taking into account conditions (5), for η = 1, the steady state E * will be stable if where with β H = α 2 being a Hopf bifurcation for system (2). If, on the other hand, we consider conditions (6), steady states will be stable if β is bounded as follows However, for η ∈ (0, 1), the stability of E * is more complexe and depends on the inequality for the roots λ in the equation (9). In order to really appreciate the validity of our predictions, it would be wise to assign some values to the constitutive parameters of Eq. (4). In that respect, we choose parameters as follows K = 4, r = 0.7, m = 1 a = 0.4, s = 0.25, h = 1. This set of parameters leads to α 1 = 0.15625, α 2 = −0.059731 and α 3 = −0.0117744. yielding max{α 1 , α 2 , α 3 } = α 1 . Let us first consider the inequality (5). We choose β = 0.17, which is in agreement with (12). With this value of β, we have one positive and stable root for η = 1, as well as for η ∈ [0, 1]. If we consider instead the inequality (6), and we choose β = 0.09 as in [10], we will have three positive equilibrium points E * 1 , E * 2 and E * 3 . The first equilibrium point is stable for η = 1, and for η ∈ [0, 1]. For the second equilibrium point, if η = 1, ∆ < 0, the discriminant is always positive, i.e., T 2 − 4∆ > 0, and both eigenvalues are real. However, one is positive and the other is negative. Trajectories approach the steady state along the eigenvector corresponding to the negative eigenvalue, but move away along the eigenvector corresponding to the positive eigenvalue. The steady state has one stable and one unstable direction. It is therefore unstable and called a saddle [13]. When η ∈ [0, 1], this equilibrium point is stable, according to inequality (15). Regarding the third equilibrium point, its stability is guaranteed for both η = 1 and η ∈ [0, 1]. While remaining in inequality (6), according to Descartes' rule of signs, it is possible to have a unique positive root. So, taking β = 0.03, for η = 1, the equlibrium point E * 4 has complex conjugate eigenvalues with a positive real part, λ = T /2. The steady state is unstable. Due to the presence of a nonzero imaginary part, perturbations grow in an oscillatory manner and spiral away from the steady state. The steady state is an unstable focus. If η ∈ [0, 1], E * 4 is found to be unstable, according to (15). These results confirm the paramount role that the toxic coefficient plays on the stability of the system. In addition, it would be interesting to note that this stability, depending on the system parameters, can either be modified or maintained.
3 Analysis of a system with toxic substances in a subdiffusive regime We now proceed to analyze the case we are interested in, namely the modifications of Turing pattern induced by toxic substances in presence of strong memory effects in the transport process. In order to appreciate these modifications, it is convenient to use the approach developped by Hernández et al [8], because it can help us to see every new results that will emerge. For this purpose, we start by linearizing around the fixed points. Eq. (1) is first written as follows where a ij are the elements of the Jacobian matrix given by Eq. (11). Applying Fourier and Laplace transforms where U and V are given by Solving Eqs. (17) yields The time evolution of u(k, t) = L −1 (U (k, s)) is given by the singularities of Eq. (19) which is the zeros of the function In order to analyze the stability of the system near the equilibrium point, we solve the equation where s 0 (k) will be studied numerically for different values of the toxic coefficient to determine when instabilities should appear in the system of Eq. (16), and if it is influenced by toxicity. In particular, the conditions for a Turing instability are and Now, let us use the approach proposed by Hernández et al. [8] to have a Turing instability in the case where α = η < 1. For this purpose, we start with the Fourier transform of Eq. (16), i.e., The eigenvalues of the stability matrix [A(k)] are given by where The stability of the system is completely determined by the nature of the eigenvalues. Real and negative eigenvalues yield a stable system, meaning that no Turing instabilities are possible. Meanwhile if both are real and one of them is positive (the other being negative), the system will no longer be stable, and the condition for this are exactly similar to the case η = 1. However, complex roots go with a critical value of η. This critical value η c only exists if T rA(k) > 0, which is given by η c is a critical value of the of the anomalous exponent or bifurcation parameter that separates a regime of stationary Turing patterns from an oscillatory cellular instability, commonly known as a Hopf−Turing bifurcation [6]. The Turing conditions are fullfiled when the stationary homogeneous state (k = 0) is stable, and the system becomes unstable under perturbations with finite wavelength. In other words, when reactions take place in the presence of subdiffusion, the condition (22) can be satisfied in two ways, either or β < α 2 , where Condition (28) cannot be fullfiled for the classical diffusion.
To fulfill the condition in Eq. (23) it is necessary that where D = 1 is the ratio of diffusion coefficients, and either or In Fig. 1, we depict the influence of the toxic substances on the parameter that constitute the anomalous exponent. It is clear that the slight variation of the efficiency of the toxicity modifies this parameter, yielding a modification of the anomalous exponent itself. However, it would be premature to assert, based solely on Fig. 1, that this modification would systematically lead to a modification in the Turing instabilities. In order to decide on this, it would be necessary that one of the two relations (22), or (23) should be satisfied. If at least one of the two is verified, then the relationship between the degree of toxicity and the Turing instabilities will be clearly established, as depicted in Fig. 2 which displays clearly the efficiency of toxicity on the Turing instabilities. These results are a direct consequence of Hernandez's work [8] on subdiffusion, where the existence of Turing structures has been demonstrated. In this work, we will take a few cases to study the phenomenon of subdiffusion under the influence of toxicity.
In the next section, we will try to numerically analyze the different possible scenarios by varying a few parameters, with the aim being to see how memory affects different ecosystems.

Numerical results and discussions
In this section, our numerical predictions are verified via direct numerical simulations of system (1). This system describes an ecosystem where preys produce a substance that is toxic to predators. In the process, suitable initial and boundary conditions have been use along with explicit and implicit schemes for time and centered difference for space derivatives. The fractional derivatives were approximated using a scheme base on the Grunwald-Letnikov definition [15,8].
We provide now some simulations as evidence of the analytical predictions derived in the previous sections. Upon this, the Turing patterns obtained for this set of parameters are depicted in Fig. 3.
In this case, we observe that before taking into account the memory effect (η = 1), there are oscillations both in predators and in preys, which means that predators indiscriminately consume toxic preys. However, if the memory effect is triggered (η < 1), we observe that after a certain time, there is stabilization both in preys and predators. This stabilization can be explained in two ways. Either the predators have known how to recognize the less toxic preys, and therefore consume only those preys whose toxicity is reduced, or the predators completely abstain from ingesting the toxic preys. This behavior has been observed in European starlings which, when their body masses have been experimentally reduced become more willing to eat prey items that have been injected with quinine, which is toxic to birds in high doses [16,17]. It is now clear that memory effects play a key role in ways in which naive predators learn to associate warningly preys with their defenses and remember to avoid them in future encounters. This field underpins an extensive body of evolutionary theory [18][19][20][21][22]. Moreover, it is obvious that memory effects regulate the effect of toxic substances on predators. In other word, predators are affected by toxic substances for a certain time, then, they become used to it, and the ecosystem becomes locally stable as shown in Figs.

(e)-(f ) and (h)-(i).
Example 2. Let a = 0.05, d 1 = 0.1, and the other parameters remain as in the first exemple. As previously, we make sure that condition (23) is satisfied. For this set of parameters, we have the structures of Fig. 4 Despite the toxicity of preys, the cohabitation between predator and prey takes place without any problems, even if the number of predators remains relatively low. However, when we take into account the memory effect, we realize that the number of preys decreases considerably, while that of predators remains relatively constant. By further reducing the fractional order of the derivative, we observe a complete stabilization of the two entities (predators and preys). Albeit the drastic fall in their respective population, no species will disappear due to the effect of fractional -order parameter or memory effect. This scenario is similar to the one described in Ref. [23] where predators deliberately choose to swallow toxic preys, choice which is a tradeoff between the benefits of obtaining nutrients and the costs of ingesting toxins. This trade-off is affected by the fact that: animals will consume more toxic preys if they are food-deprived. Indeed, animals face constant decisions about what to eat and what not to eat. While some items are never worth eating, there are many cases where the decision to eat or not should depend on the environment and the individual's current state [24][25][26]. For example, many potential preys available to wild birds are chemically defended, and so contain toxins that will be harmful in the long term or if eaten in excess [27,28]. However, such preys also contain valuable nutrients. In such cases, having lower energy reserves or poorer foraging prospects shifts the balance of costs and benefits in favour of consumption [26]. Eq. (23) is checked. The resulting structures are shown in Fig. 5. This case presents some similarities with the one depicted in example 2. The preys are not really affected by the predators, which is similar to the case discussed in Ref. [29] Here, the first conspicuous mutants would have to survive greater levels of predation than previously thought, because even when the learning process is complete, educated predators may still be prepared to eat aposematic prey. There may be another reason for this phenomenon. Indeed, a predator's ability to moderate and process toxins would be a key factor in limiting attack rates on chemically defended preys, and one that could have significant implications for the survival advantage of being aposematic [29].  (23) is also met. The resulting structures of Fig. 6. reveal that predators and prey may be simultaneously exposed to toxins. Since the natural dispersive force of the movement of each species is weak ( d 2 d 1 < 1), we can conclude that this configuration corresponds to environmental toxins exposure. Indeed, if toxin came from either preys or predators, it could not be affected as much. This means that toxin comes from the environment. The impact of environmental toxins on predator-prey dynamics has recently been investigated [30]. Since mobility is reduced, we are in a situa- tion of confinement in a toxic environment. The direct effects of toxins typically reduce organism abundance by increasing mortality or reducing fecundity. Such direct effects, therefore, alter both bottom-up food availability and top-down predatory ability. However, the indirect effects, when mediated through predator-prey interactions, may lead to counterintuitive effects. Environmental toxins also reduce population variability by preventing populations from fluctuating around a coexistence equilibrium [30] as depicted in Figs. 6 (f ) and (i).

Conclusion
The main objective of this work, was to answer the question whether toxic coefficient β, in the presence of memory can still be considered as an important parameter that can give rise to complex spatial dynamics. To this end, we have examined the influence of toxic substances on the dynamics properties of the fractionalorder ODE system (2), then we have analyzed the system describing an ecosystem where prey produces substance that is toxic for predator in the subdiffusive regime, and give general conditions for toxicity coefficient allowing Turing instability to appear. We have numerically studied in the generic model equations (1). Particularly, investigation have been conducted with emphasis on the effect of subdiffusion on pattern for-  mation in a context where preys release toxin to prevent attacks from predators. Insisting on the formation of Turing patterns. It appears clearly that the effect of toxic coefficient can either be altered by the memory effects by canceling existing Turing structures or by creating new ones. Another interesting phenomenon that we have observed is that extinction does not occur at any time, which suggests that memory allows the stabilization and preservation of ecosystems. We could not finished this work whitout mentioning the importance of the parameter in the ecosystem. Each scenario is reallistic, one has to adequately choose the parameter values. As an illustration, we have proposed four examples. In each of them, by changing the value of certain key parameters, the system configuration changed. In each configuration obtained, memory plays more of a stabilizing role while preserving the ecosystem, regardless of the value of the toxic coefficient.