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

Authors: Wendkouni Ouedraogo; Hamidou Ouedraogo; Boureima Sangaré
Journal: Journal of Mathematics
 (2018)

---

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


---

## Body

## 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 Section3 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 byP the density of the prey, by Z the density of the competitor, and by F that of the predator, according to [12], the general model is written as follows:(1)dPdT=ϕ2P-g5P,ZP-g6P,ZZdZdT=g7P,ZZ-g8P,Z,FZ-g9P,Z,FFdFdT=g10Z,FF-g11Z,FF,where (i)
ϕ2, g5, g6, g7, g8, g9, g10, and g11 are positive C∞ functions,(ii)
ϕ2(P) is the growth function of the prey,(iii)
g5(P,Z) is the function that measures the mortality due to the competition between the same individuals of the prey,(iv)
g6(P,Z) is the function that measures the mortality of the prey due to the competitor’s consumption,(v)
g7(P,Z) is the competitor’s growth function,(vi)
g8(P,Z,F) is the function that measures the mortality of the competitor from the release of the substance produced by the prey,(vii)
g9(P,Z,F) is the function that measures the mortality of the competitor from the consumption of the predator,(viii)
g10(Z,F) is the growth function of the predator,(ix)
g11(Z,F) 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.Figure 1
Typical ecological situation presented by food-chain model.According to Figure1, for any time T>0, the dynamics of the phytoplankton (prey)-zooplankton (competitor)-fish (predator) system is governed by the following ODE system:(2)dPdT=rpP-μϵP-b0P2-α0ZP+γ0P,P0=P0≥0,dZdT=-rzZ-θpZ+α1PP+γ0Z-α2FZ+γ2Z,Z0=Z0≥0,dFdT=rfF-μfF-α3F2Z+γ3,F0=F0≥0,where (i)
rp and μϵ denote, respectively, the phytoplankton growth rate and the phytoplankton external mortality rate,(ii)
b0 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)
rz denotes the external mortality rate of zooplankton,(v)
θp 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)
rf and μf 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.


### 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 δ1(x), δ2(x), and δ3(x) be, respectively, the diffusivity coefficients of P, Z, and F. Based on the work established in [12], the structured model associated with model (1) can be modeled for x∈Ω and T>0 as follows:(3)∂TP-divδ1x∇P=ϕ2P-g5P,ZP-g6P,ZZ,∂TZ-divδ2x∇Z=g7P,ZZ-g8P,Z,FZ-g9P,Z,FF,∂TF-divδ3x∇F=g10Z,FF-g11Z,FF,where Ω⊂Rn(n≥2) is the spatial domain in which species occur and zero-flux boundary condition. (4)δix∇Qx,t·νx=0,i=1,2,3,x∈∂Ω,t>0,for Q=P,Z,F, where ν is the unit normal vector to Ω along ∂Ω and nonnegative initial conditions (5)Qx,0=Q0x>0,Q=P,Z,F,x∈Ω.So, model (2) can be written by the following system:(6)∂TP=rpx-μϵx-b0xP-α0xZP+γ0xP+divδ1x∇P,∂TZ=-rzx-θpx+α1xPP+γ0x-α2xFZ+γ2xZ+divδ2x∇Z,∂TF=rfx-μfx-α3xFZ+γ3xF+divδ3x∇F.We make the following assumptions:(Ha
) All demographic parameters of system (6) are positive constants.(Hb
) The diffusivity coefficients of system (6) are independent of x.Let us considerδ1(x)=δ1, δ2(x)=δ2, and δ3(x)=δ3 in the remaining of our work.Taking into account hypotheses (Ha) and (Hb), model (6) becomes(7)∂TP=rp-μϵ-b0P-α0ZP+γ0P+δ1ΔP,Px,0=P0x>0,∂TZ=-rz-θp+α1PP+γ0-α2FZ+γ2Z+δ2ΔZ,Zx,0=Z0x>0,∂TF=rf-μf-α3FZ+γ3F+δ3ΔF,Fx,0=F0x>0.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].


## 2.1. Original Model Formulation

