Order and Chaos in a Prey-Predator Model Incorporating Refuge, Disease, and Harvesting

In this paper, a mathematical model consisting of a prey-predator system incorporating infectious disease in the prey has been proposed and analyzed. It is assumed that the predator preys upon the nonrefugees prey only according to the modified Holling type-II functional response. There is a harvesting process from the predator. The existence and uniqueness of the solution in addition to their bounded are discussed. The stability analysis of the model around all possible equilibrium points is investigated. The persistence conditions of the system are established. Local bifurcation analysis in view of the Sotomayor theorem is carried out. Numerical simulation has been applied to investigate the global dynamics and specify the effect of varying the parameters. It is observed that the system has a chaotic dynamics.


Introduction
The interaction between a predator and their prey is one of the most important interaction types in ecology as well as applied mathematics due to its wide existence in our real life. Although such interaction models may appear to be simple mathematically, they are very challenging and complicated [1]. Since the classical Lotka (1925) and Volterra (1926) model on prey-predator interaction, which was developed on the basis of chemical principle of mass action, many research works have been done in literature. Later, the simple growth and decay terms in the Lotka-Volterra model were replaced by more complex functions such as the logistic growth; prey-dependent functional responses such as Holling's I, II, and III types; ratio-dependent functional response, and predator-dependent type of functional response, see, for example, [2][3][4][5][6][7][8] and the references therein. It is well known that the two species models are very restrictive in their dynamics and can give rise to only two basic dynamics approaches to equilibrium or to a limit cycle, while ecological situations demonstrate very complex dynamics in nature. Accordingly, multispecies systems comprising of three spe-cies or more was considered to study the complex ecological community [9][10][11][12][13]. Furthermore, due to the existence of different important factors in the environment, such as stage structure, refuge, cannibalism, harvesting, toxicant, and scavenger, different prey-predator models incorporating such factors have been proposed and studied in literature, see, for example, [14][15][16][17].
On the other hand, the spread of diseases between the individuals of species is an important subject for study, which is known as an epidemiology, due to the existence of a variety of diseases in the environment. It is well known that after the SIR epidemic model of , the field of epidemiology has received a lot of attention from the researchers, see [19] and the references therein. Recently, Zhou and Yao [20] proposed and studied a host-vector epidemic model with a stage structure for the vector. Further articles have been done also in the field of epidemiology. In fact, recently, there is increasing interest from the scientists by combining ecology and epidemiology in one field that is called ecoepidemiology, which has become a major field of study, see, for example, [21][22][23][24] and the references therein. Recently, Wang et al. [25] proposed and studied a generalized ecoepidemiological system with prey refuge. The saturation incidence kinetics and a generalized functional response are used to describe the contact process and the predation process, respectively. They show that the prey refuge can control the spread of disease by the relative level of prey refuge.
Keeping the above studies in view, an ecoepidemiological model is proposed consisting of prey-predator incorporating a prey refuge, infectious disease of SI type in prey species only, and the predator preys upon the prey (susceptible and infected) according to the modified Holling type-II functional responses. It is assumed further that there is a harvesting process on the predator species due to external resources.

Model Formulation
In this section, a prey-predator system incorporating disease in a prey is formulated mathematically. It is assumed that there is a refuge for prey to protect themselves against predation and the harvesting process on predator. The following hypotheses are adopted in constricting the mathematical model that describes the dynamics of the above real-world ecoepidemiological system.
(1) In the presence of infectious disease, the prey population is divided into two compartments, namely, susceptible prey that is denoted to their density at time t by SðTÞ and infected prey, which is denoted to their density at time T by IðTÞ. Therefore, at time T, the total prey population density is XðTÞ = SðTÞ + IðTÞ (2) It is assumed that in the absence of predator that is denoted to their density at time T by YðTÞ, the susceptible prey grows logistically with carrying a capacity K > 0 and an intrinsic growth rate r > 0, while the infected prey cannot reproduce and rather than that, it competes for food and space (3) The availability of refuge for the prey in the environment keeps m ∈ ð0, 1Þ portion of them safe from predation and left ð1 − mÞ from the total density of the prey facing predator (4) The disease transfers from infected to susceptible individuals by contact with an infected rate given by β > 0 (5) There is a harvesting process on the predator, due to an external force, represented by a rate ðl 1 EÞ/ðl 2 E + l 3 YÞ, where E > 0 is the effort applied to harvest the prey; l 1 > 0 is the catchability coefficient while l 2 > 0 and l 3 > 0 are suitable constants (6) The predator consumed the prey according to the modified Holling type-II with attack rates a 1 > 0 and a 2 > 0 for susceptible prey and infected prey, respectively, half saturation constant b > 0, and infected prey's preference rate c > 0, while the conversion rate is given by e > 0 (7) Finally, there is a disease death rate given by d 1 > 0 and the natural death rate for predator d 2 > 0 According to the above hypotheses, the dynamics of the above described ecoepidemiological real-world system can be represented by the following set of differential equations.
with initial data Sð0Þ ≥ 0, Ið0Þ ≥ 0, and Yð0Þ ≥ 0. Clearly, the interaction functions in system (1) are continuously differential functions on the positive octant. Hence, they are Lipschitzian. Therefore, the solution of system (1) exits and is unique. Now, for further simplifying the above system, the following dimensionless variables and constants are used.

