Threshold analysis of an algae-zooplankton model incorporating general interaction rates and nonlinear independent stochastic components

1.   Introduction and overview

Stochasticity is inherent in ecological systems, manifesting in myriad ways that shape the destiny of their residents. Stochasticity, characterized by its inherent randomness and unpredictability, is a prime example of this variability, playing a pivotal role in ecological dynamics ^[1]. The occurrence of stochastic events, such as unpredictable fluctuations in environmental conditions or chance encounters with predators, introduces a layer of uncertainty that can profoundly impact the survival and evolution of species ^[2]. For example, consider a population that, under deterministic conditions, would persist indefinitely due to stable environmental factors and abundant resources. Despite apparent resilience, the introduction of stochastic events can dramatically alter the population's trajectory. Sudden environmental disturbances, disease outbreaks, or shifts in predator-prey dynamics can swiftly tip the balance, leading to population declines or even extinction ^[3].

Generally, the effects of stochasticity are not uniform across species or ecosystems. Some species may possess traits or behaviors that buffer them against the impacts of random events, while others may be more susceptible to stochastic shocks ^[4]. Additionally, the scale and frequency of stochastic events vary spatially and temporally, further complicating their influence on ecological systems ^[5]. Understanding the interplay between deterministic and stochastic factors is crucial for comprehending the resilience and vulnerability of ecosystems. While deterministic processes provide a framework for understanding general patterns and trends, stochastic events inject the element of chance that can shape the outcomes in unpredictable ways. Therefore, a holistic approach that considers both deterministic and stochastic forces is essential for effectively managing and conserving biodiversity in a dynamic and uncertain world ^[6].

The addition of random perturbations in the algae-zooplankton model is motivated by the inherent variability and unpredictability present in natural ecosystems ^[7]. Environmental fluctuations, such as changes in temperature, light, nutrient availability, and water quality, significantly impact the growth and interactions of both algae and zooplankton populations. Incorporating random perturbations helps capture these environmental variabilities. Additionally, immigration and emigration processes, which are stochastic in nature, can influence population sizes through currents and other physical processes. The interactions between algae and zooplankton are complex and influenced by random predation events and changes in prey availability. Stochastic elements account for the noise and uncertainty in real-world data due to measurement errors and sampling inconsistencies. Understanding how populations respond to random disturbances is crucial for assessing their stability and resilience. Therefore, incorporating random perturbations into the model provides a realistic and comprehensive framework for understanding and predicting the dynamics of these populations in fluctuating environments. In ecological modeling, various methods are employed to introduce variability into ecosystems, each serving different purposes and mimicking distinct aspects of natural dynamics. Among the most common techniques are white noise and Lévy jumps, which play crucial roles in capturing different types of variability within ecological systems.

White noise is a stochastic process characterized by random fluctuations with equal intensity across all frequencies. In ecological contexts, white noise is often used to simulate natural disturbances driven by environmental changes ^[8]. These disturbances can include fluctuations in temperature, precipitation, or resource availability, as well as disturbances caused by human activities such as logging, pollution, or habitat fragmentation ^[9]. By incorporating white noise into ecological models, researchers can replicate the unpredictability and variability inherent in natural ecosystems, providing insights into how species and communities respond to environmental fluctuations over time ^[10]. Building on this approach, certain papers have explored the dynamics between algae and zooplankton under white noise perturbations. For instance:

$\bullet$ In ^[11], the author introduced a novel approach to modeling the dynamics of algae-zooplankton systems, incorporating stochasticity through the inclusion of random environmental perturbations represented as white noise. His work yielded insights into the conditions for stochastic extinction, persistence, and the establishment of an ergodic stationary distribution within the system.

$\bullet$ In ^[12], the authors presented a nutrient-algae model enriched with the incorporation of white noise, aiming to capture the stochastic nature of ecological processes. Their study not only provided insights into the necessary conditions for extinction, persistence, and ergodicity within the model but also shed light on the dynamic interplay between nutrient availability and algae abundance in fluctuating environments.

