Equation Setup

deqs = a list of the right-hand side of the differential equations
vars = the variables, in the same order as the deqs
start = initial conditions for the variables
paramset = parameter values
tmax = maximum time to plot time series

(* selection gradient = d(dI/dt)/dv *)
D[((βmax v)/(v+βh))s i-(c+m+v)i,v]

deqs={
b(s+i)-((βmax v)/(v+βh))s i-m i+c s,
((βmax v)/(v+βh))s i-(c+m+v)i,
Var (-i-((i s v βmax)/(((v+βh))^2))+((i s βmax)/(v+βh)))
};
vars={s,i,v};
start={0.1,0.5,0.2};
paramset={b->2,m->0.5,c->0.5,βmax->1.7,βh->0.7,Var->1.0};
tmax=4;

Plot of Density vs. Time

(* find time-solution *)
addT=Module[{i},Table[vars[[i]]->Map[vars[[i]],{t}][[1]],{i,Length[vars]}]];
outset=NDSolve[Join[Thread[D[vars/.addT,t]==(deqs/.addT)],Thread[vars==start]/.addT/.t->0]/.paramset,vars,{t,0,tmax},MaxSteps->500*tmax,AccuracyGoal->Infinity];
view:=tmax;
pA=Plot[Evaluate[vars/.addT/.outset[[1]]],{t,tmax-view,tmax},PlotStyle->PlotStyleFunction,PlotRange->All,Frame->True,FrameLabel->{t}, PlotLegends->Placed[vars,Above]]

v[tmax]/.outset

(* Theoretically, this is the optimal v from R0 *)
(* the simulation will get there with a little more time *)
Sqrt[βh(c+m)]/.paramset