If we designate byP the density of the prey, by Z the density of the competitor, and by F that of the predator, according to [12], the general model is written as follows:(1)dPdT=ϕ2P-g5P,ZP-g6P,ZZdZdT=g7P,ZZ-g8P,Z,FZ-g9P,Z,FFdFdT=g10Z,FF-g11Z,FF,where (i)
ϕ2, g5, g6, g7, g8, g9, g10, and g11 are positive C∞ functions,(ii)
ϕ2(P) is the growth function of the prey,(iii)
g5(P,Z) is the function that measures the mortality due to the competition between the same individuals of the prey,(iv)
g6(P,Z) is the function that measures the mortality of the prey due to the competitor’s consumption,(v)
g7(P,Z) is the competitor’s growth function,(vi)
g8(P,Z,F) is the function that measures the mortality of the competitor from the release of the substance produced by the prey,(vii)
g9(P,Z,F) is the function that measures the mortality of the competitor from the consumption of the predator,(viii)
g10(Z,F) is the growth function of the predator,(ix)
g11(Z,F) 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.Figure 1
Typical ecological situation presented by food-chain model.According to Figure1, for any time T>0, the dynamics of the phytoplankton (prey)-zooplankton (competitor)-fish (predator) system is governed by the following ODE system:(2)dPdT=rpP-μϵP-b0P2-α0ZP+γ0P,P0=P0≥0,dZdT=-rzZ-θpZ+α1PP+γ0Z-α2FZ+γ2Z,Z0=Z0≥0,dFdT=rfF-μfF-α3F2Z+γ3,F0=F0≥0,where (i)
rp and μϵ denote, respectively, the phytoplankton growth rate and the phytoplankton external mortality rate,(ii)
b0 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)
rz denotes the external mortality rate of zooplankton,(v)
θp 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)
rf and μf 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.


## 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 δ1(x), δ2(x), and δ3(x) be, respectively, the diffusivity coefficients of P, Z, and F. Based on the work established in [12], the structured model associated with model (1) can be modeled for x∈Ω and T>0 as follows:(3)∂TP-divδ1x∇P=ϕ2P-g5P,ZP-g6P,ZZ,∂TZ-divδ2x∇Z=g7P,ZZ-g8P,Z,FZ-g9P,Z,FF,∂TF-divδ3x∇F=g10Z,FF-g11Z,FF,where Ω⊂Rn(n≥2) is the spatial domain in which species occur and zero-flux boundary condition. (4)δix∇Qx,t·νx=0,i=1,2,3,x∈∂Ω,t>0,for Q=P,Z,F, where ν is the unit normal vector to Ω along ∂Ω and nonnegative initial conditions (5)Qx,0=Q0x>0,Q=P,Z,F,x∈Ω.So, model (2) can be written by the following system:(6)∂TP=rpx-μϵx-b0xP-α0xZP+γ0xP+divδ1x∇P,∂TZ=-rzx-θpx+α1xPP+γ0x-α2xFZ+γ2xZ+divδ2x∇Z,∂TF=rfx-μfx-α3xFZ+γ3xF+divδ3x∇F.We make the following assumptions:(Ha
) All demographic parameters of system (6) are positive constants.(Hb
) The diffusivity coefficients of system (6) are independent of x.Let us considerδ1(x)=δ1, δ2(x)=δ2, and δ3(x)=δ3 in the remaining of our work.Taking into account hypotheses (Ha) and (Hb), model (6) becomes(7)∂TP=rp-μϵ-b0P-α0ZP+γ0P+δ1ΔP,Px,0=P0x>0,∂TZ=-rz-θp+α1PP+γ0-α2FZ+γ2Z+δ2ΔZ,Zx,0=Z0x>0,∂TF=rf-μf-α3FZ+γ3F+δ3ΔF,Fx,0=F0x>0.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:(8)a0=rp-μϵ,a1=rz+θp,c3=rf-μf,U1t=b0a0PT,U2t=b0α0a02ZT,t=a0T,p=c3a0,q=α3α2,y=Ya0,a=b0γ0a0,b=a1a0,x=Xa0,r=b0α0γ3a02,d=b0α0γ0a02,ϵ3=a03b0α0α1δ3,ϵ1=a0b0δ1,ϵ2=a02b0α0δ2,z=Za0,c=α1a0,U3t=b0α0α1a03FT.Thus, we obtain, respectively, systems (2) and (7) of the following forms:(9)dE1dt=1-E1+E2E1+aE1,dE2dt=-b+cE1E1+a-E3E2+dE2,dE3dt=p-qE3E2+rE3,(10)∂tU1=1-U1+U2U1+aU1+ϵ1ΔU1,∂tU2=-b+cU1U1+a-U3U2+dU2+ϵ2ΔU2,∂tU3=p-qU3U2+rU3+ϵ3ΔU3.


