Dynamics of Oxygen-Plankton Model with Variable Zooplankton Search Rate in Deterministic and Fluctuating Environments

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.


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 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][4][5], where the search rate of the predator population was constant, i.e., independent of the prey population [6][7][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 4 describes 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.

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][11][12]), we consider the following model (for details see Figure 1): ap 2 z ahp 2 + p + g − µz with initial conditions: 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
A effect 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 m rate of oxygen loss due to the biochemical reaction in a marine ecosystem B maximum phytoplankton per capita growth rate in the high oxygen limit c i , i = 0, 1, 2, 3, 4 half saturation constant of the corresponding processes γ mortality rate due to intraspecific competition among individual phytoplankton a maximally achievable search rate of zooplankton h handling time of zooplankton g half saturation constant σ natural mortality rate of phytoplankton. It is assumed that B > σ η ∈ (0, 1) maximum feeding efficiency µ mortality rate of zooplankton 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 Ac 0 c+c 0 describes the rate of oxygen production per unit of phytoplankton biomass during the daytime by photosynthetic activity; δcp c+c 2 and νcz c+c 3 indicate the respiration of phytoplankton and zooplankton, respectively, and mc is the loss of oxygen due to natural depletion in a marine ecosystem. ahp 2 +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 with an efficiency of ηc 2 c 2 +c 2 4 , which depends on the oxygen concentration (zooplankton die due to insufficient oxygen). , and H (p) p=0 = 2a g > 0. Therefore, H (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 H (p) is presented in Figure 2.

Positivity and Uniform Boundedness
Theorem 1. Solutions of (1) with (2) exist uniquely and are positive for all t ≥ 0.
From the first equation of (1), we have: From the second equation of system (1), we have: From the last equation of system (1), we have: ap 2 (θ) Therefore, c(t) > 0, p(t) > 0 and z(t) > 0 for all t ≥ 0. Hence, the theorem is proved.
Hence, every solution of (1) enters into the region:
Trivial equilibrium point E 0 (0, 0, 0) corresponding to depletion of oxygen and the extinction of plankton; − σ , and c is a positive root of the following equation: Here, Interior (coexistence) equilibrium E( c, p, z), where c, p, and z can be obtained by solving the following system of equations using the software MATHEMATICA: ap 2 ahp 2 + p + g − µ = 0.

Local Stability
Now, we will determine the stability behavior of the biologically feasible equilibrium points of system (1).
The Jacobian matrix J 1 at E 1 ( c, p, 0) is given by: Here, one eigenvalue is λ 1 = η c 2 c 2 +c 2 4 a p 2 ah p 2 + p+g − µ, and the other eigenvalues can be obtained by solving the equation: where Hence, we have the following theorem: The Jacobian matrix J at E( c, p, z) is given by: The characteristic equation corresponding to E( c, p, z) is a 23 a 32 + a 12 a 23 a 31 + a 13 (a 21 a 32 − a 22 a 31 )}. By Routh-Hurwitz's criteria [14], E( c, p, z) has three eigenvalues with negative real parts if C 1 > 0, C 3 > 0, and C 1 C 2 > C 3 . So, the local stability condition of E( c, p, z) is described in the following theorem:

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.

Transcritical Bifurcation
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 J 1 at the bifurcation point must be zero.

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.

Hopf-Bifurcation
The characteristic equation of system (1) at E( c, p, z) is given by where C i (A) for i = 1, 2, 3 were defined earlier.
To determine the Hopf-bifurcation around E( c, p, z) 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 C 1 (A), C 2 (A), and C 3 (A) are continuously differentiable functions of A in a small neighbourhood of A [H] ∈ R such that Equation (5) has: (i) a pair of imaginary eigenvalues λ = p 1 (A) ± ip 2 (A) with p 1 (A) ∈ R, p 2 (A) ∈ R, so that they become purely imaginary at A = A [H] and dp 1 dA

Proof. At A = A [H]
, the roots of the equation: , and λ 3 = −C 1 , where C 1 , C 2 and C 3 are differential functions of A. Furthermore, in the deleted neighborhood of A [H] , the roots (eigenvalues) are λ 1 Now, we will verify the transversality condition: Substituting λ(A) = p 1 (A) + ip 2 (A) into the characteristic Equation (5), we have: Differentiating with regard to A, we have: Comparing the real and imaginary parts, we obtain: and where From (8) and (9), we obtain:ṗ Hence, Theorem 7 is proved using Theorem 6.
Note: Imaginary eigenvalues are connected with any molecular process (e.g., collisions) and the reverse of that process [16].
Remark 2. Similarly, system (1) undergoes Hopf-bifurcations around E( c, p, z) taking any one of the parameters a, g, h, m, η, σ, and µ as a bifurcation parameter.