$\bullet$ In ^[13], the authors investigated the dynamics of a stochastic system involving toxin-producing algae and fish, incorporating harvesting practices. Their study yielded valuable insights into the system's behavior, offering sufficient conditions for both extinction and persistence.

$\bullet$ In ^[14], the authors explored the dynamics of a stochastic model concerning algae-zooplankton interactions, taking into account variations in algae cell size and zooplankton body size. The authors derived criteria for both extinction and ergodicity within the model.

On the other hand, pure Lévy jumps represent a type of stochastic process characterized by infrequent but intense changes in a system's state ^[15]. Unlike white noise, which generates random fluctuations of relatively constant magnitude, pure Lévy jumps simulate abrupt and extreme events that can have significant impacts on ecological dynamics. These events may include rare but severe disturbances such as wildfires, hurricanes, or disease outbreaks, as well as sudden shifts in population dynamics due to predator-prey interactions or competitive pressures ^[16]. By incorporating pure Lévy jumps into ecological models, researchers can explore the effects of these rare but influential events on the resilience and stability of ecosystems, shedding light on how ecosystems respond to and recover from extreme perturbations. Expanding upon this methodology, specific studies have delved into the interplay between algae and zooplankton dynamics, incorporating the concept of Lévy jumps. For example:

$\bullet$ In ^[17], the authors introduced a stochastic model of nutrient-plankton food chains, incorporating Lévy jumps into their framework. They established necessary conditions for both the persistence and extinction of plankton populations within the model.

$\bullet$ In ^[18], the authors introduced an ecological system subjected to high-order white noise and quadratic Lévy jumps. They delineated the conditions governing the extinction and ergodic behavior of the model, offering valuable insights into the long-term dynamics and stability of the system under stochastic perturbations.

By integrating both white noise and pure jumps into ecological modeling frameworks, researchers can capture a wide range of variability and complexity in natural systems. This comprehensive approach enables scientists to investigate how different types of disturbances and environmental changes shape the dynamics of ecological communities, helping to inform conservation and management strategies in an ever-changing world. Additionally, understanding the interplay between different sources of variability can provide valuable insights into the adaptive capacity of species and ecosystems, aiding efforts to predict and mitigate the impacts of global environmental change.

In this paper, we present a general ecological model that incorporates two distinct sources of jumps, each tailored to capture specific aspects of stochasticity within natural systems. Firstly, we integrate second-order white noises, meticulously designed to account for the subtle yet impactful fluctuations that characterize small-scale variations in environmental conditions. These fluctuations, although seemingly minor, can have profound effects on ecosystem dynamics, influencing everything from population dynamics to community structure. Secondly, we introduce second-order pure jumps, strategically employed to capture the abrupt and intense changes that typify heavy perturbations within ecological systems. These jumps represent the sudden shifts and disturbances that can arise from various factors, such as extreme weather events, sudden changes in resource availability, or outbreaks of disease.

Crucially, we introduce a parameter ranging between $0$ and $1$, offering the flexibility to personalize the perturbation characteristics and providing insight into both subtle and heavy variations. This parameterization allows us to systematically explore the effects of different levels of perturbation intensity on ecosystem dynamics, offering a nuanced understanding of how ecological systems respond to environmental variability across a spectrum of conditions.

A notable feature of our study is the establishment of a global threshold delineating the boundary between extinction and ergodicity within the model. This threshold plays a pivotal role in distinguishing between two critical scenarios: The sustained persistence of ecological communities (ergodicity) and their eventual decline and disappearance (extinction). By identifying and elucidating this threshold, we offer valuable insights into the underlying mechanisms governing the long-term stability and resilience of algae-zooplankton communities.

The remainder of this paper is organized as follows: Section 2 introduces the algae-zooplankton model and its stochastic formulation. Following this, we present some preliminary concepts and lemmas. In Section 3, we unveil the main theorems that establish the global threshold between extinction and ergodicity within the stochastic model. Section 4 provides numerical illustrations to examine the impact of noise on the dynamics of algae and zooplankton. Finally, we conclude with a discussion of our findings and a consideration of our limitations.

