SPATIOTEMPORAL COMPLEXITY OF A DIFFUSIVE PLANKTONIC SYSTEM WITH PREY-TAXIS AND TOXIC EFFECTS∗

In this paper, we propose a three-species reaction-diffusion planktonic system with prey-taxis and toxic effects, in which the zooplankton can recognize the nontoxic and toxin-producing phytoplankton and can make proper response. We first establish the existence and stability of the unique positive constant equilibrium solution by utilizing the linear stability theory for partial differential equations. Then we obtain the existence and properties of nonconstant positive solutions by detailed steady state bifurcation analysis. In addition, we obtain that change of taxis rate will result in the appearance of time-periodic solutions. Finally, we conduct some numerical simulations and give the conclusions.


Introduction
Plankton are drifting organisms that live in the surface layers of the ocean [24]. According to the trophic modes, plankton can be divided into autotrophic phytoplankton and herbivorous zooplankton. They are important in the ocean's food chain. The phytoplankton can consume carbon and produce oxygen [5]. The zooplankton are the main source of food for almost all fish larvae as they switch from their yolk sacs to catching prey. Basking sharks and blue whales feed on them directly; other large fish feed on them indirectly, by eating fish of smaller size, such as herrings. So plankton play a significant role in ecological environment, fishery economy, tourist industry, and so on [12]. And it is important to study the interactive relationship between phytoplankton and zooplankton.
Yet in another sense, some of the phytoplankton can produce toxin and contaminate seafood or kill fish [10]. From the field-collected samples and mathematical modelling method, Chattopadhayay et al. [3] concluded that the toxic substance plays one of the important role on the growth of the zooplankton population and have a great impact on phytoplankton-zooplankton interactions. They formulated the new model by means of ordinary differential equations for the first time: where P is the density of toxin-producing phytoplankton and Z is the density of zooplankton population. All the coefficients are positive constants, r and K are the intrinsic growth rate and environmental carrying capacity of toxic phytoplankton population respectively, α denotes the rate of predation of zooplankton on toxic phytoplankton population, β denotes the ratio of biomass consumed by zooplankton for its growth, µ is the mortality rate of zooplankton due to natural death as well as due to higher predation, and θ is the rate of toxin liberation by toxic phytoplankton population. Further, f (P ) represents the predation response function and g(P ) represents the distribution of toxic substances. By choosing different functional forms, it was concluded that the toxic chemicals can result in the termination of planktonic blooms. Moreover, various extended models of (1.1) have been proposed and investigated. For instance, Wang et al. [28] investigated the delay-induced fluctuation phenomenon of plankton population by Hopf bifurcation analysis; Sharma et al. [20] studied the diffusion-driven instability; Upadhyay et al. [25] explored the joint effects of toxin production and spatial heterogeneity and displayed the complex spatiotemporal patterns. More related results can be found in [1,13,19,27,30,33].
However, the previous researches only considered two planktonic populations: the phytoplankton and zooplankton. In fact, the nontoxic phytoplankton population and toxin-producing phytoplankton population have different growth habits and roles in the marine system. It is necessary to consider three interacting components consisting of the nontoxic phytoplankton, toxin-producing phytoplankton and herbivorous zooplankton. So, Chattopadhyay et al. [4] extended (1.1) to a more realistic situation of three populations as follows: (1.2) In model (1.2), P 1 (t) and P 2 (t) denote the concentrations of nontoxic phytoplankton and toxin-producing phytoplankton at time t respectively, Z(t) denotes the concentration of zooplankton at time t. The two phytoplankton populations have the same environmental carrying capacity K and don't compete with each other. The grazing pressure of zooplankton is assumed to be reduced by the presence of the toxic phytoplankton P 2 . The authors finally concluded that toxin-producing phytoplankton may be used as a bio-control agent for the harmful algal bloom's problems.
Generally, nontoxic phytoplankton and toxin-producing phytoplankton share the common survival resources, such as the carbon dioxide, inorganic salt, and so on. Combining with the assumption that the zooplankton is able to recognize the two phytoplankton populations and tries to avoid the area filled with the thick toxic phytoplankton species, Banerjee and Venturino [2] further proposed a novel mathematical system to model this situation: Here, the two phytoplankton populations have different environmental carrying capacities and compete with each other. The toxic phytoplankton can kill zooplankton and the zooplankton decreases its consumption. Banerjee and Venturino [2] analyzed the stability of equilibria and numerically showed the insurgence of brown tides. Based on model (1.3), the delayed model and reaction-diffusion model were investigated in [17] and [23], respectively. As the plankton can freely float near the surfaces of all aquatic environments and the phytoplankton have the group defense against zooplankton, it is more realistic to use the reaction-taxis-diffusion equation to describe this phenomenon. It should be noted that although many significant results on reaction-diffusion planktonic models have been derived, very few of them consider both three populations and prey-taxis effect (see, for instance, [8,11,18,32]). The spatiotemporal dynamics of predator-prey model with prey-taxis has been increasingly studied, more details can be found in [15,21,26,29,31]. Motivated by the fore-mentioned work, the aim of this paper is to investigate the effect of prey-taxis on the three-component planktonic system.
The rest of this paper is organized as follows. In Section 2, we give some basic assumptions and build our model. In Section 3, we give the sufficient condition for existence of the unique positive constant equilibrium solution, and then analyze the linear stability of the positive equilibrium solution. In Section 4, we establish the existence of nonconstant positive solutions induced by steady state bifurcation. In Section 5, we investigate the type of steady state bifurcation and stability of bifurcation solutions. In Section 6, we prove the existence of spatially inhomogeneous and time-periodic solutions induced by prey-taxis. In Section 7, we numerically study the diffusive system and illustrate the theoretical results. Finally, Section 8 is devoted to the conclusions.

