# Dynamic behavior analysis of a diffusive plankton model with defensive and offensive effects

## Dynamic behavior analysis of a diffusive plankton model with defensive and offensive effects

Chaos, Solitons and Fractals 129 (2019) 94–102 Contents lists available at ScienceDirect Chaos, Solitons and Fractals Nonlinear Science, and Nonequi...

Chaos, Solitons and Fractals 129 (2019) 94–102

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Dynamic behavior analysis of a diffusive plankton model with defensive and offensive effects Qiuyue Zhao a, Shutang Liu a,∗, Xinglong Niu b a b

The School of Control Science and Engineering, Shandong University, Jinan, 250061, China The School of Electrical and Control Engineering, North University of China, Taiyuan, 030051, China

a r t i c l e

i n f o

Article history: Received 9 April 2019 Revised 3 July 2019 Accepted 12 August 2019

Keywords: Diffusive plankton model Defense of phytoplankton Offense of zooplankton Hopf bifurcation Turing instability

a b s t r a c t This paper investigates the dynamic behavior of a diffusive plankton model with defensive and offensive effects in two cases. For the single compartment model, we ﬁrst derive the suﬃcient conditions for the stability and Hopf bifurcation of coexisting equilibrium, which implies that the changes of defense and offense can cause oscillation of planktonic population. Then the properties of Hopf bifurcation are discussed by center manifold theorem. For the spatially extended model, we obtain the suﬃcient conditions for Turing instability and Hopf bifurcation. It is observed that spatial patterns put in place, under the interaction of diffusion, defense and offense. Finally, some numerical simulations are carried out to support the analytical results.

1. Introduction Plankton, composed of phytoplankton and zooplankton, is an important indicator of the status of marine ecosystems. Phytoplankton is an indicator for measuring water quality, which links to the abundance of phytoplankton and community composition, and zooplankton, as a key component of the marine food webs, is an indicator of warm and cold currents. The rapid increase or decrease of phytoplankton populations may lead to severe ecological problems, such as the mortality of zooplankton and harmful algal blooms. It is therefore important to understand the interaction mechanism between phytoplankton and zooplankton. For this end, the dynamic behavior of plankton has received much attention from the empirical and theoretical points of view in recent decades (see [1–4] and the reference therein). The dynamic behavior of plankton is inﬂuenced by a complex interaction between chemical-physical conditions and biological processes. The former, such as pH [5], temperature [6] and turbulence [7], has a direct effect on the marine ecosystem. Moreover, many results showed that the biological processes including growth, mortality, predation and competition are affected by various factors [8–13]. For example, in Ref. [14], the authors investigate a phytoplankton-zooplankton-ﬁsh model which includes nonlinear ﬁsh harvesting and taxation. In Ref. [15], the local Hopf

bifurcation of a phytoplankton-zooplankton model with gestational delays is addressed, while in Ref. [8] the author studies a nutrientphytoplankton-zooplankton model with spatiotemporal dynamics under the effects of prey-taxis and intratrophic predation. Due to ﬂuctuations of environment, such as earthquakes, epidemics and tsunami, the birth rate, death rate, carrying capacity and other parameters of plankton ecosystem exhibit stochastic disturbances. Hence, many authors have introduced stochastic disturbances into deterministic models to reveal the effect of environmental variability on the population dynamics [16–19]. Particularly, the environmental noise causes a reduction of the Listeria monocytogenes mean value concentration [18] and [19] points that stochastic model can reproduce the phytoplankton distributions observed in real data. It should be noted that the adaptive changes of plankton can affect the population dynamics by genetic evolution, phenotypic plasticity, behavioral selection [20]. The defense of phytoplankton and offense of zooplankton are two important adaptive changes. Zooplankton often drives phytoplankton to adapt to environment by producing the effect of defense, while an increase defense capability of phytoplankton may lead to the occurrence of offense [21]. The factors of defense include physiology (e.g., toxicity, bioluminescence), morphology (e.g., silica shell, colony formation) and behavior (e.g., escape response) [22], which increase digestion time of predator, reduce predator’s motivation to search prey by producing toxic substances and lower the encounter rate with a predator [23]. As with the defense of phytoplankton, the offense of zooplankton

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

95

Table 1 Description of the parameters for model (1). Parameters

Description

r1 (u) r0 c1 K b a(u, v) a0

Intrinsic growth rate of phytoplankton Maximum growth rate Costliness of defense Carrying capacity of environment Handing time Capture rate of zooplankton Maximum capture rate Sensitivity of capture rate to the difference between the trait values Conversion eﬃciency Maximum conversion eﬃciency Costliness of offense Mortality rate of zooplankton

θ  (v ) 0 c2 r2

can enhance immunity against prey toxins, decrease handling time and increase encounter area. This paper proposes a diffusive plankton model with the effects of defense of phytoplankton and offense of zooplankton. Compared with the models in the existing literature, the intrinsic growth rate of phytoplankton and the conversion rate of zooplankton respectively depend on defense and offense in our model. Meanwhile the capture rate of zooplankton to phytoplankton is shared by defense and offense, which extends the construction conditions of biological mathematical model. Moreover, we also study the single compartment model and its spatially extended version by choosing capture rate as the control parameter. The rest of this paper is organized as follows. Section 2 describes the diffusive plankton model and Section 3 discusses the local stability and bifurcation of coexisting equilibrium of the plankton model without diffusion. In Section 4 we investigate the occurrence of Turing instability and Hopf bifurcation of the diffusive model and some numerical simulations are presented to illustrate our theoretical results in Section 5. Finally, we end this paper with some discussions in Section 6. 2. Formulation of diffusive plankton model This section considers the phytoplankton P and zooplankton Z with defense u and offense v to model the defensive and offensive effects as follows:

⎧    a(u, v )Z ⎪ P ⎪ ˙ P = P r ( u ) 1 − − , ⎨ 1 K 1 + ba(u, v )P   (v )a(u, v )P ⎪ ⎪ ⎩Z˙ = Z 1 + ba(u, v )P − r2 Z .