2.   Stochastic system and preliminaries

Drawing upon the model outlined in ^[14], we proceed to introduce a refined algae-zooplankton system that builds upon and extends previous findings in the field. Our system encapsulates the intricate interplay between algae ${y}_{0}$ and zooplankton ${y}_{1}$ dynamics, aiming to elucidate the underlying mechanisms governing their interactions within aquatic ecosystems. The model under perturbation is structured as follows:

where $r$ is the size of the algae cells, $s_1$, $s_2$, and $s_3$ are some experimental coefficients, $L$ is the significant ecological carrying capacity, $\beta_0$ is the peak consumption rate, $\gamma$ is the saturated coefficient, $\zeta_1$ and $\zeta_2$ are the usage rate coefficients, $u$ is the metric for zooplankton morphology, $\alpha$ is the transition rate of ${y}_{1}$, $\beta_1$ is the toxic emission, $d_1$ is the mortality rate of ${y}_{1}$, and $\varpi_{00}$ and $\varpi_{10}$ represent the intensity of noise for the independent standard Brownian motions ${S}_{0}(t)$ and ${S}_{1}(t)$, respectively.

In the complex ecosystems of aquatic environments, the predation dynamics between zooplankton and phytoplankton are closely tied to the size disparities among these tiny organisms. A key principle recognized by ecologists is that larger-bodied zooplankton typically prey on smaller-bodied phytoplankton ^[19,20,21]. However, as the ecosystem becomes more populated with multiple species of both phytoplankton and zooplankton, each with varying body sizes, the predation dynamics become increasingly intricate. In these multifaceted scenarios, interactions among planktonic organisms are often complicated by mutual interference, resulting in complex webs of predation and resource use. Consequently, the incidence rates that describe these interactions must take on specific forms to accurately reflect the nuanced relationships at play. Understanding these specific forms of incidence rates is crucial for elucidating the detailed mechanisms governing predator-prey interactions and species coexistence within aquatic ecosystems.

In the model described in ^[14], the authors focused on the nonlinear interaction function $\frac{\beta_0{y}_{0}(t){y}_{1}(t)}{\gamma {y}_{0}(t)+{y}_{1}(t)}$. However, real-world scenarios and certain modeling requirements may necessitate a more versatile approach. In order to accommodate a broader range of situations, this paper will explore a generalized functional interaction function $f\big({y}_{0}(t), {y}_1(t)\big)$. By doing so, we aim to encompass a wider spectrum of epidemiological dynamics and provide a more flexible framework for modeling ecological competition in various contexts. The generalized form of the system (2.1) is presented as follows:

where $f\in\mathbb{C}^2({\mathbb R}_+\times{\mathbb R}_+, {\mathbb R}_+)$ is an uniformly continuous function such that

for a given positive constant $\vartheta$. This generalized function, as described under the aforementioned conditions, includes several specific cases outlined in Table 1 (see, for example, ^{[22,23,24,25,26]}).

Having expanded the incidence rate to a more generalized form, we are now poised to introduce our novel model, enriched with second-order perturbations ^[26,27]. Consider the following system:

Here, ${S}_{k}(t)\; ({k} = 0, 1)$ denotes two mutually independent Standard Brownian motions (SBMs) of strengths $\varpi_{0i} > 0\; ({i} = 0, 1)$ and $\varpi_{1i} > 0\; ({i} = 0, 1)$, respectively. These SBMs are essentially defined on a filtered probability triple (stochastic basis) $\left(\Omega, \mathcal{F}_\Omega, (\mathcal{F}_{\{\Omega, t\}})_{t\geqslant 0}, {\mathbb P}\right)$ equipped with a filtration satisfying the usual criteria ^[4]. ${y}_{k}(t^-)\; ({k} = 0, 1)$ is the left limits of ${y}_{k}(t)\; ({k} = 0, 1)$. $\bar{D}_{k}\; ({k} = 0, 1)$ are two independent stochastic Poisson processes related to the measures $m_{k}\; ({k} = 0, 1)$ defined on a measurable set ${\mathbb R}^2\setminus\{0\}$ as follows:

