From 574c46d81e2e1ce8f2c50fb8d94dc79422566655 Mon Sep 17 00:00:00 2001 From: heqixin Date: Tue, 20 Sep 2022 16:01:53 -0400 Subject: [PATCH] v11 code --- .DS_Store | Bin 0 -> 8196 bytes GeneralizedImmunityExp1-old.R | 1247 +++++++++++++++++++++++++++++++++ batchRuns/.DS_Store | Bin 0 -> 6148 bytes 3 files changed, 1247 insertions(+) create mode 100644 .DS_Store create mode 100644 GeneralizedImmunityExp1-old.R create mode 100644 batchRuns/.DS_Store diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..42e1ce030ea6086b515dc9106a8e8e1fa1bb7d28 GIT binary patch literal 8196 zcmeHMYitx%6h3Dt^v;yWEtH{M5SFe)#7Z9&q4LP~Nf5B4+tL=0b#`ZGrNm4 zO`7+>H<^3p z-gC~Kd%pS3+?-tiz^Q;O^IfzNFw>VVKE8Od}cr=<*CX-!!@VC0HXiGk9c^hse(G9AfjDWyB0bO(%{j8TPx;nT@4 z8JYv8rHuO+ff#`s5fJI4LJkZFAb)lJ?oZ?#FGt#saubg22Q4j^P+3s5a8Wg@W;N`g z%$Qfs1V&KwyNtpSF&gwdCs(;Qa))fIoaJ@ho*NjpYZXEV+qB5!=&h_Hs0H|ZD(ii0}swFS-OmC%eA`UlI`0$#~LczW-vBt_<74U zZTG}s0?!TXzw6w$n5%Gu%_(9$i1|j~NB4zfA1n7o)xIqxtZqkw&c5YubTVh?M=IlG|>{Ch8LC zlEI9zRC=RC{~BOaF3zO}XoPmy2?t>mOen()JO^jtbvOs-;bZs=zJ`nN4P1gB;Ai*+ zeudxQGW>(psG^27ScgemkBxW_wqP4>!FJq*UD$)Y*oXZ%gu^(3qnN`-(7|zZv5eFB zI6i?-;SJgKpR%>_6mt;dkcIy-qjOEqRI;k}>B0FCwha0rJ7}1W3H=Ppg-AQD-XD|ID$ z7$pOrm!!<0&(<5_=%YLHv%$TiuJe_*AXjIxB)j~Gqw^p1v7VJ8auHY z_u~O#=Rq96BX|^#VHOQy=t<(JVCX44jgvTqPvO&e250aYJd4lc3-}_wM4WvC-@A^7 zdn*i-65pSjhowxxbv*Y3X)6rA)-fL4OavCl#p*@<-#h>B|2OlI@metgF#@+30$A0P z>FJ_P&G*q2wRV!Ohv^bc*iB0rx= O=YM?uhupiJyMF;lcN1j* literal 0 HcmV?d00001 diff --git a/GeneralizedImmunityExp1-old.R b/GeneralizedImmunityExp1-old.R new file mode 100644 index 0000000..786f2ed --- /dev/null +++ b/GeneralizedImmunityExp1-old.R @@ -0,0 +1,1247 @@ +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]) +jobid<-1351351 +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*( + eN*h_NU[1] + + eN*h_ND[1] + + eS*h_SU[1] + + eS*h_SD[1] + + eA*h_AU[1] + + eA*h_AD[1] + ) + + + B_PR = b*( + eN*h_NU[2]+ + eN*h_ND[2]+ + eS*h_SU[2]+ + eS*h_SD[2]+ + eA*h_AU[2]+ + eA*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)+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+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*( + eN*h_NU[1] + + eN*h_ND[1] + + eS*h_SU[1] + + eS*h_SD[1] + + eA*h_AU[1] + + eA*h_AD[1] + ) + + + B_PR = b*( + eN*h_NU[2]+ + eN*h_ND[2]+ + eS*h_SU[2]+ + eS*h_SD[2]+ + eA*h_AU[2]+ + eA*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)+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+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*( + eN*h_NU[1] + + eN*h_ND[1] + + eS*h_SU[1] + + eS*h_SD[1] + + eA*h_AU[1] + + eA*h_AD[1] + ) + + + B_PR = b*( + eN*h_NU[2]+ + eN*h_ND[2]+ + eS*h_SU[2]+ + eS*h_SD[2]+ + eA*h_AU[2]+ + eA*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)+1) + muWD_gen_N = 1/((1/muWD)*exp(-red*fracN)+1) + muRU_gen_N = 1/((1/muRU)*exp(-red*fracN)+1) + muRD_gen_N = 1/((1/muRD)*exp(-red*fracN)+1) + + muWU_gen_S = 1/((1/muWU)*exp(-red*fracS)+1) + muWD_gen_S = 1/((1/muWD)*exp(-red*fracS)+1) + muRU_gen_S = 1/((1/muRU)*exp(-red*fracS)+1) + muRD_gen_S = 1/((1/muRD)*exp(-red*fracS)+1) + + muWU_gen_A = 1/((1/muWU)*exp(-red*fracA)+1) + muWD_gen_A = 1/((1/muWD)*exp(-red*fracA)+1) + muRU_gen_A = 1/((1/muRU)*exp(-red*fracA)+1) + muRD_gen_A = 1/((1/muRD)*exp(-red*fracA)+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,] + +#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----- + +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/.DS_Store b/batchRuns/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..3e9528e9a780712cdd1e9378a6cfa4f471b73ea4 GIT binary patch literal 6148 zcmeHK%}T>S5Z-NTn^J@x6nb3nTCmk>DPBUYM=wV7pi&zWYA|L?lbS;*~0RhfHx64Q+B`E`Pt2U(EVYI@$t&*GuCE|Sd*?DrW8%5}8EYFf={yR%;RM@J`JF+ARM z#d>AAB9h5>%ug)?GWi$tSo|DHC5(C5l zF+dD#9|Ptruv*(!I#o;z5CcCjfct}hhG-kiG^(uwI=nukzlMkcI=&?kg+be3rV%_K zT&Du+RBoOaT&IIwm^j;DrctLeu2zP5%*y%Wg{#%UE>t+Z6& zU#9etUreD9F+dFbGX{8V;!MU+lsQ|!m4|1ofc6Lt1>-VQKtL~C0$_mqNJlxfU!V?g Yw!usz&VqK84oDXPMF@4oz%MZH1*!Q=761SM literal 0 HcmV?d00001