Abstract
This paper is aimed at the mathematical formulation, the analysis, and the numerical simulation of a prey-competitor-predator model by taking into account the toxin produced by the phytoplankton species. The mathematical study of the model leads us to have an idea on the existence of solution, the existence of equilibria, and the stability of the stationary equilibria. These results are obtained through the principle of comparison. Finally, the numerical simulations allowed us to establish a threshold of release of the toxin, above which we talk about the phytoplankton blooms.
1. Introduction
Understanding the functioning of an ecosystem is a major issue for resource and environmental management. However, this goal remains difficult to attain due to the complexity of natural systems, especially in the aquatic environment, where many processes of all types interact with living organisms [1–3]. Plankton is the basis of all aquatic food chains and phytoplankton in particular occupies the first trophic level and the fluctuations in its abundance determine the production of marine environment. Rapid increase and almost equally rapid decrease separated by periods are the common features of plankton populations. In a broad sense, planktonic blooms can be divided into two types: “spring blooms” and “red tides.” Spring blooms occur seasonally for the changes in temperature or nutrient availability which are connected with seasonal changes in thermocline depth, strength, and consequent mixing. Red tides are localized outbreaks and occur due to high water temperature [4–7]. Toxic substances produced by phytoplankton species reduce the growth of zooplankton by decreasing grazing pressure and this is one of the important common phenomena in plankton ecology [2, 6, 8, 9].
Within the broad perspective drawn above, the present paper explores and compares the coupled dynamics of the phytoplankton and the zooplankton in a number of mathematical models. The system phytoplankton-zooplankton has attracted considerable attention from various fields of research [3, 8, 10, 11]. It is an important issue in mathematical ecology. The literature abounds in models focusing on various aspects of the problem. Recently, the attention has been focused on the role of the space in explaining heterogeneity and the distribution of the species and the influence of the spatial structure on their abundance [2, 6, 8]. However, the very question of the interactions between phytoplankton, zooplankton, and fish depending on space is far from being fully elucidated. As part of our work, we will highlight a threshold value of the toxin released by phytoplankton, below which the effect of the toxin influences less the dynamics of the zooplankton-fish system. To our knowledge, this has never been addressed in the literature. The proposed model consists of three interactive components, zooplankton, herbivorous fish, and toxin-producing phytoplankton, which reduce the growth of fish and zooplankton population. The model studied here is of the reaction-diffusion type, describing the dynamics of the phytoplankton-zooplankton-fish system in the sense of the work of Courchamp et al. [12].
The paper is organized as follows. As far as Section 3 is concerned, we will establish mathematical results such as the existence of a solution, stability of equilibria, and persistence, relating to the constructed model in the Section 2. Section 4 will be devoted to numerical experiments to illustrate the mathematical results. Finally, Section 5 is devoted to the conclusion and perspectives.
2. Mathematical Model
In this section, we propose a model to describe the dynamics of the phytoplankton-zooplankton-fish system in the presence of toxin. We begin our modeling by a general model describing the dynamics of the prey-competitor-predator system, based on the equations with ordinary derivatives. And then we transform this model into a model of reaction-diffusion type while remaining in the logic of the works of Courchamp et al. [12] and Bendahmane [13]. It should be noted that the aim is to take into account the effect of the toxin on the fish population through that of zooplankton.
2.1. Original Model Formulation
If we designate by the density of the prey, by the density of the competitor, and by that of the predator, according to [12], the general model is written as follows:where (i), , , , , , , and are positive functions,(ii) is the growth function of the prey,(iii) is the function that measures the mortality due to the competition between the same individuals of the prey,(iv) is the function that measures the mortality of the prey due to the competitor’s consumption,(v) is the competitor’s growth function,(vi) is the function that measures the mortality of the competitor from the release of the substance produced by the prey,(vii) is the function that measures the mortality of the competitor from the consumption of the predator,(viii) is the growth function of the predator,(ix) is the function that measures the external mortality of the predator population.
We continue our modeling by fixing the expressions of the functions intervening in model (1). The dynamics of the system are represented by Figure 1.
According to Figure 1, for any time , the dynamics of the phytoplankton (prey)-zooplankton (competitor)-fish (predator) system is governed by the following ODE system:where (i) and denote, respectively, the phytoplankton growth rate and the phytoplankton external mortality rate,(ii) and denote, respectively, the mortality due to competition between the individuals of the phytoplankton population and the specific predation rate of the zooplankton on the phytoplankton population,(iii) is the half-saturation constant of zooplankton predation on the phytoplankton,(iv) denotes the external mortality rate of zooplankton,(v) denotes the rate of toxin production per phytoplankton species,(vi) and represent, respectively, the maximum value that the consumption rate of the phytoplankton by the zooplankton can reach and the maximum value that the reduction rate by the zooplankton can reach,(vii) is the value of zooplankton for which the elimination rate by the individual zooplankton becomes ,(viii) and represent, respectively, the growth rate and the external mortality rate of the fish,(ix) is the maximum value that the rate of reduction by the fish can reach,(x) represents the residual loss caused by high rarity of the zooplankton of individuals of the fish population.
2.2. A Spatially Structured Model
Considering the relationship between the climate and the diffusion of species and the fact of the existence of diffusion in population, system (1) is developed into a spatial system with diffusion. We expect to explore the effect of climate change on plankton population by studying the spatial dynamics of the diffusion system. We will introduce the concept of spatial structure in the model; that is to say, the population’s densities depend now on the time and space. Diffusion models are a simple and reasonable choice for modeling dispersion of populations on a spatial domain [14, 15]. Indeed, let , , and be, respectively, the diffusivity coefficients of , , and . Based on the work established in [12], the structured model associated with model (1) can be modeled for and as follows:where is the spatial domain in which species occur and zero-flux boundary condition. for , where is the unit normal vector to along and nonnegative initial conditions
So, model (2) can be written by the following system:We make the following assumptions:() All demographic parameters of system (6) are positive constants.() The diffusivity coefficients of system (6) are independent of .
Let us consider , , and in the remaining of our work.
Taking into account hypotheses () and (), model (6) becomes
Remark 1. The model presented here is a model of predator prey system with a self-diffusion but it is also possible to extend these types of system to systems with a cross-diffusion in the same logic of the works of Andreianov et al. [16, 17].
3. Mathematical Results
The mathematical results are based on the works in [3, 16, 18–20].
3.1. Reduction of Model Parameters
To simplify the writing, we will make the following changes of variables:Thus, we obtain, respectively, systems (2) and (7) of the following forms:
3.2. Existence and Boundedness of Solution
Before stating the boundedness of the solution, we give the following theorem concerning the existence and uniqueness of the solution of system (7).
Theorem 2 (see [16, 21–23]). System (7) has a unique local solution under the condition , where depends on nonnegative initial data , , and
The following theorem ensures the boundedness of the solution.
Theorem 3. The set attracts all global solutions of system (7) initiating in the positive orthant, where .
Proof. By noticing that , , and and from [10, 19], we have and Finally, we get that
The following theorem gives the existence of global solution.
Theorem 4 (see [21–23]). System (10) admits a global solution for any and for any regular initial positive functions , , and
Proof. In fact,(i)on one hand, we have , , and because is the subsolution of each equation of system (10);(ii)on the other hand, satisfies the following problem:According to the principle of comparison, we have , where is the solution of the systemIn the same order, satisfies
According to the principle of comparison, we have , where is the solution of the second equation of system (2) with initial condition .
According to [10, 19], we have Let us denote ; then Therefore, we have ,
From the third equation of system (10), we find From the comparison principle, we can deduce that , where is the solution of the third equation of system (2) with initial condition Therefore, Then, And thus,
3.3. Analysis of Stationary Solutions
The positive equilibria are the solutions of system (7), where , , and . Suppose that there exist two positive constants and such that , , and . The equilibrium of system (7) satisfiesWe make the following assumptions:
Theorem 5. Let assumptions , , , and hold; system (17) has at least one positive solution.
Proof. We will construct a pair of subsolution and oversolution of problem (17): which verifiesBy fixing and , according to Theorem 3 and under hypothesis , we have , and This implies that Thus, the right-hand term of the first inequality of (20) is satisfied. Let , and to satisfy the left term of the second inequality of system (20), we need to satisfy the following inequality: Let and According to Theorem 3 and under hypothesis , we have and
Therefore, Thus, the right-hand side of the second inequality of system (20) is satisfied. To satisfy the right term of the third inequality of system (20), we need From Theorem 3 and hypothesis , this last condition is satisfied.
3.4. Stability of Homogeneous Stationary Solutions
Let us consider a solution of system (7) in the following form:where is the positive solution of (17). Let us replace by (22) in system (10) and, by identifying all terms that are linearized in , we havewith We assume the following assumptions:
Theorem 6. Let , , and hold; then the equilibrium is locally stable.
Proof. Concerning the proof of this theorem, we use the results established in [22, 23]. In fact, we consider the following system:where is the eigenfunction of in and the eigenvalue satisfying .
So, we haveFrom (27) and (23), we have , where is the matrix defined as is stable if and only if each decreases to zero. So, has three eigenvalues with negative real parts , , such that , withFor the three eigenvalues to admit a negative real part, we must necessarily haveFor , from (30), we haveAccording to the Routh-Hurwitz criterion [24], sufficient stability conditions for positive stationary solutions can be obtained.
We will establish the stability of the homogeneous equilibrium under certain conditions. In fact, we have the following hypothesis:
Thus, we have the following theorem.
Theorem 7. Suppose that the hypotheses , , , and are satisfied; then is globally asymptotically stable.
Proof. Let us define the functionWe havewhereUnder hypothesis , the matrix is positive, and consequently .
On the other hand, we have
4. Numerical Experiments
In this section, we perform extensive numerical simulations of the spatial model system (7) in two-dimensional space using the forward finite difference method. The set of parameter values used for the numerical simulation are given in Table 1 [11, 15–17, 25]. Here, the system is studied on a spatial domain: . It is assumed that the fish, zooplankton, and phytoplankton are spread over the whole domain at the beginning.
4.1. Pattern Formation
Here, we will illustrate the mathematical predictions, by numerical simulations, concerning the behavior of the dynamics under hypotheses The qualitative results of different pattern formations due to the variation of and are shown. Here, we consider the value of released toxin These numerical results show that, for every strictly positive initial condition, under assumptions , the nonhomogeneous equilibrium is always globally asymptotically stable. Figures 2–7 show the spatial structures formation for the threes species described in (7). These numerical results confirm the mathematical results about the existence and the stability of positive equilibrium according to the values of and . In this case, we will speak of a subsistence phenomenon of the fish and zooplankton population.
Remark 8. From a biological point of view, these results (Figures 2–7) show that there is coexistence between the three populations despite the release of the toxin into the aquatic environment. This means that, despite the harmful effects of the toxin released by phytoplankton, fish and zooplankton populations persist.
4.2. Analysis of the Dynamics Behavior with Toxin Effect
In this section, we consider that and . Figures 8–11 show that, after a transitional phase, an equilibrium can be established with a strong distribution of the three populations, which is not far from the results obtained previously. As a biological interpretation, we can say that if the toxin is released below this value, the impact is not significant on the fish and zooplankton populations (Figures 8–10). A less dense distribution of fish and zooplankton than before is observed. This explains the considerable decrease of these species due to the increase in the number of toxic phytoplankton. On the other hand, we observe in Figure 11 a sudden change concerning the spatial distribution of phytoplankton (strong and bad spatial distribution).
Remark 9. Figure 12 shows the estimate of the error by using the standard norm ‖‖2 on the iterations [16, 17, 26]. These figures show a decrease in error estimation over the course of the iterations. This confirms the convergence of the numerical scheme used.
(a)
(b)
(c)
5. Conclusion and Future Work
Our interest in this paper is the construction of a reaction-diffusion model for modeling the dynamics of fish population under a diet on plankton (zooplankton and phytoplankton), taking into account the effect of the toxin. The model formulation is derived from an ODE system by considering an isotropic distribution as in [12]. It should be noted that we consider diffusion independently of the spatial variable in the construction of the reaction-diffusion model. The mathematical results allowed us to establish conditions of existence of equilibrium which depend on the demographic parameters. We have also given some results about the stability of the stationary equilibria and we have established the conditions on the inexistence of the equilibrium with strictly positive components. We continued our study through numerical experiments in order to confirm our mathematical results. Also, our numerical results have yielded interesting results on the effect of the toxin on the dynamics. This is why we are led to conclude that the release of the toxin under certain conditions in the aquatic space contributes to the regulation of the system. The phytoplankton bloom was observed during our simulations and is in perfect agreement with the biological observations.
Despite important results obtained in this paper, in order to further our study, we can consider for future work subdividing the phytoplankton population into toxic phytoplankton and nontoxic phytoplankton to extend our results to a system of cross-diffusion type.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Authors’ Contributions
All the authors read and approved the final manuscript.