Effect of water temperature on the dynamic behavior of phytoplankton–zooplankton model

Applied Mathematics and Computation 378 (2020) 125211

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Effect of water temperature on the dynamic behavior of phytoplankton–zooplankton model Qiuyue Zhao a, Shutang Liu a,∗, Xinglong Niu b a b

School of Control Science and Engineering, Shandong University, Jinan 250061, China 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 17 August 2019 Revised 26 January 2020 Accepted 8 March 2020 Available online 31 March 2020 Keywords: Phytoplankton–zooplankton model Water temperature Holling parameter Hopf bifurcation Sensitivity index

a b s t r a c t In this paper, we investigate the effect of water temperature on phytoplankton– zooplankton dynamics with Holling type III functional response. Remarkably, growth rate, capture rate, handling time, death rate and Holling parameter are considered to be temperature dependent, which makes the model different from those in the existing related literatures. First, the relationships between model parameters and water temperature are presented. Then we give the positivity and boundedness of the solutions. And parameters dependent equilibrium analysis and stability analysis are obtained, respectively. It shows that when the water temperature is lower or higher than the optimal temperature, the equilibrium density of phytoplankton increases with decreasing capture rate and increasing handling time. Meanwhile, we show that the decrease of capture rate and increase of handing time cause the density of zooplankton also increases without the optimal temperature, while the population of zooplankton tends to decrease with decrease of growth rate. In addition, sensitivity of coexisting equilibrium to the parameters of the phytoplankton– zooplankton model is analyzed. Finally, one parameter bifurcation analysis is done. © 2020 Elsevier Inc. All rights reserved.

1. Introduction Temperature is the most important factor of marine environment. Compared with terrestrial organisms, most marine organisms tend to live in relatively stable temperature environment and they are more sensitive to temperature changes. Temperature changes have important effects on the growth, reproduction, metabolism, development, distribution and food availability of marine organisms. Therefore, the question arises as to how changes in temperature affect the populations of biological communities. For this reason, the dynamics of the water temperature changes in marine ecosystems has been a focus of both experimental and theoretical researches for a few decades [1–5]. In particular, [6] examined the effects of increasing sea surface temperature and fluctuating food levels on reef-fish reproduction. The effects of increasing temperatures on the aerobic metabolism and swimming ability of coral reef fishes were investigated in [7]. In South Georgia, as the endemic and range-edge species of marine that may be at their thermal tolerance limits, the influence of temperature change to ecological diversity of marine could be severe [8]. Beisner et al. [9] investigated the effect of temperature (18 o C and 25 o C) on the stability of a common freshwater predator–prey system consisting of Daphnia pulex and phytoplankton in different types of mesocosm communities. Using a consumer-resource model, Amarasekare [10] addressed that the increase ∗

Corresponding author. E-mail addresses: [email protected] (Q. Zhao), [email protected], [email protected] (S. Liu), [email protected] (X. Niu).

https://doi.org/10.1016/j.amc.2020.125211 0 096-30 03/© 2020 Elsevier Inc. All rights reserved.

2

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

or decrease of consumer-resource oscillations by warming depends on the manner in which temperature affects intraspecific competition. Phytoplankton is the most important primary producer in marine ecosystems. To date, many studies have shown that ocean warming has a huge impact on phytoplankton. For example, Staehr and Sand-Jensen [11] investigated the extent by which temperature influences the rate of photosynthesis, respiration and growth of natural phytoplankton and tested seasonal dependence. Bestion et al. [12] developed a theoretical model to understand that the changes in temperature affect competitive interactions among phytoplankton. In [13], temperature has an indirect effect on phytoplankton diversity which mainly depends on other factors, such as grazing, subsidence or nutritional restriction depending on temperature. Toseland et al. [14] studied the impact of temperature on resource allocation and metabolism of phytoplankton. The researchers suggested that eukaryotic phytoplankton requires lower density ribosomes to produce the required amount of cellular protein at higher temperatures. Hare et al. [15] analyzed the effects of temperature on phytoplankton community structure and suggested that increasing temperature could be a primary driving force in future algal dominance shifts away from diatoms and towards smaller nanophytoplankton groups in the Bering Sea. Additionally, temperature also can act changes on zooplankton. Two major changes have been validated by experiments: (i) decreasing Cladoceran body size with rising temperature [16] and (ii) changes of critical density with water temperature [17]. However, due to the complexity of the mechanism of zooplankton, its various metabolic activities are indirectly affected by water temperature. The effects of temperature on zooplankton would contribute to modifying its predation ability and inducing shifts in the community structure and density of phytoplankton, which may change the whole marine food web. The temperature is sufficiently important so that for this paper, we explores the influence of temperature on plankton dynamics using the temperature-dependent parameters. The relationships between temperature and parameters of plankton model are modeled as hump-shaped and U-shaped in the range of physiologically relevant temperature. We consider a coupled plankton-temperature model, and it exhibits a rich dynamics for different temperatures and parameters. The rest of this paper is organized as follows. Section 2 describes the phytoplankton-zooplankton model with temperature-dependent parameters. Positivity and boundedness of solutions, parameters dependent equilibrium analysis, stability analysis and hopf bifurcation are established in Section 3. Some numerical simulations are presented to illustrate our theoretical results and Hopf bifurcation with respect to some important parameters is done in Section 4. In Section 5 we present sensitivity analysis of the coexisting equilibrium. Finally, we give some discussions and biological significance of our analytical findings in Section 6. 2. Formulation of temperature-dependence model Functional response, a significant component of predator-prey relationship, reflects the strength of predator–prey interactions. These functional response functions, which divide into various types based on densities of prey and predator, affect the stability of predator-prey dynamics and maintain the dynamic balance of community structure. Following C. S. Holling, functional responses are generally classified into four types, which are called Holling’s type I, II, III and IV. The Holling type II functional response is often modeled by a rectangular hyperbola, which follows from the assumption that prey is scarce. When predator actively seeks out large concentrations of prey, the Holling type III is more appropriate. In addition, the type III functional response can potentially stabilize fluctuations of prey and predator populations, and reduce prey extinction risk [18,19]. Therefore, we concentrate on the following phytoplankton-zooplankton model:

⎧   P ⎪ P− ⎨P˙ = r 1 −

aP n Z , 1 + ahP n n ⎪ ⎩Z˙ = −mZ + e aP Z , 1 + ahP n K

(1)

with initial conditions P(0) > 0, Z(0) > 0. n > 0 is Holling parameter. n = 1 is a type II response and n > 1 implies a type III response. P denotes the density of phytoplankton population and Z denotes the density of zooplankton population. r and K are the intrinsic growth rate and carrying capacity of phytoplankton. The parameters a is the capture rate of zooplankton on phytoplankton and e is the conversion rate of phytoplankton to zooplankton. h and m are the handling time and the natural mortality rate of zooplankton, respectively. 2.1. Temperature dependence of model parameters To explore the effect of temperature on the dynamic behavior of plankton with Holling type III function response, we consider temperature dependent parameters. (1) The natural mortality rate of zooplankton is increased exponentially with temperature. The mathematical model is given in [20] as follows:

m(T ) = m0 exp(−E/(kT )), where m0 is a constant, E is activation energy, k is Boltzmann’s constant (8.62 × 10−5 eVK−1 ), T is absolute temperature of water. The dynamics of mortality rate function with respect to T is shown in Fig. 1 in which the variation of mortality rate is similar to the experimental study [21] and the Fig. 1b of [20].

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

3

1 0.8 0.6 0.4 0.2 0 270

280

290

300

310

320

Fig. 1. Variation of natural mortality rate of zooplankton with temperature.

2.5 2 1.5 1 0.5 0 0

10

20

30

40

50

Fig. 2. Variation of intrinsic growth rate of phytoplankton with temperature.

(2) The intrinsic growth rate of phytoplankton is a hump-shaped relationship with temperature. The mathematical model is given in [20] as follows:

r (T ) = r0 exp(−(T − Topt )2 /(2S2 )), where Topt is optimal Celsius temperature at which the intrinsic growth rate r reaches its maximum, S is the function breadth, and r0 is a constant. The dynamics of intrinsic growth rate function with respect to T is shown in Fig. 2 in which the variation of intrinsic growth rate is similar to the experimental study [22] and Fig. 1a of [20]. (3) The capture rate of zooplankton is a hump-shaped relationship with temperature. The mathematical model is given in [20] as follows:

a(T ) = a0 exp(−(T − Topt )2 /(2S2 )), where Topt is optimal Celsius temperature at which the capture rate a reaches its maximum, S is the function breadth, and a0 is a constant. The dynamics of capture rate function with respect to T is shown in Fig. 3 in which the variation of capture rate is similar to the experimental study [23] and the Fig. 1c of [20]. (4) The handling time of zooplankton is a U-shaped relationship with temperature. The mathematical model is given in [20] as follows:

h(T ) = h0 exp((T − Topt )2 /(2S2 )), where Topt is Celsius optimal temperature at which the handling time h reaches its minimum, S is the function breadth, and h0 is a constant. The dynamics of handling time function with respect to T is shown in Fig. 4 in which the variation of handling time is similar to the experimental study [23] and the Fig. 1a of [20].

4

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

10 8 6 4 2 0 0

10

20

30

40

Fig. 3. Variation of capture rate of zooplankton with temperature.

3 2.5 2 1.5 1 0.5 0 5

10

15

20

25

30

35

Fig. 4. Variation of handling time of zooplankton with temperature.

(5) The Holling parameter of functional response is a U-shaped relationship with temperature. The mathematical model is given in [20] as follows:

n(T ) = n0 exp((T − Topt )2 /(2S2 )), where Topt is Celsius optimal temperature at which the Holling parameter n reaches its minimum, S is the function breadth, and n0 is a constant. The dynamics of mortality rate function with respect to T is shown in Fig. 5 in which the variation of intrinsic growth rate is similar to Fig. 1e of [20]. According to Fig. 5, it has been assumed that the functional response p(x ) = uously differentiable function defined on [0, ∞) and satisfies

p( 0 ) = 0 ,



p ( x ) > 0,

lim p(x ) =

x→∞

axn 1+ahxn

is Holling type III, which is a contin-

1 < ∞. h

It shows that p(x) is monotonic. The biologically explains that the effect of predation is low when the phytoplankton population is low, but as the number of population increases, the predation will be more intensive. However, when the nutrient concentration reaches a high level, it may have an inhibitory effect on the growth rate of the population. Due to the nonmonotonic responses, the Holling type III is not suitable for situations involving inhibition and defense. Additionally, we assume that metabolic rates and parameters change exponentially with temperature, which has enhanced our understanding of how temperature variation influences the dynamics of plankton interactions. These assumptions about the temperature dependence of model parameters also make the phytoplankton–zooplankton model closer to the actual situation and add complex of the proposed model.

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

5

1.8 1.7 1.6 1.5 1.4 1.3 1.2 0

10

20

30

40

Fig. 5. Variation of Holling parameter of functional response with temperature.

In [20], the authors described a parameterized pelagic producer-predator system for green alga Monoraphidium minutum and crustacean Daphnia hyalina. They illustrated how to derive the temperature dependence of prey carrying capacity and demonstrated that different shapes of this temperature dependence are possible under different resource supply constraints. Since model (1) is a plankton model with temperature-dependent model parameters, it can also describe the interaction between the green alga Monoraphidium minutum and the crustacean Daphnia hyalina. However, our purpose is to make a detailed analysis of this specific system from a biomathematical perspective. 3. Qualitative analysis of model (1) In this section, we first give positivity and boundedness of solutions. Then we show parameters dependent equilibrium analysis. In addition, the conditions for local and global stability of model (1), and the sufficient conditions for occurrence of Hopf bifurcation are also derived. 3.1. Positivity and boundedness Theorem 1. All solutions of model (1) are positive and uniformly bounded. Proof. Let (P(t), Z(t))  be any solution with positive initial conditions. Due to the right-hand side of model (1) is continuous and smooth on R2+ = (P, Z ) : P, Z > 0 , from the initial values P(0) > 0, Z(0) > 0, we have



