#### Abstract

In this paper, a mathematical model for the system of prey-predator with immigrant prey has been analyzed to find an approximate solution for immigrant prey population density, local prey population density, and predator population density. Furthermore, we present a novel soft computing technique named LeNN-WOA-NM algorithm for solving the mathematical model of the prey-predator system with immigrant prey. The proposed algorithm uses a function approximating ability of Legendre polynomials based on Legendre neural networks (LeNNs), global search ability of the whale optimization algorithm (WOA), and a local search mechanism of the Nelder–Mead algorithm. The LeNN-WOA-NM algorithm is applied to study the effect of variations on the growth rate, the force of interaction, and the catching rate of local prey and immigrant prey. The statistical data obtained by the proposed technique establish the effectiveness of the proposed algorithm when compared with techniques in the latest literature. The efficiency of solutions obtained by LeNN-WOA-NM is validated through performance measures including absolute errors, MAD, TIC, and ENSE.

#### 1. Introduction

A prey-predator system is one of the prevalent phenomena in nature. The interaction between predators and prey in any environment was first introduced by Lotka and Volterra in 1926 [1]. Holling in 1966 [2] elaborated prey-predictor models with different kinds of functional responses for predation. All these models are inspired by biological phenomena and presented by nonlinear ordinary and partial differential equations. Danca et al. [3] study the detailed analysis of a nonlinear prey-predator model. A connection between predators and prey has a long history and will continue as a governing theme in biomathematics because of its universal significance [4]. Researchers have made many changes by introducing different facets to the predator-prey system such as delay in predator growth, harvesting prey and predators, and providing additional food to a predator for sustaining prey population and prey diseases [5–7]. Cai et al. [8] studied the dynamical system of prey-predator with Allee effects in prey growth. A stage-structured predator-prey model with gestation delay is investigated by [9]. In [10, 11], bifurcation, chaos and dynamic behavior of the nonlinear discrete-time predator-prey system is studied. Huang et al. [12] study the stability analysis of the model with consideration of the prey refuge. Studies by Din [13], Weide et al. [14], and Gong et al. [15, 16] are presented for the discrete-time nonlinear prey-predator types of model. Hadeler and Freedman [17] provide extensive references to real-world examples of three-species ecoepidemiological systems of sound prey, infected prey, and predators.

In recent years, there has been a substantial increase in the amount of attention paid to population biology by scientists due to the significant applications it has in ecology. It bridges the gap between mathematics and biology. The dynamical systems in biology have been investigated to interpret various problems. The Lotka–Volterra population model is a well-known mathematical model that describes biological systems [18]. In real life, there is a strong correlation between size, age, and developmental stages of different populations of species. It is an important strategy to incorporate all these variables in the mathematical modeling of populations of different spices. To develop more realistic and accurate mathematical models, many scientists have suggested noise-induced models [19, 20] and spatial models [21]. In biomathematical problems, researchers utilize mathematical modeling and simulations along with biological structures to describe the phenomena by the system of nonlinear differential equations. These nonlinear models are considered stiff and unrealistic, and therefore, finding exact and semianalytical solutions for such problems is challenging because of nonlinearity. Analytical approaches for solving nonlinear problems are mostly based on Laplace or Fourier transformation, Laguerre’s integral formula, and the Grunwald–Letnikov concept [22, 23]. When tackling complex problems, these techniques may be challenging to use; in addition, the solution is provided in a closed form that necessitates the evaluation of special functions using complex expressions, such as the Mittag-Leffler function [24].

In recent years, researchers have been working to develop new techniques for finding approximate solutions to nonlinear models; e.g., the Laplace Adomian decomposition method [25], the new coupled fractional reduced differential transform method [26], the Runge–Kutta–Fehlberg method [27], the finite element method [28], the Sumudu decomposition method [29], the implicit Adams methods [30], the confidence domain technique [31], and the homotopy analysis method [32] have become much more significant to get accurate solutions. All these deterministic approaches have their own advantages, applicability, and drawbacks. With great interest, it is noted that such techniques are gradient-based and call for information about the problem beforehand. The availability of several local optima, which leads to solutions where global optimality cannot be easily ensured, is one of the fundamental limitations of gradient-based approaches. Global optimality is sought in gradient-based approaches by randomly scanning the design space from various starting points. However, this causes the technique to become sluggish and computationally inefficient for complex nonlinear optimization problems [33, 34]. In addition, the Runge–Kutta methods are self-starting and stable techniques that can be easily implemented to calculate the solution to different problems. The main drawbacks of the Runge–Kutta methods are that they take longer time to calculate solutions than other multistep methods with equivalent precision, and it is difficult for them to get accurate global estimates of the truncation error [35]. To overcome these drawbacks, a stochastic metaheuristic approach with artificial neural networks is developed, which is free of a gradient and does not require any prior information about the problem. The majority of metaheuristic techniques are inspired by natural, physical, or biological processes and make use of a variety of operators to mimic the fundamental behavior. The harmony between exploration and exploitation is a recurring subject in all metaheuristics.