Model formulation
Let P (x, t), T (x, t) and Z(x, t) denote the densities of nontoxic phytoplankton, toxin-producing phytoplankton and zooplankton at location x and time t, respectively. Throughout this paper, we make the following assumptions.
(i) For simplicity, we only consider the one-dimensional space Ω = (0, L), where L denotes the depth of the water column.
(ii) No plankton species is entering or leaving the column at the top or the bottom, which can be described by homogeneous Neumann boundary condition.
(iii) In the absence of zooplankton, the growth of both phytoplankton follows the logistic law. The two species of phytoplankton compete for a limited number of nutrients.
(iv) Zooplankton ingest both kinds of phytoplankton, the nontoxic phytoplankton can promote the growth of zooplankton, while the toxin-producing phytoplankton have a negative effect.
(v) The predation response function of zooplankton on nontoxic phytoplankton is linear, and the corresponding response function on toxin-producing phytoplankton is of Monod-Haldane.
(vi) Zooplankton mortality in plankton population models is often represented by the so-called closure term. The quadratic form has a specific rate dependent on the zooplankton biomass itself. This may be interpreted as representing either cannibalism within the zooplankton, or a predator whose biomass is proportional to that of the zooplankton [9,22]. As such, the quadratic term is adopted to describe the total mortality of zooplankton.
(vii) The zooplankton can recognize the helpfulness or harmfulness of phytoplankton and make appropriate response. More specifically, zooplankton tend to the high density regions of nontoxic phytoplankton, and keep away from the high density regions of toxin-producing phytoplankton. This kind of directional movement is the so-called prey-taxis effect.
From above assumptions, we propose the following reaction-taxis-diffusion model to describe the interactions of the three planktonic populations: (2.1) Here, ∆ is the usual Laplace operator, and ∇ is the gradient operator. In system (2.1), all the coefficients are positive constants and their meanings can be seen in Table 1.
To our knowledge, there is no result published for this novel model in existing literatures. In what follows, we discuss the stability of positive constant equilibrium solution, existence and properties of nonconstant positive steady state solutions, existence of Hopf bifurcation periodic solutions, and so on.

Symbol
Ecological interpretation The self-diffusion rates of three planktonic populations ξ The nontoxic-phytoplankton-tactic sensitivity coefficient χ The toxic-phytoplankton-tactic sensitivity coefficient The intrinsic growth rate of nontoxic phytoplankton The intrinsic growth rate of toxic phytoplankton The environmental carrying capacity of nontoxic phytoplankton The environmental carrying capacity of toxic phytoplankton The competition coefficients The maximal ingestion rate of zooplankton on nontoxic phytoplankton The maximal conversion rate into zooplankton The maximal ingestion rate of zooplankton on toxic phytoplankton The death rate of zooplankton due to toxic phytoplankton γ The semi-saturation constant d The loss rate of zooplankton due to natural mortality and higher predator

