Taxis-Driven Pattern Formation in Tri-Trophic Food Chain Model with Omnivory

: The spatiotemporal dynamics of a three-component model of a food web are considered. The model describes the interactions between populations of resources, prey, and predators that consume both species. It assumes that the predator responds to the spatial change in the resource and prey densities by occupying areas where species density is higher (prey-taxis) and that the prey population avoids areas with a high predator density (predator-taxis). This work studies the conditions for the taxis-driven instability leading to the emergence of stationary patterns resulting from Turing instability and autowaves caused by wave instability. The existence of nonconstant positive steady states for the system is assessed through a rigorous bifurcation analysis. Meanwhile, the conditions for the existence of both types of instabilities are obtained by linear stability analysis. It is shown that the presence of cross-diffusion in the system supports the formation of spatially heterogeneous patterns. For low values of the resource-tactic and predator-tactic coefficients, Turing and wave instabilities coexist. The system undergoes only Turing instability for high levels of these parameters.


Introduction
An ecological community is a part of an ecosystem containing different species of organisms interacting in space and time.The need to consider the spatial organization of a community along with its trophic structure is recognized in mathematical ecology and has been the subject of analysis in the modeling of ecological systems over the past decades.The interplay between trophic and spatial heterogeneity might induce contrasting effects depending on the internal dynamics of the system [1,2].
Among the mathematical tools for studying spatiotemporal processes in ecosystems are partial differential equations (PDEs).PDE models allow the description of numerous fundamental population processes such as dispersal, ecological invasions, dispersalmediated coexistence, and emergence of spatial patterns [3].Skellam's analogy [4] between the movement of molecules and the random movements of individual members of biological species gave impetus to the use of reaction-diffusion models in theoretical ecology.
Diffusive predator-prey models assume that species move randomly in their habitats.However, in real predator-prey communities, directed movements of predator and prey populations often occurs, which is usually shown by the predator pursuing the prey and the prey escaping from the predator.The spatial behavior of the predator moving toward the gradient direction of the prey distribution is classified as prey-taxis.Predator-taxis means that the prey moves opposite to the gradient of the predator distribution.
Starting with the classical model [5], directed movement of predators has attracted much attention in recent years and has inspired numerous works on population dynamics models of one-predator and one-prey systems with prey-taxis [6][7][8][9].These studies analyzed the prey-taxis effects on pattern formation when a spatially homogeneous equilibrium loses its stability.In [10], necessary conditions for pattern formation with diverse prey birth rates, predator mortality rates, and functional responses were investigated.The authors found that with an increase in the taxis coefficient, spatial structures disappeared, and the system stabilized.However, the taxis role as a stabilizing factor is unobvious.It is generally accepted that the ability for self-organization of a system is determined by nonlinear trophic functions.However, [11] provided a method for modeling spatial patterns emerging in trophic systems due to the directed movements of the consumer species.
The influence of taxis on system dynamics has also been studied for three-species systems.In [12], a model of a two-predator-one-prey ecosystem with cross-diffusion and a prey defense mechanism was considered.The authors obtained conditions for the emergence of stationary inhomogeneous structures with an increasing taxis coefficient.The work [13] also modeled a system in which two predators consume one prey, and the preytaxis of both predators can be either positive or negative.Negative taxis was interpreted as a predator's avoidance of areas where the prey population accumulates because prey can defend themselves as a group (see [14,15]).The authors established that the formation of stationary and periodic structures is possible when the taxis of one predator is positive and that of the other is negative.
The interactions among more than two species are modeled by food chain models.Various types of food chain models consisting of three species are available: one of which is a food chain with omnivory.A special kind of omnivory known as intraguild predation [16] is characterized by predators that feed on one or more of the resources of their own prey.This approach can be applied to model a tritrophic food chain, such as a planktivorous fish that feeds on planktonic food, including phytoplankton and zooplankton, or herbivorous and carnivorous zooplankton [17,18].A similar model can describe the interaction between phytoplankton, microzooplankton, and mesozooplankton.
This study continues the analysis of the nonspatial model [19] proposed by the author and takes into account the spatial movements of species.The presented model considers relationships between the resource, prey, and a predator that feeds on both species.This model differs from previous ones in that the predator's functional response depends on the density of one species.The use of a one-prey functional responses allows us to model the various feeding behaviors of predators that can shift from one feeding mode to another.The previous study revealed that the proposed model with self-limitation of the prey or predators and saturated consumption of them can reproduce coexistence between the three species in a resource-rich environment.Additionally, the model considers movement of the predator toward the gradient direction of both species' distributions and the anti-predator behavior of prey moving away from the direction with a high predator density.Recently, Wu et al. [20] investigated the global existence of and the boundedness in a reactiondiffusion model with predator-taxis and showed that the presence of predator-taxis can induce the disappearance of spatial patterns.Fuest [21] investigated global solutions and convergence near homogeneous steady states in a spatial predator-prey model with predator-and prey-taxis.The spatiotemporal patterns in a predator-prey model including predator-and prey-taxis were studied by Wang et al. [22].
A number of works [13, 23,24] studying the effect of taxis on the model dynamics state the stabilizing role of this factor in the case of positive values and the possibility of pattern formation in the case of negative ones.The purpose of this study is to determine the effects of positive prey-taxis and negative predator-taxis on the spatiotemporal dynamics of the system and pattern formation.The reason for the Turing pattern occurrence is the existence of nonconstant positive steady states of a model as a result of diffusion.Stationary structures can also be formed under the influence of taxis pressure.The system's ability to self-organize into spatial patterns forms a mechanism of species coexistence.Through biodiversity conservation, therefore, the influence of taxis pressure on the existence of nonconstant positive solutions is of interest as the coexistence mechanism.
The remainder of this paper is organized as follows: Section 2 introduces a model of a food chain and its dimensionless form.In Section 3, the linear bifurcation analysis is carried out and the existence of a nonconstant positive equilibrium is investigated.In Section 4, the conditions for Turing and wave instabilities are obtained.Pattern formation is numerically illustrated in Section 5.