(1)

The biological signiﬁcance of the parameters of the model are described in Table 1 and the functions r1 (u), a(u, v) and  (v ) are assumed as follows: •

A trade-off between defense u and the intrinsic growth rate 2 r1 of phytoplankton is assumed to be r1 (u ) = r0 e−c1 u [24]. Thus the impact of r1 on the dynamics decreases as prey biomass approaches the carrying capacity. When the prey biomass is low, the cost of the same defense level is higher, while when the prey biomass is high, the cost of the same defense level is lower. The capture rate of the predator to the prey is a(u, v ) = a0 θ (u−v ) , which is determined by the levels of defense and 1+e

offense [24,25]. Then, the capture rate decreases with increasing defense and increases with increasing offense. A trade-off between offense v and the conversion eﬃciency 2  is assumed to be  (v ) = 0 e−c2 v [26].

In the ocean, plankton usually distributes inhomogeneously and moves to areas with smaller population concentration because of

diffusion and current. Therefore, one can modify model (1) by introducing diffusion terms to better present the real world situation.

⎧ ∂P (u,v )Z ⎪ − d1 P = P r1 (u ) 1 − KP − 1+aba , ⎪ ( u, v ) P ⎪ ∂t ⎪ ⎨ ∂Z   (v)a(u,v)P  − d2 Z = Z 1+ba(u,v )P − r2 Z , ⎪ ∂t ⎪ ⎪ Px (0, t ) = Px (l π , t ) = 0, Zx ( 0, t ) = Zx ( l π , t ) = 0, ⎪ ⎩ P (x, 0 ) = φ (x, 0 ) > 0, Z (x, 0 ) = ϕ (x, 0 ) > 0,

x ∈ ( 0, l π ), t > 0, x ∈ ( 0, l π ), t > 0,

(2)

t > 0, x ∈ ( 0, l π ),

where d1 and d2 are diffusive coeﬃcients of phytoplankton and zooplankton, respectively. The boundary conditions are Numann boundary conditions (the region is closed), with no phytoplankton and zooplankton species entering and leaving the region at the 2 boundary. The Laplacian operator  = ∂∂x2 ∈ R models that plankton is spatially heterogeneous and migrates to regions of lower population density. 3. Dynamics analysis of single compartment model In the following, we perform the stability and bifurcation analysis for the coexisting equilibrium by calculating the eigenvalues of the Jacobian matrix of model (1), and give the stability of the bifurcated periodic solutions. First, the stability of the coexisting equilibrium is investigated. Apparently, (i) the extinction equilibrium E0 = (0, 0 ) is unstable; (ii) the boundary equilibrium E1 = (K, 0 ) is a saddle point. From the viewpoint of ecology, we are interested in the coexisting equilibrium E∗ = (P∗ , Z∗ ) which is obtained by solving the following set of equations:

⎧   ⎪ ⎨r1 (u ) 1 − KP − a(u, v )Z = 0, 1 + ba(u, v )P  (v )a(u, v )P ⎪ ⎩ − r2 Z = 0. 1 + ba(u, v )P It is easy to check that

Z∗ =

r1 ( u ) P∗ 1− (1 + ba(u, v )P∗ ), a(u, v ) K

then Z∗ > 0, provided that P∗ < K and P∗ is a positive root of the following equation:

S(P ) = b2 a(u, v )2 P 3 + P 2 ba(u, v )(2 − Kba(u, v ))



+ P 1 − 2Kba(u, v ) +

K  (v )a(u, v )2 r2 r1 ( u )

− K = 0.

2 2 Clearly, S(0 ) = −K < 0, S(K ) =  (v )ra(ru,(vu)) K > 0. Hence coexisting 2 1

equilibrium E∗ exists. Now, we will pay attention to the dynamics of model (1) with a coexisting equilibrium E∗ . The Jacobian matrix

96

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

at E∗ is given by

J=

a11 a(u, v )a22

b11 a(u, v )b22 ,

with

r1 ( u ) ba(u, v )2 Z∗ P∗ P∗ + , K (1 + ba(u, v )P∗ )2  ( v )Z∗ = , (1 + ba(u, v )P∗ )2

a11 = − a22

P∗ a(u, v ) , 1 + ba(u, v )P∗ r2 Z∗ =− . a(u, v )

b11 = − b22

The characteristic equation satisﬁes:

λ2 − Tr(J )λ + Det(J ) = 0,

Det(J ) = a(u, v )a11 b22 − a(u, v )b11 a22 .

λ1,2 = α (a(u, v )) ± iω (a(u, v )) =

Tr(J ) ±



1  2

bP∗ Z∗ −bP∗ Q

Q , bP∗ Z∗ −bP∗ Q

we obtain Tr(J ) =



= λ=iω, a(u,v )=a

baZ∗ P∗

(1 + baP∗ )3

r ( u )P

bP∗ Z∗ −bP∗ Q

E∗ is locally asymptotically stable; (iii) if a11 > 0 holds, model (1) undergoes Hopf bifurcation at E∗ when a(u, v) crosses a = √ Q . bP∗ Z∗ −bP∗ Q

Next, we give the stability of bifurcated periodic solution. To this end, we translate E∗ to the origin by the translation x¯ = P − P∗ , y¯ = Z − Z∗ . So, model (1) becomes

⎧   a(u, v )(y¯ + Z∗ ) (x¯ + P∗ ) ⎪ ⎪ − , ⎨x¯˙ =(x¯ + P∗ ) r1 (u ) 1 − K 1 + ba(u, v )(x¯ + P∗ )  (4)  (v )a(u, v )(x¯ + P∗ ) ⎪ ⎪ ⎩y¯˙ = (y¯ + Z∗ ) 1 + ba(u, v )(x¯ + P ) − r2 (y¯ + Z∗ ) . ∗ Then we rewrite (4) into the following form



x¯ y¯



+

f (x¯, y¯ , a(u, v )) h(x¯, y¯ , a(u, v ))