where $\mathbb{1}$ denotes the indicator function, $\alpha_s\in(0, 2)$, and $R_k$ is the Rosiński measure defined on ${\mathbb R}^2\setminus\{0\}$ such that $R_k(0) = 0$, $({k} = 0, 1)$, with ${ }\int_{{\mathbb R}^2\setminus\{0\}} \left(||x||^2\wedge||x||^{\alpha_s}\right)R_k({{\rm{d}}} x) < \infty, \; \; \alpha_s \in(0, 2)$ (for more details, please see ^[28]). $D_{k}\; ({k} = 0, 1)$ are two independent compensated Poisson random measures such that $D_{k}({\rm{d}}t, {\rm{d}}\xi) = \bar{D}_{k}({\rm{d}}t, {\rm{d}}\xi)-m_{k}({\rm{d}}\xi) {\rm{d}}t$. ${\nu}_{k}: {\mathbb R}^2\setminus\{0\}\to\mathbb{R}$ are the jumps size functions which are postulated to be positive and continuous on ${\mathbb R}^2\setminus\{0\}$, for all $i, k = 0, 1$. It is assumed that ${ }\int_{{\mathbb R}^2\setminus\{0\}}\nu^2_{ik }(\xi)m_{i}({{\rm{d}}} \xi)$ is finite. The variable $p$ is a constant between $0$ and $1$. This is included in order to have the possibility of personalizing the influence of the noise part (from a modeling perspective). If $p = 1$, system (2.3) simplifies to a model featuring solely second-order white noise, effectively simulating continuous and Gaussian variations commonly observed in ecological systems. Conversely, when $p = 0$, system (2.3) transforms into a model characterized by unexpected jumps with heavy tails, resulting in abrupt and drastic changes in species dynamics.

Based on the theoretical framework outlined in references ^[15,29], it becomes evident that the system described by Eq (2.3) possesses a unique positive global solution. This observation not only confirms the ecological validity of the model but also establishes its mathematical robustness and well-posedness.

As outlined in the introduction, the central focus of this paper revolves around delineating the global threshold that delineates the balance between the sustenance and extinction of zooplankton populations. To address the intricacies posed by the new stochastic system (2.3), we propose a new approach based on the ergodicity of an auxiliary equation closely resembling the equation of ${y}_{0}$. Based on the stochastic comparison theorem ^[30], we present the following equation:

Based on Lemma 2.2 in ^[31], it is evident that Eq (2.5) maintains its well-posed nature. Additionally, $\tilde{y}_0(t)$ emerges as a stochastic process endowed with the following properties:

(1) ${\mathbb P}\big(\tilde{y}_0(t)-{y}_{0}(t)\geq 0\big) = 1$.

(2) If

then there exists a unique ergodic stationary distribution $\tilde{\eta}_0$ for Eq (2.5); and

For analytical coherence, it is imperative to maintain the assumption of positivity for $\tilde{d}_0$ throughout this study. Additionally, we rigorously define the following quantities, crucial in delineating the global threshold within the framework of system (2.3):

3.   Global threshold of system (2.3)

As elucidated earlier, the crux of ecosystem examination lies in accurately demarcating the threshold that distinguishes between population persistence and decline. Therefore, the paramount aim of the forthcoming theorems is to tackle this pivotal inquiry head-on, while simultaneously enriching our comprehension of algae-zooplankton system dynamics with additional insights and refinements.

Theorem 3.1. Suppose that $\theta_0\leq\theta_1$. Under this condition, for any initial positive data $(y_0(0), y_1(0))$, the stochastic system described by Eq (2.3) asymptotically tends towards a state of exponential decline and

Proof. Utilizing Itô's formula, we readily obtain

In accordance with the stochastic comparison theorem ^[30], we deduce that

By integrating the above inequality from $0$ to $t$, we acquire