Model Formulation and Preliminary Analysis
The model consists of three coupled partial differential equations describing the changes to the densities of the resource R(x, t), prey N(x, t), and predator P(x, t) at time t and position x: ( Here, Ω = [0, L], and a homogeneous Neumann boundary condition is imposed to describe an enclosed domain; D i are the diffusion coefficients.The term ∇(Nχ P ∇P) shows the tendency of prey to move away from the high gradient of predator density with taxis rate χ P ≥ 0, and the term ∇(P∇(χ R R + χ N N)) shows the tendency of predators to move towards the direction of resource and prey gradients with rates χ R and χ N , respectively.The nonlinear functions f i (R, N, P) represent the interactions between the predators, prey, and resources as follows: ( Here, r is the maximum growth rate of the resource population, K is the resource carrying capacity, µ i and K i are the predation rates and half-saturation constants of prey and predator; α i are the assimilation coefficients, and m i and δ i are the rates of natural mortality and intraspecies predation of prey and predator, respectively.All the parameters are assumed to be positive.
Introducing the dimensionless variables and parameters and dropping the over-bars, we obtain the nondimensional model as Here, Ω = [0, 1] and The spatially homogeneous non-negative steady states of System (3) are given by: E 0 = (0, 0, 0), E 1 = (1, 0, 0), E 2 = (u 2 , v 2 , 0), and E 3 = (u 3 , 0, w 3 ), where and u 2 , u 3 are positive roots of the following equations: The positive steady state is given by E * = (u * , v * , w * ), where v * is a positive root of the equation: and u * is a positive root of the following equation: Proposition 1.We have the following results: 1.
There exist equilibria E 0 and E 1 for any parameters.

Classical Bifurcation Analysis
The main purpose of this section is to investigate the bifurcation structure of positive solutions to System (3) by regarding χ 3 as a bifurcation parameter and fixing all other parameters.

Linear Analysis
Here, W = (u, v, w) T are small spatiotemporal perturbations away from E * = (u * , v * , w * ), and The Jacobian matrix evaluated at E * is: where Since Ω = [0, 1], the solutions for System (12) can be expanded into a Fourier series [25]: where W k are the Fourier coefficients, λ k denotes the growth rate of the solution, and k is the wavenumber corresponding to mode number n. Substituting ( 14) into ( 12) leads to the following system: The characteristic equation for the above system can be written as: where ).According to the Routh-Hurwitz criteria, the equilibrium point E * is stable with respect to spatially homogeneous and nonhomogeneous perturbations if the following conditions hold: The local stability of the positive equilibrium E * for the homogeneous model is a necessary condition for Turing and wave instability.Let p 2 > 0, p 0 > 0, p 2 p 1 − p 0 > 0; then (16) are satisfied for k = 0, and E * is the stable homogeneous stationary state.Moreover, r 2 (k 2 ) > 0. The equilibrium loses its stability if one of the two remaining conditions in ( 16) is violated.If r 0 (k 2 ) < 0 for some k 2 > 0, then Equation (15) has one positive real root: this implies that Turing instability occurs, and there exist nonconstant positive steady states for System (3).If r 2 (k 2 )r 1 (k 2 ) − r 0 (k 2 ) < 0 for some k 2 > 0, then Equation (15) has a pair of complex roots with a positive real part: this means that wave instability occurs.

Existence of Nonconstant Positive Steady States
In this subsection, the existence of a nonconstant interior equilibrium is considered.In the following, the bifurcation theory of Crandall-Rabinovitz [26] and its application to the bifurcation analysis of elliptic systems developed in [27] are applied to prove the existence of a Turing bifurcation with χ 3 as the bifurcation parameter.I study a positive nonconstant solution for the following system: Here, prime denotes differentiation with respect to x. System (17) can be rewritten in the following form: where and Then, for fix χ 3 ∈ R, any solution (u, v, w) of ( 18) is a classical solution for (17).
Theorem 1. Suppose that all parameters in (17) are positive and there exists a positive constant equilibrium Here, A i , B j can be found in the proof of the theorem.If there exists some k > 0 such that χ k defined in ( 26) is positive, then there exists a constant ϵ > 0 such that System (17) with smooth functions of s.Furthermore, all nonconstant solutions of (17) must stay on the curve Proof.Let z = (u, v, w) T , z = ( ũ, ṽ, w) T .Then ( 20) is equivalent to Here, and (22) is strictly elliptic, and it satisfies Agmons's condition according to Remark 2.5 of Case 2 with N = 1 in [27].Therefore, ( 22) is a Fredholm operator with zero index thanks to Theorem 3.3 and Remark 3.4 in [27].The necessary condition for the existence of bifurcation is that the null space of D z F (z * , χ)(z) is nonempty, where z * = (u * , v * , w * ) T .By taking z = z * in (20), one has Its null space consists of solutions for the following system: Because Problem ( 24) is linear, we now consider the solution z(x) = (u(x), v(x), w(x)) of ( 24) in the form where Z k = (U k , V k , W k ) are the Fourier coefficients.Substituting these series into (24) gives rise to There exists a nontrivial solution 25) if the coefficient matrix in (25) is singular.For each k > 0, this condition holds if there exists positive χ 3 = χ k , where with where N denotes the null space, and N (D z F (z * , χ k )) = span{z k }, where Then, in accordance with bifurcation from the simple eigenvalue [26], to prove the theorem, it remains to show the validity of the following transversality condition Differentiating (23) with respect to χ 3 , we get Suppose that condition (28) fails.Then, there exists a nontrivial z = ( ũ, ṽ, w) T such that Let z = ∑ ∞ k=0 Zk cos kx, where Zk = ( Ũk , Ṽk , Wk ).Substituting these series into (29) gives rise to We observe that the coefficient matrix in (30) with χ 3 = χ k is singular, while the right-hand side is nonzero: then this leads to a contradiction.Consequently, the transversality condition (28) is verified.Then, according to the Crandall-Rabinowitz bifurcation theory [26], there exist nonconstant positive solutions (u k (s, x), v k (s, x), w k (s, x), χ k (s)) that can be parameterized as where = (0, 0, 0), and z k given in (27).Finally, we require χ k ̸ = χ j for all positive k ̸ = j in order to apply the bifurcation theory from a simple eigenvalue; therefore, which is equivalent to (21).

Conditions for Turing Instability
In this section, the effect of taxis on pattern formation is studied.According to Theorem 1, positive equilibrium E * loses its stability through bifurcation as χ 3 crosses the critical bifurcation value χ k .It is widely known that a reaction-diffusion system with two or more agents with largely different diffusion coefficients can generate spatially heterogeneous patterns.The aim of this study is to investigate the occurrence of stationary patterns as a result of taxis pressure.Therefore, I consider a special case of Model (3) with equal diffusion coefficients for all components: d 1 = d 2 = d 3 = d, and resource-tactic sensitivity depends on prey-tactic coefficient χ 2 = ηχ 3 = ηχ.Turing bifurcation occurs when r 0 (k 2 ) = 0, which is analogous to Note that ( 32) is equivalent to (26).Then System (3) undergoes Turing instability when r 0 (k 2 ) < 0 for some k 2 > 0. If the denominator in (32) is negative (positive), then The following theorem provides sufficient conditions for Turing instability.G 1 (η − η 1 ) < 0 and χ > χ T1 ; 2.

Conditions for Wave Instability
In three-component models, not only can Turing bifurcation occur, but wave bifurcation (or finite wavenumber Hopf bifurcation) can also occur.Wave bifurcation breaks spatial and temporal symmetries, and the system generates patterns that are oscillatory in space and time.This bifurcation occurs when r(k where where Then System (3) undergoes wave instability when r(k 2 ) < 0 for some k 2 > 0. If the denominator in (34) is negative (positive), then r(k 2 ) < 0 is analogous to χ > max The following theorem provides sufficient conditions for wave instability.

Numerical Scheme
I simulated System (4) using the method of lines with an n = 100-point spatial grid for a domain [0, 1], i.e., (u i , v i , w i )(t) = (u, v, w)(i∆x, t), i = 0, . . ., n, and ∆x = h = 0.01.For spatial discretization, the second-order central difference is used for the diffusive terms: where , and the first-order upwind scheme is used for the taxis terms: where ).The result of spatial discretization is the autonomous semi-discrete system: The splitting method [28,29] is used for solving this system.The function , where function G 1 contains the taxis term, and G 2 contains the diffusion and reaction terms.In the time integration, the taxis term is treated explicitly, and the remaining terms are treated linearly implicitly by approximate matrix factorization and operator splitting.Then, given an approximation y = z n at time t n and a step size τ, we compute Here Γ 1 and Γ 2 are approximate evolution operators of G 1 and G 2 , respectively.Specifically, Γ 1 (τ/2, t n )y approximates the solution of the initial-value problem As a result, we obtain ŷ (1) .For approximation of Γ 1 , the Runge-Kutta fourthorder method is used.Then, Γ 2 (τ, t n ) ŷ (1) approximates the solution of the initial-value problem 1) .The right-hand side of the ODE can be written as , where G 2a contains the terms arising from the discretization of diffusion terms, and G 2b contains the reaction terms.A linearly implicit variant of the trapezoidal splitting method was used for approximation of Γ 2 : 3) ).

