Skip to content

Commit

Permalink
new changes
Browse files Browse the repository at this point in the history
  • Loading branch information
Qixin He committed Nov 15, 2022
1 parent 984ea3a commit 9f831d1
Show file tree
Hide file tree
Showing 3 changed files with 74 additions and 79 deletions.
Binary file modified .DS_Store
View file
Binary file not shown.
153 changes: 74 additions & 79 deletions ODESImmunity_v13.R
100644 → 100755
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,12 +72,12 @@ SEAI_p<-function(t,y,p){


# 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))
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost)*(varyCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost)*(varyCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost)*(varyCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost)*(varyCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost)*(varyCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost)*(varyCost), 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
Expand Down Expand Up @@ -157,19 +157,19 @@ SEAI_p<-function(t,y,p){

#lose immunity-----
#SU back to NU due to loss of immune memory
hsu = (Bu*(1-rSU/capacity)) + rSU*SU2SD+rSU*SU2AU+att
hsu = (Bu*(1-rSU/capacity)) + att
SU2NU = hsu*exp(-hsu/ImmLost)/(1-exp(-hsu/ImmLost))

#SD back to ND due to loss of immune memory
hsd = (Bd*(1-rSD/capacity)) +1/Tau+att
hsd = (Bd*(1-rSD/capacity)) +att
SD2ND = hsd*exp(-hsd/ImmLost)/(1-exp(-hsd/ImmLost))

#AU back to SU due to loss of immune memory
hau = (Bu*(1-rAU/capacity)) + rAU*AU2AD + att
hau = (Bu*(1-rAU/capacity)) + att
AU2SU = hau*exp(-hau/ImmLost)/(1-exp(-hau/ImmLost))

#AD back to SD due to loss of immune memory
had = (Bd*(1-rAD/capacity)) + 1/Tau + att
had = (Bd*(1-rAD/capacity)) + att
AD2SD = had*exp(-had/ImmLost)/(1-exp(-had/ImmLost))

# Parasite pop ODE------
Expand Down Expand Up @@ -444,12 +444,13 @@ SEAI_resOnly<-function(t,y,p){
ctr_AD <- c((1-p0_AD[1])*p0_AD[2], (1-p0_AD[1])*(1-p0_AD[2])*PW_AD/(PW_AD+PR_AD+1e-6), (1-p0_AD[2])*p0_AD[1], (1-p0_AD[1])*(1-p0_AD[2])*PR_AD/(PW_AD+PR_AD+1e-6))

# 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))
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost)*(varyCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost)*(varyCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost)*(varyCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost)*(varyCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost)*(varyCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost)*(varyCost), 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
Expand Down Expand Up @@ -528,19 +529,19 @@ SEAI_resOnly<-function(t,y,p){

#lose immunity-----
#SU back to NU due to loss of immune memory
hsu = (Bu*(1-rSU/capacity)) + rSU*SU2SD+rSU*SU2AU+att
hsu = (Bu*(1-rSU/capacity)) + att
SU2NU = hsu*exp(-hsu/ImmLost)/(1-exp(-hsu/ImmLost))

#SD back to ND due to loss of immune memory
hsd = (Bd*(1-rSD/capacity)) +1/Tau+att
hsd = (Bd*(1-rSD/capacity)) +att
SD2ND = hsd*exp(-hsd/ImmLost)/(1-exp(-hsd/ImmLost))

#AU back to SU due to loss of immune memory
hau = (Bu*(1-rAU/capacity)) + rAU*AU2AD + att
hau = (Bu*(1-rAU/capacity)) + att
AU2SU = hau*exp(-hau/ImmLost)/(1-exp(-hau/ImmLost))

#AD back to SD due to loss of immune memory
had = (Bd*(1-rAD/capacity)) + 1/Tau + att
had = (Bd*(1-rAD/capacity)) + att
AD2SD = had*exp(-had/ImmLost)/(1-exp(-had/ImmLost))

# Parasite pop ODE------
Expand Down Expand Up @@ -763,12 +764,13 @@ SEAI_wildOnly<-function(t,y,p){
ctr_AD <- c((1-p0_AD[1])*p0_AD[2], (1-p0_AD[1])*(1-p0_AD[2])*PW_AD/(PW_AD+PR_AD+1e-6), (1-p0_AD[2])*p0_AD[1], (1-p0_AD[1])*(1-p0_AD[2])*PR_AD/(PW_AD+PR_AD+1e-6))

# 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))
h_NU = c(ctr_NU[1]+ctr_NU[2]+ctr_NU[4]*(mixCost)*(varyCost), ctr_NU[3]*(1-cloneCost)+ctr_NU[4]*(1-mixCost))
h_ND = c(ctr_ND[1]+ctr_ND[2]+ctr_ND[4]*(mixCost)*(varyCost), ctr_ND[3]*(1-cloneCost)+ctr_ND[4]*(1-mixCost))
h_SU = c(ctr_SU[1]+ctr_SU[2]+ctr_SU[4]*(mixCost)*(varyCost), ctr_SU[3]*(1-cloneCost)+ctr_SU[4]*(1-mixCost))
h_SD = c(ctr_SD[1]+ctr_SD[2]+ctr_SD[4]*(mixCost)*(varyCost), ctr_SD[3]*(1-cloneCost)+ctr_SD[4]*(1-mixCost))
h_AU = c(ctr_AU[1]+ctr_AU[2]+ctr_AU[4]*(mixCost)*(varyCost), ctr_AU[3]*(1-cloneCost)+ctr_AU[4]*(1-mixCost))
h_AD = c(ctr_AD[1]+ctr_AD[2]+ctr_AD[4]*(mixCost)*(varyCost), 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
Expand Down Expand Up @@ -848,19 +850,19 @@ SEAI_wildOnly<-function(t,y,p){

#lose immunity-----
#SU back to NU due to loss of immune memory
hsu = (Bu*(1-rSU/capacity)) + rSU*SU2SD+rSU*SU2AU+att
hsu = (Bu*(1-rSU/capacity)) + att
SU2NU = hsu*exp(-hsu/ImmLost)/(1-exp(-hsu/ImmLost))

#SD back to ND due to loss of immune memory
hsd = (Bd*(1-rSD/capacity)) +1/Tau+att
hsd = (Bd*(1-rSD/capacity)) +att
SD2ND = hsd*exp(-hsd/ImmLost)/(1-exp(-hsd/ImmLost))

#AU back to SU due to loss of immune memory
hau = (Bu*(1-rAU/capacity)) + rAU*AU2AD + att
hau = (Bu*(1-rAU/capacity)) + att
AU2SU = hau*exp(-hau/ImmLost)/(1-exp(-hau/ImmLost))

#AD back to SD due to loss of immune memory
had = (Bd*(1-rAD/capacity)) + 1/Tau + att
had = (Bd*(1-rAD/capacity)) + att
AD2SD = had*exp(-had/ImmLost)/(1-exp(-had/ImmLost))

# Parasite pop ODE------
Expand Down Expand Up @@ -1078,9 +1080,9 @@ out2<-out2%>%mutate(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(A
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")
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(No!=1){
write.table(outData,
Expand All @@ -1101,11 +1103,7 @@ currentPara$d2<-0.2
currentPara$d3<-0.2

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))
N0 = unlist(storeOut[1,])
print(N0)
out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)

Expand All @@ -1118,28 +1116,29 @@ 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))
}

#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(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD),
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")
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,
Expand All @@ -1155,16 +1154,12 @@ write.table(outData,
storeOut<-rbind(storeOut,tail(out2,1)[,2:23])
# Step 3 full treatment, wildOnly-----

currentPara$d1<-0.8
currentPara$d2<-0.8
currentPara$d3<-0.8
currentPara$d1<-0.05
currentPara$d2<-0.05
currentPara$d3<-0.05

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))
N0 = unlist(storeOut[1,])
print(N0)
out2 = ode(y=N0,times=t2,func=SEAI_wildOnly,parms=currentPara)

Expand All @@ -1177,28 +1172,28 @@ 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))
}
#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(m_N = P_acc_N/(NU+ND),m_S =P_acc_S/(SU+SD), m_A = P_acc_A/(AU+AD),
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")
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,
Expand Down Expand Up @@ -1312,9 +1307,9 @@ write.table(outData,

# Step 6 treat fully, wildOnly + resistant-----

currentPara$d1<-0.8
currentPara$d2<-0.8
currentPara$d3<-0.8
currentPara$d1<-0.05
currentPara$d2<-0.05
currentPara$d3<-0.05

t2 = seq(from=0,to=50*365,by=1)
N0 = unlist(storeOut[3,])
Expand Down
Binary file removed batchRuns/.DS_Store
View file
Binary file not shown.

0 comments on commit 9f831d1

Please sign in to comment.