### 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 (U1(·,t),U2(·,t),U3(·,t)) under the condition 0≤t<Tmax, where Tmax depends on nonnegative initial data U01(x), U02(x), and U03(x).The following theorem ensures the boundedness of the solution.Theorem 3.
The setD=0,1×0,c+c/4b×0,1/α+1/4αb+M/αb attracts all global solutions of system (7) initiating in the positive orthant, where α=4b/(b+p)2(4bc+c+4br).Proof.
By noticing thatU1≤E1, U2≤E2, and U3≤E3 and from [10, 19], we have limt→+∞¯E1(t)≤1 and limt→+∞¯E1(t)+1/cE2(t)≤1+1/4b. Finally, we get that (11)limt→+∞¯E1t+1cE2t+αE3t≤1+14b+Mαb.The following theorem gives the existence of global solution.Theorem 4 (see [21–23]).
System (10) admits a global solution (U1(·,t),U2(·,t),U3(·,t)) for any t>0 and for any regular initial positive functions U01(x)≤1, U02(x), and U03(x).Proof.
In fact,(i)
on one hand, we haveU1(x,t)≥0, U2(x,t)≥0, and U3(x,t)≥0 because 0 is the subsolution of each equation of system (10);(ii)
on the other hand,U1 satisfies the following problem:(12)∂tU1≤U11-U1+ϵ1ΔU1onΩ×0,+∞,∂νU1=0,on∂Ω×0,+∞,U1x,0=U01x≤U01≡maxΩ¯⁡U01x,x∈Ω.According to the principle of comparison, we have U1(x,t)≤U(t)≤1, where U(t)=U01et/1-U01+et is the solution of the system(13)∂U∂t=U1-U,U0=U01≤1.In the same order, U2 satisfies ∂U2/∂t=cU1U2/U1+a-bU2-U2U3/U2+d+ϵ2ΔU2.
According to the principle of comparison, we have1/c∂U2/∂t≤1/cdE2/dt, where E2 is the solution of the second equation of system (2) with initial condition E2(0)=maxΩ¯U02(x).
According to [10, 19], we have 1/c∂U2/∂t≤1/cdE2/dt+dE1/dt. Let us denote σ=1/cE2+E1; then 1/c∂U2/∂t≤dσ/dt≤1/4+b-bσ. Therefore, we have U2(x,t)≤c+c/4b, 0≤t<T.
From the third equation of system (10), we find α∂U3/∂t=αpU3-αqU32/U2+r+αϵ3ΔU3. From the comparison principle, we can deduce that α∂U3/∂t≤αdE3/dt, where E3 is the solution of the third equation of system (2) with initial condition E3(0)=maxΩ¯U3(x). Therefore, (14)α∂U3∂t≤αdE3dt+1cdE2dt+dE1dt≤E11-E1-bcE2-1cE2E3E2+d+αp-qE3E2+rE3≤14+b-αbU3+αbE3+αp-qE3E2+rE3.Then, (15)α∂U3∂t+αbU3≤14+b+M,M=14q.And thus, (16)U3x,t≤1α+14αb+Mαb,0≤t<T.


### 3.3. Analysis of Stationary Solutions

The positive equilibria are the solutions(U1(x),U2(x),andU3(x)) of system (7), where U1(x)≥0, U2(x)≥0, and U3(x)≥0. Suppose that there exist two positive constants m1<1 and m2 such that U1(x)≥m1>0, U2(x)≥m2>0, and U3(x)>0. The equilibrium of system (7) satisfies(17)-ϵ1ΔU1=U11-U1-U1U2U1+a=F∗U1,U2,U3,-ϵ2ΔU2=cU1U2U1+a-bU2-U2U3U2+d=G∗U1,U2,U3,-ϵ3ΔU3=pU3-qU32U2+r=H∗U1,U2,U3,∂U1∂ν=∂U2∂ν=∂U3∂ν=0on∂Ω.We make the following assumptions: (18)H1c+c4b≤m1+am1m21+a+m11-m1,H2b≥c1+a,1α+c4bα+Mα≤cm1m1+a-b-U3¯m2+d,H3pqm2+r≤1α+c4bα+Mα.Theorem 5.
Let assumptions(Ha), (H1), (H2), and (H3) hold; system (17) has at least one positive solution.Proof.
We will construct a pair of subsolution and oversolution of problem (17): (19)U¯1x,U¯2x,U¯3x>U1_x,U2_x,U3_x,which verifies(20)-ϵ1ΔU¯1-U¯11-U¯1-U¯1U¯2U¯1+a≥0≥-ϵ1ΔU_1-U_11-U_1-U_1U_2U_1+a,-ϵ2ΔU¯2-cU¯1U¯2U¯1+a+bU¯2-U¯2U¯3U¯2+d≥0≥-ϵ2ΔU_1-cU_1U_2U_1+a+bU_2+U_2U_3U_2+d,-ϵ3ΔU¯3-pU¯3-qU¯32U¯2+r≥0≥-ϵ3ΔU_1-pU_3-qU_32U_2+r,∂U¯1∂ν≥0≥∂U_1∂ν,∂U¯2∂ν≥0≥∂U_2∂ν,∂U¯3∂ν≥0≥∂U_3∂ν,on∂Ω.By fixing U¯1=1>U_1=m1 and U_2=m2, according to Theorem 3 and under hypothesis (H1), we have U¯2≤c+c/4b, and c+c/4b≤m1+a/m1m2/(1+a)+m1(1-m1). This implies that U¯2≤(m1+a)/m1m2/(1+a)+m1(1-m1). Thus, the right-hand term of the first inequality of (20) is satisfied. Let U¯1=1, and to satisfy the left term of the second inequality of system (20), we need U_3 to satisfy the following inequality: c/(1+a)-b-U_3/U¯2+d≤0. Let U_1=m1 and U_2=m2. According to Theorem 3 and under hypothesis (H2), we have 0≤U¯3≤1/α+c/4bα+M/α and 1/α+c/4bα+M/α≤cm1/(m1+a)-b-U¯3/m2+d.
Therefore,(21)cm1m1+a-b-U¯3m2+d≥0.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 U_3≤p/q(m2+r). From Theorem 3 and hypothesis (H3), this last condition is satisfied.