In recent times, artificial neural network (ANN)-based stochastic algorithms with global and local search optimizers have been designed to solve differential equations representing physical phenomena including flow in a circular cylindrical conduit via electrohydrodynamics [36], a model of an immobilized enzyme system that follows the Michaelis–Menten (MM) kinetics for a microdisk biosensor [37], flow of Johnson–Segalman fluid on the surface of an infinitely long vertical cylinder [38], and beam-column designs [39]. The abovementioned techniques motivate authors to design a new soft computing algorithm, the LeNN-WOA-NM algorithm, to find approximate series solutions using Legendre polynomials for the model presenting the prey-predator system with immigrant prey. The salient features of the paper are summarized as follows:(i)A mathematical model for a prey-predator system with immigrant prey is formulated and analyzed to study the influence of variations on the growth rate, force of interaction, and the catching rate of local and immigrant prey.(ii)Artificial neural network-based weighted Legendre polynomials are used to construct the model of approximate solutions for the prey-predator model. A fitness function based on mean square errors is designed to assess unknown parameters with the help of global search ability of whale optimization and a local search mechanism of the Nelder–Mead algorithm.(iii)The suggested technique can result in adequate solutions for nonlinear hard problems for which no exact algorithm exists that can solve them in a reasonable amount of time.(iv)Performance of the proposed algorithm is validated in terms of absolute errors, mean absolute deviation, Theil’s inequality coefficient, Nash–Sutcliffe efficiency, and error in Nash–Sutcliffe efficiency.(v)Approximate solutions, convergence of fitness, and performance measures obtained by the LeNN-WOA-NM algorithm for prey-predator systems with immigrant prey are shown through different graphs and tables, which shows the dominance and robustness of the proposed algorithm in solving real-world problems.(vi)Unlike most classical methods, the LeNN-WOA-NM algorithm requires no gradient information and therefore can be used with nonanalytic, black-box, or simulation-based objective functions to approximate complex problems.

#### 2. Problem Formulation

Goteti developed the mathematical model for a prey-predator system with immigrant prey to understand the interaction and communication between predators and immigrant prey. The model is assumed to follow mass action theory and consists of prey population density, which is defined as follows:where and denote the population density of local and immigrant prey, respectively. Population density of the predator is denoted by .

The following assumptions are used in formulating the mathematical model for the prey-predator system with immigrant prey:(i)In the presence of a predator, the population of prey is classified into two subcategories named local prey and immigrant prey .(ii)In absence of the predator, the population growth of local prey logistically increases with an intrinsic growth rate with environmental carrying capacity denoted by .(iii)With availability of the local and immigrant prey population, the population of the predator grows logistically with growth rates and , while suffering loss of populations is denoted by and .(iv)Immigrant and local prey can reproduce, and therefore, it is assumed that birth rates should be positive. The growth rate of immigrant prey increases at the rate with environmental carrying capacity denoted by .(v)It is assumed that immigrant prey is a natural choice of predators. and denote the positive and negative force of interaction between local and immigrant prey.(vi)It is assumed that the local prey and immigrant prey are caught by the predator at the rate of and , respectively.

A mathematical model for a prey-predator system with immigrant prey is given by the following system of differential equations:with initial populations

#### 3. Approximate Solutions and Weighted Legendre Polynomials

The Legendre polynomials are denoted by , where denotes the order of Legendre polynomials. These polynomials constitute the set of orthogonal polynomials on . The first eleven Legendre polynomials are given in Table 1. High-order Legendre polynomials are generated by the following recursive formula:

We consider an approximate series solution for equations (2)–(4) representing the prey-predator model with immigrant prey as follows:where , , and are unknown parameters.

Since order continuous derivatives of system (6) exist, we consider the first derivative of system (6) as follows:and we plug equations (6)–(7) in governing ordinary differential equations. Equations (2)–(4) will be transformed into an equivalent algebraic system of equations that can be solved for unknown parameters , , and using the LeNN-WOA-NM algorithm.

#### 4. Fitness Function Formulation

The mean square error (MSE)-based fitness function for solving the prey-predator model for equations (2)-(4) can be written as follows:where to are the fitness-based errors for equations (2)–(4) along with the initial conditions and are given as follows:

#### 5. Optimization Network

##### 5.1. Whale Optimization Algorithm

The whale optimization algorithm is a nature-inspired metaheuristic algorithm designed by Mirjalili and Lewis [40]. The working strategy of WOA is inspired by foraging behavior of humpback whales. The humpback whales chase prey or krill by swimming around them in a molded way as shown in Figure 1. The mathematical model of each phase is explained below.

###### 5.1.1. Exploration Phase

Humpback whales encircle prey for hunting. Equations (10) and (11) mathematically model this behavior as follows:where denotes the current iteration, has provided the best solution so far, and gives the location of humpback whales to prey at each step. and are coefficient vectors which are defined as follows:where is an arbitrary nominated vector. The value of reduces from two to zero during exploration as well as in the exploitation phase. During helix-shaped movement, the distance between prey and the humpback whale is given as follows: where represents the location of the whale to the prey and is defined as follows:

Furthermore, denotes the state of the logarithmic helix and is any arbitrary number. The shrinking surrounding technique of humpback whales during contracting loop is summarized as follows:

###### 5.1.2. Exploitation Phase

In the exploration phase (searching for prey), the heterogeneity of the vector is utilized. If , the position of the search agent is updated, and the entire mechanism is modeled by the following equations:where is an arbitrary position vector taken from the current population.

Figure 2 shows the flowchart of the WOA. It can be seen that the WOA creates a random initial population from candidate space and evaluates it using an error-based fitness function when the optimization process starts. After finding the best solution, the algorithm repeatedly executes the following steps until ending criterion is achieved.

##### 5.2. Nelder–Mead Algorithm

A Nelder–Mead (NM) algorithm is a direct search method also known as a downhill simplex method developed by Nelder and Mead in 1965 to solve different problems without any information about the gradient [41]. NM is a single path following a local search optimizer that can find good results if initialized with a better initial solution. A simplex consisting of vertices is set up to minimize a function with dimensions *n* [42]. The NM algorithm generates a sequence of simplices by following four basic procedures, namely, reflection, expansion, contraction, and shrink. Further details about the NM algorithm can be found in [43]. Figure 2 shows the working procedure of the NM algorithm.

In recent times, the most successful and effective trend in optimization is the action of integrating components from different methods. The foremost motivation behind the hybridization of diverse algorithmic ideas is to acquire better performing systems, which exploit and coalesce benefits of different techniques. Therefore, in this study, we have combined the global and local search optimization algorithms to achieve more efficient and robust solutions. The hybridized working procedure of the designed LeNN-WOA-NM algorithm is illustrated in Figure 2.

#### 6. Performance Measures

To examine the performance of the proposed algorithm, performance indices are defined in terms of mean absolute deviation (MAD), Theil’s inequality coefficient (TIC), and error in Nash–Sutcliffe efficiency (ENSE). Mathematical formulation for , , and of the predator and prey model in the case of MAD, TIC, and ENSE is presented as follows:where shows total input grid points. The values of performance measures of MAD, RMSE, and ENSE should be equal to zero for perfect modeling, while NSE should be equal to 1. The global versions of MAD, RMSE, and ENSE for the given mathematical model of the prey-predator system are formulated by given equations:where denotes the number of independent runs. Global fitness (GFIT) is defined as the mean of fitness values attained in independent runs. Mathematically, GFIT is given as follows:

#### 7. Result and Discussion

The numerical solutions of the prey-predator system with immigrant prey under the influence of variations in various parameters are investigated with the proposed methodology. The exact solution for these nonlinear models is not known, so the comparative study is conducted with MATLAB solver ode45 and the homotopy perturbation method [44]. Six problems of the prey-predator model are considered for different cases depending on values of coefficients, i.e, , , and . A graphical abstract of the paper is shown in Figure 3.

##### 7.1. Problem I: Influence of Variations in on the Prey-Predator Model

In this problem, the effect of variations in the intrinsic growth rate of local prey on population density is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined by

Five cases are considered, depending on the value of . Case I: Case II: Case III: Case IV: Case V:

The LeNN-WOA-NM algorithm is applied to prey-predator model equation (29) to study the influence of variations in the intrinsic growth rate of local prey. A comparison between approximate solutions obtained by the LeNN-WOA-NM algorithm for the prey-predator system with the homotopy perturbation method [44] for is illustrated in Table 2. Approximate solutions for population densities of local prey, immigrant prey, and predator are given in Tables 3–5, respectively. Absolute errors are presented in Tables 6–8 and are graphically illustrated in Figure4. Table 9 represents the statistics of global values of performance indicators during 100 independent trails. Figure 5 shows the convergence of the fitness value. It can be seen that the fitness value for each case lies around to . The percentage convergence of the fitness value and performance indicators during multiple runs is shown in Table 10. Trained neurons in the LeNN structure for obtaining best solutions are shown in Table 11. From Figure 6, the following conclusions are drawn:(i)Population density of local prey has a direct relation with the intrinsic growth rate of local prey(ii)Population density of immigrant prey has an inverse relation with the intrinsic growth rate of local prey(iii)Population density of the predator varies directly with the intrinsic growth rate of local prey

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 7.2. Problem II: Effect of Variations in on the Prey-Predator Model

In this problem, the effect of variations in the positive impact of force of interaction between local and immigrant prey on population densities of the prey-predator model is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined as follows:

Five cases are considered, depending on the value of . Case I: Case II: Case III: Case IV: Case V:

Approximate solutions for the effect of variations in on population densities of the prey-predator model are given in Tables 12–14. Absolute errors in our solution for population densities are presented in Tables 15–17, respectively. The solution of the design scheme overlaps the exact solution with absolute errors that lie between and as illustrated in Figure 7. Figure 8 shows the behavior of fitness evaluation for the prey-predator model under the influence of . Convergence analysis of performance measures is given in Tables 18 and 19. Trained neurons obtained by the LeNN-WOA-NM algorithm for different cases of Eq (31) are shown in Table 20. From the solutions (see Figure 9), the following conclusions can be drawn:(i)The negative force of interaction between local and immigrant prey has a direct relation with population density of local prey and the predator, while population density is inversely related to

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 7.3. Problem III: Effect of Variations in on the Prey-Predator Model

In this problem, the effect of variations in the increasing rate of immigrant prey on population densities of the prey-predator model is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined as follows:

Four cases are considered, depending on the value of . Case I: Case II: Case III: Case IV:

The LeNN-WOA-NM algorithm is used to optimize the population densities of equation (33). Table 21 represents the comparison between the Ranga–Kutta method and the proposed technique LeNN-WOA-NM. A comparison between population density of local prey , immigrant prey , and predator with variations in is shown in Table 22. Statistics of absolute errors are given in Table 23 and graphically presented in Figure 10. The behavior of the fitness function for different cases is shown on boxplots as demonstrated in Figure 11. The mean values of the fitness function lie around as shown in Tables 24 and 25. Unknown parameters achieved by the LeNN-WOA-NM algorithm for the given problem are given in Table 26. From solutions (see Figure 12), the following conclusion can be drawn:(i)Population density of local prey, immigrant prey, and predator varies directly with an increasing rate of immigrate prey

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 7.4. Problem IV

In this problem, the effect of variations in the negative force of interaction between local prey and immigrant prey on population densities of the prey-predator model is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined as follows:

Three cases are considered, depending on the value of . Case I: Case II: Case III:

The LeNN-WOA-NM algorithm is used to optimize the population densities of prey-predator model equation (35). Population densities of local prey, immigrant prey, and the predator are given in Tables 27-28, while absolute errors are shown in Table 29. The absolute errors in the solutions of the proposed technique lie between and as shown in Figure 13. Table 30 shows global values of performance indices and convergence of the fitness value as shown in Figure 14. The statistics shown in Table 31 illustrate that the global values of different performance functions lie between and , which highlights the robustness of the technique. The values of the weights in the LeNN structure for recreation of the approximate solutions are given in Table 32. From Figure 15, the following conclusion can be drawn:(i)Population density of local prey varies directly with variations in (ii)Population density of immigrant prey and predator varies inversely with variations in

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 7.5. Problem V: Effect of Variation in on the Prey-Predator Model

In this problem, the effect of variations in the catching rate of local prey on population densities of the prey-predator model is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined as follows:

Four cases are considered, depending on variations in . Case I: Case II: Case III: Case IV:

The LeNN-WOA-NM algorithm is used to optimize the population densities of equation (37). Table 33 represents the comparison between the Ranga–Kutta method (ode45) and the proposed technique LeNN-WOA-NM algorithm. Approximate solutions for population density of local prey , immigrant prey , and the predator with different values of are presented in Table 34. Statistics of absolute errors are given in Table 35 and graphically shown in Figure 16. The minimum and mean values of the fitness function for each case lie between and , respectively as illustrated in Figure 17. Convergence analysis of the fitness value, MAD, TIC, and ENSE is given in Tables 36-37. Unknown neurons of LeNN are shown in Table 38. From Figure 18, the following conclusion can be drawn:(i)Population density of local prey and population density of the predator vary inversely with variations in (ii)Population density of immigrant prey has a direct relation with variations in

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 7.6. Problem VI: Effect of Variations in on the Prey-Predator Model

In this problem, the effect of variations in the catching rate of immigrant prey on population densities of the prey-predator model is discussed. An error-based fitness function along with initial populations is given as follows:where to are defined as follows:

Five cases are considered, depending on the value of . Case I: Case II: Case III: Case IV:

The LeNN-WOA-NM algorithm is applied to prey-predator model equation (39) to study the influence of variations in the catching rate of immigrant prey. The results calculated by the designed algorithm are compared with those of the Runge–Kutta method using ode45 in MATLAB as shown in Table 39. Approximate solutions for population densities of local prey, immigrant prey, and the predator are given in Table 40. Absolute errors are presented in Table 41 and graphically shown in Figure 19. The minimum fitness values in Figure 20 reflect the accuracy of the solutions by the proposed algorithm. The convergence of our numerical approach is assessed by fitness values and statistical analysis of the performance indicators (MAD, TIC, and ENSE), see Table 42. Table 43 shows the correctness of the proposed algorithm to tackle real-world problems. Trained neurons in LeNN are shown in Table 44. From Figure 21, the following conclusions can be drawn:(i)Population densities of local prey, immigrant prey, and the predator vary inversely with variation in the catching rate of immigrant prey .

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

#### 8. Conclusion

In this paper, we have analyzed the system of the singular differential equation (SDE) representing the phenomena of the prey-predator model with immigrant prey. Generally, solving SDE is one of the challenging tasks, and therefore, we have designed a novel soft computing technique known as a LeNN-WOA-NM algorithm. Weighted Legendre polynomials are used to model approximate series solutions for the prey-predator model with variations in various coefficients, including the growth rate of local and immigrant prey , force interaction between local and immigrant prey , and the catching rate of local and immigrant prey . We summarize our findings as follows:(i)Variations in , and has a direct impact on population density of local prey , while population density of local prey is inversely affected by variations in and .(ii)Variations in , , and has an inverse impact on population density of immigrant prey , while population density of immigrant prey is directly affected by variations in and .(iii)Variations in and has a direct impact on population density of the predator , while population density of the predator is inversely affected by variations in , , and .(iv)Approximate solutions obtained by the LeNN-WOA-NM algorithm are compared with those obtained by the homotopy perturbation method and MATLAB solver ode45. The results show the dominance of the proposed technique.(v)Lower absolute errors in our solutions and convergence analysis of fitness evaluation, MAD, TIC, and ENSE show the accuracy and robustness of the proposed algorithm for obtaining solutions to real-world problems.

In future, this approach can be utilized to solve the complex nonlinear systems of fractional differential equations characterizing real-world problems, for instance, anomalous diffusion of contaminant from the fracture into the porous rock matrix, microbial survival and growth curves, smoking dynamics, parametric identification of Hammerstein systems with time delay, and asymmetric dead zones.

#### Abbreviations

LeNN: | Legendre neural network |

NM: | Nelder–Mead |

MAD: | Mean absolute diviation |

TIC: | Theil’s inequality coefficient |

NSE: | Nash–Sutcliffe efficiency |

ENSE: | Error in Nash–Sutcliffe efficiency |

ANNs: | Artificial neural networks |

WOA: | Whale optimization algorithm |

: | Intrinsic growth rate of local prey |

: | Increasing rate of immigrant prey |

, : | Force of interaction between local prey and immigrant prey |

: | Carrying capacity of local prey |

: | Carrying capacity of immigrant prey |

, : | Intrinsic growth rate of the predator population |

, : | Suffering loss of the predator population |

: | Catching rate of local prey |

: | Catching rate of immigrant prey. |

#### Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest regarding this study.