^{1}Doctoral Program of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Padjadjaran, Sumedang, Indonesia^{2}Department of Mathematics, Faculty of Mathematics and Natural Sciences, Universitas Padjadjaran, Sumedang, Indonesia

The predator-prey model has been extensively studied, but only studies models in a certain environment, where all parameters and initial values involved in the model are assumed to be certain. In real practice, some parameters and initial values are often uncertain. To overcome this uncertainty problem, a model can be made by using a fuzzy theoretical approach. In this paper, we develop a numerical scheme for solving two predator-prey models with a Holling type II functional response by considering fuzzy parameters and initial populations. The behavior of the model was studied qualitatively using the 5th order Runge-Kutta method of which was modified for the fuzzy system using the Zadeh extension principle. The numerical simulation results show that, when the initial populations of prey and predators are fuzzy, the behavior of the fuzzy model would be qualitatively the same as the crisp model. Finally, we conclude that the resulting fuzzy behavior represents a generalization of crisp behavior. This gives more realistic results since the solution is obtained by explicitly considering the problem of uncertainty.

## 1. Introduction

The predator-prey model is a model of the interaction between two species expressed in the form of a system of differential equations that describes the dynamic relationship between prey and predators [1]. This model was first introduced by Lotka and Volterra, so it is known as the Lotka-Volterra predator-prey model. In this model, the dynamic behavior of a simple predator-prey model is studied. Various applications of the predator-prey system, such as those in Supriatna and Possingham [2, 3] and several other modifications of the predator-prey model have been made by incorporating additional biological processes into the classic Lotka-Volterra predator-prey equation, including functional response modifications [4].

In ecological systems, the degree of predation depends on the functional response. A functional response considers the number of prey that the predator has successfully consumed per unit time. It is also introduced to describe changes in the rate of prey consumption by predators when prey density varies [4]. The most common and well-known functional response is the type-II Holling functional response. The Holling type II functional response describes the increasing rate of predator's consumption when the density of prey is low. Meanwhile, when the prey density is high, the predator's consumption is constant. In this case, it represents a phenomenon that the predator takes very little time to find the prey, and when the prey consumption rate reaches the highest level, the predator becomes easily full.

In recent years, various studies on predator-prey models with type-II Holling functional responses have been carried out, including in in Dawes and Souze [4], Jana and Kar [5], Ma et al. [6]. Those studies consider a predator-prey model in a definite environment, where all parameters affecting population size and initial values involved in the model are assumed to be crisp. However, in reality, each parameter and initial value is uncertain, unclear, or incomplete. This uncertainty is caused by inaccuracies made during the process of observation, measurement, experimentation, and so on. To overcome this problem, a model can be made using different approaches such as the stochastic approach, the fuzzy approach and the fuzzy-stochastic approach. A crisp ODE system could be more suitable to be converted into a fuzzy differential equation system whenever the parameters or the initial values are uncertain and have a degree of perseptional values.

In recent decades, the application of fuzzy theory has been widely used as a very effective tool in mathematical modeling to solve real problems that take into account uncertainty. In this approach, uncertain variables and parameters are represented by intervals and fuzzy numbers. In the study of fuzzy differential equations, the term fuzzy differential equations can be in the form of differential equations with fuzzy coefficients, differential equations with fuzzy initial values or fuzzy boundary values [7–12]. The stability of the fuzzy dynamic system in a dynamic population is studied through fuzzy differential equation and fuzzy initial value problem [13]. Various numerical solutions for systems of fuzzy equations have also been introduced in Ahmad and Hasan [10], Jayakumar et al. [14], Nayak and Chakraverty [15], Behroozpoor et al. [16], Tapaswini and Chakraverty [17], and Tapaswini and Chakraverty [18].

The fuzzy predator-prey model was first introduced in da Silva Peixoto et al. [19], where a classic deterministic predator-prey model was formulated using a fuzzy rule-based system. The development of fuzzy differential equations has resulted in new discoveries of fuzzy predator-prey models, including those made by Ahmad and Hasan [10], Pandit and Singh [20], Ak and Oru [21], Ahmad and De Baets [22], Narayanamoorthy et al. [23], Omar et al. [24], and Pal et al. [25]. The authors in Ahmad and Hasan [10], Ahmad and De Baets [22], and Omar et al. [24] used the Euler and 4th order Runge-Kutta method through the principles of Zadeh extension. While the authors in Ak and Oru [21], used the concept of generalized fuzzy derivatives. Other authors in Pandit and Singh [20], used Hukuhara derivative. Moreover, the authors in Narayanamoorthy et al. [23], used the fractional modified Euler method. On the biological perspective, there are some authors who have studied fuzzy predator-prey models with functional responses such as [20, 23, 24, 26]. They all studied a predator-prey model with fuzzy initial conditions. Fuzzy predator-prey models with functional responses have also been studied by Pal et al. [27], Yu et al. [28], Pal et al. [29], Meng and Wu [30], Mahata et al. [31], and Pal et al. [32], who presented fuzzy predator-prey harvesting models. Their work studied two species of predator-prey harvesting models by considering fuzzy parameters. Among those work, the authors in Mallak et al. [26], studied a fuzzy predator-prey model with an arctan functional response using the Hukuhara derivative approach, to describe the satiation predator's consumption.

