start

MalariaODEs.R

@@ -0,0 +1,56 @@
+require(deSolve)
+
+SEA<-function(t,y,p){
+  PSU = y[1]
+  PSD = y[2]
+  PS = y[3]
+  PRU = y[4]
+  PRD = y[5]
+  PR = y[6]
+  U = y[7]
+  D = y[8]
+  with(as.list(p),{
+    N = y[7]+y[8]
+    dPSU.dt = (b*(PSU/N)*ws*U) + (b*(PSD/N)*ws*D) - (muSU*PSU)
+    dPSD.dt = (b*(PSD/N)*ws*U) + (b*(PSU/N)*ws*D) - (muSD*PSD)
+    dPS.dt = dPSU.dt + dPSD.dt
+    dPRU.dt = (b*(PRU/N)*wr*U) + (b*(PRD/N)*D) - (muRU*PRU)
+    dPRD.dt = (b*(PRD/N)*wr*U) + (b*(PRU/N)*D) - (muRD*PRD)
+    dPR.dt = dPRU.dt + dPRD.dt
+    dU.dt = -b*((1-exp(-m))*(((PSD+PSU+((PRD+PRU)*(wr))))/(PSD+PSU+PRD+PRU))*(U)*d) + (D/Tau)
+    dD.dt = b*((1-exp(-m))*(((PSD+PSU+((PRD+PRU)*(wr))))/(PSD+PSU+PRD+PRU))*(U)*d) - (D/Tau)
+    return(list(c(dPSU.dt,dPSD.dt,dPS.dt,dPRU.dt,dPRD.dt,dPR.dt,dU.dt,dD.dt)))
+  })
+}
+
+b = 0.5
+ws = 1
+wr = 1-0.6
+muSU = 0.05
+muSD = 0.5
+muRD = 0.05
+muRU = 0.05
+Tau = 10
+d = 0.9
+m = 3
+p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,m=m)
+
+t2 = c(1,5,10,20)
+t2 = seq(from=0,to=400,by=0.01) 
+
+N0 = c(100,100,200,100,100,200,1000,1000)
+out2 = ode(y=N0,times=t2,func=SEA,parms=p2)
+
+plot(out2[,1],(out2[,2]),type="l",ylim=c(0,max(out2[,4])))
+lines(out2[,1],(out2[,3]),col="RED")
+lines(out2[,1],(out2[,4]),col="GREEN")
+legend(x=0,y=80,legend=c("PSU","PSD","PS"),bty="n",lwd = c(2,1),col=c("BLACK","RED","GREEN"))
+
+plot(out2[,1],(out2[,5]),type="l",ylim=c(0,max(out2[,7])))
+lines(out2[,1],(out2[,6]),col="RED")
+lines(out2[,1],(out2[,7]),col="GREEN")
+legend(x=0,y=80,legend=c("PRU","PRD","PR"),bty="n",lwd = c(2,1),col=c("BLACK","RED","GREEN"))
+
+plot(out2[,1],(out2[,8]),type="l",ylim=c(0,2000))
+lines(out2[,1],(out2[,9]),col="RED")
+legend(x=0,y=80,legend=c("U","D"),bty="n",lwd=c(2,1),col=c("BLACK","RED"))