A Spatiotemporal Prey-Predator Discrete Model and Optimal Controls for Environmental Sustainability in the Multifishing Areas of Morocco

In this work, we propose a multifishing area prey-predator discrete-time model which describes the interaction between the prey and middle and top predators in various areas, which are connected by their movements to their neighbors, to provide realistic description prey effects of two predators. A grid of colored cells is presented to illustrate the entire domain; each cell may represent a subdomain or area. Next, we propose two harvesting control strategies that focus on maximizing the biomass of prey, in the targeted area, and minimizing the biomass of middle and top predators coming from the neighborhood of this targeted area to ensure sustainability andmaintain a differential chain system.'eoretically, we have proved the existence of optimal controls, and we have given a characterization of controls in terms of states and adjoint functions based on a discrete version of Pontryagin’s maximum principle. To illustrate the theoretical results obtained, we propose numerical simulations for several scenarios applying the forward-backward sweep method (FBSM) to solve our optimality system in an iterative process.


Introduction
e interactions between prey and predator are considered to be one of the fundamental problems of the complex food chain and trophic network; they are also one of the relationships between the basic species in biology and the environment. To the best of our knowledge, mathematical ecology was born with the works of Lotka [1] and Volterra [2]. ey independently and almost simultaneously proposed the first mathematical model to attempt to describe the interaction between a prey population and a predator population. is model is defined by a system of two differential equations. Since the work of Volterra [2], a large part of the ecological literature has been devoted.
Many works have been done on the subject of the dynamics of fish populations where the last ones are assumed in competition [3][4][5][6][7][8][9]. In [3], the authors have proposed to define a bioeconomic model of two fish species. ey have considered that the evolution of these fish species is described by a densitydependent model taking into account the competition between species which compete with each other for space or food. In this model, they have assumed that there exist "n" fishermen who catch these two fish species. In [7], the authors have presented a bioeconomic model for several fish populations. ey have also assumed that these fish populations compete with each other for space or food. e techniques and issues associated with the bioeconomic modeling for the exploitation of marine resources have been discussed in detail in [10][11][12][13][14].
Concerning the prey-predator model, one can refer, for example, to [15,16]. In [15], the authors have considered the combined harvesting of a prey-predator system in which both the prey and the predator species obey the law of logistic growth, and some preys avoid predation by hiding. In [16], the authors have proposed to define a bioeconomic model that merges a model of competition and a model of prey-predator of three fish populations. About the fishery management model on several fishing areas, the reader can refer to the works of Mchich et al. [17][18][19][20]. In [20], Mchich et al. studied the optimal spatial distribution of the fishing effort in a multifishing zone, and they have given some efficient management measures by setting an appropriate system of tax and/or subsidies. Moreover, they have controlled the displacement of the fleets between the fishing zones in order to enlarge the complete activity. In [21], Mchich et al. studied the effects of predator-dependent migration rates of prey on the dynamics and the stability of the global system. Auger et al. [22] have proposed a predator-prey model in which predators can become infected by a disease. ey have also considered two time scales, a fast one for the disease and a comparatively slow one for predator-prey interactions and for predator mortality. In [23], the authors have looked for existence of a desirable situation such as a sustainable fishery in the sense that the resource can be exploited in a durable way but without any risk of extinction of the fish stock. Mchich et al. [19] have introduced a time-dependent control function in equations describing the fishing effort variation. is control is regarded as an investment proportion of the fishing income for each fleet.
In the aquatic ecosystem, all aquatic meal chains are primarily based on plankton, which are composed of phytoplankton and zooplankton. However, nutrients play a big role in the growth of plankton. As far as we know, the nutrient-plankton systems have been considered by a number of scholars [24][25][26][27][28][29][30][31].
As pointed out in [32], phytoplankton in particular have great contribution for our earth; for example, they furnish food for marine existence and oxygen for human life, and they additionally take in greater than half of carbon dioxide which may be contributing to global warming. Over recent years, many authors have studied phytoplankton-zooplankton models [24,30,[33][34][35][36][37][38][39]. In addition, many other prey-predator models have been developed, and the authors obtained the existence of global Hopf bifurcation in the predator-prey system with delay [40][41][42][43]. e absence of plankton has serious consequences for the survival of marine fishery resources.
ere is an increasing interest in the study and application of spatial spread [17,19,20]. Most of the models studied have been partly continuous because of their mathematical tractability.
In our work, we devise a discrete-time prey-predator model that describes the evolution of different fish populations which migrate between M areas.
In this work, we are more interested in devising a mathematical model, based on multiarea interaction, to describe the spatial diffusion of fish growth evolution that appears in a global area of interest Ω represented uniform in size, through the way of a grid of colored cells (see Figure 1). ese cells contain our different fish populations and represent subdomains of Ω, indicating that only one of these cells is targeted by our control strategy.
In epidemiology system, Zakary et al. [44], considered a spatial spread of an epidemic within a domain Ω, which is divided into different cells of uniform size and denoted by C pq .
We assume that the development of the fish-population chain is carried out by the interaction of fish from the studied area C pq to its neighbor areas. In fact, by measures, C pq can represent a small lake, dam, or maritime area. is means that these present cellular representations which can be useful for studying the evolution of various fish species.
In our case, we will assume that C pq refers to some important maritime areas of Morocco. Figure 2 illustrates an example of discrete fishing areas of Morocco. Births are included in our model. One reason of considering the discrete-time prey-predator model is that the discrete model has advantages in describing the evolution of fish-population biomass since the ecological data are usually collected in discrete time units, which would make it more convenient to use the discrete-time model [45].
ere is increasing interest in studying dynamics behavior of discrete-time predator-prey systems [46][47][48][49]. Many models were developed in continuous space and time via reaction-diffusion-based predator-prey models [50,51]. For underlying the spatiotemporal dynamic of predatorprey models, many researchers focused on two-dimensional maps [48,52,53]. In [48], the authors considered a predatorprey system that is discrete both in space and time and is described by a coupled map lattice. In their work, they have assumed that the prey affected by a weak Allee effect and the predator dynamics includes intraspecific competition.
is paper is situated in this general context. e aim of this work is to consider a multifishing area discrete-time prey-predator model for three species (prey, middle predator, and top predator) using an optimal control problem applicable to any type of species. We give a more general multiobjective optimization approach in which we consider two harvesting controls. e aim of these two optimal harvesting control strategies is to make fishing efforts by the fishing fleets during fishing in the studied area C pq while increasing one harvesting function to harvest the middle and top predators that threaten the prey over time to ensure environmental sustainability and maintain a differential chain system. ese two harvesting controls are designed in theory to address the impact of predators over time on prey growth in the C pq area. So, these controls play a key role in the ecological balance and thus in the preservation of aquatic food webs. erefore, the purpose of these two harvesting control strategies is to maximize the total number of prey and to minimize the number of middle and top predators in the C pq area by maximizing the two harvesting functions and using the minimum possible cost of applying these controls. We prove the existence of optimal control, and some techniques are used to characterize the optimal control pair in terms of states and adjoint functions. e optimality system is based on a discrete iterative scheme that converges following an appropriate test related to the forward-backward sweep method (FBSM) [54]. e numerical simulations of our strategies of control show that the control effect is effective if the strategies are used simultaneously. Other models from population dynamics and optimal control problems can be found in [55][56][57][58][59][60][61][62][63][64][65][66][67]. Basic notions and results concerning optimal control problems can be found in the books [68,69]. e paper is structured as follows: in Section 2, we present the multifishing area discrete prey-predator model without the control. In Section 3, we present the optimal control problem for the proposed model where we give some results concerning the existence of the optimal controls, and we characterize these optimal controls using Pontryagin's maximum principle in discrete time. As an application, the numerical results associated with our control problem are given in Section 4. Finally, we conclude the paper in Section 5.

The Multifishing Area Prey-Predator
Discrete-Time Model  Figure 1: [44] Cell grid of (a), the colored contour of C pq is its neighborhood V pq and (b), an example of a patch of 4 cells.  e Atlantic coasts of Morocco are characterized by their richness and biological diversity. Morocco's exclusive area is characterized by great diversity of fishery resources. ere are four fishing areas in Morocco whose relative importance in terms of activity has undergone a great change over time. e first, area 1 is the Mediterranean area between Saadia and Cap Spartel, the second, area 2 is the northern area between Cap Spartel and Cap Cantin, the third, area 3 is the center area between Cap Cantin and Cap Boujdour, and latest one, area 4 is the area between Cap Boujdour and Cap Blanc.
interest Ω, split to M × M areas or cells, uniform in size. erefore, our domain Ω can be represented as follows: C pq , with C pq p,q�1,...,M , a spatial area. (1) is model classifies the population into three compartments in the area C pq . Let x C pq i , y C pq i , and z C pq i denote the population densities of the prey, middle predator, and top predator, respectively. Let n C pq i be the number of biomass of C pq at time i (the biomass of fish population living in C pq ).
We propose the following assumptions to formulate our model which describe the evolution of the biomass of our three fish populations.
e population of prey x C pq i grows according to a logistic equation with carrying capacity k 1 and birth rate constant that of a pq , and these species are favorite foods for middle predators y C pq i (hence, the population density of middle predators y C pq i will increase).

Assumption 2.
e population of middle predators y C pq i grows according to a logistic equation with carrying capacity k 2 and birth rate constant that of b pq .
Assumption 3. e population of middle predators is the best food for the populations of top predators (hence, the population density of top predators z C pq i will increase); this later population also grows according to a logistic equation with carrying capacity k 3 and birth rate constant that of C pq .
We consider that there are interactions among those three prey-predator compartments for time unit i to time i + 1. e unit of time i can correspond to days, months, or years; it depends on the frequency of data collection and statistics. e relation between the areas C pq and C rs is defined as follows.
We define the vicinity set of C pq , the neighborhood of a cell C pq , by We assume that all these areas contain fish populations of prey and middle and top predators. We also assume that the interaction is faster than the growth of these species and no natural mortality. We present the following scenario: the fish populations of the prey and middle and top predators can move to different studied areas.
We also assume that all possible interactions between fish population R C pq i , located in C pq , and fish population S C rs i , predators coming from V pq , can be expressed mathematically as follows: i . e evolution of biomass of different fish populations with the presence of a harvesting activity is modeled by the three following equations of the dynamics within a domain Ω: Figure 3 illustrates an instance of discrete geographical domains of four major fishing areas in Morocco.
In order to show the effectiveness of the proposed model and the effect on the prey in the targeted area and middle and top predators from neighboring areas, we will present numerical simulation over a period of t � 200 days with  to ensure that the model adapt to the reality ; all simulations are performed using the parameter values in Table 1 taken from [70]. In Section 3, we present numerical values and methods. Figures 4-6 present the numerical results of biomass of the prey, middle predators, and top predators. We find that there is no significant effect until t � 50 days after the interaction of the prey and middle and top predators. e biomass of the prey sharply decreases that may lead to their absence, while there is a significant increase of biomass of middle predators until t � 150 days and just after we observe a remarkable decrease until t � 200 of these species. However, we observe that there is a big increase of top predators during these 200 days. e remarks observed doing these simulations lead us to think about the definition 4 Discrete Dynamics in Nature and Society of appropriate control strategies taking into account these remarks. e chosen strategy is to consider two harvesting functions in order to make fishing effort controls by the fishing fleets during fishing in the studied area C pq .

e Model with Controls.
In this section, we also introduce two harvesting controls H(y i ) and H(z i ) as functions of time. e aim is to make fishing effort controls by the fishing fleets during fishing in the studied area C pq while adding one harvesting function to harvest the middle and top predators y C pq i and z C pq i which threaten the prey x C pq i over time to ensure environmental sustainability and maintain a differential chain system. en, for a given area C pq , the model is given as follows:     Discrete Dynamics in Nature and Society with i � 0, . . . , N − 1 and x C pq H(y i ) and H(z i ) represent the harvesting functions associated to middle and top predators, respectively, with ξ rs and δ rs representing the catchability coefficients of fishing fleets, coming from C rs , to harvest middle predators and top predators from C pq .
And u C rs i and v C rs i represent the fishing effort controls of fishing fleets, coming from C rs , to harvest middle predators and top predators from C pq .
en, for a given area C pq , the model is given by  Discrete Dynamics in Nature and Society with i � 0, . . . , N − 1, and

An Optimal Control Approach.
We are interested in controlling the fishing area C pq . en, the problem is to minimize the objective function given by subject to system (7). Here, Ψ 1 , Ψ 2 , and Ψ 3 are positive constants to keep a balance in the size of x C pq i , y C pq i , and z C pq i , respectively.
In the objective function, A rs and B rs are the positive weight parameters, which are associated with the controls u C rs � (u Our goal is to minimize the number of middle predators y C pq and top predators z C pq while maximizing the harvesting function and increasing the number of prey x C pq in C pq .
In other words, we are seeking optimal controls u C rs * and v C rs * such that where U and V are the sets of admissible controls defined by where (u  2 . e sufficient condition for existence of an optimal control for the problem follows from the following theorem.

Theorem 1 (sufficient conditions).
For the optimal control problem given by (10), along with state equations (7), there exists a control (u C rs * ), (v C rs * ) such that J pq u C rs * , v C rs * � max J pq u C rs , v C rs u C rs ∈ U , v C rs ∈ V . (12) Proof. See eorem 1 of Dabbs [71]. At the same time, by using Pontryagin's maximum principle [72], we derive necessary conditions for our optimal control. For this purpose, we define the Hamiltonian as 8 Discrete Dynamics in Nature and Society □ Theorem 2. Given optimal controls u C rs * and v C rs * and solutions x C pq * , y C pq * , and z C pq * , there exists ζ k,i C pq , i � 0, . . . , N − 1, k � 1, 2, 3, the adjoint variables satisfying the following equations: with transversality conditions where u C rs * � (u Proof. Using Pontryagin's maximum principle [72] and setting x C pq � x C pq * , y C pq � y C pq * , z C pq � z C pq * , and 3,N � 0. To obtain the optimality conditions, we take the variation with respect to controls u C rs i , v C rs i and set it equal to zero: en, we obtain the optimal control pair By the bounds in U and V of the controls, it is easy to obtain u C * rs i and v C * rs i in the following form: □

Numerical Simulation
Now, we give the numerical results associated to system (7), and we estimate that the initial population of the prey and middle and top predators is given by x 0 � 2250, y 0 � 1100, and z 0 � 400. In order to show the importance of our work and without loss of generality, we consider a 8 × 8 grid where we focus in this grid on diagonal patch Ω 1 , Ω 2 , Ω 3 and Ω 4 with Ω 1 � C 11 , C 12 , C 21 , C 22 , Ω 2 � C 33 , C 34 , C 43 , C 44 , Ω 3 � C 55 , C 56 , C 65 , C 66 , and Ω 4 � C 77 , C 78 , C 87 , C 88 . It is Table  1: e description of parameters used for the definition of discrete-time system (7). Discrete Dynamics in Nature and Society therefore assumed that the prey distribution is not homogeneous with 600 in Ω 1 , 500 in Ω 2 , 350 in Ω 3 , and 800 in Ω 4 , where we placed 400 in Ω 1 , 200 in Ω 2 , 150 in Ω 3 , and 350 in Ω 4 for middle predators, while for the top predators, we placed 100 in Ω 1 , 100 in Ω 2 , 50 in Ω 3 , and 150 in Ω 4 for our numerical simulation cited in Table 1. All simulations are performed using the parameter values in Table 2 which is taken from [70]. Also, the upper limits of the optimality condition are considered to be u C rs max � 0.8 and v C rs max � 0.9 In this section, we present numerical results that illustrate and reinforce the effect of our control strategy. is strategy consists in applying two types of controls in order to preserve and increase the biomass of the prey and also to ensure environmental sustainability and maintain a differential chain system. Concerning the numerical method, we give numerical simulations to our system of optimality (7) which is formulated by three equations of states with initial conditions (8), adjoint equations with transversality conditions (14) and (15), and optimal control characterization (16) and (17). We apply the forward-backward scanning method (FBSM) [54] to solve our optimality system in an iterative process, and we use the explicit method of Euler of second order; the initial control variables are guessed at the beginning of the iterative method, and then the adjoint equations are resolved backwards in time. Finally, the control variables are updated with the current state. e iterative process is repeated until a tolerance criterion is reached.
In this paragraph, in order to obtain more accurate information about the impact of each control separately, we will develop our strategy by applying each control individually, and numerical analysis will show us the effectiveness of each strategy. In the first case, we will limit ourselves to controlling the four studied areas only using one of the two strategies u or v by making fishing effort control by the fishing fleets in these areas. In the second case, we will control all areas using both controls u and v.
(i) Case 1: applying only control u (ii) Case 2: applying both controls u and v 3.1. Case 1: Applying Only Control u. In this case, the control applied to the top predators (v Cpq ) is not optimized, while the control applied to the middle predators (u Cpq ) is optimized to see its impact on the evolution of the three species (prey, middle predators, and top predators). According to Figures 7-9, we investigate numerical results through a control strategy on the middle predators. In this case, we observe a decrease in the number of prey with their continuity in the studied areas, in a lesser way compared to the results obtained in Figure 4, in which we note the disappearance of these prey.
ere is also an increase in the number of middle and top predators during the 200 days. Figures 10-12, we note the effectiveness of applying simultaneously the strategies of spatial-temporal control u C pq i and v C pq i . ese controls are highly noticeable in maintaining the biomass of the prey from absence and thus maintaining a different chain system. We note that the number of prey decreases at a very low rate during the period of 200 days, while we note a slight increase in the number of middle and top predators. Table 2: List of all parameters of system (4).

Parameter
Physical interpretation a pq Intrinsic growth rate of a prey from a cell C pq b pq Intrinsic growth rate of a middle predator from a cell C pq c pq Intrinsic growth rate of a top predator from a cell C pq c pq Prey movement rate coming from its neighbor cell C rs ∈ V pq c rs Prey movement rate from a cell C pq to neighbor cell C rs ∈ V pq c C pq 1 Middle predator movement rate coming from its neighbor cell C rs ∈ V pq c C rs 1 Middle predator movement rate from a cell C pq to neighbor cell C rs ∈ V pq c C pq 2 Top predator movement rate coming from its neighbor cell C rs ∈ V pq c C rs 2 Top predator movement rate from a cell C pq to neighbor cell C rs ∈ V pq α rs Per capita predation of adequate contacts between a prey from a cell C pq and a middle predator coming from its neighbor cell C rs ∈ V pq α rs Conversion rate predation of adequate contacts between a prey from a cell C pq and a middle predator coming from its neighbor cell C rs ∈ V pq β rs Conversion rate predation of adequate contacts between a prey from a cell C pq and a top predator coming from its neighbor cell C rs ∈ V pq β rs Per capita predation of adequate contacts between a prey from a cell C pq and a top predator coming from its neighbor cell C rs ∈ V pq α 1 rs Per capita predation of adequate contacts between a middle predator from a cell C pq and a top predator coming from its neighbor cell C rs ∈ V pq α 1 rs Conversion rate predation of adequate contacts between a middle predator from a cell C pq and a top predator coming from its neighbor cell C rs ∈ V pq

Conclusion
In this paper, we have considered a multiarea prey-predator model with discrete time describing the dynamics of interaction between different species in different areas via movements. e application of a distributed optimal control pair for a spatiotemporal prey-predator model is incorporated to preserve the biomass of the prey. We have shown the existence of optimal controls to give a characterization of controls in terms of states and adjoint functions. e optimality system, which is composed by the system state, the dual system, and the characteristic of the control, is solved numerically based on the forward-backward sweep method (FBSM). Numerical simulations of the resulting optimality system have shown that the optimal strategy is more efficient if we apply simultaneously the two controls, which allow us to guarantee sustainability and maintain a differential chain system.

Data Availability
e disciplinary data used to support the findings of this study have been deposited in the Network Repository (http://www.networkrepository.com).

Conflicts of Interest
e authors declare that they have no conflicts of interest.