Our research discusses predator-prey models with Holling type II functional responses by considering fuzzy parameters and fuzzy initial populations. The motivation is that we would like to see how different is the dynamics of the sytems compared to their counterpart crisp predator-prey systems. To proceed we will present some preliminaries regarding the fuzzy number theory and fuzzy differential equation background in Sections 2 and 3 followed by the 5th order Runge-Kutta numerical scheme for fuzzy differential Section 4. In Section 5, we discuss the equilibria and their stabilities condition for two predator-prey models with fuzzy parameters and fuzzy initial values followed by the applications of the 5th order Runge-Kutta numerical scheme to those predator prey models in Section 6. Finally some discussion and conclusion are presented in Sections 7 and 8, respectivelly.

## 2. Preliminaries

Some of the basic concepts used in this paper, such as fuzzy number, the α-level of the fuzzy number, and the Zadeh extension principle, will be introduced in this section.

### 2.1. Fuzzy theory

**Definition 2.1** [33]. *Let U be a non-empty set, and A is a subset of U. The characteristic function of A is given by*

*for each* *x* ∈ *U*.

**Definition 2.2** [33]. *A fuzzy subset F of the non-empty set U is defined by a function* φ_{F}:*U* → [0, 1]*, which is called the membership function of F*.

**Definition 2.3 (α-level)** [33]. α*-level of the fuzzy subset A of U is the classical set* [*A*]^{α} *defined by*

Support of A is $supp\text{}A=\left\{x\in U:{\phi}_{A}(x)0\right\}={\left[A\text{}\right]}^{0}$.

Core of A is $core\text{}A=\left\{x\in U:{\phi}_{A}(x)=1\right\}={\left[A\text{}\right]}^{1}$.

**Definition 2.4 (Fuzzy Number)** [7]. *A fuzzy subset A is called a fuzzy number if the defined universal set is the set of all real numbers* ℝ *and satisfies the following conditions:*

(i) *All* α*-level A is not empty for* 0 ≤ α ≤ 1

(ii) *All* α*-levels of A are open intervals of* ℝ

(iii) *supp A* = {*x* ∈ ℝ: φ_{A}(*x*) > 0} *is bounded*.

The set of all fuzzy numbers is denoted by ${F}(\mathbb{R})$, and the α-level of the fuzzy number A is denoted by ${\left[A\right]}^{\alpha}=[{a}_{1}^{\alpha},$ ${a}_{2}^{\alpha}]$.

**Definition 2.5** [33]. *A fuzzy number A is called a triangular fuzzy number if its membership function has the following equation:*

*and* α*-level of A is*, [*A*]^{α} = [*a* + α (*b* − *a*), *c* − α (*c* − *b*)], *for one* α ∈ [0, 1].

### 2.2. Zadeh extension principle

Zadeh's extension principle is one of the basic ideas that encourage the expansion of non-fuzzy mathematical concepts to become fuzzy. This method was proposed by Zadeh to extend the concept from classical set theory to fuzzy set theory.

**Definition 2.6 (Zadeh Extension Principle)** [7]. *Let X and Y being two universal sets and* *f*:*X* → *Z* *are classical functions. The extension of f is a function* $\widehat{f}(A)\in {F}(Z)$, $A\in {F}(X)$ *such that*

*where* *f*^{−1} (*z*) = {*x*|*f* (*x*) = *z*}.

**Theorem 2.1** [7]. *if* *f*: ℝ ^{n} → ℝ^{n} *is a continuous function, then* $\widehat{f}:{F}({\mathbb{R}}^{n})\to {F}({\mathbb{R}}^{n})$ *is well-defined, continuous, and*

*for each* α ∈ [0, 1].

**Definition 2.7** [33]. *Suppose* *f*:ℝ × ℝ → ℝ *is a continuous function. If A and B are two fuzzy numbers, then the extension* $\widehat{f\text{}}$*via A and B, is a fuzzy subset* $\widehat{f\text{}}(A,B)$ *of* ℝ *with the membership function given by*:

*where* ${f}^{-1}\left(z\right)=\left\{\left(x,y\right)|f\left(x,y\right)=z\right\}$

*and* ${\left[\widehat{f}\left(A,B\right)\right]}^{\alpha}=f({[A]}^{\alpha},{[B]}^{\alpha})=\left\{f(x,y)|x\in \left[{a}_{1}^{\alpha},{a}_{2}^{\alpha}\right],y\in [{b}_{1}^{\alpha},{b}_{2}^{\alpha}]\hspace{0.17em}\right\}$.

## 3. Fuzzy differential equation

The initial value problem is given to be