Stability and Persistence
In this section, the stability analysis of all possible equilibrium points of system (3) is investigated. The persistence conditions of system (3) are established. Now, straightforward computation shows that the system (1) has at most five nonnegative equilibrium points. The existence conditions for each of them can be summarized as follows: The trivial equilibrium point (TEP), say E 0 = ð0, 0, 0Þ always exists. The disease-predator-free equilibrium point (DPFEP), namely, E 1 = ðx ′ , 0, 0Þ = ð1, 0, 0Þ always exists. The predator-free equilibrium point (PFEP), namely, E 2 = ð x, y, 0Þ = ðððw 4 Þ/ ðw 1 ÞÞ, ððw 1 − w 4 Þ/ðw 1 ð1 + w 1 ÞÞÞ, 0Þ exists in the xy-plane provided that The disease-free equilibrium point (DFEP) is denoted by E 3 = ðx, 0,ẑÞ wherê whilex is a positive root of the following third-order equation: with Hence, E 3 exists uniquely in the positive quadrant of xz-plane if one of the following sets of conditions holds. α 1 > 0 and α 2 > 0 α 1 > 0 and α 3 < 0 ð10Þ The coexistence equilibrium point (CEP) is denoted by while x * is a positive root of the following third-order equation.
Here, 3 Journal of Applied Mathematics with Clearly, the CEP that is denoted by E 4 exists uniquely in the interior of the positive octant provided that the following condition holds.
with one set of the following sets of conditions σ 1 > 0, σ 4 < 0, σ 2 > 0 σ 1 > 0, σ 4 < 0, σ 3 < 0 In the following, the local stability analysis of the above equilibrium points is carried out. The Jacobian matrix of system (1) at point ðx, y, zÞ can be written as Here, Therefore, the eigenvalues of the Jacobian matrix at E 0 are Hence, the TEP is a saddle point with unstable manifold in x-direction and stable manifold in the yz-plane.
The Jacobian matrix of system (3) at the DPFEP that is represented by E 1 = ð1, 0, 0Þ is given by Therefore, the eigenvalues are written as It is clear that all the eigenvalues are negative and then, the DPFEP is locally asymptotically stable provided that the following conditions are hold.
The Jacobian matrix of system (3) at the PFEP is denoted by E 2 = ð x, y, 0Þ.

4
Journal of Applied Mathematics Clearly, the characteristic equation of JðE 2 Þ can be written as Accordingly, all the eigenvalues of JðE 2 Þ have negative real parts and then, the PFEP is locally asymptotically stable provided that the following condition holds The Jacobian matrix of system (3) at the DFEP that is given by E 3 can be written as Hence, the characteristic equation of JðE 3 Þ can be written as Straightforward computations show that all the roots (eigenvalues) of the characteristic equation (28) have negative real parts, and then, the DFEP is locally asymptotically stable provided that the following conditions hold.
Now, the local stability conditions for the CEP are established in the following theorem.

Theorem 2.
Suppose that the CEP exists; then, it is locally asymptotically stable in the ℝ 3 + if the following sufficient conditions hold.
where all symbols are clearly described in the proof.

