separate P into different classes

batchRuns/ODESImmunity_v2.R

@@ -166,12 +166,12 @@ SEAI_p<-function(t,y,p){
     ic_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1)))
 
     # converted contribution of W or R with mixed infection cost
-    h_NU = c((1-ic_NU[1])*(exp(-rNU)/ic_NU[1]+(1-exp(-rNU)/ic_NU[1])*(1+ResCost/2)),(1-ic_NU[2])*(exp(-rNU)/ic_NU[2]*cloneCost+(1-exp(-rNU)/ic_NU[2])*(1-mixCost/2)))
-    h_ND = c((1-ic_ND[1])*(exp(-rND)/ic_ND[1]+(1-exp(-rND)/ic_ND[1])*(1+ResCost/2)),(1-ic_ND[2])*(exp(-rND)/ic_ND[2]*cloneCost+(1-exp(-rND)/ic_ND[2])*(1-mixCost/2)))
-    h_SU = c((1-ic_SU[1])*(exp(-rSU)/ic_SU[1]+(1-exp(-rSU)/ic_SU[1])*(1+ResCost/2)),(1-ic_SU[2])*(exp(-rSU)/ic_SU[2]*cloneCost+(1-exp(-rSU)/ic_SU[2])*(1-mixCost/2)))
-    h_SD = c((1-ic_SD[1])*(exp(-rSD)/ic_SD[1]+(1-exp(-rSD)/ic_SD[1])*(1+ResCost/2)),(1-ic_SD[2])*(exp(-rSD)/ic_SD[2]*cloneCost+(1-exp(-rSD)/ic_SD[2])*(1-mixCost/2)))
-    h_AU = c((1-ic_AU[1])*(exp(-rAU)/ic_AU[1]+(1-exp(-rAU)/ic_AU[1])*(1+ResCost/2)),(1-ic_AU[2])*(exp(-rAU)/ic_AU[2]*cloneCost+(1-exp(-rAU)/ic_AU[2])*(1-mixCost/2)))
-    h_AD = c((1-ic_AD[1])*(exp(-rAD)/ic_AD[1]+(1-exp(-rAD)/ic_AD[1])*(1+ResCost/2)),(1-ic_AD[2])*(exp(-rAD)/ic_AD[2]*cloneCost+(1-exp(-rAD)/ic_AD[2])*(1-mixCost/2)))
+    h_NU = c((1-ic_NU[1])*(exp(-rNU)/ic_NU[1]+(1-exp(-rNU)/ic_NU[1])*(1+mixCost/2)),(1-ic_NU[2])*(exp(-rNU)/ic_NU[2]*cloneCost+(1-exp(-rNU)/ic_NU[2])*(1-mixCost/2)))
+    h_ND = c((1-ic_ND[1])*(exp(-rND)/ic_ND[1]+(1-exp(-rND)/ic_ND[1])*(1+mixCost/2)),(1-ic_ND[2])*(exp(-rND)/ic_ND[2]*cloneCost+(1-exp(-rND)/ic_ND[2])*(1-mixCost/2)))
+    h_SU = c((1-ic_SU[1])*(exp(-rSU)/ic_SU[1]+(1-exp(-rSU)/ic_SU[1])*(1+mixCost/2)),(1-ic_SU[2])*(exp(-rSU)/ic_SU[2]*cloneCost+(1-exp(-rSU)/ic_SU[2])*(1-mixCost/2)))
+    h_SD = c((1-ic_SD[1])*(exp(-rSD)/ic_SD[1]+(1-exp(-rSD)/ic_SD[1])*(1+mixCost/2)),(1-ic_SD[2])*(exp(-rSD)/ic_SD[2]*cloneCost+(1-exp(-rSD)/ic_SD[2])*(1-mixCost/2)))
+    h_AU = c((1-ic_AU[1])*(exp(-rAU)/ic_AU[1]+(1-exp(-rAU)/ic_AU[1])*(1+mixCost/2)),(1-ic_AU[2])*(exp(-rAU)/ic_AU[2]*cloneCost+(1-exp(-rAU)/ic_AU[2])*(1-mixCost/2)))
+    h_AD = c((1-ic_AD[1])*(exp(-rAD)/ic_AD[1]+(1-exp(-rAD)/ic_AD[1])*(1+mixCost/2)),(1-ic_AD[2])*(exp(-rAD)/ic_AD[2]*cloneCost+(1-exp(-rAD)/ic_AD[2])*(1-mixCost/2)))
 
     #transmissibility is reduced to the number of infected hosts, 
     #receiving hosts needs to have open "resource room" for accepting new strains
@@ -388,24 +388,29 @@ allPara<-read.csv(paramFile,sep=",",header=T)
 
 currentPara<-allPara[allPara$No==No,]
 
-t2 = seq(from=0,to=100*365,by=1) 
+t2 = seq(from=0,to=30*365,by=1) 
 N0 = c(PW_NU=100,PW_ND=0,PW_SU=0,PW_SD=0,PW_AU=0,PW_AD=0,
        PR_NU=2*(1-WildStrainOnly),PR_ND=0,PR_SU=0,PR_SD=0,PR_AU=0,PR_AD=0,
        NU=10000,ND=0,SU=0,SD=0,AU=0,AD=0,
        P_acc_N=0,P_acc_S=0,P_acc_A=0,
        b=currentPara$bavgnull*(1-currentPara$seasonality))
+print(N0)
 out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
 
 out2<-as_tibble(out2)
 sumOut2<-out2%>%mutate(year = floor(time/365),PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD,
                        PR = PR_NU+PR_SU+PR_AU+PR_ND+PR_SD+PR_AD,P=PW+PR, N = NU+ND+SU+SD+AU+AD)%>%
                 group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N))
+print(tail(sumOut2))
 
-while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/sumOut2[30,"meanP"]>0.001 )|
-       (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/sumOut2[30,"meanN"]>0.001)){
-  print(out2[100*365,2:23])
+print(abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001))
+print(abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001))
+
+while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/(sumOut2[30,"meanP"]+0.0001)>0.001 )|
+       (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/(sumOut2[30,"meanN"]+0.0001)>0.001)){
+  print(out2[30*365,2:23])
   t2<-t2+30*365
-  N0 = as.numeric(out2[100*365,2:23])
+  N0 = as.numeric(out2[30*365,2:23])
   names(N0)<-colnames(out2[,2:23])
   out2 = ode(y=N0,times=t2,func=SEAI_p,parms=currentPara)
   out2<-as_tibble(out2)