1. Introduction

2. Preliminary models

2.1. One prey One predator Model

The most basic model assumes exponential growth of the prey in the absence of the predator. In the case of the presence of a predator, the growth decreases naturally and the decline is dependent on both the populations of the prey and the predator as more prey implies more food for the predators implying rise in predator population, which in turn results in decline of prey population. The same idea holds true for the predator except the sign changes as more prey implies more food for the predators including its reproduction rate and further we will have to eliminate the death rate. Putting all the parameters in the form of equations, we will have the following system of equations for prey and predator population respectively.

$$\begin{array}{c}\hfill {\displaystyle \frac{dx}{dt}}=\alpha x-\beta xy\end{array}$$

(1)

$$\begin{array}{c}\hfill {\displaystyle \frac{dy}{dt}}=-\delta y+\gamma xy\end{array}$$

(2)

The plots we get on solving the equations are periodic which further implies the closed structure of the phase plot. Depending on the initial conditions, which in this case is the initial value of y, we get concentric closed figures around the equilibrium point in the phase, which in this case is a center (The equilibrium occurs when both $x,y=0$).

Figure 1. The plot of x vs t on solving the Lotka Volterra equation.

Figure 1. The plot of x vs t on solving the Lotka Volterra equation.

Figure 2. The plot of y vs t on solving the Lotka Volterra equation.

Figure 2. The plot of y vs t on solving the Lotka Volterra equation.

Figure 3. The phase plot of y vs x on solving the Lotka Volterra equations.

Figure 3. The phase plot of y vs x on solving the Lotka Volterra equations.

Table 1. For the simulation we took the following values into consideration for equations 1 and 2.

Table 1. For the simulation we took the following values into consideration for equations 1 and 2.

Variables	Values
$\alpha$	2
$\beta$	1.1
$\gamma$	1
$\delta$	0.9

2.2. Modified Model

Volterra had given a general theory on species. For two species, the modified equations are

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dx}{dt}}=\alpha x-\beta xy-\lambda x^2\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dx}{dt}}=-\delta x+\gamma xy-\mu y^2\hfill \end{array}$$

(4)

The plots we get for the same values for the variables additionally considering $\lambda =\mu =0.1$ show damping plots. The phase plot we get on solving the equations is a spiral with focus at the equilibrium point. This gives us a damped oscillation of the populations for both the species. The equilibrium here forms a stable focus. The different spirals are for the different starting values of y.

Figure 4. The plot of x vs t on solving the modified Lotka Volterra equation.

Figure 4. The plot of x vs t on solving the modified Lotka Volterra equation.

Figure 5. The plot of y vs t on solving the modified Lotka Volterra equation.

Figure 5. The plot of y vs t on solving the modified Lotka Volterra equation.

Figure 6. The plot of x vs t on solving the modified Lotka Volterra equation.

Figure 6. The plot of x vs t on solving the modified Lotka Volterra equation.

2.3. Depth three prey predator model

The depth three food chain includes one predator, one prey and a lowest level species which is fed on by the prey. This may include producers like carrots for food chains like Fox → Rabbit → Carrots. Simple models considering functional responses for both species have been studied before. Incorporating type II functional responses for both the prey and the predator, and considering logistic growth of the lowest level in the absence of the prey, we get the following model.

$$\begin{array}{cc}& {\displaystyle \frac{dX}{dT}}=R_{0}X(1-X/K_0)-C_1{F}_1\left(X\right)Y\hfill \end{array}$$

(5)

$$\begin{array}{cc}& {\displaystyle \frac{dY}{dT}}=F_1\left(X\right)Y-F_2\left(Y\right)Z-D_{1}Y\hfill \end{array}$$

(6)

$$\begin{array}{cc}& {\displaystyle \frac{dZ}{dT}}=C_2{F}_2\left(Y\right)Z-D_{2}Z\hfill \end{array}$$

(7)

where

$$F_i\left(U\right)={\displaystyle \frac{A_{i}U}{B_i+U}}\mathrm{for}i=1,2$$

The parameters can be reduced which will further give us the equations

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dx}{dt}}=x(1-x)-f_1\left(x\right)y\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dy}{dt}}=f_1\left(x\right)y-f_2\left(y\right)z-d_{1}y\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dz}{dt}}=f_2\left(y\right)z-d_{2}z\hfill \end{array}$$

(10)

where

$$f_i\left(u\right)={\displaystyle \frac{a_{i}u}{1+b_{i}u}}\mathrm{for}i=1,2$$

considering the reductions

$$\begin{array}{cc}\hfill & x={\displaystyle \frac{X}{K_0}}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & y={\displaystyle \frac{C_{1}Y}{K_0}}\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill & z={\displaystyle \frac{C_{1}Z}{C_2{K}_0}}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & t=R_{0}T\hfill \end{array}$$

(14)

Table 2. The authors set the values, with the intention to consider biologically reasonable food chain and check if chaotic behaviors are likely to exist for the equations 8, 9 and 10.

Table 2. The authors set the values, with the intention to consider biologically reasonable food chain and check if chaotic behaviors are likely to exist for the equations 8, 9 and 10.

Variables	Values
$a_1$	5.0
$b_1$	2.0 to 6.0
$a_2$	0.1
$b_2$	2.0
$d_1$	0.4
$d_2$	0.01

Although we get to see some kind of cyclic plots, the nonuniform behavior of the plots suggests chaotic nature of the model. On varying the value of $b_1$, we start to observe chaos. The plot of the peak value of z against the values of $b_1$ gives us the bifurcation diagram as in Figure 10.

Figure 7. The plot of x vs t for the depth 3 prey predator model.

Figure 7. The plot of x vs t for the depth 3 prey predator model.

Figure 8. The plot of y vs t for the depth 3 prey predator model.

Figure 8. The plot of y vs t for the depth 3 prey predator model.

Figure 9. The plot of z vs t for the depth 3 prey predator model.

Figure 9. The plot of z vs t for the depth 3 prey predator model.

Figure 10. The plot of the peak value of z in each cycle against $b_1$ for the depth 3 prey predator model.

Figure 10. The plot of the peak value of z in each cycle against $b_1$ for the depth 3 prey predator model.

Figure 11. The phase plot for the depth 3 prey predator model considering $b_1=3$.

Figure 11. The phase plot for the depth 3 prey predator model considering $b_1=3$.

2.4. Two or more Predators and One prey model

The authors consider a modified model for the two predator one prey system, considering the Leslie Gower equations with Holling type I functional response. The equations are

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dx}{dt}}=rx(1-x/k)-{\beta }_{1}x{y}_1-{\beta }_{2}x{y}_2\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{y}_1}{dt}}=y_1{r}_1(1-y_1/x)-{\alpha }_1{y}_1{y}_2\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{y}_2}{dt}}=y_2{r}_2(1-y_2/x)-{\alpha }_2{y}_1{y}_2\hfill \end{array}$$

(17)

However the authors of [3] consider a different model ( Hanski and Korpimaki), considering three predators for the Snowshoe Hares, Lynx, Coyote and Great Horned Owl. The plots for this model are according to the given values in Table 3. The authors suggest the following model for single prey-predator relationship

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dN}{dt}}=rN\left(1-{\displaystyle \frac{N}{k}}\right)-{\displaystyle \frac{\gamma N^2}{N^2+{\eta }^2}}-{\displaystyle \frac{\alpha NP}{N+\mu }}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dP}{dt}}=sp\left(1-{\displaystyle \frac{qP}{N}}\right)\hfill \end{array}$$

(19)

which can be generalized for arbitrary number of predators as

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dN}{dt}}=rN\left(1-{\displaystyle \frac{N}{k}}\right)-{\displaystyle \frac{\gamma N^2}{N^2+{\eta }^2}}-\sum _i{\displaystyle \frac{{\alpha }_{i}N{P}_i}{N+{\mu }_i}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{P}_i}{dt}}=s_i{p}_i\left(1-{\displaystyle \frac{q_i{P}_i}{N}}\right)\hfill \end{array}$$

(21)

Figure 12. The plot of Hare population for the model in equations 20 and 21.

Figure 12. The plot of Hare population for the model in equations 20 and 21.

Figure 13. The plots of Horned owl (yellow), coyote (red) and lynx (blue) for the model in equations 20 and 21.

Figure 13. The plots of Horned owl (yellow), coyote (red) and lynx (blue) for the model in equations 20 and 21.

3. Canadian Lynx-Snowshoe Hares population data

Figure 14. The plots of Lynx and Hare population (Source: http://people.whitman.edu/~hundledr/courses/M250F03/LynxHare.txt).

Figure 14. The plots of Lynx and Hare population (Source: http://people.whitman.edu/~hundledr/courses/M250F03/LynxHare.txt).

4. Conclusion

References

1. Zhang, Zhibin & Tao, yi & Li, Zhenqing. (2007). Factors affecting hare lynx dynamics in the classic time series of the Hudson Bay Company, Canada. Climate Research. 34. 83-89. [CrossRef]
2. Hastings, A., & Powell, T. (1991). Chaos in a three-species food chain. Ecology, 72(3), 896-903. [CrossRef]
3. Tyson, R., Haines, S., & Hodges, K. E. (2010). Modelling the Canada lynx and snowshoe hare population cycle: the role of specialist predators. Theoretical Ecology, 3, 97-111. [CrossRef]
4. Savitri, D., Suryanto, A., & Kusumawinahyu, W. M. (2019, June). A dynamics behaviour of two predators and one prey interaction with competition between predators. In IOP Conference Series: Materials Science and Engineering (Vol. 546, No. 5, p. 052069). IOP Publishing.
5. harles J. Krebs, Rudy Boonstra, Stan Boutin, A.R.E. Sinclair, What Drives the 10-year Cycle of Snowshoe Hares? The ten-year cycle of snowshoe hares—one of the most striking features of the boreal forest— is a product of the interaction between predation and food supplies, as large-scale experiments in the yukon have demonstrated, BioScience, Volume 51, Issue 1, January 2001, Pages 25–35. [CrossRef]
6. Krebs, Charles & Bryant, John & Kielland, Knut & Donoghue, Mark & Doyle, Frank & Carrière, Suzanne & Difolco, Donna & Berg, Nathan & Boonstra, Rudy & Boutin, Stan & Kenney, Alice & Reid, Donald & Bodony, Karin & Putera, Judy & Timm, Henry & Burke, Toby & Maier, Julie & Golden, Howard. (2014). ARTICLE What factors determine cyclic amplitude in the snowshoe hare (Lepus americanus) cycle?. Canadian Journal of Zoology. 92. 1039-1048. [CrossRef]
7. Liu, R., Gourley, S.A., DeAngelis, D.L., and Bryant, J.P. 2012. Modeling the dynamics of woody plant–herbivore interactions with age-dependent toxicity. J. Math. Biol. 65(3): 521–552. [CrossRef] [PubMed]
8. Stenseth, N. C., Falck, W., Bjørnstad, O. N., & Krebs, C. J. (1997). Population regulation in snowshoe hare and Canadian lynx: asymmetric food web configurations between hare and lynx. Proceedings of the National Academy of Sciences, 94(10), 5147-5152. [CrossRef]