### 3.4. Stability of Homogeneous Stationary Solutions

Let us consider a solutionU(x,t) of system (7) in the following form:(22)Ux,t=U1x,t,U2x,t,U3x,t=Usx+Jx,t=U1sx,U2sx,U2sx+J1x,t,J2x,t,J3x,t,where Us(x) is the positive solution of (17). Let us replace U(x,t) by (22) in system (10) and, by identifying all terms that are linearized in J, we have(23)∂J∂t=ϵ1000ϵ2000ϵ3ΔJ+FU1∗FU2∗FU3∗GU1∗GU2∗GU3∗HU1∗HU2∗HU3∗J,with (24)FU1∗=1-2U1-aU2U1+a2,FU2∗=-U1U1+a,FU3∗=0,GU1∗=acU2U1+a2,GU2∗=cU1U1+a-b-dU3U2+a2,GU3∗=-U2U2+d,HU1∗=0,HU2∗=qU32U2+d2,HU3∗=p-2qU3U2+d.We assume the following assumptions: (25)H4FU1∗+GU2+HU3∗<0,H5FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗-GU3∗HU2∗>0,H6FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗>0.Theorem 6.
Let(H4), (H5), and (H6) hold; then the equilibrium Us 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:(26)-Δχj=βjχjonΩ,∂χj∂ν∂Ω=0,where χj is the eigenfunction of -Δ in Ω and the eigenvalue βj satisfying 0<β0<⋯<βj<⋯.
So, we have(27)J=∑j=0∞Jjtχjx,y,Jjt∈Rn.From (27) and (23), we have J˙j=CjJj, where Cj is the matrix defined as (28)Cj=FU1∗FU2∗FU3∗GU1∗GU2∗GU3∗HU1∗HU2∗HU3∗-λjϵ1000ϵ2000ϵ3.Us is stable if and only if each Jj(t) decreases to zero. So, Cj has three eigenvalues with negative real parts ξk, k=1,2,3, such that ξ3+aξ2+bξ+c=0, with(29)a=βjϵ1+ϵ2+ϵ3-FU1∗+GU2∗+HU3∗,b=βj2ϵ1ϵ2+ϵ1ϵ3+ϵ2ϵ3-βjϵ1GU2∗+ϵ1HU3∗+ϵ2FU1∗+ϵ2HU3∗+ϵ3FU1∗+ϵ3GU2∗+FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗,c=λj3ϵ1ϵ2ϵ3-βj2ϵ1ϵ2HU3∗+ϵ1ϵ3GU2∗+ϵ2ϵ3FU1∗+βjϵ1GU2HU3∗+ϵ2FU1∗HU3∗+ϵ3FU1∗GU2∗-ϵ1GU3∗HU2∗-ϵ2FU3∗HU1∗-ϵ3FU2∗GU1∗+FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗.For the three eigenvalues to admit a negative real part, we must necessarily have(30)FU1∗+GU2∗+HU3∗-βjϵ1+ϵ2+ϵ3<0,FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗-βjϵ1GU2∗+ϵ1HU3∗+ϵ2FU1∗+ϵ2HU3∗+ϵ3FU1∗+ϵ3GU2∗+βj2ϵ1ϵ2+ϵ1ϵ3+ϵ2ϵ3>0,FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗+βjϵ1GU2∗HU3∗+ϵ2FU1∗HU3∗+ϵ3FU1∗GU2∗-ϵ1GU3∗HU2∗-ϵ2FU3∗HU1∗-ϵ3FU2∗GU1∗-βj2ϵ1ϵ2HU3∗+ϵ1ϵ3GU2∗+ϵ2ϵ3FU1∗+βj3ϵ1ϵ2ϵ3>0.For j=0, from (30), we have(31)FU1∗+GU2∗+HU3∗<0,FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗>0,FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗>0.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(U1∗,U2∗,U3∗) under certain conditions. In fact, we have the following hypothesis: (32)H7b<c,1-a<U1∗<abc-b,U2∗<dc-b.Thus, we have the following theorem.Theorem 7.
Suppose that the hypotheses(Ha), (H1), (H2), and (H7) are satisfied; then (U1∗,U2∗,U3∗) is globally asymptotically stable.Proof.
Let us define the function(33)ΥU1,U2,U3=∫Ω∫U1∗U1n-U1∗ηdη+∫U2∗U2n-U2∗η+dηdη+∫U3∗U3n-U3∗ηdηdx.We have(34)dΥdt=T1+T2,T1=∫ΩU1-U1∗U1U11-U1-U1U2U1+a+U2-U2∗U2+dU2cU1U2U1+a-bU2-U2U3U2+d+U3-U3∗U3pU3-qU32U2+rdx=∫ΩU1-U1∗U1+a1-U1U1+a-U2+U2∗-1-U1∗U1∗+a+U2-U2∗U2+dU2cU1U2+d-bU1+aU2+d-U3U1+a-cU1∗U2∗+d-bU1∗+aU2∗+d-U3∗U1∗+a+U2-U2∗qU3∗U2∗+r-qU3U2+rdx=-∫ΩU1-U1∗,U2-U2∗-U3∗SU1-U1∗,U2-U2∗,U3-U3∗Tdx,where(35)
S
=
U
1
+
U
1
∗
+
a
-
1
U
1
+
a
U
3
∗
-
c
-
b
U
2
∗
+
d
U
1
+
a
0
1
U
1
+
a
a
b
-
c
-
b
U
1
∗
U
1
+
a
-
p
U
2
+
r
0
1
q
U
2
+
r
.Under hypothesis (H7), the matrix S is positive, and consequently T1<0.On the other hand, we have(36)T2=∫Ωϵ1U1-U1∗U1ΔU1+ϵ2U2-U2∗U2+dU2ΔU2+ϵ3U3-U3∗U3ΔU3dx=-ϵ1∫Ω∇U12U1∗U1dx-ϵ2∫Ω∇U22U2+dU3∗U2dx-ϵ3∫Ω∇U32U3∗U3dx<0.