Linear stability of positive constant equilibrium solution
In this section, we mainly analyze the stability of positive constant equilibrium solution. We first establish the existence of the uniqueness positive constant equilibrium solution E * = (P * , T * , Z * ) by mathematical analysis. Consider the following algebraic equations of P, T and Z: From the third equation of (3.1), we have .
Substitute (3.2) into the first two equations of (3.1) can lead to and Then, by substituting (3.4) into (3.3), we can get the following quintic equation of T : where Obviously, a 4 > 0 implies a 2 > 0. In this case, equation (3.5) has the unique positive root T * if a 0 < 0 based on the Descartes' rule of signs. Moreover, if γ+T * 2 is positive, then we can obtain the unique positive equilibrium solution E * .
For convenience, we always assume that the conditions of Proposition 3.1 are satisfied later. Next, we will conduct the stability analysis of the constant equilibrium solution E * . Linearizing system (2.1) at E * = (P * , T * , Z * ), we can obtain the following linear system: Then the characteristic equation is For convenience, some assumptions are proposed and will be useful in the following analysis.
Specifically, system (2.1) will undergo steady state bifurcation and have spatially inhomogeneous equilibrium solution if C k changes from positive to negative; system (2.1) will undergo Hopf bifurcation and have time-periodic solutions if A k > 0, B k > 0 and A k B k − C k changes from positive to negative.
For k ∈ N + , solving the equations C k = 0 and A k B k − C k = 0 respectively, we have where Obviously, H 2 and H 4 are both positive if Assumptions 3.2 and 3.1 are satisfied, respectively. Moreover, for any k k . Therefore, we get the conclusions on the linear stability of equilibrium solution E * as below.

Theorem 3.1. If Assumptions 3.1 and 3.2 hold, then the equilibrium solution E * is locally asymptotically stable when
From Theorem 3.1, it can be found that system (2.1) may undergo steady state bifurcation or Hopf bifurcation when χ > χ 0 . But the bifurcation type depends on whether χ 0 is from steady state bifurcation point set or Hopf bifurcation point set.
If χ 0 = χ S k0 < min k∈N + χ H k , then characteristic equation (3.6) has a zero root for k = k 0 and χ = χ 0 . Since E * is unstable when χ > χ 0 , characteristic equation (3.6) has at least one root with positive real part for k ̸ = k 0 and χ = χ S k , in which case Hopf bifurcation doesn't occur. If ) has a negative real root and a pair of purely imaginary roots for k = k 1 and χ = χ H k1 . And in this case, Hopf bifurcation may occur.
If χ 0 = χ S k0 = χ H k1 , then equation (3.6) has two zero roots. At this point, a Hopf/steady-state mode interaction may occur, and more details can be seen in [16]. We leave it for future research. Hence, we assume that χ S k ̸ = χ H k for k ∈ N + in the discussion that follows.

Existence of nonconstant positive steady state solutions
Next, we shall focus on the existence of nonconstant positive steady state solutions for system (2.1).
To be more precise, we will study the nonconstant solutions to the following stationary system: where P , T and Z are functions of x, ′ denotes the derivative with respect to x, and all the parameters are the same as those in (2.1). Then E * = (P * , T * , Z * ) is still the constant equilibrium solution of (4.1). To study the bifurcation problem of (4.1), we introduce the following notations: It is easy to find that F(P * , T * , Z * , χ) = 0 for any χ > 0 and F : X ×X ×X ×R → Y × Y × Y is analytic. Furthermore, for any fixed point (P ,T ,Ẑ) ∈ X × X × X , the Fréchet derivative of F is where To establish the existence of steady state bifurcation, we first need to verify the necessary condition N DF(P ,T ,Ẑ, χ) ̸ = {0}, where N represents the null space.
Proof. By the Crandall-Rabinowitz local bifurcation Theorem [6], we still need to verify where R is the range of the operator. Assume that there exists a nontrivial solution which satisfies Multiplying both sides by cos kπx L and integrating by parts over (0, L), we get  The coefficient matrix is singular when χ = χ S k , which results in contradiction. The proof is completed.