We let $M_0(t) = \int^t_0\int_{{\mathbb R}^2\setminus\{0\}}\log\big(1+\nu_{10}(\xi)\big){D}_1({\rm{d}}s, {\rm{d}}\xi)$. Then, we get

Referring to the strong law of large numbers, we derive ${ }\mathbb{P}\left(\lim\limits_{t\to\infty}\frac{1}{t}M_0(t) = 0\right) = 1$. Now, we select $0 < \varrho_1 < 1$ and consider the following quantities:

To continue, we employ the exponential inequality designed for martingales, we obtain

Thanks to the Borel-Cantelli lemma, we establish the presence of $\varrho_2(\omega)$ for all $\omega\in \Omega$, ensuring the validity of the following inequality:

verifies for any $\varrho_2$ satisfying $\varrho_2 \geq \varrho_2(\omega)$ and ${\mathbb P}\big(t\in ]\varrho_2-1, \varrho_2[\big) = 1$. Consequently,

We compute the upper limit on both sides of Eq (3.2), and then we get

Therefore,

By considering the arbitrary nature of $0 < \varrho_1 < 1$, one can deduce that as $\varrho_1$ approaches 0,

Put differently, for any arbitrarily small positive quantity $\tilde{\theta}$ satisfying $\tilde{\theta} < \min\{0.5, 0.5(\theta_1-\theta_0)\}$, we can infer the existence of $\varrho_3 = \varrho_3(\omega)$ associated with $\Omega_{\varrho_3}$ such that $\mathbb{P}(\Omega_{\varrho_3}) \geq 1 - \varrho_3$, where $\log y_1(t) \leq -0.5(\theta_1 - \theta_0)t$ for all $t \geq \varrho_3$ and $\omega \in \Omega_{\varrho_3}$. As a result, $y_1(t) \leq \exp(-0.5(\theta_1 - \theta_0)t)$ and

This, coupled with the positivity of the solution $(y_0, y_1)$, suggests

This indicates that the exponential decrease in the population of $y_1$ implies its almost certain complete eradication. With this, the proof is concluded. □

Corollary 3.1. (case of second-order Brownian motions) If $p = 1$, and

then, the extinction of the population ${y}_1$ occurs almost surely.

Corollary 3.2. (case of second-order jumps) If $p = 0$, and

then, the extinction of the population ${y}_1$ occurs almost surely.

Theorem 3.2. Suppose that $\theta_0 > \theta_1$. Under this condition, for any initial positive data $(y_0(0), y_1(0))$, the stochastic system described by Eq (2.3) asymptotically tends towards a stationary state and

where $\eta(\cdot, \cdot)$ is the unique ergodic stationary distribution associated with system (2.3).

Proof. Utilizing the Itô formula for the function $F_0(t) = \log\left(\frac{\tilde{y}_0(t)}{{y}_{0}(t)}\right)$ yields:

Based on the properties of the function $f$, it follows that:

By virtue of the positive property of the solution, we obtain:

Once more, employing the operator $\mathcal{L}$ on the second equation of (2.3) yields:

Consequently,

Consistent with the positivity of the solution, we find

Now, let us examine the following function:

By applying the operator $\mathcal{L}$ to $F_1$, we obtain

By rearranging certain terms, we achieve

Once more, let us examine this new function:

Here, $c_{\star}$ denotes a positive constant that fulfills:

Upon applying the operator $\mathcal{L}$ to $F_2$, we acquire

Then

Let $c_{\star\star}\in(0, 1)$, then we define a new function as follows:

By applying the operator $\mathcal{L}$ to $F_3$, we get

We amalgamate functions $F_2$ and $F_3$ in the following manner:

Here, $h_\star > 0$ satisfies $2+\mathcal{K}_\star-h_\star (\overbrace{\theta_0-\theta_1}^{ > 0})\leq 0$, where $\mathcal{K}_\star = \max\big\{\mathcal{A}_\star, 1\big\}$, and

To guarantee the positivity of our Lyapunov function, we introduce this new function:

where, $\underline{F}_4$ represents the lower bound of function $F_4$. When we apply the operator $\mathcal{L}$ to $F_5$, we obtain