Journal of Applied Mathematics
Proof. According to Equation (17), the Jacobian matrix at E 4 is given by with a ij = u ij ðx * , y * , z * Þ for all i, j = 1, 2, 3. Therefore, the characteristic equation of JðE 4 Þ can be written as where D 1 = − a 11 + a 22 + a 33 ½ , D 2 = a 11 a 22 − a 12 a 21 + a 11 a 33 − a 13 a 31 + a 22 a 33 − a 23 a 32 , Moreover, it is easy to verify that Recall that according to Routh-Hurwitz criterion, Equation (32) has roots (eigenvalues) with negative real parts if and only if D 1 > 0, D 3 > 0 and D 1 D 2 − D 3 > 0. Now, straightforward computation shows that all the requirements of the Routh-Hurwitz criterion are satisfied provided that the conditions (30a), (30b), (30c), (30e), and (30f) hold. Hence, the CEP is locally asymptotically stable and the proof is complete.
It is well known that the persistence of an ecological system is an important subject for both environment and ecologists. The coexistence of all the species along the entire time life is known as persistence. Mathematically, this is equivalent to unapproachable of the trajectory to any of the boundary planes of axis. It is well known that system (3) has two subsystems. The first subsystem exists in the x y-plane, in case in the absence of the predator, and can be represented in system (36). However, the second subsystem exists in the xz-plane, in the absence of disease, and can be represented in system (37).
Recall that according to Bendixson-Dulac Criteria [26], the planar system X′ = GðXÞ, where X = ðx 1 , x 2 Þ T and G = ðg, f Þ T has no periodic dynamics provided that there is a C 1 function φðx 1 , x 2 Þ (known as Dulac function) so that the expression ð∂ðφgÞÞ/ð∂xÞ + ð∂ðφf ÞÞ/ð∂yÞ has the same sign and does not equal zero almost everywhere in a simply connected region of the plane.
Consequently, by choosing the Dulac function ð1Þ/ðxyÞ, it is easy to verify that system (36) has no periodic dynamics in the interior of xy-plane, while with the Dulac function ð1Þ/ðxzÞ, system (37) has no periodic dynamics fall in the interior of xz-plane provided that the following condition holds.
Now, in order to established the persistence conditions for system (3), the average Lyapunov function will be used as shown in the following theorem.

Theorem 3.
Assuming that there is no periodic dynamics in the boundary planes; then, system (3) is persistent provided that the following conditions hold.
Proof. Consider the following function Pðx, y, zÞ = x q 1 y q 2 z q 3 , where q i , ∀i = 1, 2, 3 are positive constants to be determined. Then, we obtain 6 Journal of Applied Mathematics Now, according to the average Lyapunov function, the proof is done if we could prove that this function Λ is positive at each boundary equilibrium point. Note that at E 0 , it obtains So, for any positive constants q i , ∀i = 1, 2, 3 with q 1 sufficiently large compared with q 2 and q 3 , we get ΛðE 0 Þ > 0. At the DPFEP that is represented by Clearly, ΛðE 1 Þ > 0 under the conditions (39a) and (39b). Now, at the PFEP that is denoted by E 2 = ð x, y, 0Þ, it obtains Here, ΛðE 2 Þ > 0 under condition (39c). Finally, at the DFEP that is given by E 3 = ðx, 0,ẑÞ, we have Clearly, condition (39d) guarantees that ΛðE 3 Þ > 0. Hence, the system has no omega limit sets on the boundary planes and thus, the system is persists. Now, the global dynamics of system (3) is investigated using suitable Lyapunov functions as shown in the following theorem.
Theorem 4. Assuming that the DPFEP which is given by E 1 = ð1, 0, 0Þ is locally asymptotically stable, then, it is a globally asymptotically stable provided that Proof. Consider the positive definite function It is clear that L 1 : ℝ 3 + ⟶ ℝ and L 1 ð1, 0, 0Þ = 0 with L 1 ðx, y, zÞ > 0, for all ðx, y, zÞ ≠ ð1, 0, 0Þ. Now the derivative of L 1 with respect to time can be written, after some algebraic computation, as Therefore, by using conditions (45a) and (45b), it obtains ðdL 1 Þ/ðdtÞ which is negatively defined in ℝ 3 + . Thus, according to Lyapunov's second theorem, E 1 is globally asymptotically stable.

Theorem 5.
Assuming that the PFEP that is given by is locally asymptotically stable, then, it is globally asymptotically stable in ℝ 3 + provided that Proof. Consider the positive definite function It is clear that L 2 : ℝ 3 + ⟶ ℝ and L 2 ð x, y, 0Þ = 0 with L 2 ðx, y, zÞ > 0, for all ðx, y, zÞ ≠ ð x, y, 0Þ = ðððw 4 Þ/ðw 1 ÞÞ, ððw 1 − w 4 Þ/ðw 1 ð1 + w 1 ÞÞÞ, 0Þ. Then, the derivative of L 2 with respect to time can be written, after some algebraic steps as So by choosing the positive constants as C 1 = 1, C 2 = ðð1 + w 1 Þ/ðw 1 ÞÞ, and C 3 = ð1/w 5 Þ, the following is obtained 7 Journal of Applied Mathematics Clearly, under the condition (48), ðdL 2 /dtÞ is negatively semidefinite, because ðdL 2 /dtÞ = 0 if and only if x = x and z = 0 with any value of y. Therefore, according to Lyapunov's second theorem, E 2 is a stable point. Now, form the first equation of system (3), we have ðdx/dtÞ = x½1 − x − ð1 + w 1 Þy, which is equal to zero if and only if y = y. Hence, the set fE 2 g is the largest invariant set for which ðdL 2 /dtÞ = 0. Accordingly, the global asymptotic stability of the PFEP follows from LaSalle's invariance principle. This completes the proof.

Theorem 6.
Assuming that the DFEP that is given by E 3 = ðx, 0,ẑÞ is locally asymptotically stable, then, it is globally asymptotically stable in ℝ 3 + provided that the following sufficient conditions hold where all symbols are given in the proof.

Theorem 7.
Assuming that the positive equilibrium point E 4 = ðx * , y * , z * Þ is locally asymptotically stable, then, all the points that satisfy the following conditions belong to the basin of attraction of E 4 .
Journal of Applied Mathematics Proof. Consider the positive definite function It is clear that L 4 : ℝ 3 + ⟶ ℝ and L 4 ðx * , y * , z * Þ = 0 with L 4 ðx, y, zÞ > 0, for all ðx, y, zÞ ≠ ðx * , y * , z * Þ with x > 0, y > 0 and z > 0. Then, the derivative of L 4 with respect to time can be written, after some algebraic steps as Therefore, we obtain where with N 1 = w 2 + ð1 − mÞx + cð1 − mÞy, N 1 * = w 2 + ð1 − mÞx * + cð1 − mÞy * , and N 2 * = w 8 + z * . Clearly according to the condition (57a), then, q 11 > 0, while q 22 and q 33 are always positive. Therefore, by using the conditions (57b)-(57d), we obtain The first tow terms are clearly negative while the third term is clearly positive; therefore, the sufficient condition (57e) guarantees the negative definite of the derivative. Hence, the solution starting with any point satisfying conditions (57a)-(57e) approaches asymptotically to E 4 . So these points are belonging to the basin of attraction.

Bifurcation Analysis
Now, in order to compute the second derivative of the Jacobian matrix, system (3) is rewritten in the vector form as follows: dX/dt = FðXÞ, X = ðx, y, zÞ T , F = ðF 1 , F 2 , F 3 Þ T = ðxf 1 , yf 2 , zf 3 Þ T .(50) Therefore, according to the Jacobian matrix of system (3) at the point ðx, y, zÞ that is given by Equation (17), we obtain where V = ðv 1 , v 2 , v 3 Þ T be any vector and 9 Journal of Applied Mathematics Consequently, the conditions guarantee the occurrence of local bifurcations around the nonhyperbolic equilibrium points that are established in the following theorems.
Theorem 8. Assume that condition (22b) holds. Letting the parameter w 1 pass through the value w 1 * ≡ w 4 , then, system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occur Proof. Note that when the parameter w 1 passes through the positive value w 1 * ≡ w 4 , it is easy to verify that the equilibrium point E 1 becomes a nonhyperbolic point and the Jacobian matrix of system (3) at ðE 1 , w 1 * Þ can be written as Thus, J 1 has the following eigenvalues λ 11 ½1 = −1, λ 12
Since the partial derivative of vector field F with respect to the parameter w 1 is given by ∂F/∂w 1 = F w 1 = ð−xy, xy, 0Þ T , hence, we obtain F w 1 ð E 1 , w 1 * Þ = ð0, 0, 0Þ T . Therefore, Thus system (3) at E 1 with w 1 = w 1 * does not experience saddle-node bifurcation in view of the Sotomayor theorem [26]. Now, we have where DF w 1 represents the derivative of F w 1 with respect to X. In addition, by substituting the value of E 1 , w 1 * , V 1 , and Ψ 1 in Equation (63), then we get Accordingly, by the Sotomayor theorem, system (3) near the equilibrium point E 1 with w 1 = w 1 * possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.
Journal of Applied Mathematics Theorem 9. Assume that the parameter w 6 passes through the value w 6 * , where Then, system (3) near the PFEP undergoes a transcritical bifurcation but neither saddle-node bifurcation nor pitchfork bifurcation can occur provided that where Γ 1 , Γ 2 are given in the proof.
Proof. Note that, straightforward computation shows that as the parameter w 6 passes through the value w 6 * , the PFEP becomes a nonhyperbolic equilibrium point and the Jacobian matrix of system (3) at ðE 2 , w 6 * Þ can be written as where N 1 = w 2 + ð1 − mÞ x + cð1 − mÞ y. Hence, J 2 has three eigenvalues given by Þ/2Þ, and λ 23 ½2 = 0: Obviously, the first two eigenvalues, λ 21 ½2 and λ 22 ½2 , have a negative real part. Now, let V 2 = ðv 21 , v 22 , v 23 Þ T be the eigenvector of J 2 corresponding to the eigenvalue λ 23 ½2 = 0. Then, straightforward computation gives and v 23 represents any nonzero real number.
Since the partial derivative of vector field F with respect to the parameter w 6 is given by ð∂F/∂w 6 Þ = F w 6 = ð0, 0,−zÞ T , hence, we obtain F w 6 ð E 2 , w 6 * Þ = ð0, 0, 0Þ T . Therefore, Thus, system (3) at E 2 with w 6 = w 6 * does not experience saddle-node bifurcation in view of the Sotomayor theorem. Now, we have where DF w 6 represents the derivative of F w 6 with respect to X. In addition, by substituting the value of E 2 , w 6 * , V 2 , and Ψ 2 in Equation (63), then we get Obviously, Ψ 2 T ½D 2 Fð E 2 , w 6 * ÞðV 2 , V 2 Þ ≠ 0 under the condition (70), and hence by the Sotomayor theorem, system (3) near the equilibrium point E 2 with w 6 = w 6 * possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.

Journal of Applied Mathematics
Proof. Straightforward computation shows that as the parameter w 4 passes through the value w 4 * then, the DFEP becomes a nonhyperbolic equilibrium point and the Jacobian matrix of system (3) at ðE 3 , w 4 * Þ can be written as whereb ij = b ij for all i, j = 1, 2, 3 are given in Equation (28) withb 22 ðw 4 * Þ = 0. Hence, J 3 has three eigenvalues given by Clearly, the eigenvalues λ 31 ½3 and λ 33 ½3 have negative real parts under the conditions (29a) and (29b).  T that corresponds to the eigenvalue λ 32 ½3 = 0. Then, direct calculation shows that Ψ 3 = ð0, ψ 32 , 0Þ T , where ψ 32 is any nonzero real number.
Since the partial derivative of vector field F with respect to the parameter w 4 is given by ð∂F/∂w 4 Þ = F w 4 = ð0,−y, 0Þ T , hence, we obtain F w 4 ð E 3 , w 4 * Þ = ð0, 0, 0Þ T . Therefore, Thus, system (3) at E 3 with w 4 = w 4 * does not experience saddle-node bifurcation in view of the Sotomayor theorem. Now, we have where DF w 4 represents the derivative of F w 4 with respect to X. In addition, by substituting the value of E 3 , w 4 * , V 3 , and Ψ 3 in Equation (63), then, we get Obviously, Ψ 3 T ½D 2 Fð E 3 , w 4 * ÞðV 3 , V 3 Þ ≠ 0 under the condition (75), and hence, by the Sotomayor theorem, system (3) near the DFEP that is given by E 3 with w 4 = w 4 * possesses a transcritical bifurcation, while pitchfork bifurcation cannot occur and hence, the proof is complete.
Theorem 11. Assume that in addition to the conditions (30a), (30b), (30c), and (30e) the following conditions are satisfied a 12 a 31 a 32 < a 11 < a 12 a 21 a 22 : Then, as the parameter w 3 passes through the value w 3 * , system (3) near the CEP undergoes a saddle-node bifurcation, but neither transcritical bifurcation nor pitchfork bifurcation can occur provided that the following conditions hold where all the symbols are clearly described in the proof, while with a ij , ∀i, j = 1, 2, 3 are given in Equation (31) and N 1 * = Proof. Straightforward computation shows that the above given conditions guarantee that the parameter w 3 * is positive and the determinant (D 3 ) of the Jacobian matrix J 4 = JðE 4 , w 3 * Þ = ða ij ðw 3 * ÞÞ 3×3 that is given in Equation (32) is zero.