Numerical Simulations
Here, numerical simulations were performed to verify the analytical outcomes of the oxygen-plankton model (1). We are mainly interested in the existence and stability analysis of the interior equilibrium point E( c, p, z). For this purpose, we fixed most of the parameters as follows: (10) but varied A in a broad range. For the existence of E( c, p, z), we always take A ≥ B (otherwise, it does not exist; see Figure 3), it is also ecologically meaningful. If we choose A = 1.8, and the other parameters are selected from set (10), then, the interior equilibrium E( c, p, z) ≡ E(0.682267, 0.515313, 0.30855) exists uniquely and is locally asymptotically stable (LAS). Figure 4 depicts the stable nature of E(0.682267, 0.515313, 0.30855). If we increase A from 1.8 to 2 keeping the others fixed as in set (10), then E( c, p, z) is destabilized through Hopf-bifurcation. Figure 5 shows the oscillatory nature of system (1) around E(0.710824, 0.473597, 0.352043). If we take a very large value of A (= 10), choosing the other parameters from set (10), then system (1) again enters into a stable interior equilibrium by excluding the existence of a periodic solution. Figure 6 presents the stability nature of E( c, p, z) ≡ E(1.48635, 0.218551, 0.866809) of system (1) when A = 10. In this manner, we have found two thresholds for the parameter A: when 1.8 ≤ A < A [H] 1 = 1.966532 (threshold value) and A > A [H] 2 = 7.258206 (threshold value), a stable interior equilibrium exists; when A [H] 1 = 1.966532 < A < A [H] 2 = 7.258206, the interior equilibrium becomes unstable, and a Hopf-bifurcation occurs, leading to the occurrence of a stable periodic solution (see Figures 7 and 8). Moreover, comparing Figures 6 and 9, we found that for very large values of A, the interior equilibrium of system (1) exists, but it does not exist in the plankton-oxygen model system analyzed by Sekerci and Petrovskii [3]. This is the main difference between the proposed system (1) and the model studied by Sekerci and Petrovskii [3] (it is shown that the coexistence steady state exists unless A is too large or too small).
Again, if we increase µ (the mortality rate of the zooplankton) from 0.1 to 0.5 selecting other parameters from Figure 4, the zooplankton population can not persist in the marine ecosystem. Therefore, the coexistence steady state E( c, p, z) goes to zooplankton free equilibrium E 1 ( c, p, 0). Under the parametric values: {A = 1.8, c 0 = 1, δ = 1, c 2 = 1, m = 0.5, B = 1.8, c 1 = 1.0, γ = 0.7, σ = 0.1, ν = 0.01, µ = 0.5, h = 1.2, a = 3.0, η = 0.7, c 4 = 1, g = 0.3, and c 3 = 1}, we have obtained two planer equilibrium points E   0.958588, 1.11567, 0). Similarly, if we take A = 10 (large), but the other parameters remain the same as in Figure 10a, then we have also obtained two planer equilibria     (c) Bifurcation diagram of z Figure 7. Hopf-bifurcation diagrams of E( c, p, z) of system (1) while A varies in the interval [1.8, 2] and the others remain unchanged as in Figure 6. Here, A = A [H] 1 = 1.966532.    The qualitative nature of different steady states corresponding to bifurcation parameters σ, h, µ and m are depicted in Figures 11-14 respectively (for details see Table 2). Also, the qualitative nature of different steady states corresponding to bifurcation parameters g, η, a and γ are presented in Figures 15-18 respectively (for details see Table 2). Table 2. Nature of the steady states when the parameters are chosen from Figure 6. Here, 'H' stands for the Hopf-bifurcation point, and 'tc' stands for the transcritical bifurcation point.  Figure 14)

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, the phytoplankton growth term through parameter B, and the zooplankton mortality rate µ by turning A, B, and µ into random variables as follows: where γ 1 , γ 2 , and γ 3 are independent Gaussian white noise terms and satisfy the following conditions: Hence, our proposed stochastic system is: 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: (13) but varied A in 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).   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).

Discussion and Conclusions
A Holling type II functional response [6][7][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 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 7 and 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. 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.