# Immune HIV Dynamics
# Kirschner, D (1996)
# implementation of model I
# with the Stochastic Simulation Algorithm (Gillespie)
# Univ Lyon 1, Population Dynamics, 2014

p gV=2
p muT=0.02,muTi=0.24
p kV=2.4e-5,kT=7.4e-4
p r=0.01,N=1000,C=100
p b=10,amax=12,a1=0.25



prolif(x,y)=r*x*y/(C+y)
s(V)=5+5/(1+V)

# cumulative reaction rates
p1=s(V)
p2=p1+muT*Tn
p3=p2+prolif(Tn,V)
p4=p3+kV*Tn*V
p5=p4+muTi*Ti
p6=p5+prolif(Ti,V)
p7=p6+kT*Tn*V
p8=p7+gV*V/(b+V)

tr'=tr-log(ran(1))/p8
# choose random reaction
ch=ran(1)*p8
z1=(ch<p1)
z2=(ch<p2)&(ch>=p1)
z3=(ch<p3)&(ch>=p2)
z4=(ch<p4)&(ch>=p3)
z5=(ch<p5)&(ch>=p4)
z6=(ch<p6)&(ch>=p5)
z7=(ch<p7)&(ch>=p6)
z8=(ch<p8)&(ch>=p7)

# update the number of each species
Tn'=max(0,Tn+z1-z2+z3-z4)
Ti'=max(0,Ti+z4-z5-z6)
V'=max(0,V+N*z6-z7+z8)

# initial conditions
init Tn=1000,Ti=0,V=1;

# xpp settings
@ bound=1000000,meth=discrete,total=500000
@ xp=tr,yp=V
@ xlo=0,ylo=0,xi=500,yhi=1000
done