Stability and Hopf Bifurcation Analysis of a Nutrient-Phytoplankton Model with Delay Effect

Hindawi Publishing Corporation Abstract and Applied Analysis Volume 2014, Article ID 471507, 10 pages http://dx.doi.org/10.1155/2014/471507 Research Article Stability and Hopf Bifurcation Analysis of a Nutrient-Phytoplankton Model with Delay Effect Xinhong Pan,1,2 Min Zhao,2,3 Chuanjun Dai,1,2 and Yapei Wang1,2 1 School of Mathematics and Information Science, Wenzhou University, Wenzhou, Zhejiang 325035, China Zhejiang Provincial Key Laboratory for Water Environment and Marine Biological Resources Protection, Wenzhou University, Wenzhou, Zhejiang 325035, China 3 School of Life and Environmental Science, Wenzhou University, Wenzhou, Zhejiang 325035, China 2 Correspondence should be addressed to Min Zhao; Received 13 November 2013; Accepted 25 January 2014; Published 12 March 2014 Academic Editor: Imran Naeem Copyright © 2014 Xinhong Pan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. A delay differential system is investigated based on a previously proposed nutrient-phytoplankton model. The time delay is regarded as a bifurcation parameter. Our aim is to determine how the time delay affects the system. First, we study the existence and local stability of two equilibria using the characteristic equation and identify the condition where a Hopf bifurcation can occur. Second, the formulae that determine the direction of the Hopf bifurcation and the stability of periodic solutions are obtained using the normal form and the center manifold theory. Furthermore, our main results are illustrated using numerical simulations. 1. Introduction Phytoplankton plays a very important role as the first trophic level in aquatic ecosystems. To describe the complex dynamics of phytoplankton populations, the dynamic relationship between phytoplankton and nutrients has been investigated theoretically for a long time, as well as experimentally. Since the pioneering work of Riley et al. [1], various nutrientphytoplankton models have been proposed and analyzed [2–9]. Several of these models have been shown to predict phytoplankton dynamics successfully in specific situations. The model proposed by Taylor et al. [10] in 1986 describes the nutrient-dependent growth of a single phytoplankton population by considering sinking and variable vertical mixing: 𝜕𝑁 𝜕 𝜕𝑁 = −𝛾 (𝛼 (𝐼) 𝜙 (𝑁) − 𝜀𝑚) 𝑃 + (𝐾 (𝑧, 𝑡) ), 𝜕𝑡 𝜕𝑧 𝜕𝑧 𝜕𝑃 𝜕𝑃 𝜕𝑃 𝜕 = (𝛼 (𝐼) 𝜙 (𝑁) − 𝑚) 𝑃 − V + (𝐾 (𝑧, 𝑡) ) , 𝜕𝑡 𝜕𝑧 𝜕𝑧 𝜕𝑧 the phytoplankton is taken to be the product of two parts, a dependence (𝛼) on incident light (𝐼) and a function (𝜙) of the concentration of a single nutrient; 𝑚 is the death rate of phytoplankton; 𝜀 (0 < 𝜀 ≤ 1) is the regeneration efficiency; V is the sinking rate of phytoplankton; 𝐾(𝑧, 𝑡) is the turbulent diffusion coefficient. The diffusion coefficients for phytoplankton and nutrients are assumed to be the same for simplicity. It is well known that the abundance of the phytoplankton population is affected by many environmental factors, such as the water temperature, salinity, and sunlight intensity [11]. In system (1), 𝐼 is the light intensity and numerical simulation indicates that the light intensity can affect the result. In this study, we consider an approximated model of system (1) with delay effect as a model of a layer of phytoplankton growing over a pool of nutrients: 𝑘 𝑑𝑁 = −𝛾 (𝛽𝑁 − 𝜀𝑚) 𝑃 + (𝑁0 − 𝑁) , 𝑑𝑡 ℎ (1) where 𝑁 and 𝑃 are the concentrations of the nutrient and phytoplankton, respectively; the specific growth rate of (2) 𝑑𝑃 V 𝑘 = (𝛽𝑁 (𝑡 − 𝜏) − 𝑚 − − ) 𝑃, 𝑑𝑡 ℎ ℎ where 𝛽 is the specific growth rate of phytoplankton; ℎ is the thickness of the layer; 𝑘 is the turbulent diffusion 2 Abstract and Applied Analysis coefficient; 𝜏 is a positive delay, that is, the time required to convert nutrients into phytoplankton; 𝑁0 and 𝑃0 are the concentrations of nutrients and phytoplankton below the layer, respectively. We assume that all parameters are positive, except that 𝑃0 is negligible and equal to zero. Taylor et al. [10] described the dynamics of the above system without considering the effect of delay, mainly using numerical methods. Pardo [12] conducted a mathematical study of the local and global stability of the equilibria in the same system. He proved the positivity and boundedness of the solutions, which made sure that the model is biologically sound. He found that the interior equilibrium point of the system is locally and globally asymptotically stable if it exists and that the boundary equilibrium point is also globally asymptotically stable if the system has only one equilibrium point. Time delay is known to play important roles in biological dynamical systems, which have been studied by many researchers in recent years [13–27]. Our aim is to investigate how time delay may affect the system (2). Thus, we use the delay 𝜏 as a bifurcation parameter. The remainder of this paper is organized as follows. In Section 2, we consider the local stability of the equilibria and the condition where Hopf bifurcation can occur based on the characteristic equation. In Section 3, we derive an explicit algorithm to determine the direction of the Hopf bifurcation and the stability of the periodic solutions. The results of numerical simulations are presented to support the theoretical results in Section 4. The paper ends with a brief conclusion. In this section, we mainly consider the existence and stability of the nonnegative equilibria of system (2). The equations for the equilibria are as follows: (𝛽𝑁 − 𝑚 − 𝑘 (𝑁 − 𝑁) = 0, ℎ 0 V 𝑘 − ) 𝑃 = 0. ℎ ℎ (3) There are two solutions: 𝐸0 (𝑁0 , 0), 𝐸∗ (𝑁∗ , 𝑃∗ ) if 𝑃 is not equal to zero and 𝑁0 > 𝑁∗ exists, where 𝑁∗ = 𝑃∗ = 1 V 𝑘 (𝑚 + + ) , 𝛽 ℎ ℎ ∗ 𝑘 (𝑁0 − 𝑁 ) 𝑘 (𝛽𝑁0 − 𝑚 − V/ℎ − 𝑘/ℎ) = . 𝛾ℎ (𝛽𝑁∗ − 𝜀𝑚) 𝛽𝛾ℎ (𝑚 + V/ℎ + 𝑘/ℎ − 𝜀𝑚) Theorem 1. The equilibrium 𝐸0 (𝑁0 , 0) is stable if 𝛽𝑁0 − 𝑚 − V/ℎ − 𝑘/ℎ < 0 and unstable if 𝛽𝑁0 − 𝑚 − V/ℎ − 𝑘/ℎ > 0. We recall that 𝛽𝑁0 − 𝑚 − V/ℎ − 𝑘/ℎ > 0 is equivalent to 𝑁0 > 𝑁∗ , which is the existing condition of the unique interior equilibrium. Thus, 𝐸0 is stable only if 𝐸∗ does not exist. 2.2. Local Stability and the Hopf Bifurcation of Equilibrium 𝐸∗ . Similarly, we linearize (2) about 𝐸∗ to obtain the linear system 𝑑𝑁 𝑘 = (−𝛾𝛽𝑃∗ − ) 𝑁 − 𝛾 (𝛽𝑁∗ − 𝜀𝑚) 𝑃, 𝑑𝑡 ℎ 𝑑𝑃 = 𝛽𝑃∗ 𝑁 (𝑡 − 𝜏) . 𝑑𝑡 The corresponding characteristic equation is 𝑘 𝜆2 + (𝛾𝛽𝑃∗ + ) 𝜆 + 𝛾𝛽 (𝛽𝑁∗ − 𝜀𝑚) 𝑃∗ 𝑒−𝜆𝜏 = 0 ℎ (7) 𝜆2 + 𝐴𝜆 + 𝐵𝑒−𝜆𝜏 = 0, (8) or where 𝑘 𝐴 = (𝛾𝛽𝑃∗ + ) > 0, ℎ 𝑑𝑃 V 𝑘 = (𝛽𝑁0 − 𝑚 − − ) 𝑃. 𝑑𝑡 ℎ ℎ 𝐵 = 𝛾𝛽 (𝛽𝑁∗ − 𝜀𝑚) 𝑃∗ > 0. (9) 𝜆2 + 𝐴𝜆 + 𝐵 = 0, (10) and the two eigenvalues satisfy 𝜆 1 + 𝜆 2 = −𝐴 < 0, 𝜆 1 ⋅ 𝜆 2 = 𝐵 > 0, (11) which indicates that 𝐸∗ is locally asymptotically stable when 𝜏 = 0. When 𝜏 ≠ 0, we assume that 𝜆 = 𝜇(𝜏) + 𝑖𝑤(𝜏) is a root of (8). Substituting this in (8) we have 𝜇2 − 𝑤2 + 2𝜇𝑤𝑖 + 𝜇𝐴 + 𝐴𝑤𝑖 + 𝐵𝑒−𝜇𝜏 (cos 𝜏𝑤 − 𝑖 sin 𝜏𝑤) = 0 (12) or 𝜇2 − 𝑤2 + 𝜇𝐴 = −𝐵𝑒−𝜇𝜏 cos 𝜏𝑤, (4) 2 (...truncated)