Spatiotemporal Patterns
Now, I investigate the effect of prey-and predator-taxis on the formation of stationary and oscillatory patterns.In the simulations to be presented below, I have taken the following parameters: For the given set of parameter values, the equilibrium point E * = (0.31, 0.38, 0.05) for the spatially homogeneous model is stable.The initial condition is always the steady-state with small perturbations: u(x, 0) = u * + ξ 1 , v(x, 0) = v * + ξ 2 , w(x, 0) = w * + ξ 3 , where ξ i are small-amplitude random perturbations at one percent from the steady-state E * .
Figure 2 demonstrates the system dynamics for η = 0 and χ = 2 and three variants of χ 1 .When χ 1 = 0, the system undergoes wave instability.The top panel in the figure shows that after a long period of spatially homogeneous oscillations, a standing wave regime is established.With an increase in χ, both instabilities are possible.The system generates standing waves (middle panel) and a stationary structure (bottom panel) when χ 1 = 0.1 and χ 1 = 2, respectively.Figure 3 shows the spatial patterns for η ̸ = 0; in this case, Turing and wave bifurcations are possible.For η = 0.05 and η = 0.5, after a certain incubation period, stationary patterns emerge.When η > 0.61, only Turing bifurcation occurs.Figure 4 demonstrates the system dynamics for η = 1 and η = 2.For η = 1 (χ 2 = 2), the system generates a spatiotemporal pattern oscillating in space and time (top panel).This behavior indicates that nonlinear effects arising from the interaction between the Turing and wave modes cause the pure Turing mode to lose its stability before the wave bifurcation occurs.With an increase in η, the spatial structures achieve equilibrium (bottom panel).The initial condition for all the simulations discussed above is the stable spatially homogeneous equilibrium with small perturbations.Now, consider the system dynamics when taking an unstable spatially homogeneous equilibrium as the initial state.In the simulations to be presented below, I have taken the same set of parameter values as in the previous ones; the only exceptions are µ 1 = 0.5 and µ 2 = 0.2.For the given set of parameter values, the equilibrium point E * = (0.13, 0.34, 0.04) for the spatially homogeneous model is unstable.Oscillatory dynamics of the system without diffusion and taxis are presented in Figure 5a.As diffusion and taxis are incorporated into the system, the system behavior changes as follows.When only prey-taxis is considered, we observe a solution that is oscillatory in time and homogeneous in space (Figure 5b).When predator-taxis is considered, stationary structures are formed.Figure 6 demonstrates the solution to System (3) for χ 3 = 0.1 and three variants of prey-taxis coefficients.