where

f (x¯, y¯ , a(u, v )) =



−b11 a11

0 , −ω

x¯2 2 +

 −

− b1 11 ⎝ = − ωab11

B−1

11

,

(ba(u, v )P∗ − 1 )a(u, v )x¯y¯ 2(1 + ba(u, v )P∗ )2

 x˜˙ y˜˙

 =

−ω 0

0

ω

 x˜ y˜

 +

x˜ y˜

0 − ω1 ⎠.

 =B

x¯ , then (5) becomes y¯

f1 (x˜, y˜ ) , h1 (x˜, y˜ )

(6)

1 f (−b11 x˜, a11 x˜ − ωy˜, a(u, v )), b11 1 h1 (x˜, y˜ ) = − h(−b11 x˜, a11 x˜ − ωy˜, a(u, v )) f1 (x˜, y˜ ) = −

ω

a11

ωb11

f (−b11 x˜, a11 x˜ − ωy˜, a(u, v )),

f (−b11 x˜, a11 x˜ − ωy˜, a(u, v ) )



=

r1 ( u ) 1 ba(u, v )2 Z∗ (1 − ba(u, v )P∗ ) 2 2 − + b11 x˜ 2 K (1 + ba(u, v )P∗ )3

b11 x˜a(u, v )(ba(u, v )P∗ − 1 )(a11 x˜ − ωy˜ ) 2(1 + ba(u, v )P∗ )2

ba(u, v )3 Z∗ (2 − ba(u, v )P∗ )(−b11 x˜)3 3(1 + ba(u, v )P∗ )4

(b11 x˜)2 ba(u, v )2 (1 − ba(u, v )P∗ )(a11 x˜ − ωy˜ ) 6(1 + ba(u, v )P∗ )4 2 ˜ (b11 x ) ba(u, v )2 (3 − ba(u, v )P∗ )(a11 x˜ − ωy˜ ) − 6(1 + ba(u, v )P∗ )3   + O |(−b11 x˜, a11 x˜ − ωy˜ )|4 , −

h(−b11 x˜, a11 x˜ − ωy˜, a(u, v ) )

r1 ( u ) ba(u, v )2 Z∗ (1 − ha(u, v )P∗ ) + K (1 + ba(u, v )P∗ )3

Through the transformation

1

Theorem 1. Assume that (H1) holds and Q = ( 1 K ∗ + r2 Z∗ ) 2 . For model (1), (i) if a11 < 0 holds, then the equilibrium E∗ is locally asymptotically stable; (ii) if a11 > 0 holds, for a(u, v ) < √ Q , then the equilibrium

=J

ba(u, v )2 y¯ x¯2 + O(|(x¯, y¯ )|4 ), 3(1 + ba(u, v )P∗ )3

where

> 0.

This implies that the transversality condition is satisﬁed. Therefore, the result of model (1) follows immediately.

x¯˙ y¯˙

+



0. Since Det(J) > 0 under (H1), (3) has a pair of purely imaginary roots. Choosing a(u, v) as perturbation parameter, we obtain



ba(u, v )2 (1 − ba(u, v )P∗ )x¯y¯ x¯ 6(1 + ba(u, v )P∗ )3

(H1 )

hold. In addition, when a(u, v ) = a = √

d(Reλ(a(u, v ))) da(u, v )

+

When a(u, v ) = a, λ1,2 (a ) = ±iω, ω = a(a11 b22 − b11 a22 ), the one eigenvector corresponding to the eigenvalues λ1,2 (a) is ξ = (−b11 , a11 − iω )T . Set

B=

and

a11 b22 − b11 a22 > 0



ba(u, v )2 (1 − ba(u, v )P∗ )x¯2 y¯ 6(1 + ba(u, v )P∗ )4



Tr(J )2 − 4Det(J ) . 2

Then, it follows that the roots of (3) have negative real parts if and   r ( u )P only if a11 + a(u, v )b22 < 0, a(u, v ) < √ Q , Q= 1K ∗ + r2 Z∗

+

 (v )ba(u, v )2 Z∗ x¯2  (v )a(u, v )x¯y¯ + (1 + ba(u, v )P∗ )3 2(1 + ba(u, v )P∗ )2 r2 y¯ 2  (v )ba(u, v )2 x¯2 y¯ − − 2 3(1 + ba(u, v )P∗ )3  (v )b2 a(u, v )3 Z∗ x¯3  (v )ba(u, v )2 x¯y¯x¯ + − (1 + ba(u, v )P∗ )4 3(1 + ba(u, v )P∗ )3   + O |(x¯, y¯ )|4 .

(3)

The roots of (3) are given by

b2 a(u, v )3 Z∗ (2 − ba(u, v )P∗ )x¯3 3(1 + ba(u, v )P∗ )4

h(x¯, y¯ , a(u, v )) = −

where

Tr(J ) = a11 + a(u, v )b22 ,

(5)

 (v )ba(u, v )2 Z∗ (b11 x˜)2 b11 x˜ (v )a(u, v )(a11 x˜ − ωy˜ ) − (1 + ba(u, v )P∗ )3 2(1 + ba(u, v )P∗ )2 r2 (a11 x˜ − ωy˜ )2

=− −

2

 (v )ba(u, v )2 (b11 x˜)2 (a11 x˜ − ωy˜ )  (v )b2 a(u, v )3 Z∗ (b11 x˜)3 − − 3(1 + ba(u, v )P∗ )3 (1 + ba(u, v )P∗ )4    (v )ba(u, v )2 (b11 x˜)2 (a11 x˜ − ωy˜ ) − + O |(−b11 x˜, a11 x˜ − ωy˜ )|4 . 3(1 + ba(u, v )P∗ )3

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

We rewrite (6) into the following polar coordinate form



r˙ = α (a(u, v ))r + σ (a(u, v ))r + O(|r| ), θ˙ = ω (a(u, v )) + c(a(u, v ))r2 + O(|r|4 ). 3

4

(7)

