+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ Learn Git and GitHub without any code!
++ Using the Hello World guide, you’ll start a branch, write comments, and open a pull request. +
+ Read the guide +
+
+
+
+
+
+
+
+
+
+ Permalink
+
+
+
+
+
+
+
+
+
+
+
+
+
+ + MalariaResistance/batchRuns/ODESImmunity_hqx.R +
+ + Go to file + + +
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 
+
+
+
+
+
+
+
+
+ Qixin He
+
+
+ updated immunity.R
+
+
+
+
+
+
+
+ Latest commit
+ 0d7d35c
+ Feb 27, 2022
+
+
+
+
+
+ History
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ | + | require(deSolve) | +
| + | library(tidyverse) | +
| + | args <- commandArgs(trailingOnly = TRUE) | +
| + | print(args) | +
| + | paramFile<-args[1] | +
| + | No<-as.integer(args[2]) | +
| + | resultsFolder <-args[3] | +
| + | #paramFile<-"paramList.csv" | +
| + | #No<-2 | +
| + | #resultsFolder<-"./" | +
| + | + | +
| + | SEAI<-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 | +
| + | P = PWU+PWD+PRU+PRD | +
| + | + |
| + | # parasite/host ratio | +
| + | m = P/N | +
| + | # prevalence | +
| + | I = 1-exp(-m) | +
| + | # converted cost of resistance considering competition | +
| + | h = 1-(ResCost/(1+0.1*exp(-competition*(I-1)))) | +
| + | #adjusted parasite population, transmissibility lower in PR | +
| + | PADJ = PWU+PWD+(PRU+PRD)*h | +
| + | + |
| + | # pop size +1 is to prevent N = 0 for undefined scenario | +
| + | fracN = P_acc_N/(NU+ND+1) | +
| + | fracS = P_acc_S/(SU+SD+1) | +
| + | fracA = P_acc_A/(AU+AD+1) | +
| + | + |
| + | # probability that the infected strain has not been seen by the adaptive immune system | +
| + | infN = exp(-fracN/nstrains) | +
| + | infS = exp(-fracS/nstrains) | +
| + | infA = exp(-fracA/nstrains) | +
| + | + | +
| + | #transmissibility is reduced to the number of infected hosts, | +
| + | #receiving hosts needs to have open "resource room" for accepting new strains | +
| + | B = b*(I/m)*(1-m/capacity) | +
| + | + |
| + | # per capita biting rate for hosts who are not treated | +
| + | Bu = B*PADJ/N | +
| + | # per capita biting rate for hosts who are drug treated | +
| + | Bd = B*(PRU+PRD)*h/N | +
| + | #print(B) | +
| + | #eq b | +
| + | infD = ND*infN+SD*infS+AD*infA | +
| + | infU = NU*infN+SU*infS+AU*infA | +
| + | #print(infD) | +
| + | #print(infU) | +
| + | + |
| + | eq1 = B*h*(infU/N) | +
| + | eq2 = B*h*(infD/N) | +
| + | eq3 = B*(infU/N) | +
| + | eq4 = B*(infD/N) | +
| + | + |
| + | # Parasite pop ODE | +
| + | dPWU.dt = PWU*eq3 + PWD*eq3 - PWU*eq4 - muWU*PWU | +
| + | dPWD.dt = PWD*eq4 + PWU*eq4 - PWD*eq3 - muWD*PWD | +
| + | dPRU.dt = PRU*eq1 + PRD*eq1 - PRU*eq2 - muRU*PRU | +
| + | dPRD.dt = PRD*eq2 + PRU*eq2 - PRD*eq1 - muRD*PRD | +
| + | #sprintf("%.3f, %.3f,%.3f,%.3f",dPWU.dt,dPWD.dt,dPRU.dt,dPRD.dt) | +
| + | + | +
| + | # Generalized immunity ODE | +
| + | dNU.dt = birth + | +
| + | (ND/Tau)- # drug loose effectiveness | +
| + | (NU*Bu*infN*d1) - #naive class getting infected, and treated | +
| + | (NU*Bu*rho1) - #naive class getting enough # of infections to move to S | +
| + | (NU*Bu*muNU)- #naive class die from infections | +
| + | (NU*att)#normal death rate | +
| + | + | +
| + | + |
| + | #print(dNU.dt) | +
| + | + |
| + | dND.dt = (NU*Bu*infN*d1) - #naive class getting infected, and treated | +
| + | (ND/Tau) - # drug loose effectiveness | +
| + | (ND*Bd*rho1) -#naive class getting enough # of infections to move to S | +
| + | (ND*Bd*muND) -#naive class die from infections | +
| + | (ND*att) #normal death rate | +
| + | #print(dND.dt) | +
| + | + |
| + | dSU.dt = (NU*Bu*rho1) +#naive class getting enough # of infections to move to S | +
| + | (SD/Tau)- # drug loose effectiveness | +
| + | (SU*Bu*infS*d2)- #symptomatic class getting infected, show symptoms, and treated | +
| + | (SU*Bu*rho2) -#symptomatic class getting enough # of infections to move to AS | +
| + | (SU*att) #normal death rate | +
| + | + |
| + | dSD.dt = (ND*Bd*rho1)+ #naive class getting enough # of infections to move to S | +
| + | (SU*Bu*infS*d2)- #symptomatic class getting infected, show symptoms, and treated | +
| + | (SD/Tau) - # drug loose effectiveness | +
| + | (SD*Bd*rho2)- #symptomatic class getting enough # of infections to move to AS | +
| + | (SD*att) #normal death rate | +
| + | + |
| + | dAU.dt = (SU*Bu*rho2)+ #symptomatic class getting enough # of infections to move to AS | +
| + | (AD/Tau)- # drug loose effectiveness | +
| + | (AU*Bu*infA*omega*d3)- #asymptomatic class getting infected, show symptoms, and treated | +
| + | (AU*att) #normal death rate | +
| + | + |
| + | dAD.dt = (SD*Bd*rho2) + #symptomatic class getting enough # of infections to move to AS | +
| + | (AU*Bu*infA*omega*d3)- #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 + ND*Bd - # total number of infection received | +
| + | (NU*(Bu*muNU+att))*P_acc_N/(NU+ND+1)- #when naive U class die how many infections deducted | +
| + | (ND*(Bd*muND+att))*P_acc_N/(NU+ND+1)- #when naive D class die many infections deducted | +
| + | ((NU*Bu*rho1)+(ND*Bd*rho1))*P_acc_N/(NU+ND+1) #when naive class leave to S | +
| + | + |
| + | dP_acc_S.dt = ((NU*Bu*rho1)+(ND*Bd*rho1))*P_acc_N/(NU+ND+1)+ #when naive class leave to S | +
| + | SU*Bu + SD*Bd- #total number of new infections | +
| + | (SU+SD)*att*P_acc_S/(SU+SD+1)- # normal death | +
| + | ((SD*Bd*rho2)+ (SU*Bu*rho2))*P_acc_S/(SU+SD+1) #Symptomatic move to Asymptomatic | +
| + | + |
| + | dP_acc_A.dt = ((SD*Bd*rho2)+ (SU*Bu*rho2))*P_acc_S/(SU+SD+1)+ #Symptomatic move to Asymptomatic | +
| + | AU*Bu + AD*Bd- #total number of new infections | +
| + | (AU+AD)*att*P_acc_A/(AU+AD+1) # 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(dPWU.dt,dPWD.dt,dPRU.dt,dPRD.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(PWU=100,PWD=0,PRU=2,PRD=0, | +
| + | NU=1000,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)) | +
| + | + | +
| + | out2 = ode(y=N0,times=t2,func=SEAI,parms=currentPara) | +
| + | out2<-as_tibble(out2) | +
| + | sumOut2<-out2%>%mutate(year = floor(time/365),P = PWU+PWD+PRU+PRD, N = NU+ND+SU+SD+AU+AD)%>% | +
| + | group_by(year)%>%summarise(meanP=mean(P),meanN=mean(N)) | +
| + | + | +
| + | while((abs(sumOut2[20,"meanP"]-sumOut2[30,"meanP"])/sumOut2[30,"meanP"]>0.001 )| | +
| + | (abs(sumOut2[20,"meanN"]-sumOut2[30,"meanN"])/sumOut2[30,"meanN"]>0.001)){ | +
| + | print(out2[30*365,2:15]) | +
| + | t2<-t2+30*365 | +
| + | N0 = as.numeric(out2[30*365,2:15]) | +
| + | names(N0)<-colnames(out2[,2:15]) | +
| + | out2 = ode(y=N0,times=t2,func=SEAI,parms=currentPara) | +
| + | out2<-as_tibble(out2) | +
| + | sumOut2<-out2%>%mutate(year = floor(time/365),P = PWU+PWD+PRU+PRD, 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(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD), | +
| + | P = PWU+PWD+PRU+PRD, 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), id.vars=1) | +
| + | #ggplot(out2.melt, aes(time,value ))+geom_line()+ylim(0,NA)+ | +
| + | # facet_wrap(.~variable,scales="free_y") | +
| + | + | +
| + | outFile<-paste(resultsFolder,"/",strsplit(paramFile,split="[.]")[[1]][1],"_results.csv",sep="") | +
| + | if(file.exists(outFile)){ | +
| + | write.table(as.data.frame(cbind(currentPara,out2[30*365,])), | +
| + | file = outFile, row.names=F, col.names=F,sep=",", | +
| + | quote=F,append=T) | +
| + | }else{ | +
| + | write.csv(as.data.frame(cbind(currentPara,out2[30*365,])), | +
| + | file = outFile, row.names=F, | +
| + | quote=F) | +
| + | + |
| + | } | +
| + | + | +
| + | + | +
| + | + | +