From c35684d9bcbea5eb678550e9cfe59cec1c23b273 Mon Sep 17 00:00:00 2001 From: Qixin He Date: Sun, 27 Feb 2022 13:44:28 -0500 Subject: [PATCH] add batch run files, fix one locus running file --- ODESImmunity2-copy.R | 137 +++++++++++--------------- batchRuns/ODESImmunity_hqx.R | 180 +++++++++++++++++++++++++++++++++++ batchRuns/paramList.csv | 3 + batchRuns/submitODE.sbatch | 19 ++++ 4 files changed, 260 insertions(+), 79 deletions(-) create mode 100644 batchRuns/ODESImmunity_hqx.R create mode 100644 batchRuns/paramList.csv create mode 100755 batchRuns/submitODE.sbatch diff --git a/ODESImmunity2-copy.R b/ODESImmunity2-copy.R index 9f35222..1e731db 100644 --- a/ODESImmunity2-copy.R +++ b/ODESImmunity2-copy.R @@ -1,9 +1,8 @@ require(deSolve) #library(tidyverse) SEAI<-function(t,y,p){ - PSU = y[1] - PSD = y[2] - #PS = y[3] + PWU = y[1] + PWD = y[2] PRU = y[3] PRD = y[4] #PR = y[6] @@ -14,42 +13,22 @@ SEAI<-function(t,y,p){ AU = y[9] AD = y[10] beta = y[11] - #PMOR = y[12] - #HMOR = y[13] - #PSUACC = y[12] - #PSDACC = y[13] - #PRUACC = y[14] - #PRDACC = y[15] - #NUACC = y[18] - #NDACC = y[19] - PSUACCN2 = y[12] - PSDACCN2= y[13] + PWUACCN2 = y[12] + PWDACCN2= y[13] PRUACCN2 = y[14] PRDACCN2 = y[15] - #PSUACC210 = y[16] - #PSDACC210 = y[17] - #PRUACC210 = y[18] - #PRDACC210 = y[19] - #ACC210U = y[20] - #ACC210D = y[21] - #ACC210NU = y[22] - #ACC210ND = y[23] - #PSUACCN2 = y[24] - #PSDACCN2 = y[25] - #PRUACCN2 = y[26] - #PRDACCN2 = y[27] - PSUACC = y[16] - PSDACC = y[17] + PWUACC = y[16] + PWDACC = y[17] PRUACC = y[18] PRDACC = y[19] - PSUACCA = y[20] - PSDACCA = y[21] + PWUACCA = y[20] + PWDACCA = y[21] PRUACCA = y[22] PRDACCA = y[23] with(as.list(p),{ N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10] - m = (PSU+PSD+PRU+PRD)/N - I = 1-exp(-m)#(1-((PSU+PSD+PRU+PRD)/(N*K))) 1-exp(-m) + m = (PWU+PWD+PRU+PRD)/N + I = 1-exp(-m)#(1-((PWU+PWD+PRU+PRD)/(N*K))) 1-exp(-m) #pone = 1 - (m/K) #ptwo = m/(K-1) #pthree = m/(K-2) @@ -80,12 +59,12 @@ SEAI<-function(t,y,p){ #rho3=prob_asymptomatic #rho4=rho3 #h=wr - #totalm = (PSUACC+PSDACC+PRUACC+PRDACC)/(K*(SU+SD)) - frac1 = (PSUACC+PRUACC)/(SU) - frac2 = (PSDACC+PRDACC)/(SD) - frac3 = (PSUACCA+PRUACCA)/AU - frac4 = (PSDACCA+PRDACCA)/AD - #frac2 = PSDACC/(SD) + #totalm = (PWUACC+PWDACC+PRUACC+PRDACC)/(K*(SU+SD)) + frac1 = (PWUACC+PRUACC)/(SU) + frac2 = (PWDACC+PRDACC)/(SD) + frac3 = (PWUACCA+PRUACCA)/AU + frac4 = (PWDACCA+PRDACCA)/AD + #frac2 = PWDACC/(SD) #frac3 = PRUACC/(SU) #frac4 = PRDACC/(SD) #f1 = frac1 #*(1-rho3) @@ -97,15 +76,15 @@ SEAI<-function(t,y,p){ prob24 = min(frac2,nstrains)/nstrains rhobin1 = dbinom(nstrains,size=nstrains,prob=prob23*0.999999) rhobin2 = dbinom(nstrains,size=nstrains,prob=prob24*0.999999) - ratr = ((PRU+PRD)*h)/(PSU+PSD+((PRU+PRD)*h)) - #rats = (PSU+PSD)/(PSU+PSD+PRU+PRD) - #eps = (PSD+PSU+((PRD+PRU)*(wr)))/N + ratr = ((PRU+PRD)*h)/(PWU+PWD+((PRU+PRD)*h)) + #rats = (PWU+PWD)/(PWU+PWD+PRU+PRD) + #eps = (PWD+PWU+((PRD+PRU)*(wr)))/N #Ndeaths_demo = del*att*NU*t+del*att*ND*t #Sdeaths = del*att*SU*t+del*att*SD*t #Adeaths = (del*att*AU*t)+(del*att*AD*t) #totalAhosts = AU+AD#+(del*att*AU*t)+(del*att*AD*t) #rebirth = (PMOR)/(SU+SD) - #probbin = ((PSUACC+PSDACC+PRUACC+PRDACC))/((K)*(SU+SD)) + #probbin = ((PWUACC+PWDACC+PRUACC+PRDACC))/((K)*(SU+SD)) #mbin = pnbinom((1/rho3),size = 1,prob=probbin,lower.tail=TRUE) #meanbin = 1*(1-probbin)/probbin #mean number of treatment failures before a successful immune defense #meanstrainsS = meanbin*((1-rho3)^meanbin) @@ -115,17 +94,17 @@ SEAI<-function(t,y,p){ meanstrainsS4 = (1-rho3)^(frac2) meanstrainsS5 = (1-rho3)^(frac3) meanstrainsS6 = (1-rho3)^(frac4) - eps = (PSD+PSU+((PRD+PRU)*(h)))/(PSD+PSU+PRD+PRU) - #dPSU.dt = - (b*(PSU/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU) - #dPSU.dt = - (b*(PSU/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU) - dPSU.dt = - (b*(PSU/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU)#-(del*att*(1/N)*(NU+SU+AU)*PSU) - #dPSD.dt = - (b*(PSU/N)*ws*(ND+SD*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD) - #dPSD.dt = - (b*(PSU/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD) - dPSD.dt = - (b*(PSU/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD)#-(del*att*(1/N)*(ND+SD+AD)*PSD) - #dPRU.dt = - (b*(PRU/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU) - dPRU.dt = - (b*(PRU/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU) - #dPRD.dt = - (b*(PRU/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD) - dPRD.dt = - (b*(PRU/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD) + eps = (PWD+PWU+((PRD+PRU)*(h)))/(PWD+PWU+PRD+PRU) + #dPWU.dt = - (b*(PWU/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(NU+SU*(1-f1)*meanstrainsS1+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSU)*(PWU)#-(del*att*(1/N)*(NU+SU+AU)*PWU) + #dPWU.dt = - (b*(PWU/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(NU+SU+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSU)*(PWU)#-(del*att*(1/N)*(NU+SU+AU)*PWU) + dPWU.dt = - (b*(PWU/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(NU+SU*(meanstrainsS1+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSU)*(PWU)#-(del*att*(1/N)*(NU+SU+AU)*PWU) + #dPWD.dt = - (b*(PWU/N)*ws*(ND+SD*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ND+SD*(1-f2)*meanstrainsS2+AD*0)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSD)*(PWD)#-(del*att*(1/N)*(ND+SD+AD)*PWD) + #dPWD.dt = - (b*(PWU/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ND+SD+AD*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSD)*(PWD)#-(del*att*(1/N)*(ND+SD+AD)*PWD) + dPWD.dt = - (b*(PWU/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muSD)*(PWD)#-(del*att*(1/N)*(ND+SD+AD)*PWD) + #dPRU.dt = - (b*(PRU/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU+AU*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU) + dPRU.dt = - (b*(PRU/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU) + #dPRD.dt = - (b*(PRU/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD+AD*omega1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD) + dPRD.dt = - (b*(PRU/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND+SD*(meanstrainsS4+omega1)+AD*(omega1+meanstrainsS6))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD) #dNU.dt = (del*(AU+AD))+(ND/Tau)-(b*I*eps*NU*(d+rho1))-(muNU*NU) - (del*att*NU) dNU.dt = del+(ND/Tau)-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU) dND.dt = (b*I*eps*NU*d1)-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND) @@ -136,39 +115,39 @@ SEAI<-function(t,y,p){ #dAD.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(b*I*eps*SD*rho4)+(b*I*eps*omega1*AU*d)-(AD/Tau) - (del*att*AD) dAD.dt = (ratr*b*I*eps*SD*eta)+(b*I*eps*(omega1)*AU*d3)-(AD/Tau) - (del*att*AD) dbeta.dt = ((L*bavgnull*sin(((t/a)+ps)-kap*sin((t/a)+ps))*(kap*cos((t/a)+ps)-1))/a) - #dPMOR.dt = (muSU)*(PSUACC)+(muSD)*PSDACC+muRU*PRUACC+muRD*PRDACC + #dPMOR.dt = (muSU)*(PWUACC)+(muSD)*PWDACC+muRU*PRUACC+muRD*PRDACC #dHMOR.dt = del*att*del*t+(muNU*NUACC*I*b*eps)+(muND*ratr*NDACC*I*b*eps) - #dPSUACC.dt = - (b*(PSU/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PSUACCN2/NU)+(SD/Tau)*(PSDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PSUACC/SU)#*(1-rho3) - #dPSDACC.dt = - (b*(PSU/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PSDACCN2/ND)+(b*I*eps*SU*d2)*(PSUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PSDACC/SD)- (muSD)*(PSDACC)#death rate included because of clearance by drugs rather than immune system - #dPRUACC.dt = - (b*(PRU/N)*h*(SU*meanstrainsS3)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*meanstrainsS3)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3) - #dPRDACC.dt = - (b*(PRU/N)*h*(SD*meanstrainsS4)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*meanstrainsS4)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2)*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PRDACC/SD)#*(1-rho4) - #dPSUACCN.dt = - (b*(PSU/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ND/Tau)*(PSDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PSUACCN/NU) - (muSU)*(PSUACCN)#*(1-rho3) - #dPSDACCN.dt = - (b*(PSU/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PSUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PSDACCN/ND) - (muSD)*(PSDACCN)#*(1-rho4) - #dPRUACCN.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN/NU) - (muRU)*(PRUACCN)#*(1-rho3) - #dPRDACCN.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN/ND)- (muRD)*(PRDACCN)#*(1-rho4) - #dPSUACC210.dt = - (b*(PSU/N)*ws*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(PSDACC210/Tau)-(b*I*eps*(d1))*(PSUACC210)-(muNU*PSUACCN*(ACC210NU/NU)*I*b*eps)-(muSU*PSUACC210)-(PSUACC210*(1/(8*365))) - #dPSDACC210.dt = - (b*(PSU/N)*ws*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(PSDACC210/Tau)+(b*I*eps*(d1))*(PSUACC210)-(muND*PSDACCN*(ACC210ND/ND)*I*b*eps)-(muSD*PSDACC210)-(PSDACC210*(1/(8*365))) - #dPRUACC210.dt = - (b*(PRU/N)*h*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210U)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(PRDACC210/Tau)-(b*I*eps*(d1))*(PRUACC210)-(muNU*PRUACCN*(ACC210NU/NU)*I*b*eps)-(muRU*PRUACC210)-(PRUACC210*(1/(8*365))) - #dPRDACC210.dt = - (b*(PRU/N)*h*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210D)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(PRDACC210/Tau)+(b*I*eps*(d1))*(PRUACC210)-(muND*PSDACCN*(ACC210ND/ND)*I*b*eps)-(muRD*PRDACC210)-(PRDACC210*(1/(8*365))) + #dPWUACC.dt = - (b*(PWU/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(SU*meanstrainsS1)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PWUACCN2/NU)+(SD/Tau)*(PWDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PWUACC/SU)#*(1-rho3) + #dPWDACC.dt = - (b*(PWU/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(SD*meanstrainsS2)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PWDACCN2/ND)+(b*I*eps*SU*d2)*(PWUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PWDACC/SD)- (muSD)*(PWDACC)#death rate included because of clearance by drugs rather than immune system + #dPRUACC.dt = - (b*(PRU/N)*h*(SU*meanstrainsS3)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*meanstrainsS3)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*(d2+(rho3))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3) + #dPRDACC.dt = - (b*(PRU/N)*h*(SD*meanstrainsS4)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*meanstrainsS4)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2)*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD))*(PRDACC/SD)#*(1-rho4) + #dPWUACCN.dt = - (b*(PWU/N)*ws*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ND/Tau)*(PWDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PWUACCN/NU) - (muSU)*(PWUACCN)#*(1-rho3) + #dPWDACCN.dt = - (b*(PWU/N)*ws*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PWUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PWDACCN/ND) - (muSD)*(PWDACCN)#*(1-rho4) + #dPRUACCN.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN/NU) - (muRU)*(PRUACCN)#*(1-rho3) + #dPRDACCN.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN/ND)- (muRD)*(PRDACCN)#*(1-rho4) + #dPWUACC210.dt = - (b*(PWU/N)*ws*(ACC210U)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ACC210U)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(PWDACC210/Tau)-(b*I*eps*(d1))*(PWUACC210)-(muNU*PWUACCN*(ACC210NU/NU)*I*b*eps)-(muSU*PWUACC210)-(PWUACC210*(1/(8*365))) + #dPWDACC210.dt = - (b*(PWU/N)*ws*(ACC210D)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ACC210D)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))-(PWDACC210/Tau)+(b*I*eps*(d1))*(PWUACC210)-(muND*PWDACCN*(ACC210ND/ND)*I*b*eps)-(muSD*PWDACC210)-(PWDACC210*(1/(8*365))) + #dPRUACC210.dt = - (b*(PRU/N)*h*(ACC210U)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210U)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(PRDACC210/Tau)-(b*I*eps*(d1))*(PRUACC210)-(muNU*PRUACCN*(ACC210NU/NU)*I*b*eps)-(muRU*PRUACC210)-(PRUACC210*(1/(8*365))) + #dPRDACC210.dt = - (b*(PRU/N)*h*(ACC210D)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ACC210D)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))-(PRDACC210/Tau)+(b*I*eps*(d1))*(PRUACC210)-(muND*PWDACCN*(ACC210ND/ND)*I*b*eps)-(muRD*PRDACC210)-(PRDACC210*(1/(8*365))) #dACC210U.dt = del*(0.9)+(ACC210D/Tau)-(ACC210U*b*I*eps*d1)-(ACC210NU*muNU*I*b*eps) - (ACC210U*(1/(8*365))) #dACC210D.dt = -(ACC210D/Tau)+ (ACC210U*b*I*eps*d1)-(ACC210ND*muND*I*b*eps) - (ACC210D*(1/(8*365))) #dACC210NU.dt = del*(0.9)*(ACC210NU/ACC210U)+(ACC210ND/Tau)-(ACC210NU*b*I*eps*d1)-(ACC210NU*muNU*I*b*eps)- (ACC210NU*(1/(8*365))) #dACC210ND.dt = -(ACC210ND/Tau)+ (ACC210NU*b*I*eps*d1)-(ACC210ND*muND*I*b*eps)- (ACC210ND*(1/(8*365))) - dPSUACCN2.dt = - (b*(PSU/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ND/Tau)*(PSDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PSUACCN2/NU)#*(1-rho3) - dPSDACCN2.dt = - (b*(PSU/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PSUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PSDACCN2/ND)- (muSD)*(PSDACCN2)#*(1-rho4) - dPRUACCN2.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN2/NU)#*(1-rho3) - dPRDACCN2.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN2/ND)#*(1-rho4) - dPSUACC.dt = - (b*(PSU/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PSUACCN2/NU)+(SD/Tau)*(PSDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PSUACC/SU)#*(1-rho3) - dPSDACC.dt = - (b*(PSU/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PSDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PSUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PSDACC/SD)- (muSD)*(PSDACC)#death rate included because of clearance by drugs rather than immune system - dPRUACC.dt = - (b*(PRU/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3) - dPRDACC.dt = - (b*(PRU/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PRDACC/SD)#*(1-rho4) - dPSUACCA.dt = - (b*(PSU/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(AD/Tau)*(PSDACCA/AD)+(b*I*eps*SU*eta)*(PSUACC/SU)-(b*I*eps*(meanstrainsS5)*AU*d3)*(PSUACCA/AU) - (del*att*AU)*(PSUACCA/AU) - dPSDACCA.dt = - (b*(PSU/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PSDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PSUACCA/AU)-(AD/Tau)*(PSDACCA/AD) - (del*att*AD)*(PSDACCA/AD) - (muSD)*(PSDACCA) - dPRUACCA.dt = - (b*(PRU/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(AD/Tau)*(PRDACCA/AD)+(b*I*eps*SU*eta)*(PRUACC/SU)-(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)- (del*att*AU)*(PRUACCA/AU) - dPRDACCA.dt = - (b*(PRU/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PRDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)-(AD/Tau)*(PRDACCA/AD) - (del*att*AD)*(PRDACCA/AD) + dPWUACCN2.dt = - (b*(PWU/N)*ws*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ND/Tau)*(PWDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PWUACCN2/NU)#*(1-rho3) + dPWDACCN2.dt = - (b*(PWU/N)*ws*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) +(b*I*eps*NU*d1)*(PWUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PWDACCN2/ND)- (muSD)*(PWDACCN2)#*(1-rho4) + dPRUACCN2.dt = - (b*(PRU/N)*h*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(NU)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) +(ND/Tau)*(PRDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN2/NU)#*(1-rho3) + dPRDACCN2.dt = - (b*(PRU/N)*h*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(ND)*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (b*I*eps*NU*d1)*(PRUACCN2/NU)+ (-(ND/Tau)-(ratr*b*I*eps*ND*rho2)-(muND*ratr*ND*I*b*eps) - (del*att*ND))*(PRDACCN2/ND)#*(1-rho4) + dPWUACC.dt = - (b*(PWU/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(SU*(meanstrainsS1+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(b*I*eps*NU*rho1)*(PWUACCN2/NU)+(SD/Tau)*(PWDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PWUACC/SU)#*(1-rho3) + dPWDACC.dt = - (b*(PWU/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(SD*(meanstrainsS2+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ratr*b*I*eps*ND*rho2)*(PWDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PWUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PWDACC/SD)- (muSD)*(PWDACC)#death rate included because of clearance by drugs rather than immune system + dPRUACC.dt = - (b*(PRU/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SU*(meanstrainsS3+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (b*I*eps*NU*rho1)*(PRUACCN2/NU)+(SD/Tau)*(PRDACC/SD)+(-(b*I*eps*SU*((d2*(meanstrainsS4+omega1))+(eta))) - (del*att*SU))*(PRUACC/SU)#*(1-rho3) + dPRDACC.dt = - (b*(PRU/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(SD*(meanstrainsS4+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) + (ratr*b*I*eps*ND*rho2)*(PRDACCN2/ND)+(b*I*eps*SU*d2*(meanstrainsS4+omega1))*(PRUACC/SU)+ (-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD))*(PRDACC/SD)#*(1-rho4) + dPWUACCA.dt = - (b*(PWU/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(AU*(meanstrainsS5+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(AD/Tau)*(PWDACCA/AD)+(b*I*eps*SU*eta)*(PWUACC/SU)-(b*I*eps*(meanstrainsS5)*AU*d3)*(PWUACCA/AU) - (del*att*AU)*(PWUACCA/AU) + dPWDACCA.dt = - (b*(PWU/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*ws*(AD*(meanstrainsS6+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PWDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PWUACCA/AU)-(AD/Tau)*(PWDACCA/AD) - (del*att*AD)*(PWDACCA/AD) - (muSD)*(PWDACCA) + dPRUACCA.dt = - (b*(PRU/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PRD/N)*h*(AU*(meanstrainsS5+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(AD/Tau)*(PRDACCA/AD)+(b*I*eps*SU*eta)*(PRUACC/SU)-(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)- (del*att*AU)*(PRUACCA/AU) + dPRDACCA.dt = - (b*(PRU/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K)) - (b*(PWD/N)*h*(AD*(meanstrainsS6+omega1))*(1/K)*(((PWU+PWD+PRU+PRD)/N)-K))+(ratr*b*I*eps*SD*eta)*(PRDACC/SD)+(b*I*eps*(omega1)*AU*d3)*(PRUACCA/AU)-(AD/Tau)*(PRDACCA/AD) - (del*att*AD)*(PRDACCA/AD) #dNUACC.dt = del+(ND/Tau)#-(b*I*eps*NU*(d1+rho1)) #dNDACC.dt = (b*I*eps*NU*d1)-(ND/Tau)#-(ratr*b*I*eps*ND*rho2) - return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dbeta.dt,dPSUACCN2.dt,dPSDACCN2.dt,dPRUACCN2.dt,dPRDACCN2.dt,dPSUACC.dt,dPSDACC.dt,dPRUACC.dt,dPRDACC.dt,dPSUACCA.dt,dPSDACCA.dt,dPRUACCA.dt,dPRDACCA.dt))) + return(list(c(dPWU.dt,dPWD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dbeta.dt,dPWUACCN2.dt,dPWDACCN2.dt,dPRUACCN2.dt,dPRDACCN2.dt,dPWUACC.dt,dPWDACC.dt,dPRUACC.dt,dPRDACC.dt,dPWUACCA.dt,dPWDACCA.dt,dPRUACCA.dt,dPRDACCA.dt))) }) } @@ -226,7 +205,7 @@ plot(out2[,1],I*eps,type="l", xlab="Time", ylab="I*eps") plot(out2[,1],(out2[,2]),type="l", xlab="Time", ylab="Population", main = paste("b=",b,", d=",d1)) lines(out2[,1],(out2[,3]),col="RED") -legend("bottomright",legend=c("PSU","PSD"),bty="n",lwd = c(2,1),col=c("BLACK","RED")) +legend("bottomright",legend=c("PWU","PWD"),bty="n",lwd = c(2,1),col=c("BLACK","RED")) plot(out2[,1],(out2[,4]),type="l", ylim = c(0,max(out2[,5])),xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30") lines(out2[,1],(out2[,5]),col="RED") diff --git a/batchRuns/ODESImmunity_hqx.R b/batchRuns/ODESImmunity_hqx.R new file mode 100644 index 0000000..3473a82 --- /dev/null +++ b/batchRuns/ODESImmunity_hqx.R @@ -0,0 +1,180 @@ +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[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) + +} + + + diff --git a/batchRuns/paramList.csv b/batchRuns/paramList.csv new file mode 100644 index 0000000..a6750da --- /dev/null +++ b/batchRuns/paramList.csv @@ -0,0 +1,3 @@ +No,birth,att,muND,muNU,nstrains,rho1,rho2,omega,ResCost,capacity,muWU,muWD,muRD,muRU,Tau,d1,d2,d3,bavgnull,seasonality,dryWidth,phase,competition,modelType +1,1,9.13242E-05,0.1,0.1,20,0.5,0.05,0.01,0.6,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0,1,0,12,oneLocus_nonSeasonal +2,1,9.13242E-05,0.1,0.1,20,0.5,0.05,0.01,0.6,10,0.006666667,0.99,0.006666667,0.006666667,20,1,1,1,0.1,0.6,1,0,12,oneLocus_nonSeasonal \ No newline at end of file diff --git a/batchRuns/submitODE.sbatch b/batchRuns/submitODE.sbatch new file mode 100755 index 0000000..2960d2a --- /dev/null +++ b/batchRuns/submitODE.sbatch @@ -0,0 +1,19 @@ +#!/bin/bash +#SBATCH --job-name=oneLoucs +#SBATCH --time=2:00:00 +#SBATCH --output=oneLoucs_%A_%a.out +#SBATCH --error=oneLoucs_%A_%a.err +#SBATCH --array=1-20 +#SBATCH --tasks=1 +#SBATCH --cpus-per-task=1 +#SBATCH --mem-per-cpu=30000 +#SBATCH --partition=broadwl +#SBATCH --mail-type=ALL +#SBATCH --mail-user=YOUREMAIL + +# Print this sub-job's task ID +echo "My SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID +module load R + +Rscript ODESImmunity_hqx.R paraList.csv $SLURM_ARRAY_TASK_ID ./ +