By Taylor expansion of (7) at a(u, v ) = a, we have



r˙ = α  (a )(a(u, v ) − a )r + σ (a )r 3 + O(|r|4 ), θ˙ = ω (a ) + ω (a )(a(u, v ) − a ) + c(a )r2 + O(|r|4 ).

 1 σ (a ) = f1x˜x˜x˜ + h1x˜x˜y˜ 16   1 + f1x˜y˜ f1x˜x˜ − h1x˜y˜ (h1x˜x˜ + r2 ω ) − f1x˜x˜ h1x˜x˜ , 16ω (a ) 

f1x˜x˜x˜ = − 2

+

b11 ba(u, v )2 (1 − ba(u, v )P∗ ) 3(1 + ba(u, v )P∗ )4

4 (v )ba(u, v ) − 3(1 + ba(u, v ) )

2 2 b11 , P∗ 3

f1x˜y˜



f1x˜x˜ = −b1 −

r1 ( u ) ba(u, v )2 Z∗ (1 − ba(u, v )P∗ ) + K (1 + ba(u, v )P∗ )3

h1x˜x˜

P . Z

(8)

The characteristic equation of (8) subject to the Neumann boundary conditions is

λ2 − Trn λ + Detn = 0,

λ

n 1,2

=



Trn ±

n = 0, 1, 2, · · ·,

(9)

Tr2n − 4Detn . 2

Under (H1) and (ii) of Theorem 1, we have Trn < Tr(J) < 0 and Det(J) > 0. Therefore, Turing instability occurs when it exists n ∈ N such that Detn < 0. Denote

dm = min

Det(J ) − a(u, v )d1 b22 nl 2



1≤n≤m

a11 l 2 d1

then, when

n2

a11 − d1 l 2

n2



(m+1 )2

a(u, v )a11 x˜ ba(u, v )2 P∗ a11 − , (1 + ba(u, v )P∗ )2 (1 + ba(u, v )P∗ )2  a(u, v ) ba(u, v )2 P∗ = −a11 − 2(1 + ba(u, v )P∗ )2 2(1 + ba(u, v )P∗ )2  b11  (v )a(u, v ) − + r a 2 11 , 2(1 + ba(u, v )P∗ )2  a11 r1 ( u ) ba(u, v )2 Z∗ (1 − ba(u, v )P∗ ) =− − + b11 ω K (1 + ba(u, v )P∗ )3 a(u, v )a11 (1 − ba(u, v )P∗ ) + (1 + ba(u, v )P∗ )2  1 2 (v )ba(u, v )2 Z∗ b211 b11  (v )a(u, v )a11 2 + + + r2 a11 . ω (1 + ba(u, v )P∗ )3 (1 + ba(u, v )P∗ )2

According to Wu et al. [27], we know that the bifurcated periodic solution is stable if σ < 0, on the contrary, the bifurcated periodic solution is unstable if σ > 0.

 n2

,

l2

≤ 1 holds, we have l2

a11 d1 2 a11 − d1 nl 2

And when m2 ≤

+

h1x˜y˜



2

b11 ba(u, v )2 (1 − ba(u, v )P∗ ) 2ba(u, v )2 b11 + 3 3(1 + ba(u, v )P∗ ) 3(1 + ba(u, v )P∗ )3

ba(u, v )2 P∗ ω a(u, v )ω =− + , 2(1 + ba(u, v )P∗ )2 2(1 + ba(u, v )P∗ )2

b11 2 a(u, v )b22 + d2 ∂∂x2

and the eigenvalues are given by

2b11 ba(u, v )2 (1 − ba(u, v )P∗ )a11 2ba(u, v )2 a11 b11 + , 3 3(1 + ba(u, v )P∗ ) 3(1 + ba(u, v )P∗ )3



⎞  ∂P ∂2 ⎝ ∂ x ⎠ = a11 + d1 ∂ x2 a(u, v )a22 ∂Z

⎧ n2 n2 ⎪ ⎨Trn = −(d1 + d2 ) l2 + a11 + a(u, v )b22 = Tr(J ) − (d1 + d2 ) l2 , 4 2 Detn = d1 d2 nl 4 − (d1 b22 a(u, v ) + a11 d2 ) nl 2 + a(u, v )a11 b22 − a(u, v )a22 b11 ⎪  2  2 2 ⎩ = nl 2 − ad11 d1 d2 nl 2 − a(u, v )d1 b22 nl 2 + Det(J ), 1

b2 a(u, v )3 Z∗ (2 − ba(u, v )P∗ )b211 (1 + ba(u, v )P∗ )4

h1x˜x˜y˜ = −a11

where

b11 ba(u, v )2 (1 − ba(u, v )P∗ )a11 + 3(1 + ba(u, v )P∗ )4 +

said to be Turing unstable if it is stable in the absence of diffusion, but can becomes unstable when diffusion exists [28]. The linearization of (2) has the form:

∂x

The stability of bifurcated periodic solution is determined by the sign of σ (a), so we have

where

97

n2 l2

a11 d1

a11 d1

( n2 − 1 ) ≥ 0. 

≤ (m + 1 )2 and 1 ≤ n ≤ m hold, 0 < a11 1 −



≤ a11 1 −

n2 m2



is satisﬁed.

Based on the above analysis, we have the following theorem.

Theorem 2. Assume that (H1) and (ii) of Theorem 1 hold. a11 l 2 d1

(i) If (H2)

≤ 1 or (H3) m2 ≤

a11 l 2 d1

≤ (m + 1 )2 , d2 < dm holds,

the coexisting equilibrium E∗ of model (2) is locally asymptotically stable. a l2 (ii) If m2 ≤ 11 ≤ (m + 1 )2 and d2 > dm hold, the coexisting equid 1

librium E∗ of model (2) is Turing unstable. 4.2. Hopf bifurcation Now, we consider the occurrence of Hopf bifurcation, with a(u, v) being identiﬁed as bifurcation parameter. By the presentation of eigenvalues, we know that (9) has purely imaginary roots if

