Next Article in Journal
Identifying the Distribution and Frequency of Dust Storms in Iran Based on Long-Term Observations from over 400 Weather Stations
Next Article in Special Issue
A Harmful Algal Bloom Detection Model Combining Moderate Resolution Imaging Spectroradiometer Multi-Factor and Meteorological Heterogeneous Data
Previous Article in Journal
Phenotypic, Geological, and Climatic Spatio-Temporal Analyses of an Exotic Grevillea robusta in the Northwestern Himalayas
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Systematic Review

Harvested Predator–Prey Models Considering Marine Reserve Areas: Systematic Literature Review

by
Arjun Hasibuan
1,
Asep Kuswandi Supriatna
2,*,
Endang Rusyaman
2 and
Md. Haider Ali Biswas
3
1
Doctoral Program of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Padjadjaran, Sumedang 45363, Indonesia
2
Department of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Padjadjaran, Sumedang 45363, Indonesia
3
Mathematics Discipline, Khulna University, Khulna 9208, Bangladesh
*
Author to whom correspondence should be addressed.
Sustainability 2023, 15(16), 12291; https://doi.org/10.3390/su151612291
Submission received: 19 June 2023 / Revised: 21 July 2023 / Accepted: 27 July 2023 / Published: 11 August 2023
(This article belongs to the Special Issue Sustainable Management and Conservation of the Oceans)

Abstract

:
The United Nations has predicted the growth of the human population to reach 8.405 billion by mid-2023, which is a 70% increase in global food demand. This growth will significantly affect global food security, mainly marine resources. Most marine resources exist within complex biological food webs, including predator–prey interactions. These interactions have been researched for decades by mathematicians, who have spent their efforts developing realistic and applicable models. Therefore, this paper systematically reviews articles related to predator–prey models considering the harvesting of resources in marine protected areas. The review identifies future remodeling problems using several mathematical tools. It also proposes the use of feedback linearization consisting of both the approximation and exact methods as an alternative to Jacobian linearization. The results show that in an optimal control analysis, adding a constraint in the form of population density greater than or equal to the positive threshold value should be considered to ensure an ecologically sustainable policy. This research and future developments in this area can significantly contribute to achieving the Sustainable Development Goals (SDGs) set for 2030.

1. Introduction

