Evolutionary epidemiological model
Code for model in Duffy and Sivars-Becker, 2007, Ecology Letters. 
Written by Meghan Duffy(duffy@gatech.edu)

Clear[b,f,c,m,θ,v,n,Vc]

dsdt=b*(S+f*Q)*(1-c*(S+Q))-m*S-n*S-β*S*Q
dqdt=β*S*Q-θ*m*Q-n*Q-v*Q
dβdt1=-Vc*Q

tmax=200;

b=0.4; (* b is maximal birth rate *)
f=0.75; (* f is the proportional fecundity of infecteds, relative to uninfecteds *)
c=0.05; (* c is the strength of density dependence *)
m=0.03; (* m is the rate of selective fish predation *)
n=0.05; (* n is the rate of non-selective mortality *)
v=0.05; (* v is virulence mortality *)
θ=9; (* θ is the strength of selectivity of predation *) 
Vc=0.02; (* Vc is clonal variance; this is the parameter of greatest interest for this study *)

sol1=NDSolve[{S'[t]==b (S[t]+f Q[t]) (1-c (S[t]+Q[t]))-n S[t]-m S[t]-β1[t] S[t] Q[t],Q'[t]==β1[t] S[t] Q[t]-n Q[t]-θ m Q[t]-v Q[t],(β1)'[t]==(Piecewise[{{0, β1[t]<0}, {-Vc Q[t], 0≤β1[t]}}]),S[0]==1,Q[0]==0.001`,β1[0]==0.06`},{S,Q,β1},{t,0,tmax}];
plotsi1=Plot[Evaluate[{S[t],Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Blue,Red},PlotRange->All,ImageSize->{640,240}]
plotbeta1=Plot[Evaluate[β1[t]/.sol1],{t,0,tmax},ImageSize->{640,240},PlotStyle->{Black},PlotRange->All]
plotinf1=Plot[Evaluate[((Q[t])/(S[t]+Q[t]))/.sol1],{t,0,tmax},PlotStyle->{Black},PlotRange->All,ImageSize->{640,240}]
plotdens1=Plot[Evaluate[{S[t]+Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Black},PlotRange->All,ImageSize->{640,240}]

Vc=0.01;

sol1=NDSolve[{S'[t]==b (S[t]+f Q[t]) (1-c (S[t]+Q[t]))-n S[t]-m S[t]-β1[t] S[t] Q[t],Q'[t]==β1[t] S[t] Q[t]-n Q[t]-θ m Q[t]-v Q[t],(β1)'[t]==(Piecewise[{{0, β1[t]<0}, {-Vc Q[t], 0≤β1[t]}}]),S[0]==1,Q[0]==0.001`,β1[0]==0.06`},{S,Q,β1},{t,0,tmax}];
plotsi2=Plot[Evaluate[{S[t],Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Blue,Red},PlotRange->All,ImageSize->{640,240}]
plotbeta2=Plot[Evaluate[β1[t]/.sol1],{t,0,tmax},ImageSize->{640,240},PlotStyle->{Blue},PlotRange->All]
plotinf2=Plot[Evaluate[((Q[t])/(S[t]+Q[t]))/.sol1],{t,0,tmax},PlotStyle->{Blue},PlotRange->All,ImageSize->{640,240}]
plotdens2=Plot[Evaluate[{S[t]+Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Blue},PlotRange->All,ImageSize->{640,240}]

Vc=0.005;

sol1=NDSolve[{S'[t]==b (S[t]+f Q[t]) (1-c (S[t]+Q[t]))-n S[t]-m S[t]-β1[t] S[t] Q[t],Q'[t]==β1[t] S[t] Q[t]-n Q[t]-θ m Q[t]-v Q[t],(β1)'[t]==(Piecewise[{{0, β1[t]<0}, {-Vc Q[t], 0≤β1[t]}}]),S[0]==1,Q[0]==0.001`,β1[0]==0.06`},{S,Q,β1},{t,0,tmax}];
plotsi3=Plot[Evaluate[{S[t],Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Blue,Red},PlotRange->All,ImageSize->{640,240}]
plotbeta3=Plot[Evaluate[β1[t]/.sol1],{t,0,tmax},ImageSize->{640,240},PlotStyle->{Red},PlotRange->All]
plotinf3=Plot[Evaluate[((Q[t])/(S[t]+Q[t]))/.sol1],{t,0,tmax},PlotStyle->{Red},PlotRange->All,ImageSize->{640,240}]
plotdens3=Plot[Evaluate[{S[t]+Q[t]}/.sol1],{t,0,tmax},PlotStyle->{Red},PlotRange->All,ImageSize->{640,240}]

Vc=0.0005;

sol4=NDSolve[{S'[t]==b (S[t]+f Q[t]) (1-c (S[t]+Q[t]))-n S[t]-m S[t]-β[t] S[t] Q[t],Q'[t]==β[t] S[t] Q[t]-n Q[t]-θ m Q[t]-v Q[t],β'[t]==(Piecewise[{{0, β[t]≤0}, {-Vc Q[t], 0<β[t]}}]),S[0]==1,Q[0]==0.001`,β[0]==0.06`},{S,Q,β},{t,0,tmax}];
plotsi4=Plot[Evaluate[{S[t],Q[t]}/.sol4],{t,0,tmax},PlotStyle->{Blue,Red},PlotRange->All,ImageSize->{640,240}]
plotbeta4=Plot[Evaluate[β[t]/.sol4],{t,0,tmax},ImageSize->{640,240},PlotStyle->{Green},PlotRange->All]
plotinf4=Plot[Evaluate[((Q[t])/(S[t]+Q[t]))/.sol4],{t,0,tmax},PlotStyle->{Green},PlotRange->All,ImageSize->{640,240}]
plotdens4=Plot[Evaluate[{S[t]+Q[t]}/.sol4],{t,0,tmax},PlotStyle->{Green},PlotRange->All,ImageSize->{640,240}]

Vc=0;

sol5=NDSolve[{S'[t]==b (S[t]+f Q[t]) (1-c (S[t]+Q[t]))-n S[t]-m S[t]-β[t] S[t] Q[t],Q'[t]==β[t] S[t] Q[t]-n Q[t]-θ m Q[t]-v Q[t],β'[t]==(Piecewise[{{0, β[t]≤0}, {-Vc Q[t], 0<β[t]}}]),S[0]==1,Q[0]==0.001`,β[0]==0.06`},{S,Q,β},{t,0,tmax}];
plotsi5=Plot[Evaluate[{S[t],Q[t]}/.sol5],{t,0,tmax},PlotStyle->{Blue,Red},PlotRange->All,ImageSize->{640,240}]
plotbeta5=Plot[Evaluate[β[t]/.sol5],{t,0,tmax},ImageSize->{640,240},PlotStyle->{Purple},PlotRange->All]
plotinf5=Plot[Evaluate[((Q[t])/(S[t]+Q[t]))/.sol5],{t,0,tmax},PlotStyle->{Purple},PlotRange->All,ImageSize->{640,240}]
plotdens5=Plot[Evaluate[{S[t]+Q[t]}/.sol5],{t,0,tmax},PlotStyle->{Purple},PlotRange->All,ImageSize->{640,240}]

BetaInfplot=Show[plotinf1,plotinf2,plotinf3,plotinf4,plotinf5,AspectRatio->0.5,PlotRange->{{0,100},{0,0.5}}]

Betaplot=Show[plotbeta1,plotbeta2,plotbeta3,plotbeta4,plotbeta5,AspectRatio->0.5,PlotRange->{{0,100},{0,0.07}}]

densityplot=Show[plotdens1,plotdens2,plotdens3,plotdens4,plotdens5,AspectRatio->0.5,PlotRange->{{0,100},{0,20}}]