where *f* is continuous and ${x}_{0}\in {\mathbb{R}}^{n}$. Suppose the initial condition *x*_{0} is uncertain and is modeled by a fuzzy set, then the problem (1) converted into a fuzzy differential equation

where $f:\left[{t}_{0},\text{}T\right]\times {F}({\mathbb{R}}^{n})\to {F}({\mathbb{R}}^{n})$, ${X}_{0}\in {F}({\mathbb{R}}^{n})$.

Suppose also *L*_{t}(*x*_{0}) = *x*(*t, x*_{0}) is the solution to the problem (1), then by applying extension principle for *L*_{t}(*x*_{0}) = *x*(*t, x*_{0}) obtained ${\widehat{L}}_{t}({X}_{0})=X(t,{X}_{0})$, which is the solution of the fuzzy problem (2).

**Definition 3.1 (Equilibrium Point)** [13]. *A fuzzy number* $\overline{X}\in {F}({\mathbb{R}}^{n})$ *is the equilibrium point of* (2) *if*

*or equivalent to*

**Theorem 3.1** [13]. *If* $\overline{x}$ *is the equilibrium point of the classical system* (1), *then* ${\chi}_{\left[\overline{x}\right]}$ *is the equilibrium point of the fuzzy system* (2) *where* ${\chi}_{\left[\overline{x}\right]\text{}}$*is a characteristic function of* $\overline{x}$.

**Theorem 3.2** [13]. *Suppose* $\overline{x}$ *is the equilibrium point of the deterministic initial value problem* (2), *then*

(a) $\overline{x}\text{}$*stable if and only if* ${\chi}_{\left[\overline{x}\right]}$ *is stable for the fuzzy initial value problem* (2).

(b) $\overline{x}\text{}$*is asymptotic table, if and only if* ${\chi}_{\left[\overline{x}\right]}$ *is stable asymptotically for the fuzzy initial value problem* (2).

In many cases, a deterministic solution to problem (2) is often difficult to obtain, therefore the author considers the method introduced in [10] which modifies the 5th order Runge-Kutta method for the fuzzy model as follows.

## 4. Runge-Kutta 5th order numerical scheme

In this section, we will study a two-dimensional fuzzy differential equation system with the form:

where *f, g*:ℝ^{2} → ℝ is continuous function, and ${X}_{0},{Y}_{0}\in {F}\text{}(\mathbb{R}\text{})$

By modifying the 5th order Runge-Kutta method, the solution to the fuzzy initial value problem (3) would be:

where

Since the arguments *X*_{i} and *Y*_{i} on the right-hand side are iterative, we can define a new function as follows:

so that (4) becomes

*F*_{h} and *G*_{h}are fuzzy number-valued functions, then

and

Let ${\left[{X}_{i+1}\right]}^{\alpha}=\left[{x}_{i+1,1}^{\alpha},{x}_{i+1,2}^{\alpha}\right]$ and ${\left[{Y}_{i+1}\right]}^{\alpha}=\left[{y}_{i+1,1}^{\alpha},{y}_{i+1,2}^{\alpha}\right]$, then we get

To approximate the solution (3) at each α-level, the partition *t*_{0} < *t*_{1} < *t*_{2} < ⋯ < *t*_{N−1} < *t*_{N} = *T* is created on the interval [*t*_{0}, *T*], with *t*_{i} = *t*_{0} + *ih*, *i* = 0, 1, 2, …, *N* and the length of the partition $h=\frac{T-{t}_{0}}{N}>\text{}0.$

## 5. Predator-prey fuzzy model with type II Holling functional response

In this section, two predator-prey models with type II Holling functional response will be studied to construct a fuzzy model. The first model was built from the deterministic model introduced by Jha et al. [34]. In this model, all parameters and the initial population are assumed to be certain, whereas in reality the parameter values and the initial population number cannot be known with certainty. In the next section, the model is expressed in a fuzzy model, where the initial population and uncertain parameters are expressed in fuzzy numbers. This model is expressed in Equations (4.1), (4.2), and (4.3) which is called model I. The second model is a modification of model I by considering harvesting. First of all, the deterministic model of model I is modified by adding a harvesting factor for both predator and prey. This model is then expressed in a fuzzy model, where the initial population and parameters are considered uncertain. This model is expressed in Equations (4.4), (4.5), and (4.6) which is called model II.

**Model I:** A predator-prey model with a type II Holling functional response introduced in Jha et al. [34] is given as:

where *x* and *y* are prey and predator population density at time *t, K* is environmental carrying capacity, *a* is prey intrinsic growth rate, *b* is prey predation rate, *c* is the predator mortality rate, *d* is the predator conversion, and *A* is the constant saturation factor of the predator. All model parameters are assumed to be positive.