Discussion and Conclusions
This paper presents a model of a three-trophic food chain with an omnivorous predator and with consideration of spatial heterogeneity.The model considers the movements of the predator to be directed towards the areas with high densities of the resource and the prey, and the prey moves opposite to the gradient of the predator distribution.The proposed model provides a way to represent the dynamics of a community in which the predator can use different foraging strategies.This model, which takes into account the movement of the predator towards the gradient direction of both species' distributions, displays qualitatively different behavior of the system as the taxis coefficients change.Moreover, the model allows the assessment of the cumulative effect of prey-taxis and anti-predator mechanisms on the community dynamics.
A linear analysis of the system stability under small spatially inhomogeneous disturbances is carried out.The existence of a nonconstant positive equilibrium as a result of Turing bifurcation is proved.Sufficient conditions for Turing and wave instabilities are obtained.It is shown that the system can lose stability with a high level of prey-taxis.Analysis of the taxis coefficients' effects on the system dynamics reveals the following properties of the system.
If only the prey-taxis is considered, then either Turing or wave instability may occur depending on the resource-tactic rate.When the predator-taxis is taken into account, sta-tionary structures are formed when the resource-tactic coefficient is high.If this parameter is less than the prey-tactic coefficient, then Turing patterns and wave modes may exist.
Numerical simulations reveal that an increase in χ 1 leads to a regime shift in which oscillations in the population abundances change to equilibrium spatially heterogeneous structures.A further increase in χ 1 at a fixed η does not lead to a change in the stationary structures.As η increases, the dynamics of the system achieve equilibrium, but the structures change along with changes to the resource-tactic coefficient.The real parts of the eigenvalues plotted in Figure 7 show that with increasing η, the range of wavenumbers for unstable modes and growth rates of these modes also increase.An increase in χ 1 leads only to growth of Re(λ k ) for the most unstable mode, and the range of unstable wavenumbers practically does not change.From a biological point of view, as η increases, the predator switches its preference to the resource, and the fear-driven behavioral avoidance of predators forces prey competing for the resource and predator-free-area species to converge on a narrower suite of resource that is associated with low predation risk [30].An increase in the predator-taxis sensitivity leads to an increase in the growth rates of the most unstable modes, resulting in a decrease in the incubation time for the structures to emerge, while the spatial distribution pattern remains virtually unchanged.In real communities, such spatial heterogeneity can be caused by various reasons, both physical and biological.The heterogeneity of planktonic communities is initiated by turbulent advection [31] or the behavior of zooplankton in response to the presence of toxic phytoplankton [32].The propagation of spatial waves caused by taxis is observed in bacterial populations responding to changes in the gradient of an attractant or repellent [33].Wave solutions are demonstrated by models of the invasion of species [34] or infections [35].
The numerical simulations demonstrate a variety of system dynamics, including the formation of stable and unstable stationary and oscillatory patterns.When the parameters are selected from the Turing instability area and are close to the wave instability line, the purely stationary Turing structures become modulated (Figure 4).Similar system behavior was established in [36].The stability of Turing structures depends on whether subcritical or supercritical bifurcation occurs; the determination of the nature of bifurcation requires a weakly nonlinear analysis.This analysis can also determine the type of wave-namely, a traveling or standing wave-that arises when a homogeneous steady state goes through that bifurcation (see [37] and references therein).Such nonlinear analysis is beyond the scope of this work, which has focused exclusively on studying the existence of a nonconstant steady state to the system and the patterns appearing at the onset of Turing and wave instability, but it will form the subject of future work by the present author.
As a final comment, I would like to note one more property of the presented model.A system's ability to form patterns is considered one of the mechanisms that allows species to coexist.The anti-predator behavior of prey promotes the formation of stationary structures.The effects of refuges used by prey increase the equilibrium density of the prey population.However, the presented model demonstrates a decline in the prey population with an increase in predator-tactic sensitivity when resource-tactic sensitivity is low.This decline may be caused by the fact that prey aggregated in refuge habitats may strongly compete for limited resources [30].With an increase in the resource-tactic coefficient, the prey population begins to grow, and the predator, whose diet consists mainly of the resource, begins to decline catastrophically.Nevertheless, the system stabilizes, and all species in the community survive.

Theorem 2 .
System (3)  undergoes Turing instability around the uniform steady-state E * = (u * , v * , w * ) if one of the following sets of conditions holds:1.

Figure 7 .
Figure 7.The real parts of the eigenvalues for Equation (15) as a function of k.