## 3.1. Reduction of Model Parameters

To simplify the writing, we will make the following changes of variables:(8)a0=rp-μϵ,a1=rz+θp,c3=rf-μf,U1t=b0a0PT,U2t=b0α0a02ZT,t=a0T,p=c3a0,q=α3α2,y=Ya0,a=b0γ0a0,b=a1a0,x=Xa0,r=b0α0γ3a02,d=b0α0γ0a02,ϵ3=a03b0α0α1δ3,ϵ1=a0b0δ1,ϵ2=a02b0α0δ2,z=Za0,c=α1a0,U3t=b0α0α1a03FT.Thus, we obtain, respectively, systems (2) and (7) of the following forms:(9)dE1dt=1-E1+E2E1+aE1,dE2dt=-b+cE1E1+a-E3E2+dE2,dE3dt=p-qE3E2+rE3,(10)∂tU1=1-U1+U2U1+aU1+ϵ1ΔU1,∂tU2=-b+cU1U1+a-U3U2+dU2+ϵ2ΔU2,∂tU3=p-qU3U2+rU3+ϵ3ΔU3.


## 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 (U1(·,t),U2(·,t),U3(·,t)) under the condition 0≤t<Tmax, where Tmax depends on nonnegative initial data U01(x), U02(x), and U03(x).The following theorem ensures the boundedness of the solution.Theorem 3.
The setD=0,1×0,c+c/4b×0,1/α+1/4αb+M/αb attracts all global solutions of system (7) initiating in the positive orthant, where α=4b/(b+p)2(4bc+c+4br).Proof.
By noticing thatU1≤E1, U2≤E2, and U3≤E3 and from [10, 19], we have limt→+∞¯E1(t)≤1 and limt→+∞¯E1(t)+1/cE2(t)≤1+1/4b. Finally, we get that (11)limt→+∞¯E1t+1cE2t+αE3t≤1+14b+Mαb.The following theorem gives the existence of global solution.Theorem 4 (see [21–23]).
System (10) admits a global solution (U1(·,t),U2(·,t),U3(·,t)) for any t>0 and for any regular initial positive functions U01(x)≤1, U02(x), and U03(x).Proof.
In fact,(i)
on one hand, we haveU1(x,t)≥0, U2(x,t)≥0, and U3(x,t)≥0 because 0 is the subsolution of each equation of system (10);(ii)
on the other hand,U1 satisfies the following problem:(12)∂tU1≤U11-U1+ϵ1ΔU1onΩ×0,+∞,∂νU1=0,on∂Ω×0,+∞,U1x,0=U01x≤U01≡maxΩ¯⁡U01x,x∈Ω.According to the principle of comparison, we have U1(x,t)≤U(t)≤1, where U(t)=U01et/1-U01+et is the solution of the system(13)∂U∂t=U1-U,U0=U01≤1.In the same order, U2 satisfies ∂U2/∂t=cU1U2/U1+a-bU2-U2U3/U2+d+ϵ2ΔU2.
According to the principle of comparison, we have1/c∂U2/∂t≤1/cdE2/dt, where E2 is the solution of the second equation of system (2) with initial condition E2(0)=maxΩ¯U02(x).
According to [10, 19], we have 1/c∂U2/∂t≤1/cdE2/dt+dE1/dt. Let us denote σ=1/cE2+E1; then 1/c∂U2/∂t≤dσ/dt≤1/4+b-bσ. Therefore, we have U2(x,t)≤c+c/4b, 0≤t<T.
From the third equation of system (10), we find α∂U3/∂t=αpU3-αqU32/U2+r+αϵ3ΔU3. From the comparison principle, we can deduce that α∂U3/∂t≤αdE3/dt, where E3 is the solution of the third equation of system (2) with initial condition E3(0)=maxΩ¯U3(x). Therefore, (14)α∂U3∂t≤αdE3dt+1cdE2dt+dE1dt≤E11-E1-bcE2-1cE2E3E2+d+αp-qE3E2+rE3≤14+b-αbU3+αbE3+αp-qE3E2+rE3.Then, (15)α∂U3∂t+αbU3≤14+b+M,M=14q.And thus, (16)U3x,t≤1α+14αb+Mαb,0≤t<T.