P (t ) = P (0 )e Z (t ) = Z (0 )e

t 0

t 0





(

)

n−1

Z r 1− KP − aP 1+ahP n

n −m+e 1+aPahPn



ds

ds

> 0,

> 0.

Thus, the solutions remain in the R2+ . From first equation of model (1), we obtain





P dP ≤r 1− P, dt K Using comparison theorem, we find that lim supt→∞ P (t ) ≤ K. Let ζ (t ) = eP (t ) + Z (t ). Along the equations of model (1), we get





P dζ (t ) = eP˙ (t ) + Z˙ (t ) = er 1 − P − mZ dt K ≤ erP − mZ ≤ 2erP − θ (eP + Z ), where θ = min{r, m}. Then,

dζ (t ) + θ ζ (t ) ≤ 2erK. dt By calculation, we obtain lim supt→∞ Z (t ) ≤ 2erK θ . Therefore, the solutions of model (1) are ultimately bounded.



6

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

Theorem 2. For any given initial value (P0 , Z0 ) ∈ R2+ , there is a unique positive solution (N(t), P(t)) of model (1). Proof. Let u(t ) = ln P (t ) and v(t ) = ln Z (t ), we obtain the following equations





u˙ (t ) = r 1 −

eu(t ) K



eu(t ) −

aeu(t )n ev(t ) , 1+aheu(t )n

(2)

u (t )n v (t )

v˙ (t ) = −mev(t ) + e a1+e aheeu(t )n ,

on t ≥ 0 with initial values u(0 ) = ln P0 , v(0 ) = ln Z0 . Obviously, the coefficients of (2) satisfy local Lipschitz continuous, hence there is a unique solution of (2). Then, P (t ) = eu(t ) , Z (t ) = ev(t ) are the unique positive solution of (1) with initial values (P0 , Z0 ) ∈ R2+ .  3.2. Parameters dependent equilibrium analysis In this subsection, we investigate the influence of parameters by assessing the derivative along the coexisting equilibrium. All biological feasible equilibria of model (1) are as follows. (i) The trivial equilibrium E 0 = (0, 0 ) is always feasible. (ii) The zooplankton free equilibrium E 1 = (K, 0 ) is always feasible. (iii) The coexisting equilibrium E ∗ = (P ∗ , Z ∗ ) is feasible if (H1) e > m(T )

 P∗ =

m (T ) ea(T ) − a(T )m(T )h(T )

 n(1T )



Z ∗ = r (T ) 1 −

,



1 K n (T ) a (T )

+ h (T )

 holds, where



P ∗ 1 + a(T )h(T )P ∗n(T ) , K a(T )P ∗(n(T )−1)

which means that for existence of phytoplankton, the conversion rate of phytoplankton is greater than the product of the handling time and the natural mortality rate of zooplankton. Then, we give the effects of parameters on the coexisting equilibrium. Differentiating P∗ and Z∗ with respect to a(T) and e, we get

dP ∗ −m(T ) = da(T ) n(T )a(T )2 (e − m(T )h(T )) dP ∗ −m(T ) = de n(T )(e − m(T )h(T ))2







m (T ) ea(T ) − a(T )m(T )h(T )

m (T ) ea(T ) − a(T )m(T )h(T )

 n(1T ) −1

 n(1T ) −1

,

,



P∗ dZ ∗ −r (T ) = 1− . da(T ) K a(T )2 P ∗(n(T )−1) From above, we can observe that P∗ is a strictly decreasing function of a(T) and e, that is, increasing the capture rate of zooplankton and the conversion rate of phytoplankton lead to the decrease of coexisting equilibrium density of phytoplankton (see Fig. 6). And Z∗ is a strictly decreasing function of a(T), that is, increasing capture rate of zooplankton leads to the decrease of coexisting equilibrium density of zooplankton (see Fig. 6). Again, applying the differential equality, we obtain

dP ∗ = dm(T ) dP ∗ = dh ( T )

 

m (T ) ea(T ) − a(T )m(T )h(T ) m (T ) ea(T ) − a(T )m(T )h(T )





P∗ ∗ dZ ∗ = r (T ) 1 − P , dh ( T ) K

 n(1T ) −1  n(1T ) −1

e , n(T )(e − m(T )h(T ))2 m ( T )2 , n(T )(e − m(T )h(T ))2





P ∗ 1 + a(T )h(T )P ∗n(T ) dZ ∗ = 1− . dr (T ) K a(T )P ∗(n(T )−1)

This shows that P∗ increases as m(T) and h(T) increase, that is, increasing mortality rate, handling time of zooplankton and growth rate of phytoplankton lead to the increase of coexisting equilibrium density of phytoplankton (see Fig. 7). And Z∗ increases as h(T) and r(T) increase, that is, increasing handling time of zooplankton and growth rate of phytoplankton lead to the increase of coexisting equilibrium density of zooplankton (see Fig. 7). From the simple calculation results, we find significant effects of water temperature and parameters on coexisting equilibrium density of plankton. For both phytoplankton and zooplankton, increasing and decreasing temperature generally change the interaction between predator and prey, which may not balance their metabolic demands. 3.3. Stability analysis In this subsection, we first derive the local stability conditions by using Routh–Hurwitz criterion. Then, we show global asymptotic stability with Lyapunov method. For convenience, we drop the T in the representation of coefficients in the following subsection.

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

