This is the basic model of Tyson et al., built to reproduce the code that we assembled in InsightMaker.

Reference: Tyson, Rebecca, Sheena Haines,  Karen Hodges. Modelling the Canada lynx and snowshoe hare population cycle: The role of specialist predators. Theoretical Ecology. 3, 97–111 (2010).  https://doi.org/10.1007/s12080-009-0057-1

This is the documentation that I included in that InsightMaker model (https://insightmaker.com/insight/188321/Basic-Model-Tyson-Lynx-and-Hare#):

The basic model of Tyson, et al. demonstrates logistic growth in prey, and in predator (with prey dependence for carrying capacity).

The differential equations looks like

x’(t)=rx(1-x/K)-gamma x^2/(x^2+eta^2) - alpha y x/(x+mu)

where K is the carrying capacity of the prey, in the absence of predation.

y’(t)=sy(1- y/(x/q))

is the predator (lynx) equation. The growth term suggests exponential growth, but there is a loss term of the form sy^2/(x/q)) -- loss is proportional to population (crowding), and inversely proportional to prey density.

In this model, I scale the second equation by multiplying by q, then replace y by qy. This requires a slight change in the prey equation -- alpha replaced by the ratio of alpha/q. 

Tyson, et al. took q to be about 212 for hare and lynx -- so that it requires about 212 hare to allow for one lynx to survive at “equilibrium”. However, when alpha -- the hares/lynx/year -- gets sufficiently large (e.g. 1867 -- and that does seem like a lot of hares per lynx per year...:), limit cycles develop (rather than a stable equilibrium).

Author: Andy Long, Northern Kentucky University (2020)

Reference: Tyson, Rebecca, Sheena Haines,  Karen Hodges. Modelling the Canada lynx and snowshoe hare population cycle: The role of specialist predators. 

Theoretical Ecology. 3, 97–111 (2010). https://doi.org/10.1007/s12080-009-0057-1

Clear["Global`*"]

Manipulate[
lv=
NDSolve[
{(* x is prey: *)
x'[t]==r  x[t](1-x[t]/K) 
- gamma  x[t]^2/(x[t]^2+eta^2)
-(alpha/q) x[t]* y[t]/(x[t]+mu),
(* y is predator: *)
y'[t]==s  y[t](1- y[t]/x[t]),
x[0]==popinitial,
y[0]==popinitial
},
{x[t], y[t]},
{t,0,tend}];
popvals=Table[Evaluate[{x[t],y[t]}/.lv],{t,0,tend,.1}];
p1=ListPlot[popvals,
PlotLabel->"Phase Portrait",
AxesLabel->{"Hare","q*Lynx"}];
timeseries=Table[f,{f,{Evaluate[x[t]/.lv][[1]]/K,Evaluate[y[t]/.lv][[1]]/K}},{t,0,tend,.1}];
p2=ListLinePlot[timeseries,
DataRange->{0,tend},
PlotLegends->Placed[{"Hare","q*Lynx"},Below],
PlotLabel->"Populations, relative to K",
AxesLabel->{"Years"}
];
Show[GraphicsGrid[{{p1},{p2}}]],
{{lv,{}},None},
{{tend,28},2,100,1},
{{popinitial,1.7},.1,10,.1},
{{r,8},0.1,20,.1},
{{eta,1.25},.5,5,.1},
{{mu,0.3},.01,1,.01},
{{gamma,.1},.01,1,.01},
{{alpha,1867},100,10000,1},
{{q,212},1,600,1},
{{s,0.85},0.1,2,.1},
{{K,1.75},1,6,.1}
]