If ($\overline{x,}\overline{y}$) is an equilibrium points of (8), by setting the derivatives equal to zero, we get the equilibrium points are (0, 0), (*K*, 0), and ($\overline{x}$, $\overline{y}$), where $\overline{x}=\frac{Ac}{d-c}$, $\overline{y}=\frac{a}{b}(A+\overline{x})\left(1-\frac{\overline{x}}{K}\right)$, and ($\overline{x,}\overline{y}$) is positive when *d* > *c*. The stability at these points is

• The system unstable at (0,0)

• The system is asymptotically stable at (K,0), if $K<\frac{Ac}{d-c}$ (or $K>\frac{Ac}{d-c}$ if *d* < *c*)

• The system is asymptotically stable at ($\overline{x}$, $\overline{y}$), if $K<\frac{a{(A+\overline{x})}^{2}}{b\overline{y}}$.

Suppose the initial population of prey and predators are uncertain, i.e., X_{0} and Y_{0} become fuzzy initial populations of prey and predators, respectively, at t_{0}. Then by applying the fuzzy initial value problem, where the initial population is a fuzzy number, the fuzzy predator-prey model of the system (8) becomes

where ${X}_{0},{Y}_{0}\in {F}({\mathbb{R}}^{2})$.

Based on Theorems (3.1) and (3.2), the fuzzy equilibrium point of the system (9) is χ_{{0, 0}}, χ_{{K, 0}}, and ${\chi}_{\left\{\frac{Ac}{d-c\text{}},\text{}\frac{a}{b}(A+\overline{x})\left(1-\frac{\overline{x}}{K}\right)\right\}}$ exists if *d* > *c*. The stability at these points is:

(i) The equilibrium point χ_{{0, 0}} is unstable

(ii) The equilibrium point χ_{{K, 0}} is asymptotically stable if $K<\frac{Ac}{d-c}$ (or $K>\frac{Ac}{d-c}$ if *d* < *c*)

(iii) The equilibrium point ${\chi}_{\left\{\frac{Ac}{d-c}\text{},\text{}\frac{a}{b}\left(A+\overline{x}\right)\left(1-\frac{\overline{x}}{K}\right)\right\}}$ asymptotically stable if $K\text{}\frac{a{\left(A+\text{}\overline{x}\right)}^{2}}{b\overline{y}}$

Suppose it is assumed that the parameter *a* is also uncertain. By changing the variable *X* in (9) into *X* = (*X*_{1}, *X*_{2}) = (*X, a*), then the fuzzy model (9) becomes

where ${X}_{0},{Y}_{0},a\in {F}(\mathbb{R})$.

**Model II**. If it is assumed that both prey and predator populations in system (8) are the target of harvesting efforts, then system (8) becomes

where *E* denotes the catch effort, and *p*_{1}, *p*_{2}, respectively, shows the coefficient catching power of prey and predator, where the function *p*_{i}*Ex* adoption from Das et al. [35].

If ($\overline{x,}\overline{y}$) is an equilibrium points of (11), by setting the derivatives equal to zero, we get the equilibrium points are $(0,\text{}0),\left(-\frac{K({p}_{1}E-a)}{a},\text{}0\right),\text{}$and ($\overline{x},\overline{y})$

where $\overline{x}=-\frac{A({p}_{2}E+c)}{{p}_{2}E-d+c}$, $\overline{y}=\frac{Ad({E}^{2}K{p}_{1}{p}_{2}-EKa{p}_{2}-EK{p}_{1}d+EK{p}_{1}c-EAa{p}_{2}+Kad-Kac-Aac)}{Kb{({p}_{2}E-d+c)}^{2}}$ exists if $E<\frac{d-c}{{p}_{2}}$ and *E*(*EKp*_{1}*p*_{2} − *Kap*_{2} − *Kp*_{1} *d* + *Kp*_{1} *c* − *Aap*_{2}) > *Kac* + *Aac* − *Kad*. The stability at these points is

(i) The system is stable at (0,0) if $E>\text{}\frac{a}{{p}_{1}}$

(ii) The system is asymptotically stable at $\left(-\frac{K({p}_{1}E-a)}{a},\text{}0\right)$, if $E<\text{}\frac{a}{{p}_{1}}$

Suppose the initial population of prey and predators are uncertain, i.e., X_{0} and Y_{0} becomes the initial fuzzy population of prey and predators, respectively, at t_{0}, then the fuzzy predator-prey model obtained from the system (11) would be

where ${X}_{0},{Y}_{0}\in {F}(\mathbb{R})$.

Based on Theorems (3.1) and (3.2), the fuzzy equilibrium points of the system (12) are χ_{{0, 0}}, ${\chi}_{\left\{-\frac{K({p}_{1}E-a)}{a},0\right\}},\text{}$and ${\chi}_{\left\{-\frac{A({p}_{2}E+c)}{{p}_{2}E-d+c},\text{}\frac{Ad({E}^{2}K{p}_{1}{p}_{2}-EKa{p}_{2}-EK{p}_{1}d+EK{p}_{1}c-EAa{p}_{2}+Kad-Kac-Aac)}{Kb{({p}_{2}E-d+c)}^{2}}\text{}\right\}}$ exists if $E<\frac{d-c}{{p}_{2}}$, and *E*(*EKp*_{1}*p*_{2} − *Kap*_{2} − *Kp*_{1} *d* + *Kp*_{1} *c* − *Aap*_{2}) > *Kac* + *Aac* − *Kad*. The stability at these points is:

(i) The equilibrium point χ_{{0, 0}} is stable *if* $E>\text{}\frac{a}{{p}_{1}}$

(ii) The equilibrium point ${\chi}_{\left\{-\frac{K({p}_{1}E-a)}{a},0\right\}}$ is asymptotically stable if $E\text{}\text{}\frac{a}{{p}_{1}}$.

## 6. Numerical simulation

In this section, we will explore the solution both of the above model for some case different according to the conditions of stability at each point of equilibrium using the 5th order Runge-Kutta method. Numerical simulations were carried out to compare the behavior of the crisp system and the fuzzy system.

**Model I:** For model I, Numerical simulation is divided into two cases. It is assumed for both cases the values of the parameters *a* = 0.5, *b* = 0.254, *K* = 1, 000, and *A* = 500. Parameters *c* and *d* in this model are assumed as follows:

(i) *c* = 0.125 and *d* = 0.325

(ii) *c* = 0.325 and *d* = 0.125

Suppose the initial population of prey and predators is *X*_{0} = 1, 100 and *Y*_{0} = 900.

**Case (i):**

For the case (i) where *d* > *c* and $K<\frac{a{(A+\overline{x})}^{2}}{b\overline{y}}$ obtained three equilibrium points: (0.0), (1,000, 0), and (312.5, 1099.594). Equilibrium points (0, 0) and (1,000, 0) are unstable, and at points (312.5, 1099.594) are asymptotically stable. The results of the numerical simulation of case (i) are presented in Figures 1, 2A. Figure 1A shows a stable system toward the equilibrium point (312.5, 1099.594). The phase plane graph for case (i) is presented in Figure 2A.

Let the initial population *X*_{0} and *Y*_{0} uncertain be defined as a triangular fuzzy number with ${\left[{X}_{0}\right]}^{\alpha}=\left[1,050+50\alpha ,\text{}1,150-50\alpha \right]$ and ${\left[{Y}_{0}\right]}^{\alpha}=\left[850+50\alpha ,\text{}950-50\alpha \right]$. Based on theorems 3.1 and 3.2, the fuzzy equilibrium points are obtained: χ_{{0, 0}}, χ_{{1, 000, 0}}, and χ_{{312.5, 1099.594}}, where the equilibrium point χ_{{0, 0}} and χ_{{1, 000, 0}} is unstable, and the equilibrium point χ_{{312.5, 1099.594}} is asymptotically stable. The numerical simulation results for the fuzzy model in case (i) are presented in Figures 1, 2B. Figure 1B shows the fuzzy system is stable toward the equilibrium point χ_{{312.5, 1099.594}}. The fuzzy phase plane graph for case (i) is presented in Figure 2B.

Let other than *X*_{0} and *Y*_{0}, parameter *a* is also uncertain and is expressed in triangular fuzzy numbers with [*a*]^{α} = [0.4 + 0.1α, 0.6 − 0.1α], then the behavior of the system shown in Figure 3.

**Figure 3**. **(A)** Population growth over time. **(B)** Phase plane fuzzy for case (i) of model I with parameter a and population *X*_{0}, *Y*_{0} is fuzzy.

The simulation results can be seen that by adding a parameter as a fuzzy number, the dynamic behavior of the fuzzy system qualitatively shows the same results when only the initial populations of prey and predators are fuzzy, and this is in accordance with the behavior of the crisp system.

**Case (ii):**

For the case (ii) where *d* < *c* and $K>\frac{Ac}{d-c}$, three equilibrium points are obtained: (0.0), (1,000, 0), and (−812.5, −1114.973). The point (−812.5, −1114,973) is ignored because it has a negative value, the stability of the equilibrium point is: at point (0,0) is unstable and at point (1,000, 0) is asymptotically stable. The results of the numerical simulation of case (ii) are presented in Figures 4, 5A. Figure 4A shows the system is stable toward the equilibrium point (1,000, 0). The phase plane graph for case (ii) is presented in Figure 5A.

Suppose the initial population is uncertain and defined as a triangular fuzzy number as in case (i). Based on theorems 3.1 and 3.2, the fuzzy equilibrium points are obtained: χ_{{0, 0}}, and χ_{{1, 000, 0}}, where at the equilibrium point χ_{{0, 0}} is unstable and at the point χ_{{1, 000, 0}} asymptotically stable. The numerical simulation results for the fuzzy model in case (ii) are presented in Figures 4, 5B. Figure 4B shows the fuzzy system is stable toward the equilibrium point χ_{{1, 000, 0}}. The fuzzy phase plane graph for case (ii) is presented in Figure 5B.

Let other than *X*_{0} and *Y*_{0}, parameter *a* is also uncertain and is expressed in triangular fuzzy numbers as in case (i), then the behavior of the system is shown in Figure 6.

**Figure 6**. **(A)** Population growth over time. **(B)** Phase plane fuzzy for case (ii) of model I with parameter a and population *X*_{0}, *Y*_{0} are fuzzy.

The results of the numerical simulation are presented in Figures 1–6 shows that for both cases of model I, the solution graph with the fuzzy approach differs quantitatively, but qualitatively gives the same results as the graph from the crisp system when only the initial population is fuzzy, but slightly different when one of the parameters and the initial population is fuzzy.

**Model II:** For model II, Numerical simulation is divided into two cases. It is assumed for both cases the values of the parameters and the initial population are the same as for model I: *a* = 0.5, *b* = 0.254, *K* = 1, 000, and *A* = 500. And other parameters in this model are assumed as follows:

• *c* = 0.325, *d* = 0.125, *p*_{1} = 1, *p*_{2} = 2, and *E* = 0.275

• *c* = 0.125, *d* = 0.325, *p*_{1} = 1, *p*_{2} = 1, and *E* = 0.6

Suppose the initial population of prey and predators is *X*_{0} = 1, 100 and *Y*_{0} = 900.

**Case (i):**

For case (i) where $E<\frac{a}{{p}_{1}}$ obtained three equilibrium points: (0.0), (450,0), and (−583.3, −169.51). The point (−583.3, −169.51) is ignored because it is negative. The stability at equilibrium point (0, 0) is unstable and at (450,0) is asymptotically stable. The numerical simulation results of case (i) are presented in Figures 7, 8A. Figure 7A shows the system is stable toward the point of equilibrium (450.0). The phase plane graph for case (i) is presented in Figure 8A.

Suppose the initial population *X*_{0} and *Y*_{0} uncertain is defined as a triangular fuzzy number as in model I, then the fuzzy equilibrium point χ_{{0, 0}} is unstable, and the fuzzy equilibrium point χ_{{450, 0}} is asymptotically stable. The numerical simulation results for the fuzzy model in case (i) are presented in Figures 7, 8B. Figure 7B shows a stable fuzzy system toward the equilibrium point χ_{{450, 0}}. The fuzzy phase plane graph for case (i) is presented in Figure 8B.

Case **(ii):**

For the case (ii) where $E>\frac{a}{{p}_{1}}$ obtained three equilibrium points: (0.0), (−200,0), and (−906.25, −564.792). The points (−200,0) and (−906.25, −564.792) are ignored because they are negative. The numerical simulation results of case (ii) are presented in Figures 9, 10A. Figure 9A shows the system is stable toward the equilibrium point (0,0). The phase plane graph for case (ii) is presented in Figure 10A.

Suppose the initial population is uncertain and defined as a triangular fuzzy number as in model I, then the fuzzy equilibrium point χ_{{0, 0}} is stable. The numerical simulation results for the fuzzy model in case (ii) are presented in Figures 9, 10B. Figure 9B shows the fuzzy system is stable toward the equilibrium point χ_{{0, 0}}. The fuzzy phase plane graph for case (ii) is presented in Figure 10B.

## 7. Discussion

The research presented in this paper is an extension of the predator-prey model with a type II Holling functional response discussed by Jha et al. [34] taking into account the uncertainty in the parameters and the initial population expressed in fuzzy numbers. This model is further expanded by adding harvesting factors to both populations. In this study, the behavior of the system is only studied qualitatively by performing numerical simulations to explore the behavior of the fuzzy system and compare it with the crisp system. In conducting the simulation, we use triangular fuzzy numbers to express uncertainty in the initial population and parameters. Of the two models studied, we found the same results. In both models, the fuzzy system shows the same behavior as the crisp system when only the initial population of prey and predators is fuzzy or when one parameter is added as a fuzzy parameter. Another result obtained is, for the fuzzy model, the time required for the system to reach equilibrium is longer than the crisp model. This is due to the uncertainty in the initial value expressed in the fuzzy interval. This research only studied the behavior of the system qualitatively through numerical simulation. In this numerical simulation, all the figures of the phase planes for the fuzzy systems above are plotted for the value of the α-level equals zero. It seems that the phase planes for the crisp systems can be extracted from the fuzzy system's phase planes with the α-level equals zero. This is an interesting result that can be interpreted crisp model can be used as a special case of fuzzy model whenever the degree of uncertainty is relatively low.

## 8. Conclusion