0.5

7

0.4

0.4

0.3

0.3 0.2 0.2 0.1

0.1 0

0 0

20

40

60

80

100

20

40

60

80

100

1.4 1.2 1 0.8 0.6 0.4 0.2 0

20

40

60

80

100

Fig. 6. Decrease of coexisting equilibrium of model (1) with respect to increase of parameters.

Theorem 3. If (H1) and K totically stable.



m ea−amh

−1 n

e e−mh

−n +n<

2e e−mh

hold, then the coexisting equilibrium of model (1) is locally asymp-

Proof. The Jacobian matrix at E∗ is given by



J (a ) =

a11 (a ) a21 (a )



a12 (a ) , 0

where

a11 (a ) = −

P∗ r aP ∗n aZ ∗ P ∗(n−1) (1 + ahP ∗n − n ) aenZ ∗ P ∗(n−1) + , a12 (a ) = − , a21 (a ) = . ∗n ∗n 2 K 1 + ahP (1 + ahP ) (1 + ahP∗n )2

The characteristic equation of the Jacobian matrix is

λ2 − a11 (a )λ − a12 (a )a21 (a ) = 0.

(3)

Let λ1 and λ2 be two roots of (3). Then, we get

λ1 + λ2 = a11 (a ),

λ1 λ2 = −a12 (a )a21 (a ).

This together with Routh-Hurwitz criterion derives λ1 + λ2 < 0 and λ1 λ2 > 0 if a11 (a) < 0, which implies that model (1) is local stability around coexisting equilibrium E∗ if

aZ ∗ P ∗(n−2) (1 + ahP ∗n − n ) r < . K (1 + ahP∗n )2

(4)

Furthermore, we have

1 + ahP ∗n =

eaP ∗n . m

(5)

8

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

0.4

0.11 0.1

0.3

0.09 0.08

0.2

0.07 0.06

0.1

0.05 0.04

0 20

40

60

80

100

50

100

150

200

0.65 0.6

0.8

0.55 0.6 0.5 0.4

0.45 0.4

0.2

0.35 0

50

100

150

200

20

40

60

80

100

Fig. 7. Increase of coexisting equilibrium of model (1) with respect to increase of parameters.

Using P∗ , Z∗ and (5), inequality (4) has the following form

 K

m ea − amh

 −1n 

1+

hm −n e − mh





+n<2 1+



hm . e − mh

Thus we can say that model (1) is local stability around coexisting equilibrium E∗ .



According to the construction method of Lyapunov functional of [25,26], we can get the global asymptotic stability of coexisting equilibrium. Theorem 4. If (H1) holds, then model (1) is globally asymptotically stable around coexisting equilibrium E∗ when n ≥ 1 + ahK n . Proof. Let f (P ) =

V (t ) =



P P∗

aP n , 1+ahP n

construct a suitable Lyapunov function as follows:

e f (u ) − m du + f (u )



Z Z∗

v − Z∗ dv . v

Taking the time derivative of V(t) along the solutions of model (1) obtains

dV (t ) e f (P ) − m Z − Z∗ = P˙ + Z˙ dt f (P ) Z =





P aP n Z e f (P ) − m r 1− P− f (P ) K 1 + ahP n

= e( f (P ) − f (P ∗ ) )(h(P ) − h(P ∗ )),



 + (Z − Z ∗ ) −m + e

aP n 1 + ahP n



Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211 

9 



)(1+ahP ) +(n−1 )(P−K )−P ) where h(P ) = r (K−PKaP and h (P ) = r (ahP (K−2P )KaP . Obviously, when n ≥ 1 + ahK n , f (P ) > 0, h (P ) < 0 for all n n−1 ∗ ˙ ˙ P > 0. Then V (t ) ≤ 0 and V (t ) = 0 if and only if P = P . Therefore, E∗ is globally asymptotically stable.  n

n

3.4. Hopf bifurcation In this subsection, we discuss that there exists periodic solutions of model (1) about coexisting equilibrium E∗ , when the parameter a exceeds threshold value a∗ . e n

Theorem 5. If (H1) and mh < e < mh + a = a∗ .

hold, then model (1) undergoes a Hopf bifurcation at coexisting equilibrium when

Proof. Let

aZ ∗ P ∗(n−2) (1 + ahP ∗n − n ) r = , K (1 + ahP∗n )2 By P∗ , Z∗ and (5), we can obtain

a∗ =

m e − mh



nhm + (2 − n )e K (nhm + (1 − n )e )

n .

Then, it is found that when a = a∗ , (3) has a pair of purely imaginary roots. Again





da11  r m =  da a=a∗ Kn e − mh

 1n 1    1 +1 an

a=a∗

> 0.

Thus the transversality condition is satisfied and a simple Hopf bifurcation occurs at a = a∗ .



Next, we give an explicit formula determining the direction of Hopf bifurcation of (1) at E∗ . We translate E∗ to the origin by the translation x¯ = P − P ∗ , y¯ = Z − Z ∗ . Then model (1) is transformed into

  x¯˙ y¯˙

 

= J (a )

x¯ y¯



+



f (x¯, y¯ ) , h(x¯, y¯ )

(6)

where

f (x¯, y¯ ) = a1 x¯2 + a2 x¯y¯ + a3 x¯3 + a4 x¯2 y¯ + O(|(x¯, y¯ )|4 ), h(x¯, y¯ ) = b1 x¯2 + b2 x¯y¯ + b3 x¯3 + b4 x¯2 y¯ + O(|(x¯, y¯ )|4 ), with

r n − 1 − ahp∗n (n + 1 ) anP ∗(n−1) − anZ ∗ P ∗(n−2) , a2 = , K 2(1 + ahP ∗n )3 2(1 + ahP ∗n )2 (n − 1 )(n − 2 − 2ahP∗n (n + 1 )) + ah(1 + n )P∗n (2(1 − n ) + ahP∗n (n + 2 )) a3 = −anZ ∗ P ∗(n−3) , 6(1 + ahP ∗n )4 n − 1 − ahp∗n (n + 1 ) 2 n − 1 − ahp∗n (n + 1 ) a4 = −anP ∗(n−2) b1 = aenZ ∗ P ∗(n−2) , x¯ y¯ , ∗n 3 3(1 + ahP ) 2(1 + ahP ∗n )3 a1 = −

aenP ∗(n−1) n − 1 − ahp∗n (n + 1 ) , b4 = aenP ∗(n−2) , ∗n 2 2(1 + ahP ) 3(1 + ahP ∗n )3 (n − 1 )(n − 2 − 2ahP∗n (n + 1 )) + ah(1 + n )P∗n (2(1 − n ) + ahP∗n (n + 2 )) b3 = aenZ ∗ P ∗(n−3) . 6(1 + ahP ∗n )4 b2 =

a11 (a ) 2 ,

Set λ1,2 (a ) = α (a ) ± iω (a ) with α (a ) =

ω ( a∗ )

=



Set

−a12

 B=

( a∗ )a

21

a12 (a ) a11 (a )

(a∗ ), 0



 B−1 =

,

 

  x˜˙ y˜˙



α (a ) = ω (a )

1 2



−4a12 (a )a21 (a ) − a11 (a )2. When a = a∗ , λ1,2 (a∗ ) = ±iω (a∗ ),

the one eigenvector corresponding to the eigenvalues λ1,2 (a∗ ) is ξ = (a12 (a∗ ), a11 (a∗ ) − iω (a∗ ))T .

−ω (a )

Through the transformation

ω (a ) =

x˜ y˜

  =B

x˜ y˜

0 − ω (1a )

 .

x¯ , then (6) becomes y¯

 

−ω (a ) α (a )

1 a12 (a ) a11 (a ) ω (a )a12 (a )

 +



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

(7)

10

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

where

1 f (a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ ), a21 (a ) a11 (a ) h1 (x˜, y˜ ) = f (a (a )x˜, a11 (a )x˜ − ω (a )y˜ ) − ω (a )a12 (a ) 12 f1 (x˜, y˜ ) =

1

ω (a )

h(a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ ),

and

f (a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ ) = a1 a12 (a )2 x˜2 + a2 a12 (a )x˜(a11 (a )x˜ − ω (a )y˜ ) + a3 a12 (a )3 x˜3 +a4 a12 (a )2 x˜2 (a11 (a )x˜ − ω (a )y˜ ) + O(|(a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ )|4 ), h(a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ ) = b1 a12 (a )2 x˜2 + b2 a12 (a )x˜(a11 (a )x˜ − ω (a )y˜ ) + b3 a12 (a )3 x˜3 +b4 a12 (a )2 x˜2 (a11 (a )x˜ − ω (a )y˜ ) + O(|(a12 (a )x˜, a11 (a )x˜ − ω (a )y˜ )|4 ). We rewrite (7) into the following polar coordinate form



l˙ = α (a )l + σ (a )l 3 + O(|l |4 ), θ˙ = ω (a ) + c(a )l 2 + O(|l |4 ).

(8)

By Taylor expansion of (8) at a = a∗ , we have



l˙ = α  (a∗ )(a − a∗ )l + σ (a∗ )l 3 + O(|l |4 ), θ˙ = ω (a∗ ) + ω (a∗ )(a − a∗ ) + c(a∗ )l 2 + O(|l |4 ).

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

σ ( a∗ ) =





1 1 ( f1x˜x˜x˜ + h1x˜x˜y˜ ) + f1x˜y˜ f1x˜x˜ − h1x˜y˜ h1x˜x˜ − f1x˜x˜ h1x˜x˜ , 16 16ω (a∗ )

where

f1x˜x˜x˜ =

1 (a3 a12 (a∗ )3 + a4 a12 (a∗ )2 a11 (a∗ )), a21 (a∗ )

f1x˜y˜ = −

1 a2 a12 (a∗ )ω (a∗ ), a21 (a∗ )

h1x˜y˜ = −a11 (a∗ )a2 + b2 a12 (a∗ ), h1x˜x˜ =

h1x˜x˜y˜ = −a11 (a∗ )a4 a12 (a∗ ) + b4 a12 (a∗ )2 ,

1 (a1 a12 (a∗ )2 + a2 a12 (a∗ )a11 (a∗ )), a21 (a∗ )

f1x˜x˜ =

h1x˜x˜y˜ = −a11 (a∗ )a4 a12 (a∗ ) + b4 a12 (a∗ )2 ,

a11 (a∗ ) (a a (a∗ ) + a2 a11 (a∗ )) − ω (a∗ ) 1 12

1

ω ( a∗ )

(b1 a12 (a∗ )2 + b2 a12 (a∗ )a11 (a∗ )).

Thus, we obtain

μ2 = −

σ ( a∗ ) . α  ( a∗ )

According to Poincaré-Andronov Hopf bifurcation theorem, the following result holds. Theorem 6. If (H1) and mh < e < mh + ne are satisfied, then the direction of Hopf bifurcation of (1) at E∗ is supercritical when μ2 > 0 while the direction of Hopf bifurcation is subcritical as μ2 < 0. In addition, the bifurcated periodic solutions are unstable when σ (a∗ ) > 0 while the bifurcated periodic solutions are stable when σ (a∗ ) < 0. 4. Numerical simulations In this section, we present numerical simulations to support the analytical results by taking the parameter values in Table 1. Taking T = 25 o C, we get the stability diagram of E∗ in Fig. 8(a). If we increase the value of T, then the model changes its behavior (see Fig. 8(b), 8(c) and 8(d)). By regarding the capture rate a of zooplankton as a bifurcation parameter and using the values of the parameters mentioned in Table 1 along with T = 30 o C, it is observed that E∗ is locally asymptotically stable when a = 1.5 < a∗ (see Fig. 9) and unstable when a = 2 > a∗ (see Fig. 10), and the model undergoes a Hopf bifurcation around E∗ at a = a∗ = 1.8817. By calculation, we get μ2 > 0 and σ (a∗ ) < 0, which shows that model (1) undertakes supercritical Hopf bifurcation at capture rate a = a∗ and the bifurcated periodic solutions are stable (see Figs. 10 and 11). For a small capture rate, the plankton population stays at steady state, which will benefit the sustainable development of our plankton model, whereas at relatively high capture rate, the ecological balance of our ecosystem will be disrupted. Furthermore, the appearance of bifurcated periodic solutions of (1) causes the fluctuation of the plankton population (see Fig. 10). The stable bifurcated periodic solutions imply that the phytoplankton and zooplankton are still coexist in our plankton ecosystem with the small amplitude oscillations.

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

1

8

0.8

6

0.6

4

0.4

2

0.2

11

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0

2

4

(a)

6

8

10

6

8

10

8

10

(b)

8

8

6

6

4

4

2

2

0

0 0

2

4

6

8

10

0

2

4

(c)

(d)

Fig. 8. Phase plots of model (1) with different values of T: (a) T = 25 C; (b) T = 29.8 C; (c) T = 30 o C; (d) T = 30.5 o C. o

10

o

8

8

6

6 4 4 2

2 0

0 0

100

200

300

400

500

0

2

4

6

Fig. 9. Behavior of plankton (left panel), and the phase plot of model (1) (right panel) with a = 1.5 < a∗ and T = 30 o C.

12

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

10

8

8

6

6 4 4 2

2 0

0 0

100

200

300

400

500

0

2

4

6

Fig. 10. Behavior of plankton (left panel), and the phase plot of model (1) (right panel) with a = 2 > a∗ and T = 30 o C.

Fig. 11. Bifurcation analysis as a function of capture rate a of zooplankton with respect to population densities.

8

10

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

13

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

Description

Values

Reference

r K h a e m n T

Intrinsic growth rate of phytoplankton Carrying capacity of environment Handling time Capture rate of zooplankton Conversion efficiency Mortality rate of zooplankton Holling parameter of functional response Water temperature

r0 = 2.2, S = 12, Topt = 25 10 h0 = 0.17, S = 6.4, Topt = 20 a0 = 8.9, S = 9.4, Topt = 22 0.5 m0 = 4.4 × 108 , E = −0.55, k = 8.62 × 10−5 n0 = 1.2, S = 22, Topt = 16 -

[20] [20] [20] [20] [20] -

10

0.35 0.3

8 0.25 0.2

6

0.15

4

0.1 2 0.05 0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.2

0.65

6

0.6

5

0.55

4

0.5

3

0.45

2

0.4

1

0.35 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.3

0.4

0.5

0.6

0.7

0.8

Fig. 12. Bifurcation analysis as a function of natural mortality rate m of zooplankton with respect to population densities.

Bifurcation analysis with respect to capture rate a of zooplankton is presented in Fig. 11 for other parameters remaining the same as Table 1. From the left panels of Fig. 11, if T = 20 o C, the model remains stable. From the right panels of Fig. 11, if T = 30 o C, the model remains stable with a small range of a and the model has a Hopf bifurcation with larger a. From this figure, if the capture rate of zooplankton increases, then the ecosystem of phytoplankton and zooplankton becomes unstable with high temperature. Bifurcation analysis with respect to natural mortality rate m of zooplankton is presented in Fig. 12 for other parameters remaining the same as Table 1. From the left panels of Fig. 12, if T = 20 o C, we see that the model remains stable. From the right panels of Fig. 12, if T = 28 o C, the model remains stable with a small range of m and the model has a Hopf bifurcation with larger m. From this figure, if the natural mortality rate of zooplankton increases, then the ecosystem of phytoplankton and zooplankton becomes unstable with high temperature.

14

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

Fig. 13. Bifurcation analysis as a function of handling time h of zooplankton with respect to population densities.

Bifurcation analysis with respect to handling time h of zooplankton is presented in Fig. 13 for other parameters remaining the same as Table 1. From the left panels of Fig. 13, if T = 10oC, the model remains stable. From the right panels of Fig. 13, if T = 20 o C, the model remains stable with a small range of h and the model has a Hopf bifurcation with larger h. From this figure, if the handling time of zooplankton increases, then the ecosystem of phytoplankton and zooplankton becomes unstable with a appropriate temperature. Bifurcation analysis with respect to Holling parameter n of functional response is presented in Fig. 14 for other parameters remaining the same as Table 1. From the left panels of Fig. 14, if T = 20 o C, the model remains stable. From the right panels of Fig. 14, if T = 30 o C, the model has a Hopf bifurcation with a small range of n and the model remains stable with larger n. From this figure, if the Holling parameter of type III decreases, then the ecosystem of phytoplankton and zooplankton becomes unstable with high temperature. Bifurcation analysis with respect to intrinsic growth rate r of phytoplankton is presented in Fig. 15 for other parameters remaining the same as Table 1. From the left panels of Fig. 15, if T = 12 o C, the model remains stable. From the right panels of Fig. 15, if T = 29 o C, the model remains stable with a small range of r and the model has a Hopf bifurcation with larger r. From this figure, if the intrinsic growth rate of phytoplankton increases, then the ecosystem of phytoplankton and zooplankton becomes unstable with high temperature. Bifurcation analysis with respect to water temperature T is presented in Fig. 16 for other parameters remaining the same as Table 1. From Fig. 16, the model has two Hopf bifurcations. From this figure, if the water temperature is too low or too high, then the ecosystem of phytoplankton and zooplankton becomes unstable.

5. Sensitivity analysis To understand the dynamical behavior of phytoplankton and zooplankton, it is necessary to know the relative importance of the different factors responsible for the variety of biomass. We calculate the sensitivity indices of coexisting equilibrium

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

15

Fig. 14. Bifurcation analysis as a function of Holling parameter n of functional response with respect to population densities. Table 2 Sensitivity indices of E∗ to parameters for model (1) with temperature T = 20 o C. Parameters xi

γxPi

r K h a e m n

0 0 0.3976 −0.8197 −7.5294 7.5294 −3.9996

∗

γxZi

∗

1 0.0069 0.0528 −0.9472 0.1405 −0.1405 2.669

to the parameters of the model. Sensitivity index of a variable with respect to parameter is calculated by normalized forward sensitivity method [24]. Definition 1 [24]. The normalized forward sensitivity index of a variable, X, that depends differentiably on a parameter, b, is defined as:

γbX :=

∂X b × . ∂b X

The sensitivity index of the coexisting equilibrium E∗ is shown in Table 2 for the parametric values given in Table 1 with temperature T = 20 o C. These indices tell us how crucial each parameter is to plankton biomass. From the Table 2, we can intuitively observe the influence of almost all parameters on the coexisting equilibrium. If the sign of the sensitivity indices of E∗ is positive, then E∗ increases by increase of parameters, conversely, E∗ decreases by increase of parameters.

16

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

0.25

0.7 0.6

0.2 0.5 0.15

0.4

0.1

0.3 0.2

0.05 0.1 0

0 0.5

1

1.5

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0.5

1

1.5

0.5

1

1.5

0 0.5

1

1.5

Fig. 15. Bifurcation analysis as a function of intrinsic growth rate r of phytoplankton with respect to population densities.

Fig. 16. Bifurcation analysis as a function of water temperature T with respect to population densities.

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

17

We find that the equilibrium density of the phytoplankton is most sensitive to the conversion rate e of phytoplankton and the mortality rate m(T) of zooplankton, but it is not sensitive to the growth rate r(T) and the carrying capacity K of phytoplankton. In addition, the equilibrium density of zooplankton is again most sensitive to the Holling parameter n(T) of functional response. Thus, reducing the mortality rate of zooplankton and increasing the conversion rate would have a large effect on the population of phytoplankton at low temperature. 6. Conclusions In this paper, a interaction between phytoplankton and zooplankton in the presence of water temperature has been developed. Here, the effect of water temperature on phytoplankton and zooplankton dynamics has been considered. It is assumed that the natural mortality rate of zooplankton is proportional to the increase of water temperature. The growth rate of phytoplankton and capture rate of zooplankton are hump-shaped to the increase of water temperature, and the handling time of zooplankton and Holling parameter of functional response are U-shaped. The proposed model has been analysed theoretically as well as numerically. The positivity and boundedness of all solutions, and parameters dependent equilibrium analysis also have been discussed. It is found that, the equilibrium density of phytoplankton population increases by increase of mortality rate, increase of handling time, decrease of capture rate and decrease of conversion rate, while that of zooplankton increases with the increase of handling time and growth rate and decrease of capture rate. In fact, from Figs. 1 to 7, it is easy to show that there is an increase in the equilibrium density of both phytoplankton and zooplankton as the water temperature increase under a very restricted set of conditions. This occurs when the water temperature rises from 25 o C to 30 o C Then, by theoretical analyses, we have found that the effects of the parameters of the model play an important role in determining the stability of the coexisting equilibrium. On the one hand, the left panels of Figs. 11–13 have shown that when the water temperature is maintained at 20 o C, the model is stable. In addition, as the capture rate a(T), mortality rate m(T) and handling time h(T) of zooplankton increase, the model remains stable up to critical values of these parameters but after critical values the model becomes unstable with a higher water temperature (see the right panels of Figs. 11–13). On the other hand, it is also observed that system dynamics have oscillation as well as stable state in presence of high temperature with the change of Holling parameter n(T) of functional response and growth rate r(T) of phytoplankton (see the right panels of Figs. 14 and 15), whereas the system has only stable state in low temperature (see the left panels of Figs. 14 and 15). Finally, it point out that the water temperature plays an important role for existence of stability as well as oscillatory dynamics (see Fig. 16). Therefore, dynamics of the proposed model is much more complex due to the presence of water temperature. Our results have suggested that changes in temperature can influence the dynamics of plankton in predictable ways, given the dependence of their rates of growth, ingestion, and metabolism on temperature. The dynamics of our theoretical calculation are not as complicated as the actual situations, but the effects of temperature on density of plankton offers similar insights to [27] which pointed that the increase of sea surface temperature in the Northeast Atlantic Ocean results in the decrease of phytoplankton in warm waters and the increase of phytoplankton in cold water. [28] obtained that when the temperature increased, the population density and quantity of some warm water fishes increased, while that of cold water fishes increase. To understand the dynamics of temperature in a three species ecological models further investigations are necessary in this direction. Acknowledgment This work is partially supported by the National Natural Science Foundation of China (Nos. U1806203, 61533011, 61903343) and Science Foundation of North University of China (No. XJJ201804). References [1] B.C. Rall, O. Vucic-pestic, R.B. Ehnes, M. Emmerson, U. Brose, Temperature, predator–prey interaction strength and population stability, Glob. Chang. Biol. 16 (8) (2010) 2145–2157. [2] J.L. Lessard, D.B. Hayes, Effects of elevated water temperature on fish and macroinvertebrate communities below small dams, River Res. Appl. 19 (7) (2003) 721–732. [3] J.C. Morrill, R.C. Bales, M.H. Conklin, Estimating stream temperature from air temperature: implications for future water quality, J. Environ. Eng. 131 (1) (2005) 139–146. [4] L.S. Peck, K.E. Webb, D.M. Bailey, Extreme sensitivity of biological function to temperature in antarctic marine species, Funct. Ecol. 18 (5) (2004) 625–630. [5] A. Clarke, Temperature and the metabolic theory of ecology, Funct. Ecol. 20 (2) (2006) 405–412. [6] J.M. Donelson, P.L. Munday, M.I. McCormick, N.W. Pankhurst, P.M. Pankhurst, Effects of elevated water temperature and food availability on the reproductive performance of a coral reef fish, Mar. Ecol. Prog. Ser. 401 (2010) 233–243. [7] J.L. Johansen, G.P. Jones, Increasing ocean temperature reduces the metabolic performance and swimming ability of coral reef damselfishes, Glob. Chang Biol. 17 (9) (2011) 2971–2979. [8] O.T. Hogg, D.K.A. Barnes, H.J. Griffiths, Highly diverse, poorly studied and uniquely threatened by climate change: an assessment of marine biodiversity on south georgia’s continental shelf, PLoS One 6 (5) (2011). E19795 [9] B.E. Beisner, E. McCauley, F.J. Wrona, The influence of temperature and food chain length on plankton predator prey dynamics, Can. J. Fish. Aquat. Sci. 54 (3) (1997) 586–595. [10] P. Amarasekare, Effects of temperature on consumer–resource interactions, J. Anim. Ecol. 84 (3) (2015) 665–679. [11] P.A. Staehr, K.A.J. Sand-Jensen, Seasonal changes in temperature and nutrient control of photosynthesis, respiration and growth of natural phytoplankton communities, Freshw. Biol. 51 (2) (2006) 249–262. [12] E. Bestion, B. García-Carreras, C.E. Schaum, S. Pawar, G. Yvon-Durocher, Metabolic traits predict the effects of warming on phytoplankton competition, Ecol. Lett. 21 (5) (2018) 655–664.

18

Q. Zhao, S. Liu and X. Niu / Applied Mathematics and Computation 378 (2020) 125211

[13] A.M. Lewandowska, P. Breithaupt, H. Hillebrand, H.G. Hoppe, K. Jürgens, U. Sommer, Responses of primary productivity to increased temperature and phytoplankton diversity, J. Sea Res. 72 (2012) 87–93. [14] A. Toseland, S.J. Daines, J.R. Clark, A. Kirkham, J. Strauss, C. Uhlig, M.T. Lenton, K. Valentin, G.A. Pearson, V. Moulton, T. Mock, The impact of temperature on marine phytoplankton resource allocation and metabolis, Nat Clim Change 3 (11) (2013) 979. [15] C.E. Hare, K. Leblanc, G.R. DiTullio, R.M. Kudela, Y.H. Zhang, P.A. Lee, S. Riseman, D.A. Hutchins, Consequences of increased temperature and CO2 for phytoplankton community structure in the bering sea, Mar. Ecol. Prog. Ser. 352 (2007) 9–16. [16] K.E. Havens, R.M. Pinto-Coelho, M. Bekliogˇ lu, K.S. Christoffersen, E. Jeppesen, T. Lauridsen, A. Mazumder, G. Méthot, B.P. Alloul, U.N. Tavs¸ anogˇ lu, S¸ Erdogˇ an., J. Vijverberg, Temperature effects on body size of freshwater crustacean zooplankton from greenland to the tropics, Hydrobiologia 743 (1) (2015) 27–35. [17] A.M. Kramer, O. Sarnelle, J. Yen, The effect of mating behavior and temperature variation on the critical population density of a freshwater copepod, Limnol. Oceanogr. 56 (2) (2011) 707–715. [18] G. Kalinkat, F.D. Schneider, C. Digel, C. Guill, C.B. Rall, U. Brose, Body masses, functional responses and predator–prey stability, Ecol. Lett. 16 (9) (2013) 1126–1134. [19] O. Sarnelle, A.E. Wilson, Type III functional response in daphnia, Ecology 89 (6) (2008) 1723–1732. [20] W. Uszko, S. Diehl, G. Englund, P. Amarasekare, Effects of warming on predator–prey interactions-a resource-based approach and a theoretical synthesis, Ecol. Lett. 20 (4) (2017) 513–523. [21] M. Lindmark, M. Huss, J. Ohlberger, A. Gardmark, Temperature-dependent body size effects determine population responses to climate warming, Ecol. Lett. 21 (2) (2018) 181–189. [22] M.K. Thomas, C.T. Kremer, C.A. Klausmeier, E. Litchman, A global pattern of thermal adaptation in marine phytoplankton, Science 338 (2012) 1085–1088. [23] G. Englund, G. Öhlund, C.L. Hein, S. Diehl, Temperature dependence of the functional response, Ecol. Lett. 14 (9) (2011) 914–921. [24] N. Chitnis, J.M. Hyman, J.M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol. 70 (5) (2008) 1272. [25] Z.C. Jiang, X.H. Bi, T.Q. Zhang, S.A. Pradeep, B.G. Global, Hopf bifurcation of a delayed phytoplankton–zooplankton system considering toxin producing effect and delay dependent coefficient, Math. Biosci. Eng. 16 (5) (2019) 3807–3829. [26] Z.C. Jiang, W.Z. Zhang, J. Zhang, T.Q. Zhang, Dynamical analysis of a phytoplankton–zooplankton system with harvesting term and holling III functional response, Int. J. Bifurcat. Chaos 28 (13) (2018) 1850162. [27] A.J. Richardson, D.S. Schoeman, Climate impact on plankton ecosystems in the northeast atlantic, Science 305 (5690) (2004) 1609–1612. [28] M.A. MacNeil, N.A.J. Graham, J.E. Cinner, N.K. Dulvy, P.A. Loring, S. Jennings, N.V.C. Polunin, A.T. Fisk, T.R. McClanahan, Transitional states in marine fisheries: adapting to predicted global change, Philosop. Trans. R. Soc. B: Biol. Sci. 365 (1558) (2010) 3753–3763.