a11 − (d1 + d2 ) nl 2

2

a(u, v ) = an =

−b22

,

and Detn > 0. Thus, Trn (an ) = 0, Trj (an ) = 0, for j = n. Meanwhile,

Detn (an ) > −(d1 b22 a(u, v ) + a11 d2 )



n2 l2

(d1 − d2 )a11 − d1 (d1 + d2 )

>

n2 l2

n2 . l2

4. Dynamics analysis of spatially extended model

If it exists some n ∈ N such that Detn (an ) > 0 holds, we have

In this section, we focus on the diffusion driven instability and the bifurcation analysis for model (2) under condition (H1).

a1 >

4.1. Turing instability

αn (an ) =

2

In this subsection, we will derive suﬃcient conditions for occurrence of Turing instability for model (2). A constant solution is

d1 (d1 +d2 ) n2 d1 −d2

l

and d1 > d2 . In addition, the complex eigenvalues

α n (a(u, v)) ± iωn (a(u, v)) satisfy 

ban Z∗ P∗

(1 + ban P∗ )3

> 0.

This implies that the transversal condition is satisﬁed. Thus, we have the following result.

98

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

Theorem 3. Assume that condition (H1) and (ii) of Theorem 1 hold. If 2

there exists some n ∈ N such that the conditions a11 >

d1 (d1 +d2 ) n2 d1 −d2

l

and

d1 > d2 are satisﬁed, then model (2) undergoes the Hopf bifurcation at a(u, v ) = an . Moreover: (i) The bifurcated periodic solution from a(u, v ) = a is spatially homogeneous, which coincides with the periodic solution of model (1). (ii) The bifurcated periodic solutions from a(u, v ) = an (n ≥ 1) are spatially nonhomogeneous. Next, we shall study the stability of the bifurcated periodic solutions by applying center manifold theorem and normal form theorem of partial differential equations [29,30]. We take a(u, v ) = a (n = 0 ) as the bifurcation parameter. Let P¯ = P − P∗ , Z¯ = Z − Z∗ . For convenience, we drop the bar, then model (2) can be written as



2 d1 ∂∂x2

Ut =

0 U + JU + F (U, a ), 2 d2 ∂∂x2

0

(10)

with U = (P, Z )T and F (U, a ) = ( f (U, a ), h(U, a ))T , f, h are deﬁned in (5). By the straightforward computation, we have

fPP fZP

r1 ( u ) ba2 Z∗ (1 − haP∗ ) =− + , K (1 + baP∗ )3 a =− , (1 + baP∗ )2

fPPP = −

hPZ

fPPZ =

hPZP = −

2 (a )ba , (1 + baP∗ )3

2

Then, we have

1 ∗ q , C (q, q ) z2 + q∗ , C (q, q¯ ) zz¯ 2   1 + q∗ , C (q¯ , q¯ ) z¯2 + O |z|3 , |z|, |w|, |w|3 , 2 H (z, z¯ ) = F (U, a ) − q¯ ∗ , F (U, a ) q − q¯ ∗ , F (U, a ) q¯ ,

2

 ) =

hPPP q21 q¯1 + hPPZ q21 q¯2 + hPZP q1 q2 q¯1 + hPZZ q1 q2 q¯2 fZPP q¯1 q2 q1 + fZPZ q1 q2 q¯2 + fZ Z P q22 q¯1 + fZ Z Z q22 q¯2 hZPP q¯1 q2 q1 +hZPZ q1 q2 q¯2 +hZ Z P q22 q¯1 +hZ Z Z q22 q¯2 h11 = q∗ , C (q, q¯ ) ,

.

g02 = q∗ , C (q¯ , q¯ ) ,

1 a11 − iω , aa22 ω2 + a211

Hence q∗ , q = 1, q∗ , q¯ = 0.

T

,

T .

H02 = C (q¯ , q¯ ) − q∗ , C (q¯ , q¯ ) q − q¯ ∗ , C (q¯ , q¯ ) q¯ . Obviously, the values of H20 , H11 , H02 are 0. On the center manifold C near the origin, we have

(15)

Substituting z˙ of (14) into (15) and comparing the coeﬃcients with w˙ of (14), we obtain

−aa22 U. 2 d2 ∂∂x2 + ab22

iω + a11 1 ,− 2 b11 ω + a211

H11 = C (q, q¯ ) − q∗ , C (q, q¯ ) q − q¯ ∗ , C (q, q¯ ) q¯ ,

w˙ = wz z˙ + wz¯ z¯˙ = w20 zz¯ + w11 z¯z˙ + w11 zz¯˙ + w02 z¯z¯˙ + · · ·.

w20 = (2iω − L )−1 H20 = (0, 0 )T , w11 = −L−1 H11 = (0, 0 )T ,

It has known that the characteristic equation of (11) has a pair of purely imaginary roots  = {iω, −iω} when a(u, v ) = a. Let q and q∗ be the eigenfunctions of L and L∗ corresponding to iω and −iω:

q∗1 , q∗2 T

fPPP q21 q¯1 + fPPZ q21 q¯2 + fPZP q1 q2 q¯1 + fPZZ q1 q2 q¯2

g21 = 2 q , C (w11 , q ) + q , C (w20 , q¯ ) + q∗ , E (q, q, q¯ ) , ∗

(11)



,

H20 = C (q, q ) − q∗ , C (q, q ) q − q¯ ∗ , C (q, q ) q¯ ,

U˙ t = L U.

q = ( q1 , q2 )T =

hPP q21 + hPZ q1 q2 + hZP q2 q1 + hZZ q22

g20 = q∗ , C (q, q ) ,

b11 U, 2 d2 ∂∂x2 + ab22

d1 ∂∂x2 + a11 −b11

fPP q21 + fPZ q1 q2 + fZP q2 q1 + fZZ q22

This, together with (14), implies that 2 3