In this paper, we have developed a numerical scheme to find the solution of two predator-prey models with a Holling type II functional response by considering fuzzy parameters and fuzzy initial populations. The first model was developed from the model studied by Jha et al. [34] by replacing the initial population and one of the parameters with a fuzzy number. While the second model was developed from the first model by adding harvesting factors to both prey and predator populations. The behavior of the model was studied qualitatively using the Runge-Kutta method of order-5 which was modified for the fuzzy system using the Zadeh extension principle. The numerical simulation results show that, when the initial population prey and predators that have fuzzy values, then both fuzzy models have the same behavior as the crisp model, but the fuzzy model takes a longer time to achieve stability than the crisp model. This is due to the uncertainty in the initial population which is indicated by fuzzy intervals. Likewise, when one parameter is added with a fuzzy value, the fuzzy model has the same behavior as the crisp model. Finally, we can conclude that fuzzy behavior represents a generalization of crisp behavior, and this gives more realistic results that represent the problem of uncertainty. However, there are still much work to be done in the future, including studying the stability of the system analytically, bifurcation problems, and others. As pointed by one of the reviewer, it “would have been more interesting to place the choice of parameters on the crisp model exibiting marked sensitivity behavior to initial conditions” and this is currently under investigation by the authors.

## Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

## Author contributions

AS designed and supervised the research. IS developed the scheme and numerical simulations. EC and NA conducted the literature review. All authors contributed to the articles and approved the submitted versions.

## Funding

This research was published with funding from Universitas Padjadjaran through the Article Processing Charge (APC) Replacement Scheme.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

1. Wang Q, Zhai S, Liu Q, Liu Z. Stability and optimal harvesting of a predator-prey system combining prey refuge with fuzzy biological parameters. *Math Biosci. Eng*. (2021) 18:9094–120. doi: 10.3934/mbe.2021448

2. Supriatna AK, Possingham HP. Optimal harvesting for a predator-prey metapopulation. *Bull Math Biol*. (1998) 60:49–65. doi: 10.1006/bulm.1997.0005

3. Supriatna AK, Possingham HP. Harvesting a two-patch predator-prey metapopulation. *Nat Resour Model*. (1999) 12:481–98. doi: 10.1111/j.1939-7445.1999.tb00023.x

4. Dawes JHP, Souza MO. A derivation of Holling's type I, II and III functional responses in predator-prey systems. *J Theor Biol*. (2013) 327:11–22. doi: 10.1016/j.jtbi.2013.02.017

5. Jana S, Kar TK. Modeling and analysis of a prey-predator system with disease in the prey. *Chaos Solitons Fract.* (2013) 47:42–53. doi: 10.1016/j.chaos.2012.12.002

6. Ma Z, Wang S, Wang T, Tang H. Stability analysis of prey-predator system with Holling type functional response and prey refuge. *Adv Differ Equat.* (2017) 2017:243. doi: 10.1186/s13662-017-1301-4

7. Gomes LT, de Barros LC, Bede B. Fuzzy differential equations in various approaches (2015). doi: 10.1007/978-3-319-22575-3

8. Mizukoshi MT, Barros LC, Chalco-Cano Y, Román-Flores H, Bassanezi RC. Fuzzy differential equations and the extension principle. *Inf Sci.* (2007) 177:3627–35. doi: 10.1016/j.ins.2007.02.039

9. Bandyopadhyay A, Kar S. *System of Type-2 Fuzzy Differential Equations and Its Applications.* Vol. 31. London: Springer (2019).

10. Ahmad MZ, Hasan MK. Modeling of biological populations using fuzzy differential equations. *Int J Mod Phys Conf Ser*. (2012) 9:354–63. doi: 10.1142/S2010194512005429

11. Barros LC, Gomes LT, Tonelli PA. *Fuzzy Differential Equations: An Approach via Fuzzification of the Derivative Operator.* Vol. 230. Elsevier (2013). doi: 10.1016/j.fss.2013.03.004

12. Alamin A, Mondal SP, Alam S, Goswami A. Solution and stability analysis of non-homogeneous difference equation followed by real life application in fuzzy environment. *Sadhana Acad Proc Eng Sci*. (2020) 45:185. doi: 10.1007/s12046-020-01422-1

13. Mizukoshi MT, Barros LC, Bassanezi RC. Stability of fuzzy dynamic systems. *Int Uncertain J Fuzziness Knowl Based Syst.* (2009) 17:69–83. doi: 10.1142/S0218488509005747

14. Jayakumar T, Maheskumar D, Kanagarajan K. Numerical solution of fuzzy differential equations by Runge Kutta method of order five. *Appl Math Sci*. (2012) 6:2989–3002. doi: 10.17654/FS021020135

15. Nayak S, Chakraverty S. Numerical solution of fuzzy stochastic differential equation. *J Intell Fuzzy Syst*. (2016) 31:555–63. doi: 10.3233/IFS-162168

16. Behroozpoor AA, Vahidian Kamyad A, Mazarei MM. Numerical solution of fuzzy initial value problem (FIVP) using optimization. *Int J Adv Appl Sci.* (2016) 3:36–42. doi: 10.21833/ijaas.2016.08.007

17. Tapaswini S, Chakraverty S. Numerical solution of Fuzzy arbitrary order predator-prey equations. *Appl Appl Math Int J*. (2013) 8:647–72.