where

Utilizing the methodologies outlined in the proof of Theorem 3.2 in ^[31], along with the properties of the function $f$, we can promptly verify $F_6(t)\leq -1$, $\forall({y}_{0}, {y}_{1}) \in{\mathbb R}^{2}_{+}\setminus{H}_x = {H}^c_x$, $x > 0$, where

Alternatively, it's evident that there exists a positive constant $\kappa$ such that $F_6(t)\leq \kappa$ for all $({y}_{0}, {y}_{1})\in {H}_x$. Then, from (3.5), we derive

By virtue of the ergodicity property of the process $\tilde{y}$, we obtain

So,

Consequently,

By employing the lemma of mutually constrained possibilities ^[32], we infer the existence of a unique ergodic stationary distribution $\eta(\cdot, \cdot)$. Hence, incorporating the ergodic property into consideration, we can firmly conclude that

This implies that the system asymptotically converges to the persistent state over time. Thus, the proof is concluded. □

Corollary 3.3. (case of second-order Brownian motions) If $p = 1$, and

then, the permanence of the populations ${y}_1$ and ${y}_1$ occurs almost surely.

Corollary 3.4. (case of second-order jumps) If $p = 0$, and

then, the permanence of the populations ${y}_1$ and ${y}_1$ occurs almost surely.

Remark 3.1. From Theorem 3.1, we infer that the zooplankton population tends to extinction rapidly. This rapid decline suggests that under certain conditions, the zooplankton cannot sustain its population over time, likely due to insufficient resources, high predation pressure, or adverse environmental factors. The extinction of zooplankton can disrupt the balance of the aquatic ecosystem, affecting nutrient cycling and the food web structure. In contrast, Theorem 3.2 implies the persistence and stationarity of both algae and zooplankton populations. This means that under a different set of conditions, both populations are able to maintain stable levels over time, without the risk of extinction. The persistence of these populations indicates a balanced interaction between the species, where the growth rates of algae and zooplankton, along with other ecological factors, allow for long-term coexistence.

4.   Numerical application and discussion

In this section, we carry out a series of numerical examples to validate the findings presented in our study, utilizing the data outlined in ^[14]. We take $r = 0.3$, $u = 0.25$, $s_1 = 0.02$, $s_2 = 0.02$, $s_3 = 0.08$, $L = (2.5;3.5)$, $\zeta_1 = 2$, $\zeta_2 = 0.5$, $\alpha = 0.4$, $\beta_1 = 0.2$, and $d_1 = 0.6$. About the dual incidence functions $g_1$ and $g_2$, we consider the following general nonlinear incidence:

where $\beta_0 = 4.5$, $k_0 = 0.16$, $k_0 = 0.11$ and $k_2 = 0.012$. To simulate the stochastic component, we present a method for simulating jumps and tempered stable processes ^[28]. Firstly, we examine the sequences $\{S_{1, j}\}_{j\geqslant 1}$, $\{S_{2, j}\}_{j\geqslant 1}$, and $\{S_{3, j}\}_{j\geqslant 1}$, representing independent and identically distributed (i.i.d.) random variables taking values in the real numbers, following the distribution (2.4). Additionally, we define $\{S_{2, j}\}_{j\geqslant 1}$ and $\{S_{3, j}\}_{j\geqslant 1}$ as i.i.d. sequences of uniform random variables within the intervals $(0, 1)$ and $(0, T)$, respectively. Furthermore, consider $\{S_{4, j}\}_{j\geqslant 1}$ and $\{S_{5, j}\}_{j\geqslant 1}$ as i.i.d. sequences of random variables following an exponential distribution with a rate coefficient of $1$. It is assumed that all the aforementioned random variables are mutually independent. Then, we let ${ }\{S_{6, j}\} = \sum\limits_{k = 1}^{j} \{S_{5, k}\}$. Noticeably, $\{S_{6, j}\}$ can be regarded as a Poisson point process on the interval ${\mathbb R}_+$ with random intensity measure. In reference to the theory outlined in ^[28], if $\alpha_s\in(0, 2)$, then

