A Reaction Diffusion Model to Describe the Toxin Effect on the Fish-Plankton Population

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.


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][2][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][5][6][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 fish zooplankton phytoplankton 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.

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 reactiondiffusion 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.

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: 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  > 0, 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)  0 and  0 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)  0 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)  1 and  2 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)  2 is the value of zooplankton for which the elimination rate by the individual zooplankton becomes  2 /2, (viii)   and   represent, respectively, the growth rate and the external mortality rate of the fish, (ix)  3 is the maximum value that the rate of reduction by the fish can reach, (x)  3 represents the residual loss caused by high rarity of the zooplankton of individuals of the fish population.

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  1 (),  2 (), and  3 () 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  > 0 as follows: where Ω ⊂ R  ( ≥ 2) is the spatial domain in which species occur and zero-flux boundary condition.
So, model (2) can be written by the following system: We make the following assumptions: (H  ) All demographic parameters of system ( 6) are positive constants.
(H  ) The diffusivity coefficients of system ( 6) are independent of .
Taking into account hypotheses (H  ) and (H  ), 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].

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: (10)

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).
The following theorem ensures the boundedness of the solution.
Proof.We will construct a pair of subsolution and oversolution of problem ( 17): which verifies By fixing  1 = 1 >  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  3 ≤ (/)( 2 + ).
From Theorem 3 and hypothesis ( 3 ), this last condition is satisfied.

Stability of Homogeneous Stationary
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 have with We assume the following assumptions: Theorem 6.Let ( 4 ), ( 5 ), and ( 6 ) 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 0 <  0 < ⋅ ⋅ ⋅ <   < ⋅ ⋅ ⋅ .So, we have From ( 27) and ( 23), we have J =     , where   is the matrix defined as is stable if and only if each   () decreases to zero.So,   has three eigenvalues with negative real parts   ,  = 1, 2, 3, such that  3 +  2 +  +  = 0, with For the three eigenvalues to admit a negative real part, we must necessarily have For  = 0, from (30), we have According 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 ( * 1 ,  * 2 ,  * 3 ) under certain conditions.In fact, we have the following hypothesis: Thus, we have the following theorem.
Proof.Let us define the function We have Under hypothesis ( 7 ), the matrix  is positive, and consequently  1 < 0.
On the other hand, we have

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][16][17]25].Here, the system is studied on a spatial domain: Ω = [0, 50] × [0, 50].It is assumed that the fish, zooplankton, and phytoplankton are spread over the whole domain at the beginning.The qualitative results of different pattern formations due to the variation of  0 and  0 are shown.Here, we consider the value of released toxin   = 0.13.These numerical results show that, for every strictly positive initial condition, under assumptions ( 1 )-( 7 ), 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  0 and  0 .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.

Analysis of the Dynamics Behavior with Toxin Effect.
In this section, we consider that  0 = 0.2562 and  0 = 0.71.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.

Conclusion and Future Work
Our interest in this paper is the construction of a reactiondiffusion model for modeling the dynamics of fish population 14 Journal of Mathematics 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.

Figure 1 :
Figure 1: Typical ecological situation presented by food-chain model.

Figure 12 :
Figure 12: Corresponding convergence for error estimation in function of the iterations, where the simulated time is  = 50: (a) for  (b) for , and (c) for .
2 ,  5 ,  6 ,  7 ,  8 ,  9 ,  10 , and  11 are positive  ∞  9 (, , ) is the function that measures the mortality of the competitor from the consumption of the predator, (viii)  10 (, ) is the growth function of the predator, (ix)  11 (, ) is the function that measures the external mortality of the predator population.

Table 1 :
Numerical values for the parameters of system.
0 Half-saturation constant of the  predation on  0.7  2 The value of  (elimination rate by the individual of ) 0.4  3 The residual loss caused by high rarity of  of individuals of