This is about the wolves and moose of Isle Royale, Michigan:
http://www.isleroyalewolf.org/http%3A//www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/nao.shtml

We replicate the InsightMaker model which is available at
https://insightmaker.com/insight/2068/Isle-Royale-Predator-Prey-Interactions

Andy Long
Spring, 2018

Thanks to:
https://mathematica.stackexchange.com/questions/74214/the-lotka-volterra-predator-prey-model

Perhaps helpful stuff.... from other sources: http://reference.wolfram.com/language/tutorial/NDSolveDelayDifferentialEquations.html

Clear everything possibly in our way:

Clear["Global`*"]

Now let’s get the actual wolf/moose data, extracted from a much more extensive data set, available from 
http://www.isleroyalewolf.org/data/data/home.html

IsleRoyale=Import["http://ceadserv1.nku.edu/longa//classes/mat375/data/isleroyale_graph_data_28Dec2011.csv"];

A look at the actual data:

year=Rest[IsleRoyale[[All,1]]];
minYear=Min[year]
year=year-minYear;
wolves=Rest[IsleRoyale[[All,2]]];
moose=Rest[IsleRoyale[[All,3]]];
ListPlot[Transpose[{year+minYear,moose}]]
ListPlot[Transpose[{year+minYear,wolves}]]
(* Do you see a ring-around-the-rosey? *)
ListPlot[Transpose[{moose,wolves}]]

Now let’s make some moose and wolves! 
Model and starting parameters from the InsightMaker model:

WolfDeathRate=0.25;
WolfBirthRateFactor=0.000203;
MooseBirthRate=0.25;
MooseDeathRateFactor=0.0081;
mo=421.5;
wo=20;

I got these starting parameters by doing optimization of the fit of the ODE solution to the ratio of wolves to moose (third manipulation, below):

{WolfBirthRateFactorStart,WolfDeathRateStart,MooseBirthRateStart,MooseDeathRateFactorStart,moStart,woStart}={0.000203,0.25,0.25,0.0081,421.5,20.}; (* Fit: 0.0.0532418 *)
{WolfBirthRateFactorStart,WolfDeathRateStart,MooseBirthRateStart,MooseDeathRateFactorStart,moStart,woStart}={0.00023195,0.239829,0.269769,0.0113139,660.,20.}; (* Fit: 0.00806343 *)
{WolfBirthRateFactorStart,WolfDeathRateStart,MooseBirthRateStart,MooseDeathRateFactorStart,moStart,woStart}={0.000267409,0.239821,0.269755,0.0113679,591.,23.};(* Fit: 0.00690738*)
(* Fit: 0.00685681, from NonLinearRegression version *){WolfBirthRateFactorStart,WolfDeathRateStart,MooseBirthRateStart,MooseDeathRateFactorStart,moStart,woStart}={0.0005678933580857228`,0.47741786513634954`,0.1336612038887985`,0.005793336215934526`,656.3973146983205`,24.384517448668483`}

Manipulate[
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,tend}];
p1=ListPlot[Table[Style[{year[[t]],moose[[t]]},Hue[t/Length[moose]]],{t,1,Length[moose],1}],PlotStyle->PointSize[Medium]];
p3=Plot[Evaluate[{m[t]}/.lv],{t,0,tend},
PlotRange->All,AxesOrigin->{0,0}];
Show[p3,p1],
{{lv,{}},None},
{{tend,52},20,2000,10},
{{WolfDeathRate,WolfDeathRateStart},.01,.90,.01},
{{WolfBirthRateFactor,WolfBirthRateFactorStart},.000001,.0009,.0001},
{{MooseBirthRate,MooseBirthRateStart},.01,.90,.01},
{{MooseDeathRateFactor,MooseDeathRateFactorStart},.001,.019,.001},
{{wo,woStart},2,100,1},
{{mo,moStart},100,3000,10}
]

Manipulate[
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,tend}];
p2=ListPlot[Table[Style[{year[[t]],wolves[[t]]},Hue[t/Length[moose]]],{t,1,Length[moose],1}],PlotStyle->PointSize[Medium]];
p3=Plot[Evaluate[{w[t]}/.lv],{t,0,tend},
PlotRange->All,AxesOrigin->{0,0}];
Show[p3,p2],
{{lv,{}},None},
{{tend,52},20,2000,10},
{{WolfDeathRate,WolfDeathRateStart},.01,.90,.01},
{{WolfBirthRateFactor,WolfBirthRateFactorStart},.000001,.0009,.0001},
{{MooseBirthRate,MooseBirthRateStart},.01,.90,.01},
{{MooseDeathRateFactor,MooseDeathRateFactorStart},.001,.019,.001},
{{wo,woStart},2,100,1},
{{mo,moStart},100,3000,10}
]

Manipulate[
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,tend}];
p1=ListPlot[Table[Style[{year[[t]],wolves[[t]]/moose[[t]]},Hue[t/Length[moose]]],{t,1,Length[moose],1}],PlotStyle->PointSize[Medium]];
p3=Plot[Evaluate[{w[t]/m[t]}/.lv],{t,0,tend},
PlotRange->All,AxesOrigin->{0,0}];
Show[p3,p1],
{{lv,{}},None},
{{tend,52},20,2000,10},
{{WolfDeathRate,WolfDeathRateStart},.01,.90,.01},
{{WolfBirthRateFactor,WolfBirthRateFactorStart},.000001,.0009,.0001},
{{MooseBirthRate,MooseBirthRateStart},.01,.90,.01},
{{MooseDeathRateFactor,MooseDeathRateFactorStart},.001,.019,.001},
{{wo,woStart},2,100,1},
{{mo,moStart},100,3000,10}
]

Manipulate[
lv=
NDSolve[
{(*  moose: *)
m'[t]==MooseBirthRate m[t]-MooseDeathRateFactor m[t]w[t] ,
(* wolves: *)
w'[t]==WolfBirthRateFactor m[t] w[t]-WolfDeathRate w[t],
m[0]==mo,
w[0]==wo
},
{m[t],w[t]},
{t,0,tend}];
wolfFn[x_]:=First[First[MapThread[#1/.lv/.t->#2&,{({w[t]}),{x}}]]];
mooseFn[x_]:=First[First[MapThread[#1/.lv/.t->#2&,{({m[t]}),{x}}]]];
p1=ListPlot[Table[Style[{moose[[t]],wolves[[t]]},Hue[t/Length[moose]]],{t,1,Length[moose],1}],PlotStyle->PointSize[Medium]];
p3=ParametricPlot[{mooseFn[t],wolfFn[t]},{t,0,tend},
AspectRatio->1/2];
Show[p1,p3],
{{lv,{}},None},
{{tend,52},20,2000,10},
{{WolfDeathRate,WolfDeathRateStart},.01,.90,.01},
{{WolfBirthRateFactor,WolfBirthRateFactorStart},.000001,.0009,.0001},
{{MooseBirthRate,MooseBirthRateStart},.01,.90,.01},
{{MooseDeathRateFactor,MooseDeathRateFactorStart},.001,.019,.001},
{{wo,woStart},2,100,1},
{{mo,moStart},100,3000,10}
]