18. Tapaswini S, Chakraverty S. A New approach to fuzzy initial value problem by improved euler method. *Fuzzy Inf Eng*. (2012) 4:293–312. doi: 10.1007/s12543-012-0117-x

19. da Silva Peixoto M, de Barros LC, Bassanezi RC. Predator-prey fuzzy model. *Ecol Model*. (2008) 214:39–44. doi: 10.1016/j.ecolmodel.2008.01.009

20. Pandit P, Singh P. *Prey Predator Model With Fuzzy Initial Conditions*. International Journal of Web Engineering and Technology (IJEIT). Vol. 3. (2014). p. 65–8.

21. Ak O, Oru O. *A Prey Predator Model With Fuzzy Initial Values*. Hacettepe Journal of Mathematics and Statistics. Vol. 41. (2012). p. 387–95.

22. Ahmad MZ, De Baets B. A predator-prey model with fuzzy initial populations. In: *2009 International Fuzzy Systems Association world congress and 2009 European Society of Fuzzy Logic and Technology*. (2009). p. 1311–4. Available online at: https://www.researchgate.net/publication/221399560_A_Predator-Prey_Model_with_Fuzzy_Initial_Populations

23. Narayanamoorthy S, Baleanu D, Thangapandi K, Perera SSN. Analysis for fractional-order predator–prey model with uncertainty. *IET Syst Biol*. (2019) 13:277–89. doi: 10.1049/iet-syb.2019.0055

24. Omar AHA, Ahmed IAA, Hasan YA. The fuzzy ratio prey-predator model. *Int J Comput Sci Electron Eng*. (2015) 3:101–6.

25. Pal D, Mahapatra GS, Samanta GP. Parameter uncertainty in biomathematical model described by one-prey two-predator system with mutualism. *Int J Biomath*. (2017) 10:1750082. doi: 10.1142/S1793524517500826

26. Mallak S, Farekh D, Attili B. Numerical investigation of fuzzy predator-prey model with a functional response of the form arctan(ax). *Mathematics.* (2021) 9:17–9. doi: 10.3390/math9161919

27. Pal D, Mahapatra GS, Samanta GP. Stability and bionomic analysis of fuzzy prey–predator harvesting model in presence of toxicity: a dynamic approach. *Bull Math Biol*. (2016) 78:1493–519. doi: 10.1007/s11538-016-0192-y

28. Yu X, Yuan S, Zhang T. About the optimal harvesting of a fuzzy predator–prey system: a bioeconomic model incorporating prey refuge and predator mutual interference. *Nonlinear Dyn*. (2018) 94:2143–60. doi: 10.1007/s11071-018-4480-y

29. Pal D, Mahapatra GS, Samanta GP. Stability and bionomic analysis of fuzzy parameter based prey–predator harvesting model using UFM. *Nonlinear Dyn*. (2015) 79:1939–55. doi: 10.1007/s11071-014-1784-4

30. Meng XY, Wu YQ. Dynamical analysis of a fuzzy phytoplankton–zooplankton model with refuge, fishery protection and harvesting. *J Appl Math Comput*. (2020) 63:361–89. doi: 10.1007/s12190-020-01321-y

31. Mahata A, Mondal SP, Roy B, Alam S. Study of two species prey-predator model in imprecise environment with MSY policy under different harvesting scenario. *Environ Dev Sustain*. (2021) 23:14908–32. doi: 10.1007/s10668-021-01279-2

32. Pal D, Mahapatra GS, Samanta GP. A study of bifurcation of prey-predator model with time delay and harvesting using fuzzy parameters. *J Biol Syst*. (2018) 26:339–72. doi: 10.1142/S021833901850016X

33. de Barros LC, Bassanezi RC, Lodwick WA. *Studies in Fuzziness and Soft Computing A First Course in Fuzzy Logic, Fuzzy Dynamical Systems, and Biomathematics Theory and Applications* Berlin: Springer (2017).

34. Jha P, Ghorai S. Stability of prey-predator model with Holling type response function and selective harvesting. *J Appl Comput Math*. (2017) 6:3. doi: 10.4172/2168-9679.1000358

Keywords: predator-prey fuzzy model, Holling type II functional response, fuzzy parameter, fuzzy initial population, Zadeh extension principle, 5th order Runge-Kutta method

Citation: Sukarsih I, Supriatna AK, Carnia E and Anggriani N (2023) A Runge-Kutta numerical scheme applied in solving predator-prey fuzzy model with Holling type II functional response. *Front. Appl. Math. Stat.* 9:1096167. doi: 10.3389/fams.2023.1096167

Received: 11 November 2022; Accepted: 27 February 2023;

Published: 23 March 2023.

Edited by:

Dumitru Trucu, University of Dundee, United KingdomReviewed by:

Cyrille Bertelle, University of Le Havre, FranceShariful Alam, Indian Institute of Engineering Science and Technology, Shibpur, India

Copyright © 2023 Sukarsih, Supriatna, Carnia and Anggriani. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: I. Sukarsih, icih21002@mail.unpad.ac.id