Obviously, L∗U, V = U, LV for any U = (P1 , P2 )T , V = (Z1 , Z2 )T ∈ H2 ([0, l π ] ) × H2 ([0, l π ] ) and the inner product in H2 ([0, l π ] ) ×  lπ H2 ([0, l π ] ) is U, V = l1π 0 (P1 Z¯ 1 + P2 Z¯ 2 )dx. Thus, the linearization equation of (10) at (0,0) is

q =(

  z2 z2 z¯ z¯2 + g11 z¯z + g02 + g21 + O |z |3 , |z |, |w|, |w|3 , 2 2 2   z2 z¯2 z H (z, z¯ ) = H20 + H11 zz¯ + H02 + O |z |3 , |z |, |w|, |w|3 . 2 2

+

6 ( v )b a Z∗ , (1 + baP∗ )4

and the adjoint operator of L is



(14)

g(z, z¯ ) = g20



2

(13)

where



2

d1 ∂∂x2 + a11 aa22

On the center manifold C, we

z˙ = iωz + q∗ , F (U, a ) = iωz + g(z, z¯ ), w˙ = U˙ − z˙ q − z¯˙ q¯ = Lw + H (z, z¯ ),

C (q, q ) =

and other partial derivatives are 0. The linear operator L = L(a ) of (8) at a(u, v ) = a is



z=

Then, we rewrite (10) into the (z, w) coordinate form

2 (v )ha2 Z∗ =− , (1 + baP∗ )3

hP P P =

(12)

q∗ , U .

z2 z¯2 + w11 zz¯ + w02 + O(|z|3 ). 2 2



ba2 (1 − baP∗ ) , (1 + baP∗ )4

hZZ = −r2 ,

2 (a )ba , (1 + baP∗ )3

w = w20

E (q, q, q¯ ) = hPP

hP P Z = −

L∗ U =

where w = (w1 , w2 have

)T ,

q∗ , F (U, a ) =

ba2 P∗ = , (1 + baP∗ )2

ba (1 − baP∗ ) , (1 + baP∗ )3

2ba2 = , (1 + baP∗ )3  ( v )a = , (1 + baP∗ )2

LU =

U = zq + z¯q¯ + w(z, z¯ ),

2

fPZP = fZPP

2b2 a3 Z∗ (2 − baP∗ ) , (1 + baP∗ )4

fPZ

We compute the coordinates describing center manifold C at a(u, v ) = a. Then the solution of (10) can be written as

w02 = (−2iω − L )−1 H02 = (0, 0 )T , thus g21 = q∗ , E (q, q, q¯ ) . Consequently, each gij (0 ≤ i, j ≤ 2) is determined by the parameters of model (2) and a. Then, the standard result can be computed as

δ1 (a ) =

g20 g11 (3α (a ) + iω (a )) + 2(α 2 (a ) + ω2 (a ))

|g11 |2 α (a ) + iω (a )

|g02 |2 g21 + , 2(α (a ) + 3iω (a )) 2 Re(δ1 (a ))  β2 = Im(δ1 (a )) − ω ( a ). α  (a ) +

According to Tan et al. [31], we know the properties of Hopf bifurcation. If Re(δ 1 (a)) < 0 (Re(δ 1 (a)) > 0), then the bifurcated periodic solutions are stable (unstable). If β 2 > 0 (β 2 < 0), then the period of bifurcated periodic solutions increases (decreases).

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

99

Fig. 1. Behavior of plankton (left panel), and phase plot of model (1) (right panel) for stable conditions. The numerical results are obtained with d1 = d2 = 0, u = 0.2, v = 0.15 and initial condition (0.5, 0.6).

Fig. 2. Behavior of plankton (left panel), and phase plot of model (1) (right panel) for different values of defense and offense. The numerical results are obtained with d1 = d2 = 0, u = 0.2, v = 0.5 and initial condition (0.5, 0.6).

5. Numerical simulations In this section, we carry out some numerical simulations to verify our theoretical results by considering the following parametric values:

r0 = 0.95, c1 = 0.4, K = 6, b = 0.6, a0 = 1.85, c2 = 0.8 ,

0 = 0.5, r2 = 0.1.

θ = 5, (16)

For the above parametric values and u = 0.2, v = 0.15, it is noticed that model (1) has a coexisting equilibrium E∗ = (0.3819, 1.2814 ). Then we choose different capture rate a(u, v) of zooplankton to verify the effects of defense of phytoplankton and offense of zooplankton on the dynamics of (1). In Figs. 1 and 2, we consider a(u, v) as the bifurcation parameter and the Hopf bifurcation threshold is a(u, v ) = a = 0.9787. Hence coexisting equilibrium E∗ is locally asymptotically stable with a(u, v ) = 0.81 < a, that is u = 0.2, v = 0.15 in Fig. 1 and loses stability with a(u, v ) = 1.5125 > a, that is u = 0.2, v = 0.5 in Fig. 2. It suggests that different defense of phytoplankton and offense of zooplankton affect the stability of coexisting equilibrium. From those ﬁgures, we can obtain that stable plankton population density starts decaying with the reduced defense of phytoplankton or increased offense of zooplankton. This means that appropriate increase of phytoplankton defense and decrease zooplankton offense might be beneﬁcial for coexistence of phytoplankton and zooplankton under certain circumstances in deterministic environment. In addition, we obtain σ = −82.7657 < 0. Hence, the bifurcated periodic solution is stable. This is shown in Fig. 2. Based on the importance of stochastic disturbance, we perform numerical simulations to show the effect of Bromnian motion on dynamic behaviour of model (1). In Ref. [32], the author points that any nonlinear system x˙ (t ) = f (x(t ), t ) in Rd can be stabilized by Brownian motion provided |f(x, t)| ≤ K|x| for some K > 0. On the other hand, this system can also be destabilized by Brown-

ian motion if the dimension d > 2. According to Theorem 3.1 and 3.2 of [32], we obtain that stochastic disturbance causes extinction of phytoplankton in a short time (see the left panel of Fig. 3) by using parameter values of (16) and choosing matrix of inten0 1.6 sity of white noise with the same defense and of−1.6 0 fense (u = 0.2, v = 0.15) as Fig. 1. In addition, the white noise can transform model (1) from instability to stability, which means that white noise keeps the densities of phytoplankton and zooplankton moving up and down randomly in a small neighborhood (see

1.2 0 the right panel of Fig. 3), here we choose with the 0 1.2 same defense and offense (u = 0.2, v = 0.5) as Fig. 2. The numerical simulations further reveal that different stochastic disturbances can keep stable coexistence of both populations and may be lead to the extinctions of phytoplankton and zooplankton. Next, we obtain the combined effects of defense and offense in diffusive model (2). We keep other coeﬃcients to be the same as (16). Let l = 3, u = 0.2, v = 0.15. In Figs. 4 and 5, according to the simple calculations, the conditions (H2) and (H3) of Theorem 2 are satisﬁed by choosing d1 = 0.6, d2 = 0.2 and d1 = 0.05, d2 = 0.2, respectively. As shown in Figs. 4 and 5, E∗ is locally asymptotically stable with different diffusive coeﬃcient of phytoplankton. Especially, if the diffusion of phytoplankton is small, phytoplankton and zooplankton can coexist. In Fig. 6, we gradually increase the value of diffusive coefﬁcient d2 : the condition (ii) m2 ≤

a11 l 2 d1

= 13.9567 ≤ (m + 1 )2 and

d2 > dm = 13.4907 of Theorem 3 is satisﬁed with d1 = 0.05, d2 = 15, m = 3, n = 1, hence Turing instability occurs. This shows that the formation of the pattern is due to the accelerated diffusive movement of zooplankton. Biologically, the formation of the pattern reveals that zooplankton move quickly for the successful predation.

100

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

Fig. 3. Behavior of plankton with intensity of white noise 1.2 (left panel) and 1.6 (right panel) (model (1)). The numerical results are obtained with u = 0.2, v = 0.15 and u = 0.2, v = 0.5, respectively.

Fig. 4. Behavior of plankton with stable conditions (model (2)). The numerical results are obtained with d1 = 0.6, d2 = 0.2, u = 0.2, v = 0.15 and initial condition (0.5, 0.6).

Fig. 5. Behavior of plankton with stable conditions (model (2)). The numerical results are obtained with d1 = 0.05, d2 = 0.2, u = 0.2, v = 0.15 and initial condition (0.5, 0.6).

Fig. 6. Behavior of plankton with Turing instable conditions (model (2)). The numerical results are obtained with d1 = 0.05, d2 = 15, u = 0.2, v = 0.15 and initial condition (0.5, 0.6).

In Fig. 7, we ﬁx n = 2, u = 0.42, v = 0.65, d1 = 0.15 and d2 = 0.08. By choosing a(u, v) as the bifurcation parameter and solving numerically the equations, we obtain that a(u, v ) = 1.4051 < an = 2

1.8983 and a11 = 0.2753 >

d1 (d1 +d2 ) n2 d1 −d2

l

= 0.2190, thus Hopf bifur-

cation occurs. The result shows that increases of defense of phytoplankton and offense of zooplankton affect the stability of co-

existing equilibrium for small diffusive coeﬃcient of zooplankton. In addition, we obtain Re(δ1 (a )) = −124.5601 < 0. Hence, the bifurcated periodic solutions are stable and the period of bifurcated periodic solutions increases with β2 = 53.6071 > 0, this is shown in Fig. 7. In Fig. 8, we change u = 0.6, v = 1.5, while the other parameters remain the same as in Fig. 7. It is observed that the coexisting

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

101

Fig. 7. Behavior of plankton with unstable conditions (model (2)). The numerical results are obtained with d1 = 0.15, d2 = 0.08, u = 0.42, v = 0.65 and initial condition (0.5, 0.6).

Fig. 8. Behavior of plankton for different values of defense and offense. The numerical results are obtained with d1 = 0.15, d2 = 0.08, u = 0.6, v = 1.5 and initial condition (0.5, 0.6).

equilibrium of (2) is locally asymptotically stable. As effects of defense and offense increase, we observe transition of the bifurcated periodic solutions into locally asymptotically stable. Thus, in the plankton ecosystem, the effects of defense and offense have great impact in the development of plankton population.

CRediT authorship contribution statement

6. Discussions

Acknowledgement

In this paper, we have considered the effects of defense of phytoplankton and offense of zooplankton on the dynamics of a spatially extended plankton model. First, the stable scenario for the single compartment model has been analyzed to ﬁnd out the conditions for stable coexistence of both the phytoplankton and zooplankton populations. The defense and offense are two key parameters for the plankton model and hence the unstable behaviours of the coexisting equilibrium have been prepared by varying these two parameters. Then we have described the dynamics of the spatially extended model, for which stability, Turing instability and Hopf bifurcation have been analysed. It has been found that the diffusion could cause Turing instability, and the defense and offense can induce bifurcated periodic solution. Some numerical simulations have been also carried out. Although such a biological adaptive changing (defense and offense) can be included in the model, we disregarded a number of factors which affect the plankton community structure. In particular, the model does not include some factors (time delay, age structure, Allee effect), which are interesting elements that can be taken into account in the future.

This work is partially supported by the National Natural Science Foundation of China (No. U1806203), National Natural Science Foundation of China (No. 61533011) and Science Foundation of North University of China (No. XJJ201804).

Declaration of Competing Interest The authors declare that they have no known competing ﬁnancial interests or personal relationships that could have appeared to inﬂuence the work reported in this paper.

Qiuyue Zhao: Methodology, Writing - original draft, Writing - review & editing. Shutang Liu: Conceptualization, Supervision. Xinglong Niu: Writing - review & editing.

References [1] Turner JT, Tester PA. Toxic marine phytoplankton, zooplankton grazers, and pelagic food webs. Limnol Oceanogr 1997;42(5):1203–13. [2] Edwards AM. Adding detritus to a nutrient-phytoplankton-zooplankton model: a dynamical-systems approach. J Plankton Res 2001;23(4):389–413. [3] Yang RZ, Liu M, Zhang CR. A diffusive toxin producing phytoplankton model with maturation delay and three-dimensional patch. Comput Math Appl 2017;73(5):824–37. [4] Zhao QY, Liu ST, Tian DD. Dynamic behavior analysis of phytoplankton– zooplankton system with cell size and time delay. Chaos Solitons Fractals 2018;113:160–8. [5] Chaturvedi D, Misra OP. Modeling impact of varying pH due to carbondioxide on the dynamics of prey-predator species system. Nonlinear Anal 2019;46:374–402. [6] Uszko W, Diehl S, Englund G, Amarasekare P. Effects of warming on predator-prey interactions–a resource-based approach and a theoretical synthesis. Ecol Lett 2017;20(4):513–23. [7] Delaney MP. Effects of temperature and turbulence on the predator-prey interactions between a heterotrophic ﬂagellate and a marine bacterium. Microb Ecol 2003;45(3):218–25. [8] Giricheva E. Spatiotemporal dynamics of an NPZ model with prey-taxis and intratrophic predation. Nonlinear Dyn 2019;95(2):875–92. [9] Chen FD, Xie XD, Wang HN. Global stability in a competition model of plankton allelopathy with inﬁnite delay. J Syst Sci Complexity 2015;28(5):1070–9. [10] Kloosterman M, Campbell SA, Poulin FJ. An NPZ model with state-dependent delay due to size-structure in juvenile zooplankton. SIAM J Appl Math 2016;76(2):551–77.

102

Q. Zhao, S. Liu and X. Niu / Chaos, Solitons and Fractals 129 (2019) 94–102

[11] Ma ZP, Liu J, Yue JL. Spatiotemporal dynamics induced by delay and diffusion in a predator-prey model with mutual interference among the predator. Comput Math Appl 2018;75(10):3488–507. [12] Han RJ, Dai BX. Spatiotemporal pattern formation and selection induced by nonlinear cross-diffusion in a toxic-phytoplankton-zooplankton model with allee effect. Nonlinear Anal 2019;45:822–53. [13] Li J, Song YZ, Wan H, Zhu HP. Dynamical analysis of a toxin-producing phytoplankton-zooplankton model with refuge. Math Biosci Eng 2017;14(2):529–57. [14] Meng XY, Wu YQ. Bifurcation and control in a singular phytoplankton-zooplankton-ﬁsh model with nonlinear ﬁsh harvesting and taxation. Int J Bifurcation Chaos 2018;28(03):1850042. [15] Shi RX, Yu J. Hopf bifurcation analysis of two zooplankton-phytoplankton model with two delays. Chaos Solitons Fractals 2017;100:62–73. [16] Huang DW, Wang HL, Feng JF, Zhu ZW. Hopf bifurcation of the stochastic model on HAB nonlinear stochastic dynamics. Chaos Soliton Fract 2006;27(4):1072–9. [17] Zhao Y, Yuan SL, Zhang TH. The stationary distribution and ergodicity of a stochastic phytoplankton allelopathy model under regime switching. Commun Nonlinear Sci Numer Simul 2016;37:131–42. [18] Giuffrida A, Valenti D, Ziino G, Spagnolo B, Panebianco A. A stochastic interspeciﬁc competition model to predict the behaviour of Listeria monocytogenes in the fermentation process of a traditional sicilian salami. Eur Food Res Technol 2009;228(5):767–75. [19] Valenti D, Denaro G, La Cognata A, Spagnolo B, Bonanno A, Basilone G, Mazzola S, Zgozi S, Aronica S. Picophytoplankton dynamics in noisy marine environment. Acta Phys Pol B 2012;43(5):1227–40. [20] Mougi A, Iwasa Y. Unique coevolutionary dynamics in a predator-prey system. J Theor Biol 2011;277(1):83–9.

[21] Kiørboe T, Saiz E, Tiselius P, Andersen K. Adaptive feeding behavior and functional responses in zooplankton. Limnol Oceanogr 2018;63(1):308– 321. [22] Pancˇ ic´ M, Kiørboe T. Phytoplankton defence mechanisms: traits and trade-offs. Biol Rev 2018;93(2):1269–303. [23] Jeschke JM. Density-dependent effects of prey defenses and predator offenses. J Theor Biol 20 06;242(4):90 0–7. [24] Velzen E, Gaedke U. Reversed predator-prey cycles are driven by the amplitude of prey oscillations. Ecol Evol 2018;8(12). doi:10.1002/ece3.4184. [25] Velzen E, Gaedke U. Disentangling eco-evolutionary dynamics of predator-prey coevolution: the case of antiphase cycles. Sci Rep 2017;7(1):17125. [26] Mougi A. Predator-prey coevolution driven by size selective predation can cause anti-synchronized and cryptic population dynamics. Theor Popul Biol 2012;81(2):113–18. [27] Wu RC, Zhou Y, Shao Y, Chen LP. Bifurcation and turing patterns of reaction-diffusion activator-inhibitor model. Phys A 2017;482:597–610. [28] Wu RC, Chen MX, Liu B, Chen LP. Hopf bifurcation and turing instability in a predator-prey model with michaelis-menten functional response. Nonlinear Dyn 2018;91(3):2033–47. [29] Yi FQ, Wei JJ, Shi JP. Bifurcation and spatiotemporal patterns in a homogeneous diffusive predator-prey system. J Differ Equs 2009;246(5):1944–77. [30] Hassard BD, Kazarinoff ND, Wan YH. Theory and applications of Hopf bifurcation. Cambridge: Cambridge University Press; 1981. [31] Tan W, Yu WW, Hayat T, Alsaadi F. Turing instability and bifurcation in a diffusion predator-prey model with Beddington-DeAngelis functional response. Int J Bifurcation Chaos 2018;28(09):1830029. [32] Mao XR. Stochastic stabilization and destabilization. Syst Control Lett 1994;23(4):279–90.