## 3.3. Analysis of Stationary Solutions

The positive equilibria are the solutions(U1(x),U2(x),andU3(x)) of system (7), where U1(x)≥0, U2(x)≥0, and U3(x)≥0. Suppose that there exist two positive constants m1<1 and m2 such that U1(x)≥m1>0, U2(x)≥m2>0, and U3(x)>0. The equilibrium of system (7) satisfies(17)-ϵ1ΔU1=U11-U1-U1U2U1+a=F∗U1,U2,U3,-ϵ2ΔU2=cU1U2U1+a-bU2-U2U3U2+d=G∗U1,U2,U3,-ϵ3ΔU3=pU3-qU32U2+r=H∗U1,U2,U3,∂U1∂ν=∂U2∂ν=∂U3∂ν=0on∂Ω.We make the following assumptions: (18)H1c+c4b≤m1+am1m21+a+m11-m1,H2b≥c1+a,1α+c4bα+Mα≤cm1m1+a-b-U3¯m2+d,H3pqm2+r≤1α+c4bα+Mα.Theorem 5.
Let assumptions(Ha), (H1), (H2), and (H3) hold; system (17) has at least one positive solution.Proof.
We will construct a pair of subsolution and oversolution of problem (17): (19)U¯1x,U¯2x,U¯3x>U1_x,U2_x,U3_x,which verifies(20)-ϵ1ΔU¯1-U¯11-U¯1-U¯1U¯2U¯1+a≥0≥-ϵ1ΔU_1-U_11-U_1-U_1U_2U_1+a,-ϵ2ΔU¯2-cU¯1U¯2U¯1+a+bU¯2-U¯2U¯3U¯2+d≥0≥-ϵ2ΔU_1-cU_1U_2U_1+a+bU_2+U_2U_3U_2+d,-ϵ3ΔU¯3-pU¯3-qU¯32U¯2+r≥0≥-ϵ3ΔU_1-pU_3-qU_32U_2+r,∂U¯1∂ν≥0≥∂U_1∂ν,∂U¯2∂ν≥0≥∂U_2∂ν,∂U¯3∂ν≥0≥∂U_3∂ν,on∂Ω.By fixing U¯1=1>U_1=m1 and U_2=m2, according to Theorem 3 and under hypothesis (H1), we have U¯2≤c+c/4b, and c+c/4b≤m1+a/m1m2/(1+a)+m1(1-m1). This implies that U¯2≤(m1+a)/m1m2/(1+a)+m1(1-m1). Thus, the right-hand term of the first inequality of (20) is satisfied. Let U¯1=1, and to satisfy the left term of the second inequality of system (20), we need U_3 to satisfy the following inequality: c/(1+a)-b-U_3/U¯2+d≤0. Let U_1=m1 and U_2=m2. According to Theorem 3 and under hypothesis (H2), we have 0≤U¯3≤1/α+c/4bα+M/α and 1/α+c/4bα+M/α≤cm1/(m1+a)-b-U¯3/m2+d.
Therefore,(21)cm1m1+a-b-U¯3m2+d≥0.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 U_3≤p/q(m2+r). From Theorem 3 and hypothesis (H3), this last condition is satisfied.


## 3.4. Stability of Homogeneous Stationary Solutions