The human population is currently experiencing continuous growth and is projected to reach 8.405 billion by mid-2023, as reported by the United Nations (UN) [1]. This growth would undoubtedly lead to an increased demand for food worldwide. However, a significant portion of the global population still lacks access to sufficient food, raising concerns about potential disruptions in food security. By 2050, experts predict that the demand for food will soar, surpassing the current food demand by 70% [2].
Food security is a pressing global concern, and world leaders have recognized the significance of this issue by incorporating it into the Sustainable Development Goals (SDGs) program. Based on an annual report by the Food and Agriculture Organization (FAO), the issue of food security has yet to be fully resolved, and the target to achieve it by 2030 leaves the entire population with seven to eight years [3]. In the SDGs program, four of the seventeen goals set are closely tied to food security. Among these, Goal 14 is intended to protect marine ecosystems, which have declined substantially in recent decades, mainly due to fishing activities [4]. This decline has endangered numerous commercially harvested marine organisms, including mammals, birds, turtles, and some fish species [5]. There has been a significant surge in the consumption of aquatic food in recent years, and this trend is expected to continue [6]. With the increasing global population and changing dietary habits, the demand for aquatic food products tends to put significant pressure on marine ecosystems. Addressing food security and preserving marine ecosystems requires the coordinated efforts of all nations.
Preliminary research plays a vital role in achieving the SDGs program, specifically regarding addressing food-security issues through various scientific fields. One critical area where these investigations are expected to contribute is in managing fishery populations, as marine habitats serve as a significant source of human food. In this context, mathematics plays a significant role through mathematical modeling. Mathematical models allow for the formulation of real-world problems. These are described in mathematical terms or forms, providing valuable insights and solutions. Mathematical models are extensively used in fishery population management, particularly to understand predator–prey dynamics. For example, the problem of prudently managing Haliotis rubra abalone has driven the development of several mathematical models [7,8,9].
Predator–prey interactions are a common phenomenon in nature, where one or more species serve as prey, becoming food resources for several predators. To mathematically describe such interactions, the Lotka–Volterra model was proposed by Alfred J. Lotka and Vito Volterra. This model is widely used in ecological research to understand the dynamics of predator–prey relationships. According to Ibrahim [10], the Lotka–Volterra predator–prey model for two species is stated as follows
{ d x d t = a x b x y d y d t = c x y d y .
The predator–prey interaction model involves two variables, x and y , representing the densities of the prey and predator populations, respectively. Several parameters are essential in this model, a ,   b ,   c , and d , which depict the per capita growth rate of the prey population, the predation rate of prey, the growth rate of the predator population due to predation, and the mortality rate of predators, respectively. In the presence of selective harvesting, the model for one predator and prey affected by harvesting is stated as follows.
{ d x d t = a x b x y q 1 e 1 x d y d t = c x y d y q 2 e 2 y ,
where q i and e i   for i = 1 ,   2 are the catchability and effort level of prey ( i = 1 ) and predator ( i = 2 ) species, respectively. The functions q 1 e 1 x and q 2 e 2 y in Equation (2) are based on the CPUE (catch per unit effort) hypothesis, as stated by Ibrahim [10].
In systems involving predator–prey dynamics and harvesting, the exclusive exploitation of predators can lead to increased competition among prey and raise the risk of other species going extinct. A high predator abundance leads to a decrease in prey population density and fosters competition among these species. A decline in predator density can slow down the growth and reproduction of prey. On the other hand, when only prey are harvested, it could affect the growth of predators, potentially leading to predator extinction in circumstances where excessive harvesting of prey occurs. When both species are over-harvested, it ultimately leads to their extinction [11]. To tackle the issue of overexploitation, fishery managers must find effective solutions. One alternative policymakers can adopt is the use of a realistic and correct mathematical model of predators–prey and harvesting in fishery ecosystems. This model allows policymakers to identify sustainable harvesting strategies that maintain balanced predator–prey population densities, ensuring the long-term viability of both food and the fish population. By using this model, fishery managers can make informed decisions to preserve ecological stability and safeguard the sustainability of marine resources.
Several traditional strategies for managing fisheries’ resources have been researched, such as maximum sustainable, economic, and good yields. However, there is another strategy known as marine reserve or protected areas (MPAs). Unlike species-based methods, this strategy focuses on spatial considerations by dividing the marine zone into reserve and harvesting areas. The primary goals of an MPA are to enhance biodiversity and restore depleted fish stocks [4]. Fishing activities are restricted or prohibited by designating certain areas as reserves, allowing marine ecosystems to recover and flourish. As a result, MPAs can lead to increased catches in adjacent harvesting zones, benefiting both the fish population and the fishing industry [12]. When traditional strategies fall short, MPAs are considered a critical tool for biological conservation and sustainable resource management [4].
Based on previous research, marine reserve policy is an important solution for addressing food-security issues. Many mathematicians have conducted research on the predator–prey mathematical model considering marine reserves and harvesting areas [5,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]. This research aims to systematically review this model in the context of fishery ecosystems using the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) method. The PRISMA method uses detailed flow diagrams to enhance the suitability of collecting articles on the topic under review [30]. This process has three main stages, namely collection, selection, and review. The articles reviewed were sourced from the Scopus, Science Direct, Dimensions, and Web of Science databases. The articles were filtered through duplicate selection, relevance checks on titles and abstracts, the availability and relevance assessment of full papers, and, finally, reviewing the model and methods presented in the selected articles. The research discussed opportunities for further investigations in this area, aiming to contribute to solving food-security problems through mathematical modeling. In particular, these mathematical models aid policymakers in managing fisheries’ natural resources sustainably, preventing resource depletion and ensuring fishermen continue to gain profits (Supplementary Materials).

2. Research Methodology

This research reviewed mathematical models of predators in harvesting areas and nature reserves within fishery ecosystems using the PRISMA method. This method employed systematic flow diagram to ensure the articles selected for review were highly relevant to the topic discussed [30]. By adopting this method, the research maintained a rigorous and organized method, enhancing the credibility and accuracy of the review process.
Systematic review process in this research comprised three main stages, with the first two following the PRISMA method. The initial stage involved collecting articles, which occurred on March 7th, 2023. During this stage, a specific keyword was used to search the Scopus, Science Direct, Dimensions, and Web of Science databases. The keywords applied to each database are shown in Table 1.
In this systematic literature review, data filtering was meticulously performed for each keyword in four distinct databases.
  • Scopus database
  • Search within: Title, abstract, and keywords;
  • Document type: Article;
  • Publication stage: Final;
  • Source type: Journal;
  • Language: English.
  • Science Direct database
  • Search by: Title, abstract, or author-specified keywords;
  • Article type: Research articles.
  • Dimensions database.
  • Search in: Title and abstract;
  • Publication type: Article.
  • Web of Science database.
  • Search in: All fields;
  • Publication type: Article.
Additionally, the review considered articles published between 2010 and 2023 from all four databases, ensuring the relevance of the selected literature.
The second stage of systematic review was the selection phase. It consisted of several sub-stages to ensure the inclusion of relevant articles. This included duplicate selection, title, abstract relevance selection, full paper availability, and relevance selection. The duplicate selection sub-stage was performed semi-manually using Jabref software version 5.9. When there were similar articles among the three databases, one of them would be selected as a representative, while the rest were deleted. Furthermore, the last three sub-stages were performed manually to ensure the relevance of the articles to the topic. In the title and abstract relevance selection sub-stage, the articles that passed the duplicate selection were individually assessed based on their titles and abstracts. This step helped identify relevant articles efficiently, thereby saving time. Articles unrelated to the topic were excluded from the dataset. In the full paper availability sub-stage, articles not available in full were excluded from consideration. Finally, the full paper relevance selection stage entailed a more rigorous evaluation of articles. The entire content of each article included in the full paper availability stage was thoroughly read to ensure its direct relevance to the topic. This stage may result in a smaller number of articles compared with the title and abstract relevance selection phase. By following these systematic sub-stages, the selection process ensured that only the most pertinent articles related to predator–prey mathematical models in harvesting areas and nature reserves within fisheries’ ecosystems were included. This meticulous method guaranteed the credibility and accuracy of the research results and analysis.
In the last stage, the article review phase, models and methods researched in each article selected during the full paper relevance selection were comprehensively examined. Moreover, this stage entailed discussing possibilities for further research to develop more a realistic or applicable model to address food-security issues through mathematical modeling. The potential developments resulting from this research have significant implications for tackling food-security challenges through advanced mathematical models. These advancements would notably assist policymakers in effectively managing fisheries’ natural resources, ensuring their long-term sustainability, and averting resource depletion. From a mathematical modeling perspective, this method enables fishermen to sustain their income while promoting responsible and sustainable practices.

3. Results and Discussion

3.1. Results of Article Collection and Selection

The present research obtained a total of 1672 articles from four databases, namely 106, 318, 1248, and 1171 from Scopus, Science Direct, Dimensions, and Web of Science, respectively. The details of these results are shown in Table 1.
According to the diverse stages of the PRISMA method, the next phase involved automatic duplicate selection using Jabref software. Of the initial 1648 articles, 1195 duplicates were excluded, leaving the remaining ones for further consideration based on their titles and abstracts. The subsequent stage was the manual title and abstract relevance selection. Each of the titles and abstracts of the diverse articles were carefully reviewed individually to determine their appropriateness and relevance to the research focus, which was the deterministic mathematical modeling of predator–prey and harvesting problems, with a specific examination of marine reserve areas. As a result, 23 articles were selected for full paper availability, while 1625 irrelevant articles were eliminated. The full paper availability selection stage revealed that twenty articles were accessible, while the remaining three were unavailable. Based on full papers, the relevance selection stage confirmed that all 20 articles were pertinent to the problems associated with discussing the existence of a mathematical model, predator–prey relationships, harvesting, and marine reserve areas. All articles at this stage met the stipulated criteria. The articles, written by Khamis et al. [23], Srinivas et al. [24], Chakraborty et al. [25], Chakraborty et al. [27], Chakraborty and Kar [26], Chakraborty et al. [28], Lv et al. [29], Kar and Ghosh [5], Sharma and Gupta [13], Louartassi et al. [14], Agnihotri and Nayyer [15], Louartassi et al. [16], Pei et al. [17], Huang et al. [18], Huo et al. [19], Zhang et al. [20], Krivan and Jana [12], Huang et al. [21], Meng et al. [22], and Ibrahim [10], were thoroughly reviewed to address the research question in Section 1. The dataset for these twenty articles can be accessed at https://bit.ly/Database20Articles (accessed on 17 June 2023). Finally, the results of this collection and selection process are summarized using the PRISMA diagram shown in Figure 1.

3.2. Overview of Model from Each Article and Potential Model to Develop for Future Research

Twenty articles were selected for review based on the included stage shown in Figure 1. These articles explored the different categories of models in terms of the number of compartments used therein. Among the twenty articles, the most commonly researched model comprises three compartments, including one prey and predator species. Additionally, six out of the twenty articles examined a model with four compartments. The remaining three articles focused on a model with two compartments. Lastly, one article delved into a model with five compartments. The diverse range of models provides valuable insights into the dynamics of predator–prey relationships in fishery ecosystems.
The first category of models, i.e., those with three compartments, divided prey species into those present in the nature reserve and harvesting area. In contrast, predator species remained undivided in terms of their habitat. This concept was explored in several articles, including Khamis [23], Chakraborty et al. [25], Chakraborty and Kar [26], Chakraborty et al. [28], Lv et al. [29], Sharma and Gupta [13], Louartassi et al. [14], Agnihotri and Nayyer [15], Louartassi et al. [16], and Huo et al. [19]. In more detail, Sharma and Gupta [13], Louartassi et al. [14], and Agnihotri and Nayyer [15] reported that the prey and predator species in this model are fish and birds, respectively. Khamis [23], Chakraborty et al. [25], Chakraborty et al. [27], Chakraborty and Kar [26], Chakraborty et al. [28], and Lv et al. [29] stated that both the prey and predator species in this model are fish. In the three-compartment model, the variables r , u , and U represent prey species in the marine reserve, harvesting, and predator area, respectively. The general model for the three-compartment analysis is stated as follows:
{ d r d t = G r ( r ) M r ( r ) + M u ( u ) P r ( r , U ) d u d t = G u ( u ) + M r ( r ) M u ( u ) P u ( u , U ) H u ( u ) d U d t = G U ( U ) + N r ( r , U ) + N u ( u , U ) H U ( U ) .
The functions G r ( r ) ,   G u ( u ) , and G U ( U ) represent the intrinsic growth of prey species in the marine reserve ( r ) and harvesting area ( u ) as well as that of the predator species ( U ) , respectively. Additionally, the functions M r ( r ) and M u ( u ) denote the migration functions of prey species in the marine reserve and harvesting area. The functions H u ( u ) and H U ( U ) correspond to the harvesting functions of the prey and predator species. P r ( r , U ) and P u ( u , U ) describe the predation functions between predator and prey in the marine reserve and harvesting area. Finally, the functions N r ( r , U ) and N u ( u , U ) are the numerical response functions of predator species that depend on the prey in the marine reserve and harvesting area. A numerical response refers to how predation affects either the reproduction or mortality rate of predators based on the population density of prey, predators, or both [31]. These functions were widely used in all the articles that discussed a predator–prey model of fisheries’ ecosystems in marine reserves and harvesting areas.
Khamis [23], in 2011, researched a marine reserve model and proposed the following hypothesis:
{ d r d t = g r r ( 1 r K r ) + m u u m r r p r r U d u d t = g u u ( 1 u K u ) m u u + m r r p u u U h u E u u d U d t = n u u U + n r r U d U c U 2 .
Based on Equation (4), the functions G r ( r ) and G u ( u )   used by Khamis [23] are logistic growth functions, namely
G r ( r ) = g r r ( 1 r K r )   and   G u ( u ) = g u u ( 1 u K u ) .
The parameters g r > 0 and g u > 0 depict the intrinsic growth rates of prey species in the marine reserve and harvesting area, respectively. Then, K r > 0 and K u > 0 represent the carrying capacity of prey species in the marine reserve and harvesting area. In addition, the researched migration functions M r ( r ) and M u ( u ) exhibit a linear equation, stated as follows.
M r ( r ) = m r r   dan   M u ( u ) = m u u .
The migration rates of prey species in the marine reserve and harvesting area are represented by the parameters 0 < m r 1 and 0 < m u 1 , respectively. Furthermore, predation on prey species in the marine reserve and harvesting area follows the Holling Type (HT) I predation function form. This form is stated as follows.
P r ( r , U ) = p r r U   and   P u ( u , U ) = p u u U .
The parameters 0 < p r 1 and 0 < p u 1 represent the predation rate for prey species in the marine reserve and harvesting area, respectively. In Equation (4), harvesting only occurs for prey species in the harvesting area and does not affect predators. The harvesting function used in Equation (4) is stated as follows
H u ( u ) = h u E u u .
The parameter 0 < h u 1 represents the harvesting rate of prey species in the harvesting area. The model does not include a natural reproductive function regarding predator species. Instead, in Equation (4), predator reproductive function is directly related to the predation effect expressed by the numerical response function. As predation occurs on both types of prey species, there are two numerical response functions, stated as follows.
N r ( r , U ) = n r r U   and   N u ( u , U ) = n u u U .
The parameters n r > 0 and n u > 0 represent the predation conversion rate with respect to predator species growth or reproduction. Khamis introduced two additional terms not listed in Equation (3), namely the competition effect c U 2 and natural mortality d U . The parameters c and d represent the level of competition and natural mortality of predator species, respectively.
In 2011, Chakraborty et al. [25] conducted unique research on predation which specifically targeted species within the harvesting area. The influence of HT II or the Michaelis–Menten predation function in this particular context was investigated and stated as follows.
P u ( u , U ) = p u u U a u + U .
The parameter a u denotes the half-capturing saturation constant. The numerical response function considers the time delay caused by predator gestation. In other words, predators need a delay in pregnancy after preying on a prey population, which directly influences the effect of predation. The numerical response function is stated as follows.
N u ( u , U ) = n u p u u ( t τ ) U ( t τ ) a u + u ( t τ ) ,
where 0 < n u 1 and 0 < p u 1 denotes the predation conversion and predation rates, respectively. In addition, Chakraborty et al. [25] added a constraint in the form of an algebraic equation to the designed model. The algebraic equation is stated as follows:
( p h u u c ) E u m = 0 ,
The parameters p , c , and m represent the constant cost of fishing per unit of effort, the constant price per unit of harvested fish biomass, and the total profit earned from fishing, respectively. The model was designed to ensure that the fishing profit remains constant at the value of m . Chakraborty et al. [25] developed an algebraic differential equation model which is stated as follows:
{ d r d t = g r r ( 1 r K r ) + m u u m r r d u d t = r u u ( 1 u K u ) m u u + m r r p u u U a u + u h u E u u d U d t = n u p u u ( t τ ) U ( t τ ) a u + u ( t τ ) d U c U 2 ( p h u u c ) E u m = 0 .
The models in (4) and (5) are modified variations of the Lotka–Volterra predator–prey model, building upon its basic structure. There is another form of predator–prey model, known as the Leslie–Gower predator–prey model. This simpler model considers one prey ( x ) and predator ( y ). The Leslie–Gower model is equivalent to the one presented in (1), which is stated as follows:
{ d x d t = a x b x y d y d t = y ( c d y x ) .
The difference between the Lotka–Volterra and Leslie–Gower model is evident in the second equation of (6). In the Leslie–Gower model, this equation includes the term d y x . This term simply means that as prey population increases and x , the change in predator density reaches a maximum value, depicted by d . Conversely, a assuming predator becomes extinct, x would also face extinction. The Leslie–Gower model states that the carrying capacity of a predator environment is directly proportional to the availability of prey [32]. The ability of a predator population to thrive depends on the abundance of prey in its ecosystem.
In 2012, Chakraborty and Kar [26] developed a model based on the Leslie–Gower research. Meanwhile, in 2013, Chakraborty et al. [28] also designed another model in accordance with the same basis. The first and second equations of Chakraborty and Kar [26] and Chakraborty et al. [28] share certain similarities and are stated as follows
{ d r d t = g r r ( 1 r K α ) m ( r α K u ( 1 α ) K ) p r r U a r + r d u d t = r u u ( 1 u K ( 1 α ) ) + m ( r α K u ( 1 α ) K ) p u u U a u + u H u ( u ) .
The primary difference between the two research efforts lies in the harvesting function used in each article. Chakraborty and Kar [26] applied the same harvesting function used by Khamis [23] and Chakraborty et al. [25]. Unlike the research conducted by Chakraborty et al. [28], the harvesting function used is stated as follows
H u ( u ) = h u E u u γ u E u + v u u
Equation (7) employs the same predation function as Chakraborty et al. [25], with γ u and v u as positive constants. Furthermore, the functions G r ( r ) and G u ( u ) in (7) exhibit minor differences. According to the research by Khamis [23] and Chakraborty et al. [25], Chakraborty and Kar [26], as well as Chakraborty et al. [28], set K r = α K and K u = ( 1 α ) K . Additionally, the migration function used differs from that of Khamis [23] and Chakraborty et al. [25]. The migration rate is proportional to the carrying capacity of each marine reserve and harvesting area. As a result, the carrying capacity plays a critical role in influencing the migration of each prey species.
The part yet to be presented in (7) is the equation of predator species. The predator equation formulated by Chakraborty and Kar [26] is stated as follows
d U d t = s U ( 1 γ U r + u ) .
Then, the predator equation formulated by Chakraborty et al. [28] is stated as follows
d U d t = s U ( 1 γ U r + β u ) d U .
The difference between the two predator equations lies in the parameter β and the inclusion of natural mortality, denoted by the term d U . However, they share certain common parameters, such as s ,   γ , and β , which represent the intrinsic growth rate of predator species, the number of prey required to support one predator at equilibrium, and a prey–prey conversion factor, respectively.
The research conducted by Lv et al. [29] in 2013 focuses on a model that closely resembles the one in (3) with time τ = 0 (no time delay in the numerical response function). The fourth equation of (3) was not included in the research conducted by Lv et al. [29]. Therefore, the Lv et al. [29] model is simpler than the one Chakraborty et al. researched [25]. The models proposed by Lv et al. [29] and Sharma and Gupta [13] in 2014, as well as Khamis [23], exhibited significant similarities. Both Lv et al. [29] and Sharma and Gupta [13] did not assume any form of competition in predators in their model. The model can be stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u P r ( r , U ) d u d t = g u u ( 1 u K u ) + m r r m u u P u ( u , U ) h u E u d U d t = d U + N r ( r , U ) + N u ( u , U ) h U E U U u .
The two models differ in their respective predation functions and functional and numerical responses. Specifically, Lv et al. [29] used a different predation function and numerical response, stated as follows
P r ( r , U ) = 0 ,   P u ( u , U ) = p u u U a u + u   ( HT   II ) ,   N r ( r , U ) = 0 ,   and   N u ( u , U ) = n u u U   ( HT   I )
Meanwhile, Sharma and Gupta [13] used distinct predation and numerical response functions, stated as follows:
P r ( r , U ) = p r r U ,   P u ( u , U ) = p u u U ,   N r ( r , U ) = n r r U ,   and   N u ( u , U ) = n u u U .
The functions P r , P u , N r , and N u correspond to the predation function ( P ) and numerical response function ( N )   targeting prey species within the marine reserve ( r ) and harvesting area ( u ) , respectively.
Agnihotri and Nayyer [15] stated that the model developed is a modification of Sharma and Gupta [13]. The primary difference is that Agnihotri and Nayyer [15] used the intrinsic growth function G u = g u u , a linear function. This difference lies in the second equation of (10), which is stated as follows:
d u d t = g u u + m r r m u u p u u U h u E u u
The research studies conducted by Lv et al. [29], Sharma and Gupta [13], and Agnihotri and Nayyer [15] are all three-compartment evaluations that included the harvesting of predator species.
Louartassi et al. [16] conducted research in 2019 examining a model that closely resembled the one proposed by Agnihotri and Nayyer [15], particularly regarding the G u function used. The model formulated by Louartassi et al. [16] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u t r r 2 c r r u d u d t = g u u + m r r m u u t u u 2 c u r u p u u U a u + u h u E u u d U d t = n u u U a u + u d U w U h U E U U .
Louartassi et al. [16] assumed predator harvesting similar to that used by Lv et al. [29], Sharma and Gupta [13], and Agnihotri and Nayyer [15]. In contrast to Khamis [23], Louartassi et al. [16] assumed that there was competition among prey species rather than predators. The competition among prey takes the form of interactions between species in the nature reserve and those in the harvesting area. The competition function is depicted by C r ( r , u ) = c r r u , and C u ( r , u ) = c u r u for species in the nature reserve and harvesting area, respectively. In addition, Louartassi et al. [16] introduced extra terms, t r r 2 and t u u 2 , to account for changes in prey density and w U to represent changes in predator density. These terms illustrate the impact of external toxic substances infecting prey species in the marine reserve and harvesting area, as well as predators. The use of quadratic and linear functions for prey and predator species, respectively, suggested that the effect of infection by external toxicity is more significant on prey compared with predators.
The research by Louartassi et al. [14] differs significantly from all the investigations on the three compartments. Louartassi et al. [14] stated that the developed model was a modification of the research by Sharma and Gupta [13]. The modification lies in the capture effort, denoted as E . The research on the other three compartments assumed E as a constant value that remains unchanged over time, while Louartassi et al. [14] assumed the function changes with time t , denoted as E ( t ) . They found that assuming E as a constant value does not align with real-world observations. According to their results, the required fishing effort should be reduced as the fish population increases. In order to address this fact, Louartassi et al. [14] referred to the fishing effort function proposed by Idels et al. [33]. The fishing effort functions proposed by Idels et al. [33] for prey species in the harvesting area and predators are stated as follows, respectively:
E u ( u ) = α u ( t ) β u ( t ) 1 u d u d t   and   E U ( U ) = α U ( t ) β U ( t ) 1 U d U d t .
The parameters α i 0 and β i 0 for i { u ,   U } are continuous functions of time t . Louartassi et al. [14] modified Sharma and Gupta [13], incorporating the fishing effort function proposed by Idels et al. [33]. The resulting model proposed by Louartassi et al. [14] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u p r r U d u d t = g u 1 h u β u u ( 1 u K u ) + m r 1 h r β r r m u 1 h u β u u p u 1 h u β u u U h u α u 1 h u β u d U d t = d 1 h U β U U + n r 1 h U β U r U + n u 1 h U β U u U h U α U 1 h U β U U U .
In the research model proposed by Louartassi et al. [14], it was assumed that n r = n u = n and p r = p u = p in Equation (13). This assumption resulted in Equation (13) being equivalent to the model represented in (14).
{ d r d t = g r r ( 1 r K r ) m r r + m u u p r U d u d t = g u 1 h u β u u ( 1 u K u ) + m r 1 h u β u r m u 1 h u β u u p 1 h u β u u U h u α u 1 h u β u u d U d t = d 1 h U β U U + n 1 h U β U ( r U + u U ) h U α U 1 h U β U U .
The research of Huo et al. [19] introduced a unique research model that sets it apart from all other three-compartment models. Furthermore, Huo et al. [19] examined compartments comprising prey species in both the harvesting and marine reserve area, as well as the harvesting effort. The model proposed by Huo et al. [19] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u d u d t = g u u ( 1 u K u ) + m r r m u u h u E u u d E u d t = [ α β [ q ( p τ ) u c ] γ ] E u .
The parameters α and β in the model represent coefficients of proportionality. Then, α β is a measure of the reaction between effort and perceived rent, known as the stiffness parameter. The term [ q ( p τ ) u c ] E u in the equation represents the net economic income earned by fishermen, with p , τ , and c denoting the constant price per unit of fish landed, applied tax per unit of fish landed, and constant cost per unit of harvesting effort, respectively. Huo et al. [19] introduced an additional assumption in Equation (15), where it was presumed that capital depreciates at a rate of γ . Therefore, the change in effort is influenced by the term [ α β [ q ( p τ ) u c ] γ ] E u . This specific term is a unique aspect that is not present in all models with three compartments.
Among the 20 articles on models with four compartments, 1 of them is the research conducted by Srinivas et al. [24] in 2011. Srinivas et al. [24] examined the model by dividing prey species into two areas, namely those in the marine reserve ( r ) and harvesting area ( u ) . Then, predator species were categorized into two groups, immature ( I ) and mature ( M ) predators. The research model proposed by Srinivas et al. [24] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u d u d t = g u u ( 1 u K u ) + m r r m u u p u u M h u E u u d I d t = α M p u u M 2 w I + M α I p u u I M w w I + M d M d t = α I p u u I M w w I + M d M h M E M M .
The parameters α M and α I depict the birth rate of mature predator species and the coefficient from immature to mature predator, respectively. Srinivas et al. [24] assumed that mature predators consume prey in the harvesting area in the ratio w : 1 . This ratio reflects the relative consumption between one immature and a mature predator. Additionally, Srinivas et al. [24] assumed that all weighted individuals received the same amount of prey in the harvesting area. As a result, a fraction of prey is consumed by mature predators, and this is expressed in terms of β u M M w I + M . Then, a small portion is consumed by immature predators, expressed as β u M w I w I + M . In Equation (16), harvesting occurs for both prey species in the harvesting area ( u ) and mature predator ( M ), as expressed by the following function.
H u ( u ) = h u E u u ,   and   H M ( M ) = h M E M M .
In 2012, Chakraborty et al. [27] researched a model that was perceived as the one developed by Chakraborty et al. [25] in 2011. While referencing the work of Chakraborty et al. [25], Chakraborty et al. [27] did not explicitly state that their model was a modification of the earlier one proposed by Chakraborty et al. [25]. The main difference between the two model lies in the fourth equation. Chakraborty et al. [27] added a constant, σ , and removed another, denoted as m , from the right-hand segment of Equation (5). The modified equation is stated as follows
[ ( p σ ) h u u c ] E u .
The addition of the parameter σ serves as a regulatory measure to control exploitation by imposing a tax σ per unit biomass of the landed fish population. This tax influences the net economic income earned by fishermen, as stated in Equation (17). Furthermore, Chakraborty et al. [27] assumed that the level of harvesting effort could vary over time, depending on the perceived rent, whether positive or negative. Chakraborty et al. [27] represented this dynamic reaction model by multiplying the stiffness parameter ( λ ) in Equation (17). This parameter measures the intensity of the reaction between harvesting effort and perceived rent. Under this assumption, the harvesting effort E u becomes a time-dependent dynamic variable, as stated in the following equation:
d E u d t = λ [ ( p σ ) h u u c ] E u .
Equation (18) is a modified version of the fourth equation of (5), which was formulated by Chakraborty et al. [25]. Therefore, the model researched by Chakraborty et al. [27] focuses solely on dividing prey species into the marine reserve and harvesting area without making further distinctions between predators.
There are similarities between the model formulated by Huang et al. [18] and Chakraborty et al. [27]. They both examined a model with one prey species, which was further divided into those present in the marine reserve and harvesting area. Additionally, a single undivided predator species was considered. Both models introduced a dynamic capture effort variable that changes with time t . The dynamic capture effort variable equation from Huang et al. [18] is similar to the one in (18). There are distinctions in the equations used for prey species in the marine reserve and harvesting area as well as predator species. The model proposed by Huang et al. [18] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r r + ρ r p r r U r + u d u d t = g u u ( 1 u K u ) + m r r r + ρ r p u u U r + u h u E u u d U d t = n r r U r + u + n u u U r + u d U d E u d t = λ [ ( p σ ) h u u c ] E u .
In Equation (19), migration occurs exclusively from species in the marine reserve and harvesting area. The migration rate used is not a constant m r but a function dependent on prey density in the marine reserve area, given by m r ( r ) = m r r r + ρ , where ρ is the half-saturation level of migration. This function shows that the natural migration rate m r is influenced by the saturation function of prey. When the species density in the marine reserve area goes to infinity, the migration rate goes to m r . Conversely, when the species density in the marine reserve area goes to zero, the migration rate goes to zero.
The model researched by Zhang et al. [20] in 2015 is similar to the one analyzed by Chakraborty et al. [27]. However, there is a difference in the harvesting function used for prey species, represented as H u ( u ) = h u E u u . The equation for predator species is markedly distinct from that researched by Chakraborty et al. [27]. The equation for predator species is stated as follows:
d U d t = g U ( 1 U η u ) .
Based on Equation (20), the predator–prey model base researched by Zhang et al. [20] is similar to that of the Leslie–Gower predator–prey model. Meanwhile, the model base researched by Chakraborty et al. [27] is rooted in the Lotka–Volterra predator–prey model. A significant difference from the model researched by Zhang et al. [20] lies in the predator–prey model base used.
In 2013, Kar and Ghosh [5] introduced a model that differed from those explored by Srinivas et al. [24], Chakraborty et al. [27], and Huang et al. [18]. This model divided both prey and predator species into two distinct areas, namely the prey species in the marine reserve ( r ) and harvesting area ( u ) , as well as the predator species in the marine reserve ( R ) and harvesting area ( U ). The Kar and Ghosh [5] model appeared to be an extension of the one formulated by Chakraborty and Kar [26] and Chakraborty et al. [28] in 2012 and 2013, respectively. Interestingly, both research studies were not cited in the Kar and Ghosh reference section [5]. Using the development process for the models by Chakraborty and Kar [26] and Chakraborty et al. [28], it could be inferred that the predator–prey model proposed by Kar and Ghosh [5] was also based on Leslie–Gower. In addition to the division of predator species into two areas, another notable difference between Kar and Ghosh [5], as well as the research by Chakraborty and Kar [26] and Chakraborty et al. [28], is the predation function used. Kar and Ghosh [5] adopted the HT I function, which is similar to that used by Chakraborty and Kar [26]. The model proposed by Kar and Ghosh [5] is stated as follows:
{ d r d t = g r r ( 1 r K α 1 ) m 1 ( r α 1 K u ( 1 α 1 ) K ) p r r U d u d t = g u u ( 1 u K ( 1 α 1 ) ) + m 1 ( r α 1 K u ( 1 α 1 ) K ) p u u U h u E u u d R d t = g R R ( 1 R α 2 r ) m 2 ( R α 2 r U α 2 u ) d U d t = g U U ( 1 U ( 1 α 2 ) u ) + m 2 ( R α 2 r U α 2 u ) .
Parameter α 1 is a constant proportion of the marine reserve area. Meanwhile, α 2 denotes a constant proportion of the carrying capacity of predator species. g R and g U are the intrinsic growth rates of predators. There is migration since predators are divided into two distinct areas (21). Then, Kar and Ghosh [5] assumed that the migration rate of prey species from the marine reserve to the harvesting area and vice versa is equal to m 1 . The same also applies to predator species with a value of m 2 .
The next model with four compartments to be reviewed is that of Pei et al. [17]. Interestingly, Pei et al. [17] modified the model proposed by Lv et al. [29] in (10). The original model formulated by Lv et al. [29] suggested that prey species have higher mobility than predators [17]. Pei et al. [17] expanded on this concept by stating that predator species consist of larger animals such as fish and exhibit even greater mobility compared with prey. In addition, the model highlighted that predator species are more dominant in being harvested [17], making them more susceptible to exploitation [34]. Based on these facts, the model modified by Pei et al. [17] is stated as follows:
{ d r d t = g r r ( 1 r K r ) m r r + m u u p r r R a r + r d u d t = g u u ( 1 u K u ) + m r r m u u p u u U a u + u h u E u u d R d t = n r r R α r + r d R R m R R + m U U d U d t = n u u U a u + u d U R + m R R m U U h U E U U .
The next group of models discussed are the two-compartment types. Among this group, Ibrahim [10] conducted research specifically focusing on a model with an implicit marine reserve. The model does not divide species into two compartments for marine reserve and harvesting area in this unique method. Instead, each species is represented by a single compartment. The marine reserve model is implicitly expressed as a percentage, indicating the proportion that could be harvested and the one designated for protection, all within the same compartment. Ibrahim [10] sheds more light on this novel model through comprehensive research and analysis of its implications.
{ d u d t = g u u ( 1 u K u ) s h u E u u d E u d t = s h u E u ( u e c ) .
The variable u in Equation (23) represents the size of prey species and is subject to a logistic growth function. Ibrahim [10] adopted a unique method in which prey species were not partitioned into two areas (marine reserve and harvesting). The model limited the harvestable portion of prey population to s u , with s indicating the proportion available for harvesting within the range of 0 < s 1 . Meanwhile, ( 1 s ) u is the amount of stock that cannot be harvested. In order to account for human harvesting activities, Ibrahim [10] assumed that the fishing effort E u ( t ) is represented by the combined population size of fishers and predators. This perspective treats fish harvesting as a form of predation, where human fishers and natural predators contribute to the pressure on the prey population. Integrating these elements into the model, Ibrahim [10] successfully captured the intricate interactions between prey, predators, and human harvesting practices within the ecosystem.
The next two-compartment model has an explicit marine reserve. This model was researched by Křivan and Jana [12] and Huang et al. [21]. The explicit marine reserve divides a species into two distinct compartments, one for protected areas and another for harvesting zones. Křivan and Jana [12] and Huang et al. [21] focused on a single-species model, which differs from the explicit marine reserve frameworks in the three- and four-compartment models. Specifically, Křivan and Jana [12] examined the model formulated by Kar and Matsuda [35]. The model formulated by Kar and Matsuda [35], cited in the research by Křivan and Jana [12], is stated as follows:
{ d u d t = g u u ( 1 u K u ) δ ( d u r u + d r u r ) E u u d r d t = g r r ( 1 r K r ) + δ ( d u r u + d r u r ) .
The parameter δ represents the degree of tendency of individuals to disperse or migrate. Meanwhile, d u r and d r u are the dispersal (migration) probability from the marine reserve to the harvesting area and vice versa. Using the model in (24), Křivan and Jana [12] examined the effects of dispersal modes between harvesting and marine reserve areas on maximum sustainable yield and balanced population abundance. Their research offers valuable insights into the role of dispersal in the management of harvesting and marine reserve strategies, contributing to more effective and sustainable resource conservation methods.
In contrast to the research by Křivan and Jana [12], Huang et al. [21] researched a distinct model that features parameters within specific intervals. The model researched by Huang et al. [21] is stated as follows.
{ d u d t = [ g u l , g u u ] u F u ( u ) + m u 2 u 2 + a 2 u h u E u u d r d t = [ g r l , g r u ] r F r ( r ) m u 2 u 2 + a 2 u h r E r r .
In Equation (25), the model incorporated parameter intervals for the growth rates of species in the marine reserve and harvesting area, expressed as [ g u l , g u u ] and [ g r l , g r u ] . The level 1 subscript denotes the area of the species, i.e., the marine reserve ( r ) and harvesting ( u ), while the level 2 subscript denotes the lower ( l ) and upper limits ( u ) of the growth rates. In addition, [ g u l , g u u ] and [ g r l , g r u ] are represented as exponent functions g u l 1 γ 1 g u u γ 1 and g r l 1 γ 2 g r u γ 2 for γ 1 , γ 2 [ 0 , 1 ] , respectively. The expressions [ g u l , g u u ] u F u ( u ) and [ g r l , g r u ] r F r ( r ) are used to describe the growth of a species in the harvesting and marine reserve areas, respectively. F u ( u ) and F r ( r ) are functions that fulfill the following criteria.
  • F u ( u ) < 0   F r ( r ) < 0 ;
  • There exists a constant, K u ,   K r > 0 , which is the carrying capacity of the species in the harvesting area and marine reserve, respectively; hence, F u ( K u ) = 0 and F r ( K r ) = 0 ;
  • { F u ( u ) > 0 u < K u F u ( u ) < 0 u > K u F r ( r ) > 0 r < K r F r ( r ) < 0 r > K r .
Based on Equation (25), migration occurs solely from the marine reserve to the harvesting area, and the migration rate follows the function m u 2 u 2 + a 2 . When u = a , then m 2 population of the marine reserve migrates to the harvesting area. Interestingly, Huang et al. [21] assumed that harvesting continued within the marine reserve area, as stated in Equation (25) and expressed in terms of h r E r r .
The last group of models discussed comprised five compartments, and was exclusively researched by Meng et al. [22]. This model includes density compartments for nutrients ( N ) , phytoplankton ( P ), zooplankton ( Z ), fish in the marine reserve ( F R ), and fish in the harvesting area ( F U ). Meng et al. [22] revealed the existence of both super-predators and top predators in this model, with fish serving as the top predator. The following is the model researched by Meng et al. [22].
{ d N d t = D ( n 0 N ) p N N P + c 1 d 1 P + c 2 d 2 Z d P d t = n N N P ( D + d 1 ) P c P 2 p P P Z d Z d t = n P P Z ( D + d 2 ) Z p Z F U Z F U p Z F R Z F R d F U d t = n Z F U Z F U m U F U + m R F R ( D + d 3 ) F U h F U E F U F U d F R d t = n Z F R Z F R + m U F U m R F R ( D + d 4 ) F R .
Meng et al. [22] also explored the concept of delay time. A detailed description of the parameters in (26) is shown in Table 2.
A Phylogenetic Tree diagram is shown in Figure 2, presenting the relationships between the 20 articles reviewed in this research. The diagram effectively captures the intersections and distinctions of each article, providing a summarized view of their content and graphical representations. In the diagram, a mathematical model of predator–prey dynamics with marine reserve and harvesting areas can be divided based on the number of compartments of each. Among the 20 articles, there were models with two, three, four, and five compartments, represented by 3, 10, 6, and 1 article(s), respectively, as shown in Figure 2. Furthermore, articles 10 and 6 were further divided based on similar characteristics, namely the predator–prey model base used in the form of the Lotka–Volterra [13,14,15,16,17,18,19,23,24,25,27,29] or Leslie–Gower [5,20,26,28] model. The model with two compartments was subdivided based on the presence of an implicit or explicit marine reserve. Articles with two, three, four, and five compartments were categorized until their distinguishing characteristics were obtained. Through this structured organization, the diagram effectively highlights the significant characteristics of each article and offers an insightful overview of the various predator–prey models employed in the context of marine reserves and harvesting areas.
Table 3 shows a comprehensive classification of the 20 articles based on various criteria, including the presence of an implicit or explicit marine reserve, the number of compartments, Lotka–Volterra or Leslie–Gower model, algebraic equation, stage structure, time delay, the harvesting function of prey and predators, and the predation function. For example, Khamis et al. [23] researched a model with an explicit marine reserve comprising three compartments, namely r , u , and U , representing prey species in marine reserves and harvesting areas and predators. The base model adopted by Khamis et al. [23] is Lotka–Volterra. It fails to examine the algebraic equation, stage structure, and time delay. Harvesting exclusively occurs on prey species with a harvesting function of h u E u u . Finally, the predation function used is the Holling Type (HT) I function. In order to delve further into the parameters used in Table 3, a detailed description is shown in Table 4.
Figure 3 shows the interconnectedness of the twenty articles classified into three distinct relations based on the connecting lines. The connections represent the relationship between two model, namely, Model A (published earlier) and Model B (published later). A solid line indicates that Model A is a development of Model B. Model A is also explicitly mentioned as a development of Model B. A dotted line signifies that Model B is a modification or development of Model A. However, Model B is not explicitly mentioned as a development of Model A. The article from Model A is cited as a reference in Model B. A dashed line shows that Model B only modifies Model A. In this case, the article from Model A is not listed as a reference in the article from Model B. For instance, the model in Chakraborty et al. [27] is a development of the model in Chakraborty et al. [25]. The line connecting the two articles is a dotted line, showing that Chakraborty et al. [27] do not explicitly mention it as a development of Chakraborty et al. [25]. Chakraborty et al. [25] is referenced in the article by Chakraborty et al. [27]. The modifications made to the model of Chakraborty et al. [25] involve adding a time-dependent capture effort variable equation. The model in Lv et al. [29] is solely a development of the model in Khamis [23]. This means that the model is not listed as a reference in the research of Lv et al. [29]. The connection between the research of Lv et al. [29] and Khamis [23] is represented with a dashed line. Lv et al. [29] expanded the Khamis [23] model by incorporating Holling Type II predation and harvesting functions into predators. Louartassi et al. [14] explicitly mentioned that their model was a development of the research of Sharma and Gupta [13]. Louartassi et al. [14] modified the capture effort variable into a non-constant effort variable based on the work of Sharma and Gupta [13].
E u ( u ) = α u ( t ) β u ( t ) 1 u d u d t   and   E U ( U ) = α U ( t ) β U ( t ) 1 U d U d t .
The analysis of the twenty models shown in Figure 2, Table 3, and Figure 3 presents numerous opportunities for future research and model development. All models described in the articles are formulated as systems of ordinary differential equations. However, various mathematical model options can be explored for future investigations. Potential directions for model development include a difference model [36,37,38], fractional difference [39,40], fractional differential [41,42], and partial differential [43,44,45]. Given that the data obtained in applications are discrete in time, the inclusion of a difference-equation-based model becomes particularly relevant and should be considered during development. Ghanbari dan Djilali [46] emphasized the importance of considering memory effects on mutual interactions between species when developing a differential or difference equation model. The intrinsic nature of memory effects in dynamical systems significantly influences ecological dynamics. In specific contexts, such as marine reserves, incorporating partial differential equations into model development can provide valuable insights into the dynamics of each population based on spatial or time factors. These developments can lead to a deterministic or stochastic model, as shown in [47], where ordinary differential equations were employed. The ultimate objective of such model enhancements is to provide more realistic representations, bringing models closer to real-world complexities and incorporating ecological and economic sustainability considerations. The advancement of these models can potentially address food-security issues effectively.
Srinivas et al. [24] developed a mathematical model for predator–prey dynamics in marine reserves and harvesting areas considering both immature and mature predator species. However, Srinivas et al. [24] did not account for dividing predator species into two-stage structures within each area. To address this limitation, a stage structure model was added to the predator–prey model for the reserve area. This incorporation was expected to provide an additional method to address food-security concerns related to fisheries’ resources. This means that two policies can be implemented to effectively tackle food-security challenges, namely, the establishment of marine reserves and the adoption of stage-selective harvesting. Stage-selective harvesting involves determining the most suitable stage class for harvesting, whether it is one or all of the stage classes. When all stage classes are harvested, identifying their thresholds is vital to ensure ecological and economic sustainability. For instance, in a fish population with two-stage classes, namely, juvenile and adult, targeting adult-stage fish for harvesting leads to a small or non-existent reproduction population, thereby impacting future regeneration. Conversely, focusing solely on juvenile-stage fish for harvesting might lead to a shortage of the adult-stage breeding population, further affecting future regeneration, making stage-selective harvesting crucial. The envisaged model, integrating stage-selective harvesting and marine reserves, explores whether this combination can effectively address food-security issues. This integrated method should take into account ecological and economic sustainability aspects.

3.3. Overview of the Methods of Each Model for Each Article and Methods That Can Be Used for Future Research

This section reviews the twenty articles, focusing on the methods to achieve their respective research objectives. The primary goals of most of this research were to develop mathematical models, test their biologically meaningful solution properties, and analyze their dynamic behavior. Several articles investigated biologically meaningful solution properties, including positivity [28] and boundedness [10,14,16,17,18,27,28,29]. The positivity test ensured that model solutions remained positive, as populations cannot have negative values. The boundedness test confirmed that model solutions remained finite, acknowledging the limited growth of the populations under research due to resource constraints [48]. These and other biologically meaningful solution tests were crucial in making the models more realistic in addressing ecology-related problems. However, some other biologically meaningful model solution properties, such as existence, uniqueness existence, and uniqueness [48]; nonnegativity [48]; persistence [49]; and permanence [49], were not explored in the twenty articles.
Several methods were employed across the investigations in researching the dynamics of the models. These methods included local stability analyses using Jacobian linearization and Routh–Hurwitz criterion [10,13,14,15,16,18,19,20,21,23,24,27,28,29], global stability analysis using the Lyapunov method [10,13,14,15,16,19,20,23,24,27,28,29] or geometric method [28], bionomic equilibrium [10,13,21,23,24,29], Hopf bifurcation analysis with singularity-induced bifurcation [25], the theory of normal forms and center manifold [27], or integral trace of the Jacobian matrix [10]. Additionally, sensitivity analyses [23,27] and examination of the existence of limit cycles using the Bendixson–Dulac criterion [10] were also used.
The Jacobian linearization method is commonly used for analyzing local stability. However, there is a different method for linearizing nonlinear models known as feedback linearization. This method primarily aims to simplify the original model into an equivalent form which is distinct from Jacobian linearization using the Taylor series. Feedback linearization encompasses two types, namely, approximation linearization and exact linearization. The latter considers higher-order terms, making it a global linearization applicable to all state and output spaces [50], while Jacobian linearization is valid only in specific regions around an equilibrium point [51]. Feedback linearization also ensures local and global stability [51], whereas Jacobian linearization can only guarantee local stability. Approximation feedback linearization was introduced by [52] in 1984 as an alternative to exact feedback linearization due to the limitations of exact linearization for some systems [53]. Cardoso and Schnitman [54] compared approximation feedback linearization with Jacobian linearization in the context of a pendulum system, and the approximation feedback linearization was found to be superior. Consequently, feedback linearization has found broader applications and has been successfully used in various fields, such as tracking problems, control of robot arms, helicopters, airplanes, medical equipment, etc. [50]. Its use in mathematical model for ecological and biological problems, particularly in predator–prey scenarios, still needs to be improved. Significant instances have been where approximation and exact linearization were applied to predator–prey problems. Singh and Gakkhar in 2013 [55] and Singh in 2016 [51] used these techniques in their research. dos Reis et al. [50] also employed exact feedback linearization to model population growth of the Aedes aegypti mosquito, which transmits the Dengue virus.
Another research goal is to determine the optimal control of the analyzed model using the Pontryagin maximum principle. Specifically, certain research studies focused on identifying the optimal harvesting strategies for the constructed model [13,17,19,21,23,24,26,28,29], while others investigated tax policies to determine the optimal tax for harvesting [10,18,27]. The time frames considered in the optimal control research varied, with some research efforts using finite time intervals [10,17,28] and others exploring infinite time horizons [13,18,19,21,23,24,26,27,29]. For sustainability research, infinite time frames are more suitable, as sustainability problems often span long-term periods. However, this does not exclude the possibility of periodic and continuous control evaluations, which are still valuable in practice, as control strategies may need to adapt over time. Control strategies that can be used long-term are needed as a solution in sustainability research. The research of optimal control in this research is limited to the control variables and models researched. The focus should be on the population density of each species to ensure biological significance, ensuring it remains above a positive threshold value. Previous works by Biswas et al. [56] and dos Reis et al. [50] using the Pontryagin maximum principle and Variable Neighborhood Search (VNS) methods serve as a valuable references to enhance the assessment of optimal control in constructed models.
In addition to the mentioned issues, several other ecological factors deserve consideration in developing realistic models. These factors include (i) source and sink subpopulation [7,8], (ii) different levels of intraspecific competition [57,58], and (iii) coupled growth function [59]. These elements play a critical role in determining optimal harvesting strategies and the choice between MSY or Quasi MSY as harvesting strategies, as shown in [60]. Incorporating industrial and engineering aspects is also important and can enhance the analysis of harvesting problems, as shown in [61].

4. Conclusions

In conclusion, this research systematically reviewed articles related to predator–prey mathematical models in marine reserves and harvesting areas within fishery ecosystems. The review process consisted of three stages, namely, collection, selection, and review. Data collection and selection were based on the PRISMA process diagram, which led to the identification of twenty relevant articles for the review stage. During this stage, each model and method used in the articles were thoroughly examined, and future research opportunities were discussed.
The review stage revealed numerous modeling opportunities for future research. All twenty articles relied on ordinary differential equations for their model, but there is potential for exploring alternative mathematical modeling options, such as difference equations, fractional differences, partial differentials, and fractional differentials, which may be more suitable for addressing specific problems. Introducing stage structure is another valuable option to enhance existing models, as individuals in different stage classes exhibit distinct biological properties and capabilities. Incorporating stage structure into predator–prey models that consider marine reserves can positively impact the resolution of food-security issues by providing more realistic models.
In analyzing a nonlinear dynamical model, such as a fishery model, determining the stability of a steady-state solution is regarded as one of the most important steps. Although most research primarily utilizes the Jacobian linearization method, further research needs to explore other linearization methods, such as approximation and exact feedback linearization. These methods are rarely used for linearizing nonlinear models in biological and ecological mathematical research. Moreover, optimal control research on the models from the twenty relevant articles is currently limited to the constructed and applied control model. To enhance biological relevance, adding an additional condition to the optimal control problem is recommended, requiring population density to be greater than or equal to a sufficiently large positive threshold value. The results from this research and the suggested development provide valuable insights and knowledge in Mathematical Biology and Ecology. Additionally, they offer alternative suggestions to policymakers for effective management of the natural resources of fisheries, ensuring sustainability and profitability.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/su151612291/s1, bit.ly/PRISMA-Materials.

Author Contributions

Conceptualization, methodology, formal analysis, writing—original draft preparation, A.H. and A.K.S.; writing—review and editing, investigation, interpretation, A.H.; visualization, M.H.A.B.; supervision, E.R.; project administration, A.H. and A.K.S.; funding acquisition, A.K.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Universitas Padjadjaran with grant number 1549/UN6.3.1/PT.00/2023 through the Padjadjaran Doctoral Program Scholarship scheme.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are grateful to the Universitas Padjadjaran for providing the Padjadjaran Doctoral Program Scholarship and Academic Leadership Grant 2023.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Mcfarlane, I.; Black, K.; Madonia, K.; Jensen, J.; Kollodge, R.; Daldin, J.; Jayaram, T.; Ratcliffe, L.; Trautwein, C.; Baker, D.; et al. UNFPA State of World Population: 8 Billion Lives, INFINITE Possibilities the Case for Rights and Choices; UNPFA: New York, NY, USA, 2023; Available online: https://www.unfpa.org/sites/default/files/swop23/SWOP2023-ENGLISH-230329-web.pdf (accessed on 20 July 2023).
  2. FAO; IFAD; UNICEF; WFP; WHO. Food Security and Nutrition in the World the State of Building Climate Resilience for Food Security and Nutrition; FAO: Rome, Italy, 2018; ISBN 9789251305713. [Google Scholar]
  3. FAO; IFAD; UNICEF; WFP; WHO. The State of Food Security and Nutrition in the World 2022. Repurposing Food and Agricultural Policies to Make Healthy Diets More Affordable; FAO: Rome, Italy, 2022; ISBN 978-92-5-136499-4. [Google Scholar]
  4. Ghosh, B.; Pal, D.; Kar, T.K.; Valverde, J.C. Biological Conservation through Marine Protected Areas in the Presence of Alternative Stable States. Math. Biosci. 2017, 286, 49–57. [Google Scholar] [CrossRef] [PubMed]
  5. Kar, T.K.; Ghosh, B. Sustainability and Economic Consequences of Creating Marine Protected Areas in Multispecies Multiactivity Context. J. Theor. Biol. 2013, 318, 81–90. [Google Scholar] [CrossRef] [PubMed]
  6. FAO. In Brief to the State of World Fisheries and Aquaculture 2022; FAO: Rome, Italy, 2022; ISBN 978-92-5-136367-6. [Google Scholar]
  7. Supriatna, A.K.; Possingham, H.P. Optimal Harvesting for a Predator-Prey Metapopulation. Bull. Math. Biol. 1998, 60, 49–65. [Google Scholar] [CrossRef] [Green Version]
  8. Supriatna, A.K.; Possingham, H.P. Harvesting a Two-Patch Predator-Prey Metapopulation. Nat. Resour. Model. 1999, 12, 481–498. [Google Scholar] [CrossRef]
  9. Supriatna, A.K.; Tuck, G.N.; Possingham, H. On the Exploitation of a Two-Patch Metapopulation with Delayed Juvenile Recruitment and Predation. J. Indones. Math. Soc. 2003, 8, 139–150. [Google Scholar]
  10. Ibrahim, M. Optimal Harvesting of a Predator-Prey System with Marine Reserve. Sci. Afr. 2021, 14, e01048. [Google Scholar] [CrossRef]
  11. Bellier, E.; Sæther, B.E.; Engen, S. Sustainable Strategies for Harvesting Predators and Prey in a Fluctuating Environment. Ecol. Modell. 2021, 440, 109350. [Google Scholar] [CrossRef]
  12. Křivan, V.; Jana, D. Effects of Animal Dispersal on Harvesting with Protected Areas. J. Theor. Biol. 2015, 364, 131–138. [Google Scholar] [CrossRef]
  13. Sharma, A.; Gupta, B. Harvesting Model for Fishery Resource with Reserve Area and Bird Predator. J. Mar. Biol. 2014, 2014, 218451. [Google Scholar] [CrossRef] [Green Version]
  14. Louartassi, Y.; El Alami, J.; Elalami, N. Harvesting Model for Fishery Resource with Reserve Area of Bird Predator and Modified Effort Function. Malaya J. Mat. 2017, 5, 660–666. [Google Scholar] [CrossRef] [Green Version]
  15. Agnihotri, K.; Nayyer, S. Stability Analysis of a Predator (Bird)–Prey (Fish) Harvesting Model in the Reserved and Unreserved Area. Malaya J. Mat. 2018, 6, 678–684. [Google Scholar] [CrossRef] [PubMed]
  16. Louartassi, Y.; Alla, A.; Hattaf, K.; Nabil, A. Dynamics of a Predator–Prey Model with Harvesting and Reserve Area for Prey in the Presence of Competition and Toxicity. J. Appl. Math. Comput. 2019, 59, 305–321. [Google Scholar] [CrossRef]
  17. Pei, Y.; Chen, M.; Liang, X.; Li, C. Model-Based on Fishery Management Systems with Selective Harvest Policies. Math. Comput. Simul. 2019, 156, 377–395. [Google Scholar] [CrossRef]
  18. Huang, L.; Cai, D.; Liu, W. Optimal Tax Policy of a One-Predator–Two-Prey System with a Marine Protected Area. Math. Methods Appl. Sci. 2021, 44, 6876–6895. [Google Scholar] [CrossRef]
  19. Huo, H.F.; Jiang, H.M.; Meng, X.Y. A Dynamic Model for Fishery Resource with Reserve Area and Taxation. J. Appl. Math. 2012, 2012, 794719. [Google Scholar] [CrossRef] [Green Version]
  20. Zhang, Y.; Li, J.; Jie, Y.; Yan, X. Optimal Taxation Policy for a Prey-Predator Fishery Model with Reserves. Pac. J. Optim. 2015, 11, 137–155. [Google Scholar]
  21. Huang, L.; Cai, D.; Liu, W. Optimal Harvesting of an Abstract Population Model with Interval Biological Parameters. Adv. Differ. Equ. 2020, 2020, 285. [Google Scholar] [CrossRef]
  22. Meng, X.Y.; Wu, Y.Q.; Li, J. Bifurcation Analysis of a Singular Nutrient-Plankton-Fish Model with Taxation, Protected Zone and Multiple Delays. Numer. Algebra Control. Optim. 2020, 10, 391. [Google Scholar] [CrossRef] [Green Version]
  23. Khamis, S.A.; Tchuenche, J.M.; Lukka, M.; Heiliö, M. Dynamics of Fisheries with Prey Reserve and Harvesting. Int. J. Comput. Math. 2011, 88, 1776–1802. [Google Scholar] [CrossRef]
  24. Srinivas, M.N.; Srinivas, M.A.S.; Das, K.; Gazi, N.H. Prey Predator Fishery Model with Stage Structure in Two Patchy Marine Aquatic Environment. Appl. Math. 2011, 2, 1405–1416. [Google Scholar] [CrossRef] [Green Version]
  25. Chakraborty, K.; Chakraborty, M.; Kar, T.K. Bifurcation and Control of a Bioeconomic Model of a Prey-Predator System with a Time Delay. Nonlinear Anal. Hybrid Syst. 2011, 5, 613–625. [Google Scholar] [CrossRef]
  26. Chakraborty, K.; Kar, T.K. Economic Perspective of Marine Reserves in Fisheries: A Bioeconomic Model. Math. Biosci. 2012, 240, 212–222. [Google Scholar] [CrossRef] [PubMed]
  27. Chakraborty, K.; Jana, S.; Kar, T.K. Effort Dynamics of a Delay-Induced Prey-Predator System with Reserve. Nonlinear Dyn. 2012, 70, 1805–1829. [Google Scholar] [CrossRef]
  28. Chakraborty, K.; Das, K.; Kar, T.K. An Ecological Perspective on Marine Reserves in Prey-Predator Dynamics. J. Biol. Phys. 2013, 39, 749–776. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Lv, Y.; Yuan, R.; Pei, Y. A Prey-Predator Model with Harvesting for Fishery Resource with Reserve Area. Appl. Math. Model. 2013, 37, 3048–3062. [Google Scholar] [CrossRef]
  30. Page, M.J.; McKenzie, J.E.; Bossuyt, P.M.; Boutron, I.; Hoffmann, T.C.; Mulrow, C.D.; Shamseer, L.; Tetzlaff, J.M.; Akl, E.A.; Brennan, S.E.; et al. The PRISMA 2020 Statement: An Updated Guideline for Reporting Systematic Reviews. J. Clin. Epidemiol. 2021, 134, 178–189. [Google Scholar] [CrossRef] [PubMed]
  31. Berardo, C.; Geritz, S.; Gyllenberg, M.; Raoul, G. Interactions between Different Predator–Prey States: A Method for the Derivation of the Functional and Numerical Response. J. Math. Biol. 2020, 80, 2431–2468. [Google Scholar] [CrossRef]
  32. Huo, H.F.; Wang, X.; Castillo-Chavez, C. Dynamics of a Stage-Structured Leslie-Gower Predator-Prey Model. Math. Probl. Eng. 2011, 2011, 149341. [Google Scholar] [CrossRef] [Green Version]
  33. Idels, L.V.; Wang, M. Harvesting Fisheries Management Strategies with Modified Effort Function. Int. J. Model. Identif. Control 2008, 3, 83. [Google Scholar] [CrossRef] [Green Version]
  34. Pauly, D.; Christensen, V.; Dalsgaard, J.; Froese, R.; Torres, F. Fishing down Marine Food Webs. Science 1998, 279, 860–863. [Google Scholar] [CrossRef]
  35. Kar, T.K.; Matsuda, H. A Bioeconomic Model of a Single-Species Fishery with a Marine Reserve. J. Environ. Manag. 2008, 86, 171–180. [Google Scholar] [CrossRef] [PubMed]
  36. Hong, B.; Zhang, C. Neimark–Sacker Bifurcation of a Discrete-Time Predator–Prey Model with Prey Refuge Effect. Mathematics 2023, 11, 1399. [Google Scholar] [CrossRef]
  37. Khaliq, A.; Ibrahim, T.F.; Alotaibi, A.M.; Shoaib, M.; El-moneam, M.A. Dynamical Analysis of Discrete-Time Two-Predators One-Prey Lotka—Volterra Model. Mathematics 2022, 10, 4015. [Google Scholar] [CrossRef]
  38. Du, X.; Han, X.; Lei, C. Behavior Analysis of a Class of Discrete-Time Dynamical System with Capture Rate. Mathematics 2022, 10, 2410. [Google Scholar] [CrossRef]
  39. Saadeh, R.; Abbes, A.; Al-Husban, A.; Ouannas, A.; Grassi, G. The Fractional Discrete Predator–Prey Model: Chaos, Control and Synchronization. Fractal Fract. 2023, 7, 120. [Google Scholar] [CrossRef]
  40. Elsadany, A.A.; Matouk, A.E. Dynamical Behaviors of Fractional-Order Lotka–Volterra Predator–Prey Model and Its Discretization. J. Appl. Math. Comput. 2015, 49, 269–283. [Google Scholar] [CrossRef]
  41. Areshi, M.; Seadawy, A.R.; Ali, A.; Alharbi, A.F.; Aljohani, A.F. Analytical Solutions of the Predator–Prey Model with Fractional Derivative Order via Applications of Three Modified Mathematical Methods. Fractal Fract. 2023, 7, 128. [Google Scholar] [CrossRef]
  42. Bentout, S.; Djilali, S.; Kumar, S. Mathematical Analysis of the Influence of Prey Escaping from Prey Herd on Three Species Fractional Predator-Prey Interaction Model. Phys. A Stat. Mech. Its Appl. 2021, 572, 125840. [Google Scholar] [CrossRef]
  43. Xie, Y.; Zhao, J.; Yang, R. Stability Analysis and Hopf Bifurcation of a Delayed Diffusive Predator–Prey Model with a Strong Allee Effect on the Prey and the Effect of Fear on the Predator. Mathematics 2023, 11, 1996. [Google Scholar] [CrossRef]
  44. Wang, M. Diffusion-Induced Instability of the Periodic Solutions in a Reaction-Diffusion Predator-Prey Model with Dormancy of Predators. Mathematics 2023, 11, 1875. [Google Scholar] [CrossRef]
  45. Jin, D.; Yang, R. Hopf Bifurcation in a Predator-Prey Model with Memory Effect and Intra-Species Competition in Predator. J. Appl. Anal. Comput. 2023, 13, 1321–1335. [Google Scholar] [CrossRef] [PubMed]
  46. Ghanbari, B.; Djilali, S. Mathematical Analysis of a Fractional-Order Predator-Prey Model with Prey Social Behavior and Infection Developed in Predator Population. Chaos Solitons Fractals 2020, 138, 109960. [Google Scholar] [CrossRef]
  47. Shao, Y.; Kong, W. A Predator–Prey Model with Beddington–DeAngelis Functional Response and Multiple Delays in Deterministic and Stochastic Environments. Mathematics 2022, 10, 3378. [Google Scholar] [CrossRef]
  48. Beay, L.K.; Suryanto, A.; Darti, I. Trisilowati Hopf Bifurcation and Stability Analysis of the Rosenzweig-MacArthur Predator-Prey Model with Stage-Structure in Prey. Math. Biosci. Eng. 2020, 17, 4080–4097. [Google Scholar] [CrossRef]
  49. Ackleh, A.S.; Hossain, M.I.; Veprauskas, A.; Zhang, A. Persistence and Stability Analysis of Discrete-Time Predator–Prey Models: A Study of Population and Evolutionary Dynamics. J. Differ. Equ. Appl. 2019, 25, 1568–1603. [Google Scholar] [CrossRef]
  50. dos Reis, C.A.; Florentino, H.d.O.; Cólon, D.; Rosa, S.R.F.; Cantane, D.R. An Approach of the Exact Linearization Techniques to Analysis of Population Dynamics of the Mosquito Aedes Aegypti. Math. Biosci. 2018, 299, 51–57. [Google Scholar] [CrossRef] [PubMed]
  51. Singh, A. Stabilization of Prey Predator Model via Feedback Control. In Applied Analysis in Biological and Physical Sciences; Springer: New Delhi, India, 2016; Volume 186. [Google Scholar]
  52. Krener, A.J. Approximate Linearization by State Feedback and Coordinate Change. Syst. Control Lett. 1984, 5, 181–185. [Google Scholar] [CrossRef]
  53. Krener, A.J. Feedback Linearization. In Mathematical Control Theory; Springer: New York, NY, USA, 1999; pp. 66–98. [Google Scholar]
  54. Cardoso, G.S.; Schnitman, L. Analysis of Exact Linearization and Aproximate Feedback Linearization Techniques. Math. Probl. Eng. 2011, 2011, 205939. [Google Scholar] [CrossRef]
  55. Singh, A.; Gakkhar, S. Stabilization of Modified Leslie-Gower Prey-Predator Model. Differ. Equ. Dyn. Syst. 2014, 22, 239–249. [Google Scholar] [CrossRef]
  56. Biswas, M.H.A.; Haque, M.M.; Mallick, U.K. Optimal Control Strategy for the Immunotherapeutic Treatment of HIV Infection with State Constraint. Optim. Control Appl. Methods 2019, 40, 807–818. [Google Scholar] [CrossRef]
  57. Supriatna, A.K. Maximum Sustainable Yield for Marine Metapopulation Governed by Coupled Generalised Logistic Equations. J. Sustain. Sci. Manag. 2012, 7, 201–206. [Google Scholar]
  58. Husniah, H.; Anggriani, N.; Supriatna, A.K. System Dynamics Approach in Managing Complex Biological Resources. ARPN J. Eng. Appl. Sci. 2015, 10, 1685–1690. [Google Scholar]
  59. Supriatna, A.K.; Husniah, H. Sustainable Harvesting Strategy for Natural Resources Having a Coupled Gompertz Production Function. In Proceedings of the 3rd International Congress on Interdisciplinary Behavior and Social Sciences, ICIBSoS, Bali, Indonesia, 1–2 November 2014; CRC Press: London, UK, 2015. [Google Scholar]
  60. Husniah, H.; Supriatna, A.K. Marine Biological Metapopulation with Coupled Logistic Growth Functions: The MSY and Quasi MSY. AIP Conf. Proc. 2014, 1587, 51–56. [Google Scholar]
  61. Husniah, H.; Supriatna, A.K. Optimal Number of Fishing Fleet for a Sustainable Fishery Industry. In Proceedings of the 2015 International Conference on Technology, Informatics, Management, Engineering and Environment, TIME-E 2015, Samosir Island, Indonesia, 7–9 September 2015; pp. 18–24. [Google Scholar]
Figure 1. Diagram of PRISMA for the article collection and selection stage.
Figure 1. Diagram of PRISMA for the article collection and selection stage.
Sustainability 15 12291 g001
Figure 2. Phylogenetic Tree diagram mathematical model of predators–prey in marine reserves and harvesting areas [5,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].
Figure 2. Phylogenetic Tree diagram mathematical model of predators–prey in marine reserves and harvesting areas [5,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].
Sustainability 15 12291 g002
Figure 3. Diagram of connectedness between each article’s mathematical model of predator–prey in marine reserves and harvesting areas [5,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].
Figure 3. Diagram of connectedness between each article’s mathematical model of predator–prey in marine reserves and harvesting areas [5,10,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].
Sustainability 15 12291 g003
Table 1. A summary of the article collection results for the four databases, namely Scopus, Science Direct, Dimensions, and Web of Science.
Table 1. A summary of the article collection results for the four databases, namely Scopus, Science Direct, Dimensions, and Web of Science.
KeywordsThe Amount of Data from (Articles)Total
ScopusScience DirectDimensionsWeb of Science
(“harvest” OR “harvesting”) AND (“predator” OR “prey”) AND (“model” OR “system”)106318124811712843
Table 2. Description of the parameters in (26).
Table 2. Description of the parameters in (26).
Parameter(s)Description
D The rate of washout for nutrients.
n 0 The rate of constant input.
p N and p P The rate of predation of phytoplankton on nutrients and zooplankton on phytoplankton.
c i The decomposition rate of phytoplankton ( i = 1 ) and zooplankton ( i = 2 ) mortality.
d i The rate of natural mortality of phytoplankton ( i = 1 ), zooplankton ( i = 2 ), fish in the harvesting area ( i = 3 ) , and fish in the marine reserve ( i = 4 ).
n N   and n P The conversion rate of predation on nutrients and phytoplankton, respectively.
p Z F U and p Z F R The predation rate of zooplankton by fish in the harvesting and marine reserve area, respectively.
n Z F U   and n Z F R The conversion rate of fish predation within the harvesting and marine reserve area on zooplankton, respectively.
m U and m R The rate of fish migration from the marine reserve to the harvesting area and vice versa, respectively.
c The level of interspecific competition of phytoplankton.
h U The catchability of fish in the harvesting area.
E U The constant harvesting effort of the species within the harvesting area.
Table 3. Classification of the twenty articles used as review material.
Table 3. Classification of the twenty articles used as review material.
(Authors, Year, Citation Number)Implicit/
Explicit Marine Reserve Area
CompartmentsLeslie–Gower (LG)/Lotka–Volterra (LV)Algebraic EquationStage StructureTime DelayHarvesting
Function
Functional Response
PreyPredator
(Khamis et al., 2011, [23])Explicit r , u , U LV××× h u E u u ×HT 1 I
(Srinivas et al., 2011, [24]) Explicit r , u ,   I , M LV×× h u E u u h M E M M HT I
(Chakraborty et al., 2011, [25])Explicit r , u , U LV× h u E u u ×HT II
(Chakraborty and Kar, 2012, [26]) Explicit r , u , U LG××× h u E u u ×HT II
(Chakraborty et al., 2012, [27])Explicit r , u , U , E u LV×× h u E u u ×HT II
(Huo et al., 2012, [19])Explicit r , u , E u LV××× h u E u u ××
(Chakraborty et al., 2013, [28])Explicit r , u , U LG××× h u E u u γ u E u + v u u ×HT II
(Kar and Ghosh, 2013, [5])Explicit r , u , R , U LG××× h u E u u ×HT I
(Lv et al., 2013, [29])Explicit r , u , U LV××× h u E u u h U E U U HT II
(Sharma and Gupta, 2014, [13])Explicit r , u , U LV××× h u E u u h U E U U HT I
(Křivan and Jana, 2015, [12])Explicit r , u LV××× E u u ××
(Zhang et al., 2015, [20])Explicit r , u , U ,   E u LG××× h u E u u ×HT I
(Louartassi et al., 2017, [14]) Explicit r , u , U LV××× h u E u u h U E U U HT I
(Agnihotri and Nayyer, 2018, [15])Explicit r , u , U LV××× h u E u u h U E U U HT I
(Louartassi et al., 2019, [16])Explicit r , u , U LV××× h u E u u h U E U U HT II
(Pei et al., 2019, [17]) Explicit r , u , R , U LV××× h u E u u h U E U U HT II
(Huang et al., 2020, [21])Explicit r , u LV××× h u E u u ××
(Meng et al., 2020, [22])Explicit N , P , Z , F U , F R LV××× h F U E F U F U HT I
(Huang et al., 2021, [18])Explicit r , u , U , E u LV××× h u E u u ×HT II
(Ibrahim, 2021, [10])Implicit u , E u LV××× s h u E u u ××
1 HT: Holling Type.
Table 4. The descriptions of the parameters are in Table 3.
Table 4. The descriptions of the parameters are in Table 3.
Parameter(s)Description
u Prey species in the harvesting area.
r Prey species in the marine reserve area.
U Predator species in the harvesting area.
R Predator species in the marine reserve area.
I Immature predator.
M Mature predator.
N Nutrients.
P Phytoplankton.
Z Zooplankton.
F R Fish in the marine reserve area.
F U Fish in the harvesting area.
h u ,   h U , h M ,   and   h F U The catchability coefficients are u , U , M , and F U , respectively.
E u ,   E U ,   E M ,   and   E F U Effort harvesting on u , U , M , and F U , respectively.
s Fraction of implicit marine reserve.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Hasibuan, A.; Supriatna, A.K.; Rusyaman, E.; Biswas, M.H.A. Harvested Predator–Prey Models Considering Marine Reserve Areas: Systematic Literature Review. Sustainability 2023, 15, 12291. https://doi.org/10.3390/su151612291

AMA Style

Hasibuan A, Supriatna AK, Rusyaman E, Biswas MHA. Harvested Predator–Prey Models Considering Marine Reserve Areas: Systematic Literature Review. Sustainability. 2023; 15(16):12291. https://doi.org/10.3390/su151612291

Chicago/Turabian Style

Hasibuan, Arjun, Asep Kuswandi Supriatna, Endang Rusyaman, and Md. Haider Ali Biswas. 2023. "Harvested Predator–Prey Models Considering Marine Reserve Areas: Systematic Literature Review" Sustainability 15, no. 16: 12291. https://doi.org/10.3390/su151612291

APA Style

Hasibuan, A., Supriatna, A. K., Rusyaman, E., & Biswas, M. H. A. (2023). Harvested Predator–Prey Models Considering Marine Reserve Areas: Systematic Literature Review. Sustainability, 15(16), 12291. https://doi.org/10.3390/su151612291

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop