Exponential Time Differencing Method for Studying Prey-Predator Dynamic during Mating Period

This paper addresses a first numerical simulation to the nonlinear dynamic system of equations that describes the prey-predator model at the predator mating period. Some male species accompany the females during the mating period. In this case, both male and female feed on the same prey. The presented work shows the numerical solution for this specific case of the prey-predator mathematical model via an exponential time differencing method. In addition, the paper provides the biological implication of the solution.


Introduction
Many applications and phenomena in physics, biology, the medical field, and other fields can be described as a mathematical model in order to understand or predict the dynamic of these phenomena or applications [1][2][3]. The prey-predator model is a famous biomathematical model that has been investigated widely to study the dynamic of the ecology system and an intersection between prey and predator [4]. The Lotka-Volterra model is known as the predator-prey model for the study population process [5][6][7]. The model has been studied when considering the influence of diseases [8], prey refuge [9], food chain [10], or competitive populations [11]. The model also was developed to study the intersection between two different populations of predators and one prey [12,13]. Recently, the model was modified to study the intersection between two predators from the same species (female and male) and one prey. In some species, the male stays with the female all the time, and they feed together at the mating period. This model is given as follows [14]: where u and v are the total numbers of prey and predator, respectively. a 1 is the prey's growth rate. a 2 , a 3 , and a 4 are the prey's decay rate due to competition on food supply between the male and the female, between one predator and one prey, and between two predators and one prey, respectively. a 5 is the mortality rate for the predator. a 6 is the decay rate of the predator due to competition on the food supply between the male and the female. D 2 is the predator's diffusion term. In Reference [14], the system (1) was rescaled and became a dimensionless nonlinear system as follows: where k 1 = ffiffiffiffiffiffiffiffi ffi a 4 a 1 p /a 2 = a 3 /a 2 and k 2 = a 5 /ða 1 k 1 Þ. The mathematicians put forth more effort in the computational field to improve different techniques to be able to find the consistent solutions for realistic applications, for example, the geometric-qualitative method, semiapproximate method, qualitative method, numerical method, and analytical method. Analytical methods generate different solutions based on many factors such as used ansatz. There are many methods to obtain the analytical solutions such as the ðG ′ /G 2 Þ-expansion approach [15] and the direct algebraic method [16][17][18][19]. They are able to obtain the solutions without an initial condition. Each method can introduce new solutions to the system. The solutions of system (2) were found by two analytical methods, the ðG ′ /G 2 Þ-expansion approach and the direct algebraic method [15]. Since the analytical methods give different solutions, it is difficult to determine which solution predicts the applications. On the other hand, the numerical method for the initial value problem is the method that finds the approximate solution and converges to only one solution based on the initial condition that is obtained from the application, for example, the finite difference method [20][21][22], Adomian decomposition method (ADM) [23], Laplacian Adomian decomposition method (LADM), multistage differential transform method [24,25], and exponential time differencing method [26]. The numerical solution for system (2) will be introduced in this paper via the exponential time differencing method since the method is fourth-order accurate and is not local convergent comparing to ADM or LADM. The results either determine which analytical solution that predicts the problem or find a new solution.
The paper is organized as follows: the next section presents the numerical results and the scheme of the used method, the third section discusses the numerical solution from the biological perspective, the fourth section is a comparison with the previous studies for the aforementioned model, and the last section is the summary of our analysis and results.

Numerical Simulation
In this section, we give the approximate numerical solution of system (2). First, the system is discretized in space utilizing the Fourier spectral method. Hence, we obtain the following system of ordinary differential equations (ODEs) that depends only on time: where the parameter c = k 2 − k 2 k 1 is associated to the linear part, k represents the wave numbers, both termŝ represent the nonlinear forcing term, and fft is the MATLAB command that represents the fast Fourier transform (FFT). Second, we integrate the obtained uncoupled system of ODEs (3) in time with the aid of the fourth-order exponential time differencing method ETD4RK [27][28][29]. Applying the ETD4RK method to (3b) takes the formula a k,n =v k,n e cΔt/2 + e cΔt/2 − 1 À ÁĜ k,n c , b k,n =v k,n e cΔt/2 + e cΔt/2 − 1 À ÁĜ a k,n , t n + Δt/2 ð Þ c , C k,n = a k,n e cΔt/2 + e cΔt/2 − 1 where Δt represents the time step (a similar formula is obtained when applying the ETD4RK method to (3a) with c = 1).

Computational and Mathematical Methods in Medicine
The time integration starts with initial conditions (taken from [30]) given bŷ The results of our experiment are shown in figures. The solution is stable at 0 ≤ k 2 ≤ ð1 + ffiffi ffi 5 p Þ/2 based on the stability analysis in Reference [14]. A numerical solution is produced in Figure 1 for different value of k 2 . We realized that solutions are changed with k 2 . When k 2 = 1:5,U increases while V decreases until they reach the steady state. For k 2 = 0:1, the solution is periodic. The parameter k 2 = ða 2 a 5 Þ/ða 1 ffiffiffiffiffiffiffiffi ffi a 4 a 1 p Þ = ða 5 a 2 Þ/ða 1 a 3 Þ has direct relation with a 5 and a 2 and negative relation with a 1 , a 3 , and a 4 .

Biological Relevance of the Solution
A new term in the studied model is a 4 v 2 u. Thus, the numerical solution is used to study the effect of this term on a traditional prey-predator model. The solution of the traditional model is an undamped oscillation, and Lotka explained this behavior because of the law of mass action [5], while the presented model (2) has an oscillation solution with small k 2 at a bounded space and the solution converges to the steady state out of this space, because a 4 v 2 u is related to the interaction between prey and two predators at the mating period which disappears after this period and also happens in a specific location where the predators give birth.
The large rate of a 4 indicates that the predators feed on a large amount of prey in this period, which lead to a large size of the predator population than the prey population. This is not an optimal dynamic of population to live normally (see Figure 2 when a 4 = 10). Therefore, the predator population will decrease because of absence of a food source, while the prey population will grow again. The new term a 4 v 2 u has more effect than a 3 uv to decay the prey population as we see by comparing Figures 3 and 2.

Connection with the Previous Studies
The analytical solutions of Equations (2) was found by Aljahdaly and Alqudah [31] using ðG ′ /GÞ and the generalized auxiliary equation [15]. The solutions are Uðx, tÞ and Vðx, tÞ which denote the population of the prey and predator, respectively, in a single patch. The solutions predict the two following behaviors 4.1. The ðG ′ /GÞ-Expansion Method (Solution (1)). As a result of two predators' sexual interaction and the high prey consumption during the predator's mating period, V increases  (2)). The predators graze in the mating period where the prey density (U) is large. Then, U decrease and V increase. At a certain time, U will be very low, and the predators will either die or move to different patches. Then, the prey population will recover and grow again in the same location.

Exponential Time Differencing
Method. For a large value of k 2 , we obtain the numerical solution that converges to solution (1), while for small k 2 , the numerical solution approximates solution (2). The numerical solution gives the important reasons behind the two different solutions by the analytical methods.

Conclusion
The prey-predator system with the extra term a 4 v 2 u (1) is a new nonlinear dynamic system of equations in the literature. It has been solved analytically with two different methods which give two different solutions. The importance of this paper is that the numerical solution gives connection between the obtained analytical solutions from the biological perspective. The numerical solution predicts the dynamic of the studied model as well as proves the effect of the new extra term on the model.

Data Availability
The data supporting this research are from previously reported studies, which have been cited.