converges almost surely and uniformly for t within the interval $[0, T]$ to a Lévy process, where ${ } ||\rho|| = \int_{{\mathbb R}^2\setminus\{0\}}|x|^{\alpha_s} R_k({{\rm{d}}} x)$. In practice, we can devise a method for generating a jump process with specified parameters at discrete time instances $t_i$. Here, ${t_i}_{i\in[0, I]}$ denotes a partition of the interval $[0, T]$, characterized by uniformly sized sub-intervals and a mesh size $\Delta t = T/I$, where $I\in \mathbb{N}$. Subsequently, we employ the following algorithm:

(1) Choose a specific time duration, denoted as T, and partition the time interval $[0, T]$ into $I$ equal-sized segments.

(2) Specify and fix a numerical parameter, denoted as $N$.

(3) Generate numerical replications or simulations of independent quantities ${S_{i, j}}$, where $i$ ranges from $1$ to $6$, with a sample size of $N$.

(4) Calculate the value of $m_t$.

Let us now direct our attention to the following system:

By incorporating the ergodic property of Eq (2.5), we define the parameters $\theta_0$ and $\theta_1$ as follows:

4.1.   Case 1: Almost sure extinction of $y_1$

Firstly, we consider $L = 2.5$, $p = 0.5$ and we take the noise parameters as follows: $\varpi_{00} = 0.02$, $\varpi_{01} = 0.01$, $\varpi_{10} = 0.021$, $\varpi_{11} = 0.012$, ${\nu}_{00} = 0.01$, ${\nu}_{01} = 0.009$, ${\nu}_{10} = 0.02$ and ${\nu}_{11} = 0.007$. In this case, we get $\theta_1 = 0.0526$ and $\theta_2 = 0.07492$. According to Theorem 3.1, it is highly probable that the population $y_1$ will become extinct. This outcome is visually exemplified in Figure 1, which provides a numerical illustration of the aforementioned result.

Secondly, we take $p = 1$ and we establish specific values for the white noise parameters $\varpi_{00} = 0.02$, $\varpi_{01} = 0.01$, $\varpi_{10} = 0.021$ and $\varpi_{11} = 0.012$. In this case, we get $\theta_1 = 0.0526$ and $\theta_2 = 0.1894$. As a result, we have verified that the prerequisite outlined in Corollary 3.1 has been met and the extinction of $y_1$ occurs almost surely. This result is visually depicted in Figure 2.

Now, we set $p = 0$. The jump intensities are defined as follows: ${\nu}_{00} = 0.01$, ${\nu}_{01} = 0.009$, ${\nu}_{10} = 0.02$ and ${\nu}_{11} = 0.007$. Here, we get $\theta_1 = 0.0526$ and $\theta_2 = 0.118$. As a result, we have verified that the prerequisite outlined in Corollary 3.2 has been met. This result is visually depicted in Figure 3.

Upon scrutinizing the aforementioned figures, a discernible pattern emerges: The extinction time is contingent upon the nature of the noise. Notably, a significant presence of noise hastens the extinction process considerably. This observation underscores the critical influence of noise intensity on the dynamics of extinction within the system.

4.2.   Case 2: Permanence and ergodicity of $y_0$ and $y_1$

We consider $L = 3.5$, $p = 0.5$ and we take the noise parameters as follows: $\varpi_{00} = 0.001$, $\varpi_{01} = 0.001$, $\varpi_{10} = 0.004$, $\varpi_{11} = 0.0012$, ${\nu}_{00} = 0.001$, ${\nu}_{01} = 0.001$, ${\nu}_{10} = 0.01$ and ${\nu}_{11} = 0.003$. In this case, we get $\theta_1 = 2.1982$ and $\theta_2 = 0.9862$. As per Theorem 3.2, there is a strong likelihood of the population $y_0$ and $y_1$ persisting, alongside the existence of a stationary distribution. This assertion is vividly illustrated in Figure 4.