Let us consider a solutionU(x,t) of system (7) in the following form:(22)Ux,t=U1x,t,U2x,t,U3x,t=Usx+Jx,t=U1sx,U2sx,U2sx+J1x,t,J2x,t,J3x,t,where Us(x) is the positive solution of (17). Let us replace U(x,t) by (22) in system (10) and, by identifying all terms that are linearized in J, we have(23)∂J∂t=ϵ1000ϵ2000ϵ3ΔJ+FU1∗FU2∗FU3∗GU1∗GU2∗GU3∗HU1∗HU2∗HU3∗J,with (24)FU1∗=1-2U1-aU2U1+a2,FU2∗=-U1U1+a,FU3∗=0,GU1∗=acU2U1+a2,GU2∗=cU1U1+a-b-dU3U2+a2,GU3∗=-U2U2+d,HU1∗=0,HU2∗=qU32U2+d2,HU3∗=p-2qU3U2+d.We assume the following assumptions: (25)H4FU1∗+GU2+HU3∗<0,H5FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗-GU3∗HU2∗>0,H6FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗>0.Theorem 6.
Let(H4), (H5), and (H6) hold; then the equilibrium Us 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:(26)-Δχj=βjχjonΩ,∂χj∂ν∂Ω=0,where χj is the eigenfunction of -Δ in Ω and the eigenvalue βj satisfying 0<β0<⋯<βj<⋯.
So, we have(27)J=∑j=0∞Jjtχjx,y,Jjt∈Rn.From (27) and (23), we have J˙j=CjJj, where Cj is the matrix defined as (28)Cj=FU1∗FU2∗FU3∗GU1∗GU2∗GU3∗HU1∗HU2∗HU3∗-λjϵ1000ϵ2000ϵ3.Us is stable if and only if each Jj(t) decreases to zero. So, Cj has three eigenvalues with negative real parts ξk, k=1,2,3, such that ξ3+aξ2+bξ+c=0, with(29)a=βjϵ1+ϵ2+ϵ3-FU1∗+GU2∗+HU3∗,b=βj2ϵ1ϵ2+ϵ1ϵ3+ϵ2ϵ3-βjϵ1GU2∗+ϵ1HU3∗+ϵ2FU1∗+ϵ2HU3∗+ϵ3FU1∗+ϵ3GU2∗+FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗,c=λj3ϵ1ϵ2ϵ3-βj2ϵ1ϵ2HU3∗+ϵ1ϵ3GU2∗+ϵ2ϵ3FU1∗+βjϵ1GU2HU3∗+ϵ2FU1∗HU3∗+ϵ3FU1∗GU2∗-ϵ1GU3∗HU2∗-ϵ2FU3∗HU1∗-ϵ3FU2∗GU1∗+FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗.For the three eigenvalues to admit a negative real part, we must necessarily have(30)FU1∗+GU2∗+HU3∗-βjϵ1+ϵ2+ϵ3<0,FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗-βjϵ1GU2∗+ϵ1HU3∗+ϵ2FU1∗+ϵ2HU3∗+ϵ3FU1∗+ϵ3GU2∗+βj2ϵ1ϵ2+ϵ1ϵ3+ϵ2ϵ3>0,FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗+βjϵ1GU2∗HU3∗+ϵ2FU1∗HU3∗+ϵ3FU1∗GU2∗-ϵ1GU3∗HU2∗-ϵ2FU3∗HU1∗-ϵ3FU2∗GU1∗-βj2ϵ1ϵ2HU3∗+ϵ1ϵ3GU2∗+ϵ2ϵ3FU1∗+βj3ϵ1ϵ2ϵ3>0.For j=0, from (30), we have(31)FU1∗+GU2∗+HU3∗<0,FU1∗GU2∗+FU1∗HU3∗+GU2∗HU3∗-FU2∗GU1∗-FU3∗HU1∗+GU3∗HU2∗>0,FU1∗GU3∗HU2∗+FU2∗GU1∗HU3∗+FU3∗GU2∗HU1∗-FU1∗GU2∗HU3∗-FU2∗GU3∗HU1∗-FU3∗GU1∗HU2∗>0.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(U1∗,U2∗,U3∗) under certain conditions. In fact, we have the following hypothesis: (32)H7b<c,1-a<U1∗<abc-b,U2∗<dc-b.Thus, we have the following theorem.Theorem 7.
Suppose that the hypotheses(Ha), (H1), (H2), and (H7) are satisfied; then (U1∗,U2∗,U3∗) is globally asymptotically stable.Proof.
Let us define the function(33)ΥU1,U2,U3=∫Ω∫U1∗U1n-U1∗ηdη+∫U2∗U2n-U2∗η+dηdη+∫U3∗U3n-U3∗ηdηdx.We have(34)dΥdt=T1+T2,T1=∫ΩU1-U1∗U1U11-U1-U1U2U1+a+U2-U2∗U2+dU2cU1U2U1+a-bU2-U2U3U2+d+U3-U3∗U3pU3-qU32U2+rdx=∫ΩU1-U1∗U1+a1-U1U1+a-U2+U2∗-1-U1∗U1∗+a+U2-U2∗U2+dU2cU1U2+d-bU1+aU2+d-U3U1+a-cU1∗U2∗+d-bU1∗+aU2∗+d-U3∗U1∗+a+U2-U2∗qU3∗U2∗+r-qU3U2+rdx=-∫ΩU1-U1∗,U2-U2∗-U3∗SU1-U1∗,U2-U2∗,U3-U3∗Tdx,where(35)
S
=
U
1
+
U
1
∗
+
a
-
1
U
1
+
a
U
3
∗
-
c
-
b
U
2
∗
+
d
U
1
+
a
0
1
U
1
+
a
a
b
-
c
-
b
U
1
∗
U
1
+
a
-
p
U
2
+
r
0
1
q
U
2
+
r
.Under hypothesis (H7), the matrix S is positive, and consequently T1<0.On the other hand, we have(36)T2=∫Ωϵ1U1-U1∗U1ΔU1+ϵ2U2-U2∗U2+dU2ΔU2+ϵ3U3-U3∗U3ΔU3dx=-ϵ1∫Ω∇U12U1∗U1dx-ϵ2∫Ω∇U22U2+dU3∗U2dx-ϵ3∫Ω∇U32U3∗U3dx<0.


## 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: Ω=[0,50]×[0,50]. It is assumed that the fish, zooplankton, and phytoplankton are spread over the whole domain at the beginning.Table 1
Numerical values for the parameters of system.
Parameter Description Value γ 0 Half-saturation constant of theZ predation on P 0.7 γ 2 The value ofZ (elimination rate by the individual of Z) 0.4 γ 3 The residual loss caused by high rarity ofZ of individuals of F 0.9 r p Growth rate of the phytoplankton population 1.6 r z External mortality rate of the zooplankton population 0.4 δ 1 Diffusivity coefficient of the phytoplankton population 0.03 δ 2 Diffusivity coefficient of the zooplankton population 0.03 μ ϵ External mortality rate of phytoplankton population 0.9 δ 3 Diffusion coefficient of the fish population 0.05 α 1 Consumption rate of the phytoplankton population byZ 0.75 α 2 Consumption rate of the zooplankton population byF 0.4 α 3 Mortality rate of the fish population due to the competition 0.122 μ f External mortality rate of the fish population 0.02 r f Growth rate of the fish population 0.4


### 4.1. Pattern Formation

Here, we will illustrate the mathematical predictions, by numerical simulations, concerning the behavior of the dynamics under hypotheses(H1)–(H7). The qualitative results of different pattern formations due to the variation of b0 and α0 are shown. Here, we consider the value of released toxin θp=0.13. These numerical results show that, for every strictly positive initial condition, under assumptions (H1)–(H7), 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 b0 and α0. In this case, we will speak of a subsistence phenomenon of the fish and zooplankton population.Figure 2
Spatial distribution of the system withb0=0.25 and α0=0.6.Figure 3
Spatial distribution of the system withb0=0.38 and α0=0.54.Figure 4
Spatial distribution of the system withb0=0.4032 and α0=0.342.Figure 5
Spatial distribution of the system withb0=0.2562 and α0=0.71.Figure 6
Spatial distribution of the system withb0=0.61 and α0=0.69.Figure 7
Spatial distribution of the system withb0=0.69 and α0=0.595.Remark 8.
From a biological point of view, these results (Figures2–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 thatb0=0.2562 and α0=0.71. 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).Figure 8
Dynamics behavior forθp=0.18.Figure 9
Dynamics behavior forθp=0.29.Figure 10
Dynamics behavior forθp=0.46.Figure 11
Dynamics behavior forθp=0.935.Remark 9.
Figure12 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.Figure 12
Corresponding convergence for error estimation in function of the iterations, where the simulated time ist=50: (a) for P (b) for Z, and (c) for F.
(a) (b) (c)


## 4.1. Pattern Formation

Here, we will illustrate the mathematical predictions, by numerical simulations, concerning the behavior of the dynamics under hypotheses(H1)–(H7). The qualitative results of different pattern formations due to the variation of b0 and α0 are shown. Here, we consider the value of released toxin θp=0.13. These numerical results show that, for every strictly positive initial condition, under assumptions (H1)–(H7), 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 b0 and α0. In this case, we will speak of a subsistence phenomenon of the fish and zooplankton population.Figure 2
Spatial distribution of the system withb0=0.25 and α0=0.6.Figure 3
Spatial distribution of the system withb0=0.38 and α0=0.54.Figure 4
Spatial distribution of the system withb0=0.4032 and α0=0.342.Figure 5
Spatial distribution of the system withb0=0.2562 and α0=0.71.Figure 6
Spatial distribution of the system withb0=0.61 and α0=0.69.Figure 7
Spatial distribution of the system withb0=0.69 and α0=0.595.Remark 8.
From a biological point of view, these results (Figures2–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 thatb0=0.2562 and α0=0.71. 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).Figure 8
Dynamics behavior forθp=0.18.Figure 9
Dynamics behavior forθp=0.29.Figure 10
Dynamics behavior forθp=0.46.Figure 11
Dynamics behavior forθp=0.935.Remark 9.
Figure12 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.Figure 12
Corresponding convergence for error estimation in function of the iterations, where the simulated time ist=50: (a) for P (b) for Z, and (c) for F.
(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.



---

*Source: 2037093-2018-05-10.xml*