diff --git a/GeneralizedImmunityExp1.R b/GeneralizedImmunityExp1.R index 4b633fc..2f596b1 100644 --- a/GeneralizedImmunityExp1.R +++ b/GeneralizedImmunityExp1.R @@ -6,9 +6,9 @@ paramFile<-args[1] No<-as.integer(args[2]) resultsFolder <-args[3] jobid<-as.integer(args[4]) -#paramFile<-"paramList.csv" -#No<-1 -#resultsFolder<-"./" +paramFile<-"ParamListGen.csv" +No<-3 +resultsFolder<-"./" outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="") SEAI_p<-function(t,y,p){ @@ -124,12 +124,15 @@ SEAI_p<-function(t,y,p){ NU2ND = (Bu*(1-rNU/capacity)*infN*d1) #naive U class getting enough # of infections to move to SU - NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + NU2SU = (rNU*muRU_gen_N)*rho1 + #naive U class die from infections NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) #naive D class getting enough # of infections to move to SD - ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + ND2SD = (rND*muRD_gen_N)*rho1 #naive D class die from infections NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) @@ -138,10 +141,12 @@ SEAI_p<-function(t,y,p){ SU2SD = (Bu*(1-rSU/capacity)*infS*d2) #SU class getting enough # of infections to move to AU - SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + SU2AU = (rSU*muWU_gen_S)*rho2 #SD class getting enough # of infections to move to AD - SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + SD2AD = (rSD*muWD_gen_S)*rho2 #AU class getting infected, show symptoms, and treated AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) @@ -296,23 +301,24 @@ SEAI_p<-function(t,y,p){ # track the mean number of times infected per immunity class------ - dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ - NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received - (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted - (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted - (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received + (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S + (ImmLost+att)*P_acc_N dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S - SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + - SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections - att*P_acc_S- # normal death + SU*SU2AU/rho2 + + SD*SD2AD/rho2 - #total number of new infections + (ImmLost+att)*P_acc_S- # normal death (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic - AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + - AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + + AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections (ImmLost+att)*P_acc_A # normal death + a <- 365/(2*pi) db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a)) #db.dt =0 @@ -425,6 +431,8 @@ SEAI_resOnly<-function(t,y,p){ # per capita biting rate for hosts who are drug treated Bd = B_PR #print(B) + # probability that the infected strain has not been seen by the adaptive immune system + # probability that the infected strain has not been seen by the adaptive immune system #infN = (1-1/nstrains)^fracN infN = 1 #for generalized immunity model, infectivity values are fixed over time @@ -437,12 +445,15 @@ SEAI_resOnly<-function(t,y,p){ NU2ND = (Bu*(1-rNU/capacity)*infN*d1) #naive U class getting enough # of infections to move to SU - NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + NU2SU = (rNU*muRU_gen_N)*rho1 + #naive U class die from infections NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) #naive D class getting enough # of infections to move to SD - ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + ND2SD = (rND*muRD_gen_N)*rho1 #naive D class die from infections NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) @@ -451,14 +462,15 @@ SEAI_resOnly<-function(t,y,p){ SU2SD = (Bu*(1-rSU/capacity)*infS*d2) #SU class getting enough # of infections to move to AU - SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + SU2AU = (rSU*muWU_gen_S)*rho2 #SD class getting enough # of infections to move to AD - SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + SD2AD = (rSD*muWD_gen_S)*rho2 #AU class getting infected, show symptoms, and treated - AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) - + AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) # Parasite pop ODE------ ## sensitive------- @@ -573,24 +585,25 @@ SEAI_resOnly<-function(t,y,p){ # track the mean number of times infected per immunity class------ - dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ - NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received - (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted - (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted - (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received + (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S + (ImmLost+att)*P_acc_N dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S - SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + - SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections - att*P_acc_S- # normal death + SU*SU2AU/rho2 + + SD*SD2AD/rho2 - #total number of new infections + (ImmLost+att)*P_acc_S- # normal death (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic - AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + - AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + + AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections (ImmLost+att)*P_acc_A # normal death + a <- 365/(2*pi) db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a)) #db.dt =0 @@ -702,6 +715,8 @@ SEAI_wildOnly<-function(t,y,p){ # probability that the infected strain has not been seen by the adaptive immune system # probability that the infected strain has not been seen by the adaptive immune system + # probability that the infected strain has not been seen by the adaptive immune system + # probability that the infected strain has not been seen by the adaptive immune system #infN = (1-1/nstrains)^fracN infN = 1 #for generalized immunity model, infectivity values are fixed over time @@ -714,12 +729,15 @@ SEAI_wildOnly<-function(t,y,p){ NU2ND = (Bu*(1-rNU/capacity)*infN*d1) #naive U class getting enough # of infections to move to SU - NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + NU2SU = (rNU*muRU_gen_N)*rho1 + #naive U class die from infections NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) #naive D class getting enough # of infections to move to SD - ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + #ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + ND2SD = (rND*muRD_gen_N)*rho1 #naive D class die from infections NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) @@ -728,14 +746,15 @@ SEAI_wildOnly<-function(t,y,p){ SU2SD = (Bu*(1-rSU/capacity)*infS*d2) #SU class getting enough # of infections to move to AU - SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + #SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + SU2AU = (rSU*muWU_gen_S)*rho2 #SD class getting enough # of infections to move to AD - SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + #SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + SD2AD = (rSD*muWD_gen_S)*rho2 #AU class getting infected, show symptoms, and treated - AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) - + AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) # Parasite pop ODE------ ## sensitive------- @@ -850,21 +869,21 @@ SEAI_wildOnly<-function(t,y,p){ # track the mean number of times infected per immunity class------ - dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ - NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received - (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted - (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted - (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + dP_acc_N.dt = NU*NU2SU/rho1 + ND*ND2SD/rho1- # total number of infections received + (NU*(NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN- #when naive class leave to S + (ImmLost+att)*P_acc_N dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S - SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + - SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections - att*P_acc_S- # normal death + SU*SU2AU/rho2 + + SD*SD2AD/rho2 - #total number of new infections + (ImmLost+att)*P_acc_S- # normal death (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic - AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + - AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + AU*(Bu*(1-rAU/capacity)*(1-infA)+rAU*muWU) + + AD*(Bd*(1-rAD/capacity)*(1-infA)+rAD*muRD)- #total number of new infections (ImmLost+att)*P_acc_A # normal death @@ -897,7 +916,6 @@ currentPara<-allPara[allPara$No==No,] - # Step 1 get no drug treating intensity, wildOnly----- currentPara$d1<-0 @@ -942,8 +960,8 @@ out2<-out2%>%mutate(PW = PW_NU+PW_SU+PW_AU+PW_ND+PW_SD+PW_AD, ##you can check the output through this graph #out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1) -#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ -# facet_wrap(.~variable,scales="free_y") +ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ + facet_wrap(.~variable,scales="free_y") outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1))) if(file.exists(outFile)){ diff --git a/batchRuns/GeneralizedImmunityExp.R b/batchRuns/GeneralizedImmunityExp.R new file mode 100644 index 0000000..d7bdea5 --- /dev/null +++ b/batchRuns/GeneralizedImmunityExp.R @@ -0,0 +1,1242 @@ +require(deSolve) +library(tidyverse) +#args <- commandArgs(trailingOnly = TRUE) +#print(args) +#paramFile<-args[1] +#No<-as.integer(args[2]) +#resultsFolder <-args[3] +#jobid<-as.integer(args[4]) +paramFile<-"ParamListGen.csv" +No<-3 +#resultsFolder<-"./" +outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="") + +SEAI_p<-function(t,y,p){ + + with(as.list(c(y,p)),{ + # total host population size + #print(t) + N = NU+ND+SU+SD+AU+AD + # total parasite population size + #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD + #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD + #P = PW+PR + + # parasite/host ratio + rNU <- (PW_NU+PR_NU)/(NU+1e-6) + rND <- (PW_ND+PR_ND)/(ND+1e-6) + rSU <- (PW_SU+PR_SU)/(SU+1e-6) + rSD <- (PW_SD+PR_SD)/(SD+1e-6) + rAU <- (PW_AU+PR_AU)/(AU+1e-6) + rAD <- (PW_AD+PW_AD)/(AD+1e-6) + + # prevalence + I_NU = 1-exp(-rNU) + I_ND = 1-exp(-rND) + I_SU = 1-exp(-rSU) + I_SD = 1-exp(-rSD) + I_AU = 1-exp(-rAU) + I_AD = 1-exp(-rAD) + + # proportion of hosts that have no parasites from W or + p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) + p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) + p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) + p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) + p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) + p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) + + # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite + ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2])) + ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2])) + ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2])) + ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2])) + ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2])) + ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2])) + + # converted contribution of W or R with mixed infection cos---- + h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost)) + h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost)) + h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost)) + h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost)) + h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost)) + h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost)) + + #transmissibility is reduced to the number of infected hosts, + #receiving hosts needs to have open "resource room" for accepting new strains + B_PW = b*( + h_NU[1] + + h_ND[1] + + h_SU[1] + + h_SD[1] + + h_AU[1] + + h_AD[1] + ) + + + B_PR = b*( + h_NU[2]+ + h_ND[2]+ + h_SU[2]+ + h_SD[2]+ + h_AU[2]+ + h_AD[2] + ) + + + + # per capita biting rate for hosts who are not treated + Bu = (B_PW+B_PR) + # per capita biting rate for hosts who are drug treated + Bd = B_PR + #print(B) + + # pop size +1 is to prevent N = 0 for undefined scenario + fracN = P_acc_N/(NU+ND+1e-6) + fracS = P_acc_S/(SU+SD+1e-6) + fracA = P_acc_A/(AU+AD+1e-6) + + muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1) + + + # probability that the infected strain has not been seen by the adaptive immune system + #infN = (1-1/nstrains)^fracN + infN = 1 #for generalized immunity model, infectivity values are fixed over time + #infS = (1-1/nstrains)^fracS + infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5 + #infA = (1-1/nstrains)^fracA + infA = 1 + #print(infN) + #naive class getting infected, and treated + NU2ND = (Bu*(1-rNU/capacity)*infN*d1) + + #naive U class getting enough # of infections to move to SU + NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #naive U class die from infections + NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) + + #naive D class getting enough # of infections to move to SD + ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + + #naive D class die from infections + NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) + + #SU class getting infected, show symptoms, and treated + SU2SD = (Bu*(1-rSU/capacity)*infS*d2) + + #SU class getting enough # of infections to move to AU + SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + + #SD class getting enough # of infections to move to AD + SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + + #AU class getting infected, show symptoms, and treated + AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) + + + + # Parasite pop ODE------ + ## sensitive------- + dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections + (PW_ND/Tau)- # drug loose effectiveness + (PW_NU*NU2ND) - #naive class getting infected, and treated + (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S + (PW_NU*NUInfDie)- #naive class die from infections + (PW_NU*att)- #normal death rate + (PW_NU*muWU_gen_N)#clearance of infection + + + #print(dNU.dt) + + dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections + (PW_NU*NU2ND) - #naive class getting infected, and treated + (PW_ND/Tau) - # drug loose effectiveness + (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S + (PW_ND*NDInfDie) -#naive class die from infections + (PW_ND*att)- #normal death rate + (PW_ND*muWD_gen_N)#clearance from drug + + #print(dND.dt) + + dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections + (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S + (PW_SD/Tau)- # drug loose effectiveness + (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU + (PW_SU*att)- #normal death rate + (PW_SU*muWU_gen_S)#clearance of infection + + dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections + (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S + (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PW_SD/Tau) - # drug loose effectiveness + (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD + (PW_SD*att)- #normal death rate + (PW_SD*muWD_gen_S)#clearance from drug + + dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections + (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU + (PW_AD/Tau)- # drug loose effectiveness + (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PW_AU*att)- #normal death rate + (PW_AU*muWU_gen_A)#clearance from infection + + dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) + #new infections + (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD + (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PW_AD/Tau)- # drug loose effectiveness + (PW_AD*att)- #normal death rate + (PW_AD*muWD_gen_A)#clearance from drug + + ## resistant------- + dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections + (PR_ND/Tau)- # drug loose effectiveness + (PR_NU*NU2ND) - #naive class getting infected, and treated + (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S + (PR_NU*NUInfDie)- #naive class die from infections + (PR_NU*att)- #normal death rate + (PR_NU*muRU_gen_N)#clearance of infection + + #print(dNU.dt) + + dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections + (PR_NU*NU2ND) - #naive class getting infected, and treated + (PR_ND/Tau) - # drug loose effectiveness + (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S + (PR_ND*NDInfDie) -#naive class die from infections + (PR_ND*att)- #normal death rate + (PR_ND*muRD_gen_N)#clearance from drug + + #print(dND.dt) + + dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections + (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S + (PR_SD/Tau)- # drug loose effectiveness + (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU + (PR_SU*att)- #normal death rate + (PR_SU*muRU_gen_S)#clearance of infection + + dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections + (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S + (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PR_SD/Tau) - # drug loose effectiveness + (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD + (PR_SD*att)- #normal death rate + (PR_SD*muRD_gen_S)#clearance from drug + + dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections + (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU + (PR_AD/Tau)- # drug loose effectiveness + (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PR_AU*att)- #normal death rate + (PR_AU*muRU_gen_A)#clearance from infection + + dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections + (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD + (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PR_AD/Tau)- # drug loose effectiveness + (PR_AD*att)- #normal death rate + (PR_AD*muRD_gen_A)#clearance from drug + + # Generalized immunity ODE------ + dNU.dt = birth + + (ND/Tau)- # drug loose effectiveness + NU*NU2ND - #naive class getting infected, and treated + NU*NU2SU - #naive class getting enough # of infections to move to S + NU*NUInfDie- #naive class die from infections + (NU*att)#normal death rate + + + #print(dNU.dt) + + dND.dt = NU*NU2ND - #naive class getting infected, and treated + (ND/Tau) - # drug loose effectiveness + ND*ND2SD -#naive class getting enough # of infections to move to S + ND*NDInfDie -#naive class die from infections + (ND*att) #normal death rate + #print(dND.dt) + + dSU.dt = NU*NU2SU +#naive class getting enough # of infections to move to S + (SD/Tau)- # drug loose effectiveness + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + SU*SU2AU -#symptomatic class getting enough # of infections to move to AU + (SU*att) #normal death rate + + dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + (SD/Tau) - # drug loose effectiveness + SD*SD2AD- #symptomatic class getting enough # of infections to move to AD + (SD*att) #normal death rate + + dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS + (AD/Tau)- # drug loose effectiveness + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AU*att) #normal death rate + + dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AD/Tau)- # drug loose effectiveness + (AD*att) #normal death rate + + #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt) + + + # track the mean number of times infected per immunity class------ + dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ + NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received + (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + + dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S + SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + + SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections + att*P_acc_S- # normal death + (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic + + dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic + AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + + AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + (ImmLost+att)*P_acc_A # normal death + + a <- 365/(2*pi) + db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a)) + #db.dt =0 + #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt) + #out1<- list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt, + # dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt, + # dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt, + # dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt)) + #print(out1) + + return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt, + dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt, + dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt, + dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))) + }) +} + +SEAI_resOnly<-function(t,y,p){ + + with(as.list(c(y,p)),{ + # total host population size + N = NU+ND+SU+SD+AU+AD + # total parasite population size + #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD + #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD + #P = PW+PR + + # parasite/host ratio + rNU <- (PW_NU+PR_NU)/(NU+1e-6) + rND <- (PW_ND+PR_ND)/(ND+1e-6) + rSU <- (PW_SU+PR_SU)/(SU+1e-6) + rSD <- (PW_SD+PR_SD)/(SD+1e-6) + rAU <- (PW_AU+PR_AU)/(AU+1e-6) + rAD <- (PW_AD+PW_AD)/(AD+1e-6) + + # prevalence + I_NU = 1-exp(-rNU) + I_ND = 1-exp(-rND) + I_SU = 1-exp(-rSU) + I_SD = 1-exp(-rSD) + I_AU = 1-exp(-rAU) + I_AD = 1-exp(-rAD) + + # proportion of hosts that have no parasites from W or + p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) + p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) + p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) + p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) + p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) + p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) + + # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite + ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2])) + ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2])) + ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2])) + ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2])) + ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2])) + ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2])) + + # converted contribution of W or R with mixed infection cos---- + h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost)) + h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost)) + h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost)) + h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost)) + h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost)) + h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost)) + + #transmissibility is reduced to the number of infected hosts, + #receiving hosts needs to have open "resource room" for accepting new strains + B_PW = b*( + h_NU[1] + + h_ND[1] + + h_SU[1] + + h_SD[1] + + h_AU[1] + + h_AD[1] + ) + + + B_PR = b*( + h_NU[2]+ + h_ND[2]+ + h_SU[2]+ + h_SD[2]+ + h_AU[2]+ + h_AD[2] + ) + + fracN = P_acc_N/(NU+ND+1e-6) + fracS = P_acc_S/(SU+SD+1e-6) + fracA = P_acc_A/(AU+AD+1e-6) + muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1) + + # per capita biting rate for hosts who are not treated + Bu = (B_PW+B_PR) + # per capita biting rate for hosts who are drug treated + Bd = B_PR + #print(B) + # probability that the infected strain has not been seen by the adaptive immune system + #infN = (1-1/nstrains)^fracN + infN = 1 #for generalized immunity model, infectivity values are fixed over time + #infS = (1-1/nstrains)^fracS + infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5 + #infA = (1-1/nstrains)^fracA + infA = 1 + #print(infN) + #naive class getting infected, and treated + NU2ND = (Bu*(1-rNU/capacity)*infN*d1) + + #naive U class getting enough # of infections to move to SU + NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #naive U class die from infections + NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) + + #naive D class getting enough # of infections to move to SD + ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + + #naive D class die from infections + NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) + + #SU class getting infected, show symptoms, and treated + SU2SD = (Bu*(1-rSU/capacity)*infS*d2) + + #SU class getting enough # of infections to move to AU + SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + + #SD class getting enough # of infections to move to AD + SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + + #AU class getting infected, show symptoms, and treated + AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) + + + # Parasite pop ODE------ + ## sensitive------- + dPW_NU.dt = 0 + + #print(dNU.dt) + + dPW_ND.dt = 0 + + #print(dND.dt) + + dPW_SU.dt = 0 + + dPW_SD.dt = 0 + + dPW_AU.dt = 0 + + dPW_AD.dt = 0 + + #resistant-- + dPR_NU.dt = B_PR*NU*infN*(1-rNU/capacity) + #new infections + (PR_ND/Tau)- # drug loose effectiveness + (PR_NU*NU2ND) - #naive class getting infected, and treated + (PR_NU*NU2SU) - #naive class getting enough # of infections to move to S + (PR_NU*NUInfDie)- #naive class die from infections + (PR_NU*att)- #normal death rate + (PR_NU*muRU_gen_N)#clearance of infection + + #print(dNU.dt) + + dPR_ND.dt = B_PR*ND*infN*(1-rND/capacity) + #new infections + (PR_NU*NU2ND) - #naive class getting infected, and treated + (PR_ND/Tau) - # drug loose effectiveness + (PR_ND*ND2SD) -#naive class getting enough # of infections to move to S + (PR_ND*NDInfDie) -#naive class die from infections + (PR_ND*att)- #normal death rate + (PR_ND*muRD_gen_N)#clearance from drug + + #print(dND.dt) + + dPR_SU.dt = B_PR*SU*infS*(1-rSU/capacity) + #new infections + (PR_NU*NU2SU) +#naive class getting enough # of infections to move to S + (PR_SD/Tau)- # drug loose effectiveness + (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PR_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU + (PR_SU*att)- #normal death rate + (PR_SU*muRU_gen_S)#clearance of infection + + dPR_SD.dt = B_PR*SD*infS*(1-rSD/capacity) + #new infections + (PR_ND*ND2SD)+ #naive class getting enough # of infections to move to S + (PR_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PR_SD/Tau) - # drug loose effectiveness + (PR_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD + (PR_SD*att)- #normal death rate + (PR_SD*muRD_gen_S)#clearance from drug + + dPR_AU.dt = B_PR*AU*infA*(1-rAU/capacity) + #new infections + (PR_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU + (PR_AD/Tau)- # drug loose effectiveness + (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PR_AU*att)- #normal death rate + (PR_AU*muRU_gen_A)#clearance from infection + + dPR_AD.dt = B_PR*AD*infA*(1-rAD/capacity) + #new infections + (PR_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD + (PR_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PR_AD/Tau)- # drug loose effectiveness + (PR_AD*att)- #normal death rate + (PR_AD*muRD_gen_A)#clearance from drug + + # Generalized immunity ODE------ + dNU.dt = birth + + (ND/Tau)- # drug loose effectiveness + NU*NU2ND - #naive class getting infected, and treated + NU*NU2SU - #naive class getting enough # of infections to move to S + NU*NUInfDie- #naive class die from infections + (NU*att)#normal death rate + + + #print(dNU.dt) + + dND.dt = NU*NU2ND - #naive class getting infected, and treated + (ND/Tau) - # drug loose effectiveness + ND*ND2SD -#naive class getting enough # of infections to move to S + ND*NDInfDie -#naive class die from infections + (ND*att) #normal death rate + #print(dND.dt) + + dSU.dt = NU*NU2SU +#naive class getting enough # of infections to move to S + (SD/Tau)- # drug loose effectiveness + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + SU*SU2AU -#symptomatic class getting enough # of infections to move to AU + (SU*att) #normal death rate + + dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + (SD/Tau) - # drug loose effectiveness + SD*SD2AD- #symptomatic class getting enough # of infections to move to AD + (SD*att) #normal death rate + + dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS + (AD/Tau)- # drug loose effectiveness + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AU*att) #normal death rate + + dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AD/Tau)- # drug loose effectiveness + (AD*att) #normal death rate + + #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt) + + + # track the mean number of times infected per immunity class------ + dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ + NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received + (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + + dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S + SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + + SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections + att*P_acc_S- # normal death + (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic + + dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic + AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + + AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + (ImmLost+att)*P_acc_A # normal death + + + a <- 365/(2*pi) + db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a)) + #db.dt =0 + #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt) + + return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt, + dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt, + dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt, + dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))) + }) +} + +SEAI_wildOnly<-function(t,y,p){ + + with(as.list(c(y,p)),{ + # total host population size + N = NU+ND+SU+SD+AU+AD + # total parasite population size + #PW = PW_NU+PW_ND+PW_SU+PW_SD+PW_AU+PW_AD + #PR = PR_NU+PR_ND+PR_SU+PR_SD+PR_AU+PR_AD + #P = PW+PR + # parasite/host ratio + rNU <- (PW_NU+PR_NU)/(NU+1e-6) + rND <- (PW_ND+PR_ND)/(ND+1e-6) + rSU <- (PW_SU+PR_SU)/(SU+1e-6) + rSD <- (PW_SD+PR_SD)/(SD+1e-6) + rAU <- (PW_AU+PR_AU)/(AU+1e-6) + rAD <- (PW_AD+PW_AD)/(AD+1e-6) + + # prevalence + I_NU = 1-exp(-rNU) + I_ND = 1-exp(-rND) + I_SU = 1-exp(-rSU) + I_SD = 1-exp(-rSD) + I_AU = 1-exp(-rAU) + I_AD = 1-exp(-rAD) + + # proportion of hosts that have no parasites from W or + p0_NU <- c(exp(-rNU*PW_NU/(PW_NU+PR_NU+1e-6)), exp(-rNU*PR_NU/(PW_NU+PR_NU+1e-6))) + p0_ND <- c(exp(-rND*PW_ND/(PW_ND+PR_ND+1e-6)), exp(-rND*PR_ND/(PW_ND+PR_ND+1e-6))) + p0_SU <- c(exp(-rSU*PW_SU/(PW_SU+PR_SU+1e-6)), exp(-rSU*PR_SU/(PW_SU+PR_SU+1e-6))) + p0_SD <- c(exp(-rSD*PW_SD/(PW_SD+PR_SD+1e-6)), exp(-rSD*PR_SD/(PW_SD+PR_SD+1e-6))) + p0_AU <- c(exp(-rAU*PW_AU/(PW_AU+PR_AU+1e-6)), exp(-rAU*PR_AU/(PW_AU+PR_AU+1e-6))) + p0_AD <- c(exp(-rAD*PW_AD/(PW_AD+PR_AD+1e-6)), exp(-rAD*PR_AD/(PW_AD+PR_AD+1e-6))) + + # contributions of wild strain clones, wild strain in mixed, resistant clone and resistant mix in a bite + ctr_NU <- c(exp(-rNU)*(1/p0_NU[1]-1),(1-p0_NU[1])*(1-exp(-rNU)/p0_NU[1]),exp(-rNU)*(1/p0_NU[2]-1),(1-p0_NU[2])*(1-exp(-rNU)/p0_NU[2])) + ctr_ND <- c(exp(-rND)*(1/p0_ND[1]-1),(1-p0_ND[1])*(1-exp(-rND)/p0_ND[1]),exp(-rND)*(1/p0_ND[2]-1),(1-p0_ND[2])*(1-exp(-rND)/p0_ND[2])) + ctr_SU <- c(exp(-rSU)*(1/p0_SU[1]-1),(1-p0_SU[1])*(1-exp(-rSU)/p0_SU[1]),exp(-rSU)*(1/p0_SU[2]-1),(1-p0_SU[2])*(1-exp(-rSU)/p0_SU[2])) + ctr_SD <- c(exp(-rSD)*(1/p0_SD[1]-1),(1-p0_SD[1])*(1-exp(-rSD)/p0_SD[1]),exp(-rSD)*(1/p0_SD[2]-1),(1-p0_SD[2])*(1-exp(-rSD)/p0_SD[2])) + ctr_AU <- c(exp(-rAU)*(1/p0_AU[1]-1),(1-p0_AU[1])*(1-exp(-rAU)/p0_AU[1]),exp(-rAU)*(1/p0_AU[2]-1),(1-p0_AU[2])*(1-exp(-rAU)/p0_AU[2])) + ctr_AD <- c(exp(-rAD)*(1/p0_AD[1]-1),(1-p0_AD[1])*(1-exp(-rAD)/p0_AD[1]),exp(-rAD)*(1/p0_AD[2]-1),(1-p0_AD[2])*(1-exp(-rAD)/p0_AD[2])) + + # converted contribution of W or R with mixed infection cos---- + h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost)) + h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost)) + h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost)) + h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost)) + h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost)) + h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost), ctr_AD[3]*(1-cloneCost)+ctr_AD[4]*(1-mixCost)) + + #transmissibility is reduced to the number of infected hosts, + #receiving hosts needs to have open "resource room" for accepting new strains + B_PW = b*( + h_NU[1] + + h_ND[1] + + h_SU[1] + + h_SD[1] + + h_AU[1] + + h_AD[1] + ) + + + B_PR = b*( + h_NU[2]+ + h_ND[2]+ + h_SU[2]+ + h_SD[2]+ + h_AU[2]+ + h_AD[2] + ) + + + fracN = P_acc_N/(NU+ND+1e-6) + fracS = P_acc_S/(SU+SD+1e-6) + fracA = P_acc_A/(AU+AD+1e-6) + + muWU_gen_N = 1/((1/muWU)*exp(-red*(fracN-((2/3)*muWU)))+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*(fracN-((2/3)*muWD)))+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*(fracN-((2/3)*muRU)))+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*(fracN-((2/3)*muRD)))+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*(fracS-((2/3)*muWU)))+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*(fracS-((2/3)*muWD)))+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*(fracS-((2/3)*muRU)))+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*(fracS-((2/3)*muRD)))+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*(fracA-((2/3)*muWU)))+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*(fracA-((2/3)*muWD)))+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*(fracA-((2/3)*muRU)))+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*(fracA-((2/3)*muRD)))+1) + + # per capita biting rate for hosts who are not treated + Bu = (B_PW+B_PR) + # per capita biting rate for hosts who are drug treated + Bd = B_PR + #print(B) + + + # probability that the infected strain has not been seen by the adaptive immune system + # probability that the infected strain has not been seen by the adaptive immune system + # probability that the infected strain has not been seen by the adaptive immune system + #infN = (1-1/nstrains)^fracN + infN = 1 #for generalized immunity model, infectivity values are fixed over time + #infS = (1-1/nstrains)^fracS + infS = 1#red is the reduction in infectivity due to generalized immunity and is between 0 and 0.5 + #infA = (1-1/nstrains)^fracA + infA = 1 + #print(infN) + #naive class getting infected, and treated + NU2ND = (Bu*(1-rNU/capacity)*infN*d1) + + #naive U class getting enough # of infections to move to SU + NU2SU = (Bu*(1-rNU/capacity)*rho1*1/((infN/(muWU_gen_N))+((1-infN)))) + #naive U class die from infections + NUInfDie = (Bu*(1-rNU/capacity)*infN*muWU_gen_N*muNU) + + #naive D class getting enough # of infections to move to SD + ND2SD = (Bd*(1-rND/capacity)*rho1*1/((infN/(muRU_gen_N))+((1-infN)))) + + #naive D class die from infections + NDInfDie = (Bd*(1-rND/capacity)*infN*muRD_gen_N*muND) + + #SU class getting infected, show symptoms, and treated + SU2SD = (Bu*(1-rSU/capacity)*infS*d2) + + #SU class getting enough # of infections to move to AU + SU2AU = (Bu*(1-rSU/capacity)*rho2*1/((infS/(muWU_gen_S))+((1-infS)))) + + #SD class getting enough # of infections to move to AD + SD2AD = (Bd*(1-rSD/capacity)*rho2*1/((infS/(muRD_gen_S))+((1-infS)))) + + #AU class getting infected, show symptoms, and treated + AU2AD = (Bu*(1-rAU/capacity)*infA*omega*d3) + + + # Parasite pop ODE------ + ## sensitive------- + dPW_NU.dt = B_PW*NU*infN*(1-rNU/capacity) + #new infections + (PW_ND/Tau)- # drug loose effectiveness + (PW_NU*NU2ND) - #naive class getting infected, and treated + (PW_NU*NU2SU) - #naive class getting enough # of infections to move to S + (PW_NU*NUInfDie)- #naive class die from infections + (PW_NU*att)- #normal death rate + (PW_NU*muWU_gen_N)#clearance of infection + + + #print(dNU.dt) + + dPW_ND.dt = B_PW*ND*infN*(1-rND/capacity) + #new infections + (PW_NU*NU2ND) - #naive class getting infected, and treated + (PW_ND/Tau) - # drug loose effectiveness + (PW_ND*ND2SD) -#naive class getting enough # of infections to move to S + (PW_ND*NDInfDie) -#naive class die from infections + (PW_ND*att)- #normal death rate + (PW_ND*muWD_gen_N)#clearance from drug + + #print(dND.dt) + + dPW_SU.dt = B_PW*SU*infS*(1-rSU/capacity) + #new infections + (PW_NU*NU2SU) +#naive class getting enough # of infections to move to S + (PW_SD/Tau)- # drug loose effectiveness + (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PW_SU*SU2AU) -#symptomatic class getting enough # of infections to move to AU + (PW_SU*att)- #normal death rate + (PW_SU*muWU_gen_S)#clearance of infection + + dPW_SD.dt = B_PW*SD*infS*(1-rSD/capacity) + #new infections + (PW_ND*ND2SD)+ #naive class getting enough # of infections to move to S + (PW_SU*SU2SD)- #symptomatic class getting infected, show symptoms, and treated + (PW_SD/Tau) - # drug loose effectiveness + (PW_SD*SD2AD)- #symptomatic class getting enough # of infections to move to AD + (PW_SD*att)- #normal death rate + (PW_SD*muWD_gen_S)#clearance from drug + + dPW_AU.dt = B_PW*AU*infA*(1-rAU/capacity) + #new infections + (PW_SU*SU2AU)+ #symptomatic class getting enough # of infections to move to AU + (PW_AD/Tau)- # drug loose effectiveness + (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PW_AU*att)- #normal death rate + (PW_AU*muWU_gen_A)#clearance from infection + + dPW_AD.dt = B_PW*AD*infA*(1-rAD/capacity) + #new infections + (PW_SD*SD2AD) + #symptomatic class getting enough # of infections to move to AD + (PW_AU*AU2AD)- #asymptomatic class getting infected, show symptoms, and treated + (PW_AD/Tau)- # drug loose effectiveness + (PW_AD*att)- #normal death rate + (PW_AD*muWD_gen_A)#clearance from drug + + ## resistant------- + dPR_NU.dt = 0 + + #print(dNU.dt) + + dPR_ND.dt = 0 + + #print(dND.dt) + + dPR_SU.dt = 0 + + dPR_SD.dt = 0 + + dPR_AU.dt = 0 + + dPR_AD.dt = 0 + # Generalized immunity ODE------ + dNU.dt = birth + + (ND/Tau)- # drug loose effectiveness + NU*NU2ND - #naive class getting infected, and treated + NU*NU2SU - #naive class getting enough # of infections to move to S + NU*NUInfDie- #naive class die from infections + (NU*att)#normal death rate + + + #print(dNU.dt) + + dND.dt = NU*NU2ND - #naive class getting infected, and treated + (ND/Tau) - # drug loose effectiveness + ND*ND2SD -#naive class getting enough # of infections to move to S + ND*NDInfDie -#naive class die from infections + (ND*att) #normal death rate + #print(dND.dt) + + dSU.dt = NU*NU2SU +#naive class getting enough # of infections to move to S + (SD/Tau)- # drug loose effectiveness + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + SU*SU2AU -#symptomatic class getting enough # of infections to move to AU + (SU*att) #normal death rate + + dSD.dt = ND*ND2SD+ #naive class getting enough # of infections to move to S + SU*SU2SD- #symptomatic class getting infected, show symptoms, and treated + (SD/Tau) - # drug loose effectiveness + SD*SD2AD- #symptomatic class getting enough # of infections to move to AD + (SD*att) #normal death rate + + dAU.dt = SU*SU2AU+ #symptomatic class getting enough # of infections to move to AS + (AD/Tau)- # drug loose effectiveness + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AU*att) #normal death rate + + dAD.dt = SD*SD2AD + #symptomatic class getting enough # of infections to move to AS + AU*AU2AD- #asymptomatic class getting infected, show symptoms, and treated + (AD/Tau)- # drug loose effectiveness + (AD*att) #normal death rate + + #sprintf("%.0f, %.0f,%.0f,%.0f,%.0f,%.0f",dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt) + + + # track the mean number of times infected per immunity class------ + dP_acc_N.dt = NU*Bu*(1-rNU/capacity)*(1-infN) + ND*Bd*(1-rND/capacity)*(1-infN)+ + NU*Bu*(1-rNU/capacity)*infN*muWU_gen_N + ND*Bd*(1-rND/capacity)*infN*muRD_gen_N- # total number of infections received + (NU*(att+NUInfDie))*fracN- #when naive U class die how many infections deducted + (ND*(att+NDInfDie))*fracN- #when naive D class die many infections deducted + (NU*NU2SU+ND*ND2SD)*fracN #when naive class leave to S + + dP_acc_S.dt = (NU*NU2SU+ND*ND2SD)*fracN+ #when naive class leave to S + SU*Bu*(1-rSU/capacity)*(infS*muWU_gen_S+(1-infS)) + + SD*Bd*(1-rSD/capacity)*(infS*muRD_gen_S+(1-infS))- #total number of new infections + att*P_acc_S- # normal death + (SU*SU2AU+SD*SD2AD)*fracS #Symptomatic move to Asymptomatic + + dP_acc_A.dt = (SU*SU2AU+SD*SD2AD)*fracS+ #Symptomatic move to Asymptomatic + AU*Bu*(1-rAU/capacity)*(infA*muWU_gen_A+(1-infA)) + + AD*Bd*(1-rAD/capacity)*(infA*muRD_gen_A+(1-infA))- #total number of new infections + (ImmLost+att)*P_acc_A # normal death + + + a <- 365/(2*pi) + db.dt = bavgnull * seasonality / a *(1- dryWidth *cos(phase+t/a))*sin(phase+t/a-dryWidth*sin(phase+t/a)) + #db.dt =0 + #sprintf("%.0f, %.0f,%.0f,%.0f",dP_acc_N.dt,dP_acc_S.dt,dP_acc_A.dt,dbeta.dt) + + return(list(c(dPW_NU.dt,dPW_ND.dt,dPW_SU.dt,dPW_SD.dt,dPW_AU.dt,dPW_AD.dt, + dPR_NU.dt,dPR_ND.dt,dPR_SU.dt,dPR_SD.dt,dPR_AU.dt,dPR_AD.dt, + dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt, + dP_acc_N.dt, dP_acc_S.dt, dP_acc_A.dt, db.dt))) + }) +} + +allPara<-read.csv(paramFile,sep=",",header=T) + +currentPara<-allPara[allPara$No==No,] +setwd("/home/chaillet/MalariaResistance/batchRuns") +#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=0,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, +# b=allPara[1,]$bavgnull*(1-allPara[1,]$seasonality)) +#for (i in 1:42){ +# out9 = ode(y=N0,times=t2,func=SEAI_p,parms=allPara[i,]) +# out9_freqRes[i] = (out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13])/(out9[,2]+out9[,3]+out9[,4]+out9[,5]+out9[,6]+out9[,7]+out9[,8]+out9[,9]+out9[,10]+out9[,11]+out9[,12]+out9[,13]) +#} + + + + + +# Step 1 get no drug treating intensity, wildOnly----- +require(deSolve) +currentPara$d1<-0 +currentPara$d2<-0 +currentPara$d3<-0 + +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=0,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_wildOnly,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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>1000), id.vars=1) +ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ + facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1))) + +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} + +storeOut<-tail(out2,1)[,2:23] +# Step 2 treat occasionally, wildOnly----- + +currentPara$d1<-0.5 +currentPara$d2<-0.5 +currentPara$d3<-0.5 + +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=0,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_wildOnly,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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1) +ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ + facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1))) +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} + +storeOut<-rbind(storeOut,tail(out2,1)[,2:23]) +# Step 3 full treatment, wildOnly----- + +currentPara$d1<-1 +currentPara$d2<-1 +currentPara$d3<-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=0,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_wildOnly,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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1) +#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ +# facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=1,currentPara,tail(out2,1))) +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} + +storeOut<-rbind(storeOut,tail(out2,1)[,2:23]) + +# Step 4 resistant strains only, no treatment ------ +currentPara$d1<-0 +currentPara$d2<-0 +currentPara$d3<-0 + +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=0,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_resOnly,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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + out2 = ode(y=N0,times=t2,func=SEAI_resOnly,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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1) +#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ +# facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1))) +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} + +# Step 5 treat occasionally, wildOnly + resistant----- + +currentPara$d1<-0.5 +currentPara$d2<-0.5 +currentPara$d3<-0.5 + +t2 = seq(from=0,to=30*365,by=1) +N0 = unlist(storeOut[2,]) +N0["PR_NU"]<-2 +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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + 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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1) +#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ +# facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1))) +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} + +# Step 6 treat fully, wildOnly + resistant----- + +currentPara$d1<-1 +currentPara$d2<-1 +currentPara$d3<-1 + +t2 = seq(from=0,to=30*365,by=1) +N0 = unlist(storeOut[3,]) +N0["PR_NU"]<-2 +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)) + +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[30*365+1,2:23]) + names(N0)<-colnames(out2[,2:23]) + 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)) +} +out2<-as.data.frame(out2) +out2<-out2%>%mutate(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, prev = 1-exp(-P/N)) + +##you can check the output through this graph +#out2.melt<-reshape2::melt(out2%>%filter(time%%10==0,time>10), id.vars=1) +#ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ +# facet_wrap(.~variable,scales="free_y") +outData<-as.data.frame(cbind(jobid=jobid,wildOnly=0,currentPara,tail(out2,1))) +if(file.exists(outFile)){ + write.table(outData, + file = outFile, row.names=F, col.names=F,sep=",", + quote=F,append=T) +}else{ + write.csv(outData, + file = outFile, row.names=F, + quote=F) + +} diff --git a/batchRuns/ParamListGen.csv b/batchRuns/ParamListGen.csv index 82219b7..090299b 100644 --- a/batchRuns/ParamListGen.csv +++ b/batchRuns/ParamListGen.csv @@ -1,43 +1,43 @@ -No,birth,att,muND,muNU,nstrains,rho1,rho2,omega,cloneCost,mixCost,ImmLost,capacity,muWU,muWD,muRD,muRU,Tau,d1,d2,d3,bavgnull,seasonality,dryWidth,phase,modelType,red -1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0 -2,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.01 -3,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.2 -4,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.3 -5,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.4 -6,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.5 -7,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0 -8,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.1 -9,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.2 -10,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.3 -11,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.4 -12,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.5 -13,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0 -14,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.1 -15,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.2 -16,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.3 -17,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.4 -18,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.5 -19,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0 -20,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.1 -21,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.2 -22,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.3 -23,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.4 -24,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.5 -25,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0 -26,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.1 -27,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.2 -28,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.3 -29,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.4 -30,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.5 -31,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0 -32,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.1 -33,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.2 -34,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.3 -35,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.4 -36,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.5 -37,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0 -38,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.1 -39,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.2 -40,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.3 -41,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.4 -42,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.5 +No,birth,att,muND,muNU,nstrains,rho1,rho2,omega,cloneCost,mixCost,ImmLost,capacity,muWU,muWD,muRD,muRU,Tau,d1,d2,d3,bavgnull,seasonality,dryWidth,phase,modelType,red,eN,eS,eA +1,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0,,, +2,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.1,,, +3,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.2,1,0.7,0.5 +4,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.3,,, +5,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.4,,, +6,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,Generalized,0.5,,, +7,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0,,, +8,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.1,,, +9,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.2,,, +10,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.3,,, +11,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.4,,, +12,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.2,0,1,0,Generalized,0.5,,, +13,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0,,, +14,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.1,,, +15,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.2,,, +16,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.3,,, +17,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.4,,, +18,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.3,0,1,0,Generalized,0.5,,, +19,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0,,, +20,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.1,,, +21,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.2,,, +22,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.3,,, +23,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.4,,, +24,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.4,0,1,0,Generalized,0.5,,, +25,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0,,, +26,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.1,,, +27,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.2,,, +28,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.3,,, +29,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.4,,, +30,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.5,0,1,0,Generalized,0.5,,, +31,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0,,, +32,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.1,,, +33,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.2,,, +34,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.3,,, +35,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.4,,, +36,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.6,0,1,0,Generalized,0.5,,, +37,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0,,, +38,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.1,,, +39,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.2,,, +40,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.3,,, +41,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.4,,, +42,1,9.13E-05,0.1,0.1,20,0.5,0.05,0.01,0.1,0.9,0.000273973,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.7,0,1,0,Generalized,0.5,,, diff --git a/batchRuns/analyze_codes.r b/batchRuns/analyze_codes.r index b837b6d..1c72a65 100755 --- a/batchRuns/analyze_codes.r +++ b/batchRuns/analyze_codes.r @@ -1,9 +1,12 @@ setwd("~/MalariaResistance/batchRuns/") +getwd() +setwd("~/Downloads") library(tidyverse) library(cowplot) library(dplyr) #oneLR<-read.csv("v4_param_results.csv") -oneLR<-read.csv("v5_param_results.csv") +#oneLR<-read.csv("v5_param_results.csv") +oneLR<- read.csv("ParamListGen_results(8).csv") straRef<-data.frame(nstrains=ceiling(10^(seq(0.7,2.3,0.15))),rawStrains = 10^(seq(0.7,2.3,0.15))) oneLR<-oneLR%>%left_join(straRef,by="nstrains") oneLR<-oneLR%>%mutate(resistP = PR/P, @@ -23,6 +26,15 @@ oneLR<-oneLR%>%mutate(resistP = PR/P, adjPrev = (I_NU*NU+I_ND*ND+I_SU*SU+I_SD*SD+I_AU*AU+I_AD*AD)/N ) +ggplot(data=oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(log10(adjB),exp(-red*P_acc_N/(NU+ND))))+geom_point()+geom_line(aes(color=as.factor(red))) + +ggplot(data=oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(log10(adjB),P_acc_N))+geom_point()+geom_line(aes(color=as.factor(red))) + +ggplot(data=oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(log10(adjB),NU+ND))+geom_point()+geom_line(aes(color=as.factor(red))) + +ggplot(data=oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(log10(adjB),P_acc_S))+geom_point()+geom_line(aes(color=as.factor(red))) + +ggplot(data=oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(log10(adjB),exp(-red*P_acc_S/(SU+SD))))+geom_point()+geom_line(aes(color=as.factor(red))) # resistant prevalence as a function of biting rate and strains @@ -30,10 +42,16 @@ ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),log10(rawStrains)))+ geom_raster(aes(fill=resistP))+scale_fill_viridis_c()+ facet_wrap(d1~seasonality, labeller="label_both") +ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),red))+ + geom_raster(aes(fill=resistP))+scale_fill_viridis_c()+ + facet_wrap(d1~seasonality, labeller="label_both") + +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1), aes(log10(adjB),resistP))+geom_point()+geom_line(aes(color=as.factor(red)))#+facet_wrap(d1~cloneCost, labeller="label_both") + # host populations ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),log10(rawStrains)))+ geom_raster(aes(fill=N))+scale_fill_viridis_c()+ - facet_wrap(d1~seasonality, labeller="label_both") + facet_wrap(d1~seasonality~, labeller="label_both") ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),N)) + geom_point() + geom_line(aes(color=as.factor(log10(rawStrains)))) @@ -48,11 +66,11 @@ ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),(AU+AD)/N)) + geom_poin #host-specific prevalence -ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),(I_SU*SU+I_SD*SD)/(SU+SD))) + geom_point() + geom_line(aes(color=as.factor(log10(rawStrains)))) +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1), aes(log10(adjB),(I_SU*SU+I_SD*SD)/(SU+SD))) + geom_point() + geom_line(aes(color=as.factor(red))) -ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),(I_NU*NU+I_ND*ND)/(NU+ND))) + geom_point() + geom_line(aes(color=as.factor(log10(rawStrains)))) +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1), aes(log10(adjB),(I_NU*NU+I_ND*ND)/(NU+ND))) + geom_point() + geom_line(aes(color=as.factor(red))) -ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),(I_AU*AU+I_AD*AD)/(AU+AD))) + geom_point() + geom_line(aes(color=as.factor(log10(rawStrains)))) +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1), aes(log10(adjB),(I_AU*AU+I_AD*AD)/(AU+AD))) + geom_point() + geom_line(aes(color=as.factor(red))) ggplot(oneLR%>%filter(wildOnly==0,d1>0), aes(log10(adjB),(PR_AU+PR_AD)/(PR_AU+PR_AD+PW_AU+PW_AD))) + geom_point() + geom_line(aes(color=as.factor(log10(rawStrains)))) @@ -72,6 +90,14 @@ ggplot(oneLR%>%drop_na()%>%filter(wildOnly==0,d1>0), aes(log10(adjB*adjPrev*365) scale_color_viridis_c(option="turbo")+xlab("log10(EIR per year)")+ ylab("percentage of resistance")+facet_grid(d1~seasonality,labeller="label_both")+ theme_cowplot() + +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5), aes(adjB*adjPrev*365,resistP))+ + geom_point(aes(color=red))+ + scale_color_viridis_c(option="turbo")+xlab("log10(EIR per year)")+ + ylab("percentage of resistance")+facet_grid(d1~seasonality,labeller="label_both")+ + theme_cowplot() + + # AU prevalence vs resistance ggplot(oneLR%>%filter(wildOnly==0,d1>0),aes(I_AU,resistP))+ geom_point(aes(color=log10(nstrains/adjB),shape=as.factor(seasonality)))+scale_color_viridis_c(option="plasma")+xlab("Prevalance in Asymptomatic")+ylab("percentage of resistance")+ @@ -111,14 +137,14 @@ ggplot(oneLR%>%filter(wildOnly==0,d1>0),aes(log10(adjB),(AU+AD)/N))+ # prevalence vs resistance -ggplot(oneLR%>%filter(wildOnly==0,d1>0),aes(adjPrev,resistP))+ - geom_point(aes(color=log10(nstrains/adjB),shape=as.factor(seasonality)))+scale_color_viridis_c(option="plasma")+ +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1),aes(adjPrev,resistP))+ + geom_point(aes(color=red,shape=as.factor(seasonality)))+scale_color_viridis_c(option="plasma")+ xlab("Prevalance")+ylab("percentage of resistance")+ facet_grid(d1~.,labeller="label_both")+ theme_cowplot()+ geom_smooth() - +ggplot(oneLR%>%filter(wildOnly==0,d1>0.5,cloneCost==0.1,red!=0),aes(log10(adjB),adjPrev))+geom_point()+geom_line(aes(color=as.factor(red)))