We take $p = 1$ and we establish specific values for the white noise parameters $\varpi_{00} = 0.001$, $\varpi_{01} = 0.001$, $\varpi_{10} = 0.004$, $\varpi_{11} = 0.0012$. In this case, we get $\theta_1 = 2.1982$ and $\theta_2 = 1.3412$. Consequently, we have confirmed the fulfillment of the condition outlined in Corollary 3.3, ensuring the almost certain persistence of $y_0$ and $y_1$. This outcome is graphically illustrated in Figure 5, providing a visual representation of the verified result.

Now, we set $p = 0$. The jump intensities are defined as follows: ${\nu}_{00} = 0.001$, ${\nu}_{01} = 0.001$, ${\nu}_{10} = 0.01$ and ${\nu}_{11} = 0.003$. Here, we get $\theta_1 = 2.1982$ and $\theta_2 = 1.2981$. Henceforth, we have substantiated the fulfillment of the prerequisite delineated in Corollary 3.4. This validation is visually elucidated through the representation provided in Figure 6.

Expanding upon the insights drawn from the aforementioned figure, it becomes apparent that the characteristics of the noise exert a discernible influence on the magnitude of variation. Specifically, we observe a pronounced effect: as the intensity of the noise increases, there is a notable augmentation in the amplitude of variation. This correlation underscores the pivotal role played by noise levels in amplifying the variability within the system dynamics, thereby emphasizing the significance of noise intensity in shaping the observed outcomes.

Plankton size is a critical variable because it influences various biological and ecological processes. For example, smaller plankton may have faster reproduction rates but are more susceptible to predation, while larger plankton might be more resilient but have slower population growth. By studying different plankton sizes, we can observe how size-specific traits affect the overall population dynamics and stationary distribution. In instances of persistence, we explore the influence of plankton size and the parameter $p$ on the morphology of the stationary distribution, as illustrated in Figure 7. This detailed analysis illuminates the interplay between plankton size and system parameters, providing insights into how these factors shape the stationary distribution under varying conditions. By examining different plankton sizes and adjusting $p$, we can observe the impact on population dynamics and structure, ultimately gaining a deeper understanding of the conditions that favor particular plankton sizes and how populations adapt over time. This comprehensive examination is crucial for predicting plankton responses to environmental changes and for effective marine ecosystem management.

5.   Conclusions

The stochastic dynamics inherent in marine ecosystems are of paramount importance, shaping the stability and resilience of these intricate environments. In this paper, we presented a novel nonlinear algae-zooplankton model, that offers a comprehensive framework for understanding the dynamics influenced by general interaction rates and second-order noise. Specifically, we distinguish between white noise and jump components, allowing for a detailed examination of their respective impacts on population dynamics.

Central to our investigation is the delineation of the global threshold that governs the delicate balance between species extinction and persistence. Through rigorous analysis, in Theorems 3.1 and 3.2, we probed the intricate interplay of various ecological factors, shedding light on the fundamental mechanisms underlying ecosystem dynamics.

In the numerical simulation section, we investigate the repercussions of noise on the long-term behavior of algae and zooplankton populations. By systematically varying parameters and plankton sizes, we unraveled the nuanced effects of noise intensity on population dynamics and stationary distributions. Through visual representations and quantitative analyses, we provided insights into how different types of noise shape the ecological trajectories of these crucial marine species, offering valuable implications for ecosystem management and conservation efforts.

The analysis described in this paper can be extended to describe the dynamics of a mycoloop model in aquatic food webs. This model involves interactions between phytoplankton, chytrids, and zooplankton ^[33]. By incorporating additional hypotheses and exploring more complex interactions, the model can provide deeper insights into the ecological processes governing mycoloop dynamics. We aim to pursue this idea in our future work, further enhancing our understanding of aquatic ecosystem dynamics.

Author contributions

Yassine Sabbar: Conceptualization, Writing Original Draft, Software, Formal Analysis; Aeshah A. Raezah: Conceptualization, Software, Validation. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through Large Research Project under grant number RGP 2/174/45.

Conflict of interest

The authors declare no conflict of interest.