I am also looking for the same , may be you can use XPPAUT :

For my case the code looks like this :

g(t)=heav(t-tm)*((t-tm)^m)*exp(-c*(t-tm))

intx=int{g(t-t')*x}

m=3

c=1.0

#m=1*flr(x/10.0)

#c=1+0.01*x

x'=(-k*x+mi*am*(theta^n)/(theta^n+intx^p))

par tm=3.4

par a=5,p=2,tau=0.2,k=2.2, mi=1.0,am=128,theta=1.33,n=5.5

@ delay=10,total=200,xlo=100,xhi=200,ylo=0,yhi=5

@ maxstor=10000000

done