Stability of steady state bifurcation
In this section, we further study the stability of spatially inhomogeneous solution (P k (s, x), T k (s, x), Z k (s, x)) obtained in Theorem 4.1. We have the following expansions: where (φ i , ψ i , γ i ) ∈ Z and κ i (i = 1, 2) are constants. If κ 1 ̸ = 0, then the sign of κ 1 determines the stability of bifurcating solutions; if κ 1 = 0, then the bifurcation branch Γ k is of pitch-fork type and the sign of κ 2 determines the stability of bifurcating solutions for κ 2 ̸ = 0, and so on. Then we compute the values of κ 1 and κ 2 . For convenience, we first give some useful formulae: For x ∈ (0, L), substituting (5.1) into stationary system (4.1) and collecting the terms of s 2 , we have Then multiplying both sides of (5.2)-(5.4) by cos kπx L and integrating them over (0, L) respectively, we can obtain Moreover, from (φ 1 , ψ 1 , γ 1 ) ∈ Z, we also get We combine (5.5), (5.6) and (5.8) in the following system Then the determinant of the coefficient matrix of (5.9) is Hence the coefficient matrix of (5.9) is singular and the equations only have zero solution. It follows that κ 1 = 0 and the bifurcation branch obtained in Theorem 4.1 is of pitch-fork type.
Similar to the previous steps, substituting (5.1) into (4.1) and collecting the terms of s 3 , we have We also have φ 2 Integrating both sides of (5.10) and (5.11) over (0, L), we obtain Integrating both sides of (5.2)-(5.4) over (0, L) may lead to Again multiplying both sides of (5.2)-(5.4) by cos 2kπx L and integrating them over (0, L), we get Proof. For every k ∈ N + , linearizing (4.1) at (P k (s, x), T k (s, x), Z k (s, x), χ S k ), we obtain the following eigenvalue problem: where (P, T, Z) ∈ X × X × X . Then P k (s, x), T k (s, x), Z k (s, x), χ S k is asymptotically stable if and only if each eigenvalue λ(s) has negative real part. If s → 0, then λ = λ(0) is the simple eigenvalue of DF(P * , T * , Z * , χ S k ) = λ(P, T, Z), i.e., 21) has one-dimensional eigenspace N (DF(P * , T * , Z * , χ S k ))= {(D k , F k , 1) cos kπx/L}. Multiplying both sides of (5.21) by cos kπx/L and integrating them over (0, L), we can obtain that λ = 0 is an eigenvalue of the following matrix: From the proof of Theorem 3.1, for any positive integer k, the above matrix always has an eigenvalue with positive real part for Based on the standard perturbation theory in [14], for s being small, there exists an eigenvalue λ(s) to the linearized problem above that has a positive real part. Therefore the bifurcation branch P k (s, x), T k (s, x), Z k (s, x), χ S k is unstable for s ∈ (−δ, δ).

Existence of time-periodic solutions
Here, we establish the existence of time-periodic solutions to system (2.1) induced by Hopf bifurcation. In this section, we assume that χ H k < χ S k for any k ∈ N + . If χ=χ H k and B k > 0, then A k B k = C k and equation (3.6) can be reduced to whose roots are Then we claim that d Re λ dχ χ=χ H k > 0. In fact, regarding λ = λ(χ) and differentiating both sides of (3.6) with respect to χ, we get Hence, Therefore, we can draw the conclusion on the existence of time-periodic solutions to system (2.1).

Numerical simulations
In this section, we mainly verify the theoretical results obtained in the previous sections by some numerical examples.

Conclusion
In this paper, we have proposed a nontoxic phytoplankton-toxic phytoplanktonzooplankton reaction-diffusion system with prey-taxis effect. The model is based on the fact that the zooplankton can recognize the two kinds of phytoplankton and make proper response, which has generalized some existing mathematical models. By regarding the tactic coefficient as bifurcation parameter, we find the crucial threshold value χ 0 , below which the constant equilibrium solution is stable, and above which the equilibrium solution is unstable. We also find some more complex spatiotemporal dynamics when the tactic coefficient is larger than the threshold value, such as the spatially homogeneous or inhomogeneous time-periodic oscillations, spatially heterogeneous patter formation, and so on. Main results reveal that prey-taxis effect plays a vital role in the spatiotemporal behavior of the planktonic system.