Dynamics of Oxygen-Plankton Model with Variable Zooplankton Search Rate in Deterministic and Fluctuating Environments
Abstract
This research was funded by the Spanish Government and European Commission for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19
Full text
Citation: Mondal, S.; Samanta, G.; De la Sen, M. Dynamics of Oxygen-Plankton Model with Variable Zooplankton Search Rate in Deterministic and Fluctuating Environments. Mathematics 2022,10, 1641. https://doi.org/10.3390/ math10101641 Academic Editors: Sophia Jang and Jui-Ling Yu Received: 3 April 2022 Accepted: 5 May 2022 Published: 11 May 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). mathematics Article Dynamics of Oxygen-Plankton Model with Variable Zooplankton Search Rate in Deterministic and Fluctuating Environments Sudeshna Mondal 1, Guruprasad Samanta 1and Manuel De la Sen 2,* 1Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India; [email protected] (S.M.); [email protected] or [email protected] (G.S.) 2Institute of Research and Development of Processes, University of the Basque Country, 48940 Leioa, Bizkaia, Spain *Correspondence: [email protected] Abstract: It is estimated by scientists that 50–80% of the oxygen production on the planet comes from the oceans due to the photosynthetic activity of phytoplankton. Some of this production is consumed by both phytoplankton and zooplankton for cellular respiration. In this article, we have analyzed the dynamics of the oxygen-plankton model with a modified Holling type II functional response, based on the premise that zooplankton has a variable search rate, rather than constant, which is ecologically meaningful. The positivity and uniform boundedness of the studied system prove that the model is well-behaved. The feasibility conditions and stability criteria of each equilibrium point are discussed. Next, the occurrence of local bifurcations are exhibited taking each of the vital system parameters as a bifurcation parameter. Numerical simulations are illustrated to verify the analytical outcomes. Our findings show that (i) the system dynamics change abruptly for a low oxygen production rate, resulting in depletion of oxygen and plankton extinction; (ii) the proposed system has oscillatory behavior in an intermediate range of oxygen production rates; (iii) it always has a stable coexistence steady state for a high oxygen production rate, which is dissimilar to the outcome of the model of a coupled oxygen-plankton dynamics where zooplankton consumes phytoplankton with classical Holling type II functional response. Lastly, the effect of environmental stochasticity is studied numerically, corresponding to our proposed system. Keywords: oxygen-plankton model; modified Holling type II; stability analysis; local bifurcations MSC: 37M05; 92D25; 92D40 1. Introduction Plankton are the numerous series of organisms observed in water or air that are not able to propel themselves against water currents or wind, respectively. The individual organisms constituting plankton are known as plankters. In the ocean, they offer a vital source of meals to many small and massive aquatic organisms, including bivalves, fish and whales. The plant types of the plankton community are referred to as phytoplankton, they acquire their strength through photosynthesis, as do trees and different plants on land. This means phytoplankton need to have solar light, so they live within the properly-lit floor layers of oceans and lakes. Zooplankton are the animal components of the planktonic network, and they are the principle food supply for fish and other aquatic animals. Phytoplankton are not the best meal source for zooplankton; however, they offer a massive quantity of oxygen for human and different dwelling animals after soaking up carbon dioxide via photosynthesis from the environment. Some of this oxygen production is consumed by both phytoplankton and zooplankton because of respiration [ 1 , 2 ]. Furthermore, a decrease in the Mathematics 2022,10, 1641. https://doi.org/10.3390/math10101641 https://www.mdpi.com/journal/mathematics
Mathematics 2022,10, 1641 2 of 24 oxygen production rate by phytoplankton may have a disastrous effect for living animals, including humankind. Therefore, the study of the possible range of oxygen production rates is important to sustain system dynamics. Mathematical modeling is a research tool that can reveal the dynamic properties of the oxygen-plankton model. Recently, researchers have analyzed a mathematical model of oxygen-plankton interactions witha Holling type II functional response [ 3 – 5 ], where the search rate of the predator population was constant, i.e., independent of the prey population [ 6 – 8 ]. However, it seems reasonable that predators can vary their search rates based on the availability of prey. In 1977, Hassel et al. [ 9 ] experimentally observed that the search rate of various invertebrate predators, specifically zooplankton, depended on the biomass of the prey (phytoplankton) population. In 2020, Dalziel et al. [ 10 ] analyzed the dynamics of a predator–prey model with a variable predator search rate. In 2021, Mondal and Samanta [ 11 ] studied the dynamic nature of a predator–prey model with the impact of a predator’s fear, where the search rate of the predator depended on the biomass of the prey species. Recently, they also investigated the dynamic behavior of a toxin-producing plankton model where the zooplankton’s search rate depended on the biomass of the phytoplankton population, rather than being assumed constant [ 12 ]. Motivated from the above discussions, we proposed and analyzed the dynamic behavior of the oxygen-plankton model with a variable zooplankton search rate, rather than constant, where oxygen is produced by the photosynthetic activity of phytoplankton during the daytime and consumed by phyto and zooplankton for their respiration. This article is organized as follows: we have focused on the construction of the basic model in Section 2. The derivation of positivity and uniform boundedness is shown in Section 3. Section 4describes the feasibility criteria and stability conditions of all the equilibria. Furthermore, the occurrence of local bifurcations are exhibited in Section 5. In Section 6, we conduct numerical simulations using MATLAB to validate the analytical findings. The impact of the oxygen production rate on the existence of the interior equilibrium point as well as the main qualitative difference between the proposed model and the system analyzed by Sekerci and Petrovskii [ 3 ] are discussed. This section also consists of the effect of environmental stochasticity on the proposed oxygen-plankton model by perturbing some parameters of the system with Gaussian white noise terms. This work ends with a discussion and the outcomes of the analytical consequences. 2. Construction of Basic Model A marine ecosystem is a complicated system with many nonlinearly interacting species, organic substances, and inorganic chemical components. Correspondingly, a "realistic” ecosystem model can consist of many equations. In this article, we are mostly interested in the dynamics of the oxygen-plankton model, where oxygen is produced by the photosynthetic activity of phytoplankton. Revisiting an oxygen-plankton model system given in [ 3 , 5 ] and taking a modified Holling type II functional response, where the search rate of the predator (zooplankton) depends on the biomass of the prey (phytoplankton), rather than being constant (for details, see [10–12]), we consider the following model (for details see Figure 1): dc dt =Ac0p c+c0−δcp c+c2−νcz c+c3−mc dp dt =Bc c+c1−γpp−ap2z ahp2+p+g−σp(1) dz dt = ηc2 c2+c2 4!.ap2z ahp2+p+g−µz with initial conditions: c(0)>0, p(0)>0, z(0)>0. (2)
Mathematics 2022,10, 1641 3 of 24 Here, c is the amount of oxygen, and p and z are the biomass of phytoplankton and zooplankton, respectively. All the parameters are positive due to their biological meaning and are described in Table 1: Table 1. Description of biologically meaningful parameters. Parameters Descriptions Aeffect of environmental factors on the rate of oxygen production due to the photosynthesis of phytoplankton δmaximum per capita phytoplankton respiration rate νmaximum per capita zooplankton respiration rate mrate of oxygen loss due to the biochemical reaction in a marine ecosystem Bmaximum phytoplankton per capita growth rate in the high oxygen limit ci,i=0, 1, 2, 3, 4 half saturation constant of the corresponding processes γmortality rate due to intraspecific competition among individual phytoplankton amaximally achievable search rate of zooplankton hhandling time of zooplankton ghalf saturation constant σnatural mortality rate of phytoplankton. It is assumed that B>σ η∈(0, 1)maximum feeding efficiency µmortality rate of zooplankton μz ap2z ahp2+p+g νcz c+c3 δcp c+c2 mc γp2+σ p Ac0p c+c0 zooplankton phytoplankton oxygen Figure 1. Graphical scheme representing the interactions among oxygen, phytoplankton, and zooplankton, where phytoplankton produce oxygen through photosynthetic activity in sunlight and consume it during the night for their respiration; zooplankton depend on phytoplankton for their growth and consume oxygen for their respiration. Description of system (1): • The term Ac0 c+c0 describes the rate of oxygen production per unit of phytoplankton biomass during the daytime by photosynthetic activity; δcp c+c2 and νcz c+c3 indicate the respiration of phytoplankton and zooplankton, respectively, and mc is the loss of oxygen due to natural depletion in a marine ecosystem. • The term Bcp c+c1 describes the growth of phytoplankton depending on the amount of available oxygen. The function ap2 ahp2+p+g is named as a modified Holling type II functional response, based on the premise that the zooplankton’s search rate is dependent on the biomass of phytoplankton, rather than being constant (for details, see [10,11] ). Again, the consumed phytoplankton biomass is transformed into zooplankton biomass
Mathematics 2022,10, 1641 4 of 24 with an efficiency of ηc2 c2+c2 4 , which depends on the oxygen concentration (zooplankton die due to insufficient oxygen). The following are properties of a modified Holling type II functional response H(p) = ap2 ahp2+p+g 1. H(p)is a smooth function, and H(p) = 0 for p=0. 2. H0(p) = ap(p+2g) (ahp2+p+g)2> 0, i.e., H increases with p and lim p→∞H(p) = 1 h , i.e., H(p) saturates at 1 hfor a large prey population. 3. H00(p) = −2a2hp3−6a2ghp2+2ag2 (ahp2+p+g)3 , and H00(p)p=0=2a g> 0. Therefore, H00(p) has a unique positive root, and it changes sign from positive to negative at the unique inflection point. A graphical representation of H(p) and H00(p) is presented in Figure 2. 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 p H(p) (a)pverses H(p) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −4 −2 0 2 4 6 8 10 p H"(p) (b)pverses H00(p) Figure 2. Graphical representations of ( a ) H(p) and ( b ) H00(p) for the parametric set {a= 3, h= 1.2, g=0.3}. 3. Positivity and Uniform Boundedness Theorem 1. Solutions of (1) with (2) exist uniquely and are positive for all t ≥0. Proof. Since the right hand sides of (1) are completely continuous functions and locally Lipschitzian in the domain R3 + , solutions of (1) with (2) exist uniquely in [ 0, ξ) , where 0<ξ≤∞[13]. From the first equation of (1), we have: c(t) = c(0)exp−Zt 0δp(θ) c(θ) + c2 +νz(θ) c(θ) + c3 +md(θ) +Zt 0 Ac0p(u) c(u) + c0expZu tδp(θ) c(θ) + c2 +νz(θ) c(θ) + c3 +md(θ)du >0, since c(0)>0. From the second equation of system (1), we have: p(t) = p(0)expZt 0Bc(θ) c(θ) + c1−γp(θ)−ap(θ)z(θ) ahp2(θ) + p(θ) + g−σdθ>0, since p(0)>0. From the last equation of system (1), we have: z(t) = z(0)exp"Zt 0( ηc2(θ) c2(θ) + c2 4!.ap2(θ) ahp2(θ) + p(θ) + g−µ)dθ#>0, since z(0)>0.
Mathematics 2022,10, 1641 5 of 24 Therefore, c(t)>0, p(t)>0 and z(t)>0 for all t≥0. Hence, the theorem is proved. Theorem 2. Solutions of (1) with (2) are uniformly bounded. Proof. From the second equation of system (1), we obtain: dp dt ≤Bp −γp2−σp = (B−σ)p(1−p B−σ γ) ⇒lim sup t→∞ p(t)≤B−σ γ. Let Ω=c+p+z. Then, dΩ dt =dc dt +dp dt +dz dt =Ac0p c+c0−δcp c+c2−νcz c+c3−mc +Bc c+c1−γpp−ap2z ahp2+p+g−σp + ηc2 c2+c2 4!ap2z ahp2+p+g−µz ≤Ac0p c+c0 +Bcp c+c1 +ap2z ahp2+p+g ηc2 c2+c2 4−1!−γp2−{mc +σp+µz} ≤Ac0p c+c0 +Bcp c+c1−γp2−{mc +σp+µz}, since 0 <η<1 ≤(A+B)p−γp2−{mc +σp+µz} ≤(A+B)2 4γ−{mc +σp+µz}. (3) Let κ=min{m,σ,µ}. Then, from (3), we obtain: dΩ dt +κΩ≤(A+B)2 4γ. Using the differential inequality: 0<Ω(c(t),p(t),z(t))≤(A+B)2 4γκ 1−e−κt+e−κtΩ(c(0),p(0),z(0)). ∴0<Ω(c(t),p(t),z(t))≤(A+B)2 4γκ +e, for any e>0, as t→∞. Hence, every solution of (1) enters into the region: W=(c,p,z)∈R3 +: 0 <p(t)≤B−σ γ; 0 <c(t) + p(t) + z(t)≤(A+B)2 4γκ +e,e>0.
Mathematics 2022,10, 1641 6 of 24 4. Existence of Equilibria of (1) with Stability Analysis 4.1. Equilibrium Points System (1) has the following equilibrium points (steady states): 1. Trivial equilibrium point E0( 0, 0, 0 ) corresponding to depletion of oxygen and the extinction of plankton; 2. Planer equilibrium point E1(e c,e p, 0) (zooplankton free), where e p=1 γBe c e c+c1−σ , and e cis a positive root of the following equation: X1c4+X2c3+X3c2+X4c+X5=0. Here, X1=−mγ , X2=−(c1+c2)−c0+ (B−γ)δ , X3=−c1c2−c0(c1+c2) + (B− γ)(A−δ)c0+δγ,X4=−c0c1c2+ (B−γ)Ac0c2−γc0(A−δ),X5=−γAc0c1c2. 3. Interior (coexistence) equilibrium b E(b c , b p , bz) , where b c , b p , and bz can be obtained by solving the following system of equations using the software MATHEMATICA: Ac0p c+c0−δcp c+c2−νcz c+c3−mc =0, Bc c+c1−γp−apz ahp2+p+g−σ=0, ηc2 c2+c2 4!.ap2 ahp2+p+g−µ=0. 4.2. Local Stability Now, we will determine the stability behavior of the biologically feasible equilibrium points of system (1). The Jacobian matrix J0at E0(0, 0, 0)is given by: J0= −m A 0 0−σ0 0 0 −µ . Here, the eigenvalues are λ1=−m< 0, λ2=−σ< 0, and λ3=−µ< 0. Since all eigenvalues are negative, so, E0(0, 0, 0)is always locally asymptotically stable (LAS). The Jacobian matrix J1at E1(e c,e p, 0)is given by: J1= −Ac0e p (e c+c0)2−δc2e p (e c+c2)2−mme c e p−νe c e c+c3 Bc1e p (e c+c1)2−γe p−ae p2 ahe p2+e p+g 0 0 ηe c2 e c2+c2 4ae p2 ahe p2+e p+g−µ . Here, one eigenvalue is λ1=ηe c2 e c2+c2 4ae p2 ahe p2+e p+g−µ , and the other eigenvalues can be obtained by solving the equation: λ2−Q1λ+Q2=0, (4) where Q1=−Ac0e p (e c+c0)2−δc2e p (e c+c2)2−m−γe p< 0 and Q2=γe phAc0e p (e c+c0)2+δc2e p (e c+c2)2+mi− Bc1me c (e c+c1)2>0. Hence, we have the following theorem: Theorem 3. E1(e c,e p, 0)is LAS if ηe c2 e c2+c2 4ae p2 ahe p2+e p2+g−µ<0.
Mathematics 2022,10, 1641 7 of 24 The Jacobian matrix bJat b E(b c,b p,bz)is given by: bJ= a11 a12 a13 a21 a22 a23 a31 a32 a33 where a11 =−Ac0b p (b c+c0)2−δc2b p (b c+c2)2−νc3bz (b c+c3)2−m< 0, a12 =−δb c b c+c2+Ac0 b c+c0=b c b pnνbz b c+c3+mo> 0, a13 =−νb c b c+c3 < 0, a21 =Bc1b p (b c+c1)2> 0, a22 =Bb c b c+c1− 2 γb p−ab pbz(b p+2g) (ahb p2+b p+g)2−σ=−γb p− ab pbz(g−ahb p2) (ahb p2+b p+g)2 , a23 =−ab p2 ahb p2+b p+g< 0, a31 =2ηc2 4b c (b c2+c2 4)2 ab p2bz ahb p2+b p+g> 0, a32 =ηb c2 (b c2+c2 4) ab pbz(b p+2g) (ahb p2+b p+g)2>0 and a33 =0. The characteristic equation corresponding to b E(b c,b p,bz)is λ3+C1λ2+C2λ+C3=0 where C1=−(a11 +a22) , C2=−a23a32 −a13a31 +a11a22 −a12a21 , and C3=−{−a11a23a32 + a12a23a31 +a13(a21a32 −a22a31)}. By Routh-Hurwitz’s criteria [ 14 ], b E(b c , b p , bz) has three eigenvalues with negative real parts if C1> 0, C3> 0, and C1C2>C3 . So, the local stability condition of b E(b c , b p , bz) is described in the following theorem: Theorem 4. b E(b c,b p,bz)is LAS if a22 <0and a11a22 >a12a21. 5. Local Bifurcations A local bifurcation occurs when a parameter change causes the stability (or instability) of an equilibrium (or fixed point) to change. In continuous systems, this corresponds to the real part of an eigenvalue of an equilibrium passing through zero. 5.1. Transcritical Bifurcation Theorem 5. System (1) undergoes a transcritical bifurcation if µ[tc]=ηe c2 e c2+c2 4ae p2 ahe p2+e p+g. Proof. To prove a transcritical bifurcation, we apply Sotomayor’s theorem [ 14 ] by considering µ as the bifurcation parameter. According to this theorem, one eigenvalue of J1 at the bifurcation point must be zero. The eigenvectors of J1= [pij] and (J1)T corresponding to the zero eigenvalue are obtained as: V=(0, v2, 1)T and W=(0, 0, 1)T , respectively, where v2=−p13 p12 and p11 =−Ac0e p (e c+c0)2−δc2e p (e c+c2)2−m , p12 =me c e p , p13 =−νe c e c+c3 , p21 =Bc1e p (e c+c1)2 , p22 =−γe p , p23 =−ae p2 ahe p2+e p+g, and p31 =p32 =p33 =0. Compute ∆1,∆2, and ∆3as follows: ∆1=WT·Fµe c,e p, 0; µ[tc]= (0, 0, 1)· ∂F1 ∂µ ∂F2 ∂µ ∂F3 ∂µ (E1(e c,e p,0);µ[tc]) ⇒∆1= (0, 0, 1)· 0 0 −z (E1(e c,e p,0);µ[tc]) =0, where F=(F1,F2,F3)T, and F1,F2, and F3are given by: F1=Ac0p c+c0−δcp c+c2−νcz c+c3−mc, F2=Bc c+c1−γpp−ap2z ahp2+p+g−σp, F3=ηc2 c2+c2 4·ap2z ahp2+p+g−µz.
Mathematics 2022,10, 1641 8 of 24 ∆2=WT·hDFµe c,e p, 0; µ[tc]Vi= (0, 0, 1)· ∂2F1 ∂c∂µ ∂2F1 ∂p∂µ ∂2F1 ∂z∂µ ∂2F2 ∂c∂µ ∂2F2 ∂p∂µ ∂2F2 ∂z∂µ ∂2F3 ∂c∂µ ∂2F3 ∂p∂µ ∂2F3 ∂z∂µ (E1(e c,e p,0);µ[tc]) · 0 v2 1 ⇒∆2= (0, 0, 1)· 0 0 0 0 0 0 0 0 −1 (E1(e c,e p,0);µ[tc]) . 0 v2 1 =−16=0. ∆3=WT·hD2Fe c,e p, 0; µ[tc](V,V)i= (0, 0, 1)·D ∂F1 ∂cv1+∂F1 ∂pv2+∂F1 ∂zv3 ∂F2 ∂cv1+∂F2 ∂pv2+∂F2 ∂zv3 ∂F3 ∂cv1+∂F3 ∂pv2+∂F3 ∂zv3 (E1(e c,e p,0);µ[tc]) . v1 v2 v3 ⇒∆3= (0, 0, 1)· ∂2F1 ∂2cv2 1+∂2F1 ∂2pv2 2+∂2F1 ∂2zv2 3+2∂2F1 ∂c∂pv1v2+2∂2F1 ∂c∂zv1v3+2∂2F1 ∂p∂zv2v3 ∂2F2 ∂2xv2 1+∂2F2 ∂2yv2 2+∂2F2 ∂2zv2 3+2∂2F2 ∂x∂yv1v2+2∂2F2 ∂x∂zv1v3+2∂2F2 ∂y∂zv2v3 ∂2F3 ∂2xv2 1+∂2F3 ∂2yv2 2+∂2F3 ∂2zv2 3+2∂2F3 ∂x∂yv1v2+2∂2F3 ∂x∂zv1v3+2∂2F3 ∂y∂zv2v3 (E1(e c,e p,0);µ[tc]) ⇒∆3=2ae p(e p+2g) (ahe p2+e p+g)2×ηe c2 (e c2+c2 4)v26=0. Thus, by Sotomayor’s theorem [ 14 ], system (1) exhibits a trancritical bifurcation at µ=µ[tc]. Remark 1. Similarly, it can be proved that system (1) exhibits transcritical bifurcations taking any one of the parameters h, σ, m, η, a, and γas a bifurcation parameter. 5.2. Hopf-Bifurcation The characteristic equation of system (1) at b E(b c,b p,bz)is given by λ3+C1(A)λ2+C2(A)λ+C3(A) = 0, (5) where Ci(A)for i=1, 2, 3 were defined earlier. To determine the Hopf-bifurcation around b E(b c , b p , bz) of system (1), let us consider A as the bifurcation parameter. For this purpose, let us first state the following Theorem: Theorem 6 (Hopf-Bifurcation Theorem [ 15 ]) . If C1(A) , C2(A) , and C3(A) are continuously differentiable functions of A in a small neighbourhood of A[H]∈Rsuch that Equation (5) has: (i) a pair of imaginary eigenvalues λ=p1(A)±ip2(A) with p1(A)∈R , p2(A)∈R , so that they become purely imaginary at A =A[H]and dp1 dA |A=A[H]6=0, (ii) the other eigenvalue is negative at A=A[H] , then a Hopf-bifurcation occurs around b E(b c , b p , bz) at A=A[H] (i.e., a stability change of b E(b c , b p , bz) accompanied by the creation of a limit cycle at A =A[H]). Theorem 7. System (1) possesses a Hopf-bifurcation around b E(b c , b p , bz) when A passes through A[H], provided C1(A[H])>0, C3(A[H])>0, and C1(A[H])C2(A[H]) = C3(A[H]). Proof. At A=A[H], the roots of the equation: λ2+C2(λ+C1)=0 are λ1=i√C2 , λ2=−i√C2 , and λ3=−C1 , where C1 , C2 and C3 are differential functions of A . Furthermore, in the deleted neighborhood of A[H] , the roots (eigenvalues) are λ1(A) = p1(A) + ip2(A) , λ2(A) = p1(A)−ip2(A) , and λ3=p3(A) ( p3(A) = −C1 ), where pi(A) are real for i=1, 2, 3.
Mathematics 2022,10, 1641 9 of 24 Now, we will verify the transversality condition: d dA (Re λi(A))A=A[H]6=0, i=1, 2. Substituting λ(A) = p1(A) + ip2(A)into the characteristic Equation (5), we have: (p1+ip2)3+C1(A)(p1+ip2)2+C2(A)(p1+ip2)+C3(A) = 0 (6) Differentiating with regard to A, we have: 3(p1+ip2)2(˙ p1+i˙ p2)+2C1(p1+ip2)( ˙ p1+i˙ p2)+˙ C1(p1+ip2)2 +C2(˙ p1+i˙ p2)+˙ C2(p1+ip2)+˙ C3=0 (7) Comparing the real and imaginary parts, we obtain: X1˙ p1−X2˙ p2+X3=0 (8) and X2˙ p1+X1˙ p2+X4=0 (9) where X1=3p2 1−p2 2+2C1p1+C2 X2=6p1p2+2C1p2 X3=˙ C1p2 1−p2 2+˙ C2˙ p1+˙ C3 X4=2˙ C1p1p2+˙ C2p2. From (8) and (9), we obtain: ˙ p1=−(X1X3+X2X4) X2 1+X2 2 . Now, X3=˙ C1p2 1−p2 2+˙ C2p1+˙ C36=˙ C1p2 1−p2 2+˙ C2p1+˙ C1C2+C1˙ C2 [since C36=C1C2in a deleted neighborhood of A[H]] At A=A[H], •Case 1: p1=0, p2=√C2 X1=−2C2,X2=2C1√C2,X36=C1˙ C2,X4=√C2˙ C2 Therefore, X2X4+X1X36=2C1C2˙ C2−2C1C2˙ C2=0 So, X2X4+X1X36=0 at A=A[H], when p1=0, p2=√C2. •Case 2: p1=0, p2=−√C2 X1=−2C2,X2=−2C1√C2,X36=C1˙ C2,X4=−√C2˙ C2 So, X2X4+x1X36=2C1C2˙ C2−2C1C2˙ C2=0 So, X2X4+X1X36=0 at A=A[H], when p1=0, p2=−√C2. ∴d dA (Re λi(A))|A=A[H]6=0, for i=1, 2 and p3(A[H]) = −C1(A[H])<0. Hence, Theorem 7is proved using Theorem 6. Note: Imaginary eigenvalues are connected with any molecular process (e.g., collisions) and the reverse of that process [16].
Mathematics 2022,10, 1641 16 of 24 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5 2 2.5 3 3.5 4 4.5 µ Oxygen (a) Bifurcation diagram of c 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 µ Phytoplankton (b) Bifurcation diagram of p 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 µ Zooplankton (c) Bifurcation diagram of z Figure 13. Bifurcation diagram of system (1) taking µ as the bifurcation parameter, while the others remain unchanged, as in Figure 6. b E(b c , b p , bz) is unstable with a periodic solution when µ∈(µ[H]1=0.131580, µ[H]2= 0.354676 ) and stable when µ∈( 0, 0.131580 )∪( 0.354676, 0.474246 ) . Again, b E(b c,b p,bz)goes to stable zooplankton free equilibrium E1(e c,e p, 0)when µ>µ[tc]=0.474246. 0.5 1 1.5 2 2.5 −0.5 0 0.5 1 1.5 2 2.5 m Oxygen (a) Bifurcation diagram of c 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 m Phytoplankton (b) Bifurcation diagram of p 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 m Zooplankton (c) Bifurcation diagram of z Figure 14. Bifurcation diagrams of system (1) taking m as the bifurcation parameter, while the others remain unchanged, as in Figure 6. Here, b E(b c , b p , bz) is stable when m∈( 0.0, m[H]= 0.6533 ) and unstable with a periodic solution when m∈( 0.6533, 2.287 ) . When m>m[tc]= 2.287, trivial equilibrium E0= ( 0, 0, 0 ) exists corresponding to the depletion of oxygen and the extinction of plankton.
Mathematics 2022,10, 1641 17 of 24 0.2 0.4 0.6 0.8 0 1 2 3 4 5 g c H (a) Bifurcation diagram of c 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 g p H (b) Bifurcation diagram of p 0.2 0.4 0.6 0.8 0 1 2 3 g z H (c) Bifurcation diagram of z Figure 15. Bifurcation diagrams of system (1) while g varies from [ 0.09, 1 ] and the others remain unchanged, as in Figure 6. Here, the interior equilibrium b E(b c , b p , bz) is unstable with a periodic solution when g∈[0.09, g[H]=0.226067)and stable when g>g[H]=0.226067. 0 0.5 1 0 2 4 6 η c H H H H BP (a) Bifurcation diagram of c 0 0.5 1 0 1 2 3 4 η p H H BP H H (b) Bifurcation diagram of p 0 0.5 1 0 1 2 3 η z BP H H (c) Bifurcation diagram of z Figure 16. Bifurcation diagrams of system (1) taking η as the bifurcation parameter, while the others remain unchanged, as in Figure 6. Here, the zooplankton free equilibrium E1(e c , e p , 0 ) is stable when 0 <η<η[tc]= 0.147603, and the interior equilibrium b E(b c , b p , bz) is stable when η∈ ( 0.147603, η[H]1=0.198177)∪(η[H]2= 0.530864, 1 ) and unstable with a periodic solution when η∈(0.198177, 0.530864). Here, ‘BP’ stands for the transcritical bifurcation point.
Mathematics 2022,10, 1641 18 of 24 0 5 10 0 2 4 a c H H H BP (a) Bifurcation diagram of c 0 5 10 0 1 2 3 a p BP H (b) Bifurcation diagram of p 0 5 10 0 1 2 3 4 a z H H H BP (c) Bifurcation diagram of z Figure 17. Bifurcation diagrams of system (1) taking a as the bifurcation parameter, while the others remain unchanged, as in Figure 6. Here, the zooplankton free equilibrium E1(e c , e p , 0 ) is stable when 0 <a<a[tc]= 0.109643, and the interior equilibrium b E(b c , b p , bz) is stable when a∈( 0.109643, a[H]=5.325675)and unstable with periodic solution when a>a[H]=5.325675. 2 4 6 8 10 0 0.5 1 1.5 2 2.5 3 γ c stable interior equilibrium BP stable zooplankton free equilibrium Figure 18. Bifurcation diagram of system (1) taking γ as the bifurcation parameter, while the others remain unchanged, as in Figure 6. Here, ‘BP’ appears at γ=γ[tc]=4.479066. Effect of Environmental Noise on System (1) In a marine ecosystem, the oxygen-plankton model is affected by the environmental noise due to the inherent stochasticity of the weather conditions. For environmental noise, some of the parameters of system (1) change randomly over time. In this study, we have assumed that the stochasticity affects the oxygen production term through parameter A ,
Mathematics 2022,10, 1641 19 of 24 the phytoplankton growth term through parameter B , and the zooplankton mortality rate µby turning A,B, and µinto random variables as follows: A→A+γ1(t) B→B+γ2(t)(11) µ→µ+γ3(t) where γ1 , γ2 , and γ3 are independent Gaussian white noise terms and satisfy the following conditions: <γj(t)>=0 and <γj(t1),γj(t2)>=α2 jδj(t1−t2), for j=1, 2, 3 where αj are the intensities or strengths of the random perturbations, δ is the Dirac delta function defined by: (δ(x) = 0, for x6=0 R∞ −∞δ(x)dx =1 and <·>is the ensemble average of the considered stochastic process. Introducing Gaussian white noises, system (1) can be formulated as: dc dt =(A+γ1(t))c0p c+c0−δcp c+c2−νcz c+c3−mc dp dt =(B+γ2(t))cp c+c1−γp2−ap2z ahp2+p+g−σp dz dt = ηc2 c2+c2 4!.ap2z ahp2+p+g−(µ+γ3(t))z i.e., dc dt =Ac0p c+c0−δcp c+c2−νcz c+c3−mc +γ1(t)c0p c+c0 dp dt =Bcp c+c1−γp2−ap2z ahp2+p+g−σp+γ2(t)cp c+c1 dz dt = ηc2 c2+c2 4!.ap2z ahp2+p+g−µz−γ3(t)z dc dt =Ac0p c+c0−δcp c+c2−νcz c+c3−mc +c0p c+c0·α1 dw1 dt dp dt =Bcp c+c1−γp2−ap2z ahp2+p+g−σp+cp c+c1·α2 dw2 dt dz dt = ηc2 c2+c2 4!·ap2z ahp2+p+g−µz−α3zdw3 dt where γ1=α1dw1 dt , γ2=α2dw2 dt , and γ3=α3dw3 dt . Here, w=w1(t),w2(t),w3(t)t≥0 represents three-dimensional standard Brownian motion.
Mathematics 2022,10, 1641 20 of 24 Hence, our proposed stochastic system is: dc =Ac0p c+c0−δcp c+c2−νcz c+c3−mc +c0 c+c0α1pdw1 dp =Bcp c+c1−γp2−ap2z ahp2+p+g−σp+p c+c1α2cdw2(12) dz = ηc2 c2+c2 4!.ap2z ahp2+p+g−µz−α3zdw3. The effect of environmental noise on the dynamics of system (12) is analyzed numerically by the Euler Maruyama method in MATLAB. For this purpose, we chose the parametric set as follows: {c0=1, δ=1, c2=1, m=0.5, B=1.8, c1=1.0, γ=0.7, σ=0.1, ν=0.01, µ=0.1, h=1.2, a=3.0, η=0.7, c4=1, g=0.3, c3=1, α1=α2=α3=0.001}, (13) but varied Ain a broad range. When we took A= 10, while the other parameters remained the same as in set (13), then the effect of the Gaussian white noises on the stochastic system (12) were as depicted in Figure 19. Furthermore, Figure 19 shows that the oxygen, phytoplankton, and zooplankton varied around the deterministic coexistence steady-state values 1.48635, 0.218551, and 0.866809, respectively. Hence, system (12) is persistent. In this context, we repeated the stochastic simulations 20000 times, and the numerical results are depicted in Figure 20, which shows the stationary distribution of c(t) , p(t) , and z(t) at time t= 600. Moreover, when we chose A= 1.8, while the remaining parameters remained the same as in set (13), then system (12) was also persistent (see Figure 21). 0 20 40 60 80 100 120 140 160 180 200 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 Time t Oxygen 0 20 40 60 80 100 120 140 160 180 200 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 Time t Phytoplankton 0 20 40 60 80 100 120 140 160 180 200 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 Time t Zooplankton Figure 19. Stochastic trajectories of system (12) when A= 10 and the remaining parameters are same as in set (13). Initial conditions are c(0) = 1, p(0) = 0.3 and z(0) = 1.
Mathematics 2022,10, 1641 21 of 24 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 0 500 1000 1500 c Frequency 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0 200 400 600 800 1000 1200 p Frequency 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 0 100 200 300 400 500 600 700 800 z Frequency Figure 20. Histograms of system (12) with the parameters chosen from Figure 19. 0 100 200 300 400 500 600 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time t c, p, z c p z Figure 21. Persistence of system (12) when A=1.8 and the remaining parameters stay unaltered as in Figure 19. Again, if we take µ= 0.5, while the other parameters remain the same as in set (13), then, it is noted from Figure 22 that the zooplankton population can not persist in system (12) for any of the following choices: (a) A=1.8 and (b) A=10. Furthermore, it is observed from Figure 23 that system (12) becomes extinct for any of the following choices: (a) A= 1.5, (b) σ= 1.0, and (c) m= 2.9, while the other parameters remain the same as in set (13).
Mathematics 2022,10, 1641 22 of 24 0 20 40 60 80 100 120 140 160 180 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time t c, p, z c p z (a)A=1.8 and µ=0.5 0 20 40 60 80 100 120 140 160 180 200 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Time t c, p, z c p z (b)A=10.0 and µ=0.5 Figure 22. Extinction of the zooplankton in system (12) when ( a ) A= 1.8 and µ= 0.5, ( b ) A= 10.0 and µ=0.5 and remaining parameters are chosen from set (13). 0 20 40 60 80 100 120 140 160 180 200 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time t c, p, z c p z (a)A=1.5 0 10 20 30 40 50 60 70 80 90 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time t c, p, z c p z (b)σ=1.0 0 10 20 30 40 50 60 70 80 90 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time t c, p, z c p z (c)m=2.9 Figure 23. Depletion of oxygen and extinction of plankton corresponding to system (12) when (a)A=1.5, (b)σ=1.0, (c)m=2.9 and the remaining parameters are chosen from set (13). 7. Discussion and Conclusions A Holling type II functional response [ 6 – 8 ] is predicated on the assumption that the search rate of a predator is constant, i.e., independent of the prey population. However, it seems reasonable that the predator can vary their search rate based on the availability of prey. In particular, it is estimated that 50–80% of the oxygen production on Earth comes from the oceans due to the photosynthetic activity of phytoplankton. Some of this production is
Mathematics 2022,10, 1641 23 of 24 consumed by both phytoplankton and zooplankton for cellular respiration. Furthermore, zooplankton consume phytoplankton with a modified Holling type II functional response, based on the premise that the zooplankton search rate is dependent on phytoplankton (for details, see [ 10 , 11 ]). The goal of this article was to investigate the behavior of the oxygen-plankton model with a modified Holling type II functional response. The following summarizes our findings: • The coexistence steady state is stable when 1.8 ≤A< 1.966532, and it loses its stable nature through Hopf-bifurcation when 1.966532 <A< 7.258206 (see Figures 7and 8). • The dynamic behavior of system (1) changes abruptly for a low oxygen production rate (0 <A< 1.8), resulting in the depletion of oxygen and plankton extinction (see Figure 3). This depletion of oxygen production will be a consequence of the global ecological disaster. • System (1) always has a stable coexistence steady state for a high oxygen production rate (see Figure 6), i.e., the sustainability of oxygen production is possible when A is large ( A> 7.258206). This result is opposite to the outcome shown by Sekerci and Petrovskii [ 3 ] because they observed that the system dynamics were not sustainable for a high oxygen production rate. This is the main qualitative difference between the modified Holling type II (variable search rate, as mentioned in the proposed model) and the Holling type II functional responses. Therefore, the study of the modified Holling type II functional response is ecologically meaningful for the sustainability of the dynamics of system (1), if the net oxygen production rate is above a certain critical valve (A≥1.8). Moreover, the effect of environmental noise has a strong impact due to the inherent stochasticity of weather conditions. So, our proposed deterministic system (1) was compared with a corresponding stochastic model (12) incorporating Gaussian white noises in the system parameters A,B, and µ, as mentioned in (11). In the future, a realistic model can be proposed to explore the effects of spatial diffusion in the pattern formation through diffusion-driven instability. Author Contributions: Conceptualization, S.M., G.S. and M.D.l.S.; Methodology, S.M., G.S. and M.D.l.S.; Investigation, S.M., G.S. and M.D.l.S.; Formal analysis, S.M., G.S. and M.D.l.S.; Writing— original draft preparation, S.M., G.S. and M.D.l.S.; Writing—review and editing, S.M., G.S. and M.D.l.S. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Spanish Government and European Commission for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data used to support the findings of the study are available within the article. Acknowledgments: The authors are grateful to the anonymous referees, for their careful reading, valuable comments, and helpful suggestions, which have helped them to improve the presentation of this work significantly. The third author (Manuel De la Sen) is grateful to the Spanish Government and European Commission for its support through grant RTI2018-094336-B-I00 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19. Conflicts of Interest: The authors declare that they have no conflict of interest regarding this work. References 1. Harris, G.P. Phytoplankton Ecology: Structure, Function and Fluctuation; Springer: Dordrecht, The Netherlands, 1986. 2. Moss, B.R. Ecology of Fresh Waters: Man and Medium, Past to Future; Wiley: London, UK, 2009. 3. Sekerci, Y.; Petrovskii, S. Mathematical modelling of plankton-oxygen dynamics under the climate change. Bull. Math. Biol. 2015 , 77, 2325–2353. [CrossRef] [PubMed]
Mathematics 2022,10, 1641 24 of 24 4. Gokce, A.; Yazar, S.; Sekerci, Y. Delay induced nonlinear dynamics of oxygen-plankton interactions. Chaos Solitons Fractals 2020 , 141, 110327. [CrossRef] 5. Sekerci, Y.; Petrovskii, S. Global Warming Can Lead to Depletion of Oxygen by Disrupting Phytoplankton Photosynthesis: A Mathematical Modelling Approach. Geosciences 2018,8, 201. [CrossRef] 6. Holling, C.S. The components of predation as revealed by a study of small-mammal predation of the european pine sawfly. Can. Entomol. 1959,91, 293–320. [CrossRef] 7. Holling, C.S. Some characteristics of simple types of predation and parasitism. Can. Entomol. 1959,91, 385–398. [CrossRef] 8. Holling, C.S. The functional response of predators to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Can. 1965,97 (Suppl. S45), 5–60. [CrossRef] 9. Hassell, M.P.; Lawton, J.H.; Beddington, J.R. Sigmoid functional responses by invertebrate predators and parasitoids. J. Anim. Ecol. 1977,46, 249–262. [CrossRef] 10. Dalziel, B.D.; Thomann, E.; Medlock, J.; Leenheer, P.D. Global analysis of a predator-prey model with variable predator search rate. J. Math. Biol. 2020,81, 159–183. [CrossRef] [PubMed] 11. Mondal, S.; Samanta, G.P. Impact of fear on a predator-prey system with prey-dependent search rate in deterministic and stochastic environment. Nonlinear Dyn. 2021,104, 2931–2959. [CrossRef] 12. Mondal, S.; Samanta, G. Dynamics of a delayed toxin producing plankton model with variable search rate of zooplankton. Math. Comput. Simul. 2022,196, 166–191. [CrossRef] 13. Hale, J.K. Theory of Functional Differential Equations; Springer: New York, NY, USA, 1977. 14. Perko, L. Differential Equations and Dynamical Systems; Springer: New York, NY, USA, 2001. 15. Murray, J.D. Mathematical Biology; Springer: New York, NY, USA, 1993. 16. Summers, D.; Scott, J.M.W. Systems of first-order chemical reactions. Math. Comput. Model. 1988,10, 901–909. [CrossRef]