Numerical Simulation
In this section, system (3) is simulated numerically. The objectives understand the dynamical behavior of the solution of the system as the time increasing with the parameters varied and confirm our obtained analytical results. A variety of mathematical tools is used in this investigation including the bifurcation diagrams, chaotic attractors, phase portrait, and time series. Four steps of the predictor-corrector method are used to solve system (3), and then, the obtained trajectories are drawn using MATLAB version 6. Now, the following set of hypothetical parameters are used to solve system (3) numerically.
In order to investigate the behavior of the system (3) as a function of c, bifurcation diagram is drawn in Figure 1.
According to the bifurcation diagram given in Figure 1, system (3) has a complex dynamics including chaos. Moreover, the system enters to chaos through cascades of periodic doubling and then, after a specific value (c ≥ 0:58), the system loses the persistence, and the trajectory approaches asymptotically to PFEP in the xy-plane as shown in Figure 2 for typical values of c.
The dynamical behavior of system (3) as a function of refuge rate m in the range ð0,0:15Þ is investigated with the help of bifurcation diagram as presented in Figure 3.
Clearly, system (3) has chaotic dynamics where in between there are small periodic regions; moreover, the system approaches asymptotically to PFEP in the xy-plane as m increases above a specific value (m ≥ 0:12), see Figure 4 for typical values of m. The dynamics of system (3) as a function of w 1 in the range (0:3 ≤ w 1 ≤ 0:55) is shown in the bifurcation diagram given in Figure 5 and the typical attractors given in Figure 6.
Further investigation for the effect of varying w 2 in the range ð0:15, 0:21Þ on the dynamical behavior of system (3) is carried out using the bifurcation diagram and phase portraits with their time series at typical values of w 2 as shown in Figures 7 and 8, respectively.
According to the above figures, system (3) has a range of complex dynamics including chaos given by ð0:16, 0:21Þ, while it approaches asymptotically to PFEP otherwise.
Other parameters have been investigated too, and the obtained results are summarized in Table 1.

Discussion and Conclusion
In this paper, an ecoepidemiological model consisting of the prey-predator system with infectious disease in prey has been proposed and analyzed. As a defensive property against predation, it is assumed that the prey has a refuge and the predator preys upon the available prey individuals only according to the modified Holling type-II functional response. There is a harvesting process applied on the predator species only with a nonlinear harvesting rate. The existence, uniqueness, and boundedness of the solution of the system are studied. The stability analyses around all possible equilibrium points are investigated. The persistence conditions of the system are derived analytically and shown numerically. The occurrences of local bifurcation near the nonhyperbolic equilibrium points are investigated. Finally, numerical simulations are carried out to confirm our obtained analytical results and specify the control set of parameters using biologically feasible set of parameters given in (87). It is observed that the system is very sensitive to variation in the parameter set and it has a complex dynamics for a wide range of parameters including chaos. Different types of attracting sets and bifurcation diagrams are drawn to present the behavior of the system. Finally, in the following, we summarize the obtained numerical results depending on the set of data given by Equation (87).
(1) The system (3) is sensitive for varying in the preference rate of predation as shown in the bifurcation diagram represented by Figure 1, so that it has periodic dynamics in the interior of the positive octant and then periodic doubling leading to chaos, as shown in Figure 2, for the range c < 0:58. Otherwise,   (3) using the data (87) with varying a specific parameter.

The parameter
The range The dynamical behavior of system (3)