diff --git a/ODESImmunity2-copy.R b/ODESImmunity2-copy.R deleted file mode 100644 index 1e731db..0000000 --- a/ODESImmunity2-copy.R +++ /dev/null @@ -1,500 +0,0 @@ -require(deSolve) -#library(tidyverse) -SEAI<-function(t,y,p){ - PWU = y[1] - PWD = y[2] - PRU = y[3] - PRD = y[4] - #PR = y[6] - NU = y[5] - ND = y[6] - SU = y[7] - SD = y[8] - AU = y[9] - AD = y[10] - beta = y[11] - PWUACCN2 = y[12] - PWDACCN2= y[13] - PRUACCN2 = y[14] - PRDACCN2 = y[15] - PWUACC = y[16] - PWDACC = y[17] - PRUACC = y[18] - PRDACC = y[19] - 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 = (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) - #pfour = m/(K-3) - #pfive = m/(K-4) - #psix = m/(K-5) - #pseven = m/(K-6) - #peight = m/(K-7) - #pnine = m/(K-8) - #pten = m/(K-9) - ##I = pnbinom(K,size=K,prob=p) - pnbinom(0,size=K,prob=p) - #I = pnbinom(K,size=K,prob=pone) + pnbinom(K-1,size=K,prob=pone) + pnbinom(K-2,size=K,prob=pone) + pnbinom(K-3,size=K,prob=pone) + pnbinom(K-4,size=K,prob=pone) + pnbinom(K-5,size=K,prob=pone)+ pnbinom(K-6,size=K,prob=pone)+ pnbinom(K-7,size=K,prob=pone)+ pnbinom(K-8,size=K,prob=pone)+pnbinom(K-9,size=K,prob=pone) - pnbinom(0,size=K,prob=pone) - pnbinom(0,size=K-1,prob=ptwo) - pnbinom(0,size=K-2,prob=pthree) - pnbinom(0,size=K-3,prob=pfour) - pnbinom(0,size=K-4,prob=pfive) - pnbinom(0,size=K-5,prob=psix)- pnbinom(0,size=K-6,prob=pseven)- pnbinom(0,size=K-7,prob=peight)- pnbinom(0,size=K-8,prob=pnine)#-pnbinom(0,size=K-9,prob=pten) - #I = dnbinom(K,size=1,prob=pone) - kap = -1 - a =365/(2*pi) - ps = 0 - b = abs(beta) - S = 1-wr - k=12 - h = 1-(S/(1+0.1*exp(-k*(I-1)))) - # rho_inverse = vector(mode="numeric",length=nstrains) - # for (i in 1:nstrains){ - # rho_inverse[i] = 1/i - # } - #n_infections = nstrains*sum(rho_inverse) - #prob_infection = nstrains/n_infections - #prob_asymptomatic = 1/n_infections - #rho3=prob_asymptomatic - #rho4=rho3 - #h=wr - #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) - #f2 = frac2 #* - #f3 = frac3 - #f4 = frac4 - #totalm = m - prob23 = min(frac1,nstrains)/nstrains - 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)/(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 = ((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) - meanstrainsS1 = (1-rho3)^(frac1) - meanstrainsS2 = (1-rho3)^(frac2) - meanstrainsS3 = (1-rho3)^(frac1) - meanstrainsS4 = (1-rho3)^(frac2) - meanstrainsS5 = (1-rho3)^(frac3) - meanstrainsS6 = (1-rho3)^(frac4) - 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) - dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d2*(meanstrainsS1+omega1)+(eta))) - (del*att*SU) - dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2*(meanstrainsS1+omega1))-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD) - #dAU.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d) - (del*att*AU) - dAU.dt = (AD/Tau)+(b*I*eps*SU*eta)-(b*I*eps*(omega1)*AU*d3) - (del*att*AU) - #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)*(PWUACC)+(muSD)*PWDACC+muRU*PRUACC+muRD*PRDACC - #dHMOR.dt = del*att*del*t+(muNU*NUACC*I*b*eps)+(muND*ratr*NDACC*I*b*eps) - #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))) - 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(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))) - }) -} - -bavgnull = 1 - -nstrains = 20 - -del = 1 -att = 1/365/30#.000000000001 -muND = 0.1 -muNU = 0.1 -rho1 = 0.5#prob of immunity gain from NU to SU -rho2 = 0.5#prob of immunity gain from ND to SD -rho3 = 1/nstrains #prob of immunity gain from SU to AU -rho4 = rho3#=1/nstrains #prob of immunity gain from SD to AD -omega1 = 1/100 #prob of appearance of new antigen -omega2 = omega1 #prob of appearance of new antigen only in resistant strains - - -L = 1 - -ws = 1 -wr = 1-0.6 -#muSU = 1 -muSU=1/150#0.001 -#muSD = 20 -muSD=0.99 -#muRD = 1 -muRD = muSU#0.001 -#muRU = 1 -muRU = muSU#0.001 -Tau = 20#/365 -d1 = 1 -d2 = d1#/6 -d3 = d1#/10 -#m = 3 -K=10 -p2 = list(bavgnull=bavgnull,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K) - - -t2 = seq(from=0,to=40000,by=1) - - -N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200) -out2 = ode(y=N0,times=t2,func=SEAI,parms=p2) -P<-rowSums(out2[,2:5]) -N<-rowSums(out2[,6:11]) -#plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5])) -I<-1-P/N/K -PS<-rowSums(out2[,2:3]) -PR<-rowSums(out2[,4:5]) -eps<-(PS+wr*PR)/N -plot(out2[,1],I,type="l", xlab="Time", ylab="I") -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("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") -legend("bottomright",legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED")) - -plot(out2[,1],(out2[,12]),type="l", xlab = "Time", ylab="I", main = "b=0.2239*365,d=0.30") -plot(out2[,1],(out2[,13]),type="l", xlab = "Time", ylab="eps", main = "b=0.2239*365,d=0.30") - -#parasite host ratio -plot(out2[,1],apply(out2[,2:5],1,sum)/apply(out2[,6:11],1,sum),type="l", xlab="Time", ylab="parasite/hosts", main = paste("b=",b,", d=",d1)) - -plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,6])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30") -lines(out2[,1],(out2[,7]),col="RED") -legend("bottomright",legend=c("NU","ND"),bty="n",lwd=c(2,1),col=c("BLACK","RED")) - -plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,8])),type="l", xlab="Time",ylab="Population", main = "b=0.2239*365,d=0.30") -lines(out2[,1],(out2[,9]),col="RED") -legend("bottomright",legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED")) - -plot(out2[,1],out2[,10],type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30") -lines(out2[,1],out2[,11],col="RED") - -plot(out3[,1],out3[,12],xlim=c(0,2000)) - -plot(out2[,1],out2[,17],type="l") -plot(out2[,1],out2[,18],type="l") -plot(out2[,1],out2[,19],type="l") -plot(out2[,1],out2[,20],type="l") - -plot(out2[,1],out2[,21],type="l") -plot(out2[,1],out2[,22],type="l") -plot(out2[,1],out2[,23],type="l") -plot(out2[,1],out2[,24],type="l") - - - - -I1 = (1/(K*(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11])))*(-(((out3[,3]+out3[,2]+out3[,4] + out3[,5])/(out3[,6]+out3[,7]+out3[,8]+out3[,9]+out3[,10]+out3[,11]))-K)) - -plot(out2[,1],I1) - -out3eq = matrix(data=0,nrow=6,ncol=366) -out3eq2 = matrix(data=0,nrow=6,ncol=366) -#out3eq3 = vector(mode = "numeric", length = 22) -#out3eq4 = vector(mode = "numeric", length = 22) -#out3eq5 = vector(mode="numeric", length=22) -#out3eq6 = vector(mode="numeric", length=22) -out3eq7 = matrix(data=0,nrow=6,ncol=366) -m = matrix(data=0,nrow=6,ncol=366) -Ivec = vector(mode="numeric",length = 6) -mavg = vector(mode="numeric",length=6) -favg = vector(mode="numeric",length=6) -mnavg = vector(mode="numeric",length=19) -#K = c(0,1,2,3,4,5,6,7,8,9,10) -#out4 = vector(mode="numeric", length=22) -#bvalues = c(0,0.0001,0.0005,0.05,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.70,0.80,0.90,1) -#bvalues = c(-4,-3.5,-3,-2.8,-2.5,-2.2,-2,-1.9,-1.7,-1.5,-1.2,-1.1,-1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,2,4) -#bvalues = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24) -bvalues = seq(from=-1.95,to=-1.45,by=0.025) -wrvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1) -dvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1) -muSUvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1) -delvalues = c(10,20,40,50,60,70,90,100,150,200,250,300,350,400,450,500,550,600,700,800,900,1000) -rho3values = c(0.0005,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11,0.12) -omegavalues = c(0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.024,0.026,0.028,0.030,0.032,0.034,0.036,0.038,0.040,0.042,0.044,0.046,0.048,0.050,0.052) -for (i in 1:6){ - #del=delvalues[i] - bavgnull=10^(bvalues[i]) - #wr=wrvalues[i] - #d1=dvalues[i] - #muSU = muSUvalues[i]*0.1 - #rho3 = rho3values[i] - #omega1 = omegavalues[i] - wr=1-0.9 - K = 10 - p2 = list(bavgnull=bavgnull,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att) - N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,100,100,70,70)#,500,500,500,500,150,150,150,150,200,200,200,200) - out3 = ode(y=N0,times=t2,func=SEAI,parms=p2) - out3eq[i,] = out3[39635:40000,3]+out3[39635:40000,2] +out3[39635:40000,4] + out3[39635:40000,5] - #out3eq2[i,] = (out3[39635:40000,4]+out3[39635:40000,5])/(out3[39635:40000,2] + out3[39635:40000,3] + out3[39635:40000,4] + out3[39635:40000,5]) - #out3eq3[i] = out3[399999,6]+out3[399999,7] - #out3eq4[i] = out3[399999,8]+out3[399999,9] - #out3eq5[i] = out3[399999,10]+out3[399999,11] - #out4 [i] = (1/(K*(out3[399999,6]+out3[399999,7]+out3[399999,8]+out3[399999,9]+out3[399999,10]+out3[399999,11])))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K)) - #out4[i] = (1/(K*(K)))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K)) - #out3eq6[i] = (out3[3999,7]+out3[3999,9]+out3[3999,11])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]) - out3eq7[i,] = (out3[39635:40000,6]+out3[39635:40000,7]+out3[39635:40000,8]+out3[39635:40000,9]+out3[39635:40000,10]+out3[39635:40000,11]) - m[i,] = out3eq[i,]/out3eq7[i,] - mavg[i] = mean(1 - exp(-m[i,])) - #mnavg[i] = mean(1-exp(-(out3[39635:40000,17]+out3[39635:40000,18] +out3[39635:40000,19] + out3[39635:40000,20])/(out3[39635:40000,21]+out3[39635:40000,22]))) - #ptwo = m[i]/(K-1) - #pthree = m[i]/(K-2) - #pfour = m[i]/(K-3) - #pfive = m[i]/(K-4) - #psix = m[i]/(K-5) - #pseven = m[i]/(K-6) - #peight = m[i]/(K-7) - #pnine = m[i]/(K-8) - #pten = m[i]/(K-9) - #Ivec[i] = dnbinom(K,size=1,prob=pone) - #favg[i] = mean(out3eq2[i,]) -} - -#plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") -#plot(bvalues,1-out4)#,ylab="Proportion Resistant") -plot(bvalues,mavg)#,type="l") -#plot(bvalues,mavg,type="l") -#print(favg) -print(mavg) -plot(mavg,favg) - -#rho3values = c(1/20,1/50,1/100,1/250,1/400,1/500,1/750,1/1000,1/1250,1/1500) -rho3values = seq(from = 0.005, to = 0.05, by = 0.005) -mat6 = matrix(data=0,nrow=10,ncol=21) -etavalues = seq(from=0.01,to=0.21,by=0.04) -#mat77 = matrix(data=0,nrow=10,ncol=19) -#mat7 = matrix(data=0,nrow=10,ncol=6) -#mat8 = matrix(data=0,nrow=22,ncol=22) -#mat9 = matrix(data=0,nrow=22,ncol=22) -#mat10 = matrix(data=0,nrow=22,ncol=22) -#mat11 = matrix(data=0,nrow=22,ncol=22) -#mat12 = matrix(data=0,nrow=22,ncol=22) -for (i in 5:6){ - for (j in 1:21){ - #d = dvalues[i] - # A = Avalues[j] - #rho3=rho3values[i] - #rho4=rho3 - rho3 = rho3values[i] - rho4=rho3 - eta = 0.05 - #wr = wrvalues[i] - wr=1-0.9 - #muSU = muSUvalues[i]*0.01 - bavgnull = 10^(bvalues[j]) - #omega1 = omegavalues[j] - p2 = list(eta=eta,bavgnull=bavgnull,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att) - N0 = c(500,500,500,500,333,333,333,333,333,333,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200) - out5 = ode(y=N0,times=t2,func=SEAI,parms=p2) - # PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5] - # mat[i,j] = (out5[3999,4] + out5[3999,5])/PTotal - #mat2[i,j] = (out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])/6000 - #mat3[i,j] = out5[3999,6]/out5[3999,7] - mat6[i,j] = mean((out5[39635:40000,4]+out5[39635:40000,5])/(out5[39635:40000,2] + out5[39635:40000,3] + out5[39635:40000,4] + out5[39635:40000,5])) - #mat7[i,j] = mean(1-exp(-(out5[39635:40000,2]+out5[39635:40000,3]+out5[39635:40000,4]+out5[39635:40000,5])/((out5[39635:40000,6]+out5[39635:40000,7]+out5[39635:40000,8]+out5[39635:40000,9]+out5[39635:40000,10]+out5[39635:40000,11])))) - #mat77[i,j] = mean(1-exp(-(out5[39635:40000,17]+out5[39635:40000,18] +out5[39635:40000,19] + out5[39635:40000,20])/(out5[39635:40000,21]+out5[39635:40000,22]))) - #mat8[i,j]=(out5[3999,10]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - #mat9[i,j] = (out5[3999,8]+out5[3999,9])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - #mat10[i,j] = (out5[3999,6]+out5[3999,7])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - #mat11[i,j] = (out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - #mat12[i,j]=(out5[3999,7]+out5[3999,9]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - } -} -print(mat6) -#print(mat7) -#print(mat77) -#print(mat8) -#print(mat9) -#print(mat10) -#print(mat11) -#print(mat12) - -plot(bvalues,mat6[1,],type="l",ylim=c(0,1)) -lines(bvalues,mat6[2,]) -lines(bvalues,mat6[3,]) -lines(bvalues,mat6[4,]) -lines(bvalues,mat6[5,]) -lines(bvalues,mat6[6,]) -lines(bvalues,mat6[7,]) -lines(bvalues,mat6[8,]) -lines(bvalues,mat6[9,]) -lines(bvalues,mat6[10,]) - -plot(bvalues,mat77[9,],type="l",ylim=c(0,1)) -lines(bvalues,mat77[2,],col="RED") -lines(bvalues,mat77[3,],col="GREEN") -lines(bvalues,mat77[4,],col="BLUE") -lines(bvalues,mat77[5,]) -lines(bvalues,mat77[6,]) -lines(bvalues,mat77[7,]) -lines(bvalues,mat77[8,]) -lines(bvalues,mat77[9,]) -lines(bvalues,mat77[10,]) - -plot(mat77[1,],mat6[1,]) -plot(mat77[2,],mat6[2,]) -plot(mat77[3,],mat6[3,]) -plot(mat77[4,],mat6[4,]) -plot(mat77[5,],mat6[5,]) -plot(mat77[6,],mat6[6,]) - -plot(mat77[9,2:7],mat6[9,]) - -plot(mat77[1,9:12],mat6[1,9:12],type="l",ylim=c(0,1),xlim=c(0,1)) -lines(mat77[2,1:4],mat6[2,1:4]) -lines(mat77[3,5:8],mat6[3,5:8]) - -plot(bvalues,mat7[1,],type="l",ylim=c(0,1)) -lines(bvalues,mat7[2,]) -lines(bvalues,mat7[3,]) -lines(bvalues,mat7[4,]) -lines(bvalues,mat7[5,]) -lines(bvalues,mat7[6,]) -lines(bvalues,mat7[7,]) -lines(bvalues,mat7[8,]) -lines(bvalues,mat7[9,]) -lines(bvalues,mat7[10,]) - -plot(mat7[1,],mat6[1,],type="l",ylim=c(0,1)) -lines(mat7[2,],mat6[2,]) -lines(mat7[3,],mat6[3,]) -lines(mat7[4,],mat6[4,]) -lines(mat7[5,],mat6[5,]) -lines(mat7[6,],mat6[6,]) -lines(mat7[7,],mat6[7,]) -lines(mat7[8,],mat6[8,]) -lines(mat7[9,],mat6[9,]) -lines(mat7[10,],mat6[10,]) - - - - - - - -mat7vec = c(mat7) -mat7.df <- as.data.frame(c(t(mat7vec))) -mat73.df <-cbind(mat7.df,vector50,vector62) - -mat77vec = c(mat77) -mat77.df <-as.data.frame(c(t(mat77vec))) -mat773.df <-cbind(mat77.df,vector50,vector62) -#mat6[3,12] = 0 -#mat6[1,12] = 9.395332e-02 -#mat6[5,11] = 0 -#mat6[6,11] = 0 -mat6vec <- c(mat6) -mat6.df <- as.data.frame(c(t(mat6vec))) -vector50 <- c(rho3values,rho3values,rho3values,rho3values,rho3values,rho3values)#,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values,rho3values) -mat62.df <- cbind(mat6.df,vector50) - -vector60 = matrix(data=0,nrow=6,ncol=10) -for(i in 1:6){ - for(j in 1:10){ - vector60[i,] <- rep(bvalues[i],10) - } -} -vector62 <- c(t(vector60)) -mat63.df <- cbind(mat62.df,vector62) - -rows <- rownames(mat63.df) - -columns <- colnames(mat63.df) -log_transmission_intensity <- vector62 -nstrains_inverse <- vector50 -prevalence <- mat7vec -frequency_resistant <- mat6vec -naive_prevalence <- mat77vec -ggplot(data=mat63.df,mapping=aes(x=log_transmission_intensity,y=frequency_resistant))+geom_line(mapping=aes(color=as.factor(nstrains_inverse))) - -ggplot(data=mat73.df,mapping=aes(x=log_transmission_intensity,y=prevalence))+geom_line(mapping=aes(color=as.factor(nstrains_inverse))) - -mixed.df <-cbind(mat63.df,mat7.df) - -ggplot(data=mixed.df,mapping=aes(x=prevalence,y=frequency_resistant))+geom_line(mapping=aes(color=as.factor(nstrains_inverse)))+geom_smooth()+geom_point() - - -ggplot(data=mat773.df,mapping=aes(x=log_transmission_intensity,y=naive_prevalence))+geom_line(mapping=aes(color=as.factor(nstrains_inverse))) - -ggplot(data=mixed.df,mapping=aes(x=prevalence,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity)))+geom_smooth() - -mixed2.df <- cbind(mat63.df,mat77.df) - -ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity))) - -ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=frequency_resistant)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity))) - -ggplot(data=mixed.df,mapping=aes(x=nstrains_inverse,y=prevalence)) +geom_line(mapping=aes(color=as.factor(log_transmission_intensity))) - -ggplot(data=mixed2.df,mapping=aes(x=naive_prevalence,y=frequency_resistant)) + geom_point()+geom_smooth() - -vec600 = vector(mode="numeric",length=6) -vec700 = vector(mode="numeric",length=6) -vec770 = vector(mode="numeric",length=6) -for (i in 1:6){ - vec600[i] = mat6[i,i] - vec700[i] = mat7[i,i] - vec770[i] = mat77[i,i] -} -plot(vec700,vec600) - -require(ggplot2) \ No newline at end of file diff --git a/ODESSeasonality.R b/ODESSeasonality.R index a5aa9b9..054690e 100644 --- a/ODESSeasonality.R +++ b/ODESSeasonality.R @@ -15,6 +15,18 @@ SEAI<-function(t,y,p){ AD = y[10] x = y[11] beta=y[12] + PSUACCN2 = y[13] + PSDACCN2= y[14] + PRUACCN2 = y[15] + PRDACCN2 = y[16] + PSUACC = y[17] + PSDACC = y[18] + PRUACC = y[19] + PRDACC = y[20] + PSUACCA = y[21] + PSDACCA = y[22] + PRUACCA = y[23] + PRDACCA = y[24] with(as.list(p),{ #gammainit = vector() @@ -150,16 +162,35 @@ SEAI<-function(t,y,p){ #L = 10 #kap = 1 a =365/(2*pi) - b = abs(beta)#+100-(L) + ps =0 + b = abs(beta)#+(bavgnull*L)+bavgnull) #b =cos(kap*sin(t/a)-(t/a))-L + frac1 = (PSUACC+PRUACC)/(SU) + frac2 = (PSDACC+PRDACC)/(SD) + frac3 = (PSUACCA+PRUACCA)/AU + frac4 = (PSDACCA+PRDACCA)/AD + + + meanstrainsS1 = (1-rho3)^(frac1) + meanstrainsS2 = (1-rho3)^(frac2) + meanstrainsS3 = (1-rho3)^(frac1) + meanstrainsS4 = (1-rho3)^(frac2) + meanstrainsS5 = (1-rho3)^(frac3) + meanstrainsS6 = (1-rho3)^(frac4) 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) - ratr = ((PRU+PRD)*wr)/(PSU+PSD+((PRU+PRD)*wr)) + + S = 1-wr + k=12 + h = 1-(S/(1+0.1*exp(-k*(I-1)))) + + + ratr = ((PRU+PRD)*h)/(PSU+PSD+((PRU+PRD)*h)) ratr1 = (PRU+PRD)/(PSU+PSD+PRU+PRD) - rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*wr)) + rats = (PSU+PSD)/(PSU+PSD+((PRU+PRD)*h)) rats1 = (PSU+PSD)/(PSU+PSD+PRU+PRD) #P1 = ratr1 #ps = 0.34 @@ -173,14 +204,25 @@ SEAI<-function(t,y,p){ #B = cr*(gamma^2) #sigma1 = (ratr*A)+(rats*B) #sigma2 = (rats*(1-A)) - - eps = (PSD+PSU+((PRD+PRU)*(wr)))/(PSD+PSU+PRD+PRU) + n=2 + gammavec = c(((n-1)/n),x) + gammavec2 = c(0,x) + if (x<0){ + gamma = max(gammavec2) + }else{ + gamma = min(gammavec) + } + sigma = 0.1 + plus = (1-sigma)*n*gamma+(1-n*gamma) + eps = (PSD+PSU+((PRD+PRU)*(h)))/(PSD+PSU+PRD+PRU) r = 0.5*(1-(((m*exp(-m)))/(1-exp(-m)))) + #r=0 upsilon = (1/K)*(K-((PSU+PSD+PRU+PRD)/N)) - delR = b*upsilon*(PRU+PRD)*wr#*I + delR = b*upsilon*(PRU+PRD)*h#*I delS = b*upsilon*(PSU+PSD)*ws#*I - delP = delS +delR - n=2 + delP = delS+delR + delPx = delS*plus+delR + #TauRN1=((P1*((1-d1)+(d1*TauRD1)))-(d1*TauRD1))/(1-d1) #calculate TauRN1 from equation A4 #FUNC2<-function(smallp,E,S,Tau1,x1,f1,c1,h1,freqc1,freqh1){ # for(k in 1:cmax){ @@ -222,42 +264,50 @@ SEAI<-function(t,y,p){ #dPRU.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU)#-(del*att*(1/N)*(NU+SU+AU)*PRU) #dPRD.dt = -((1-sigma2)*b*((PRU+PRD)/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))-(sigma1*b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K))- (b*(PRU/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(PRD)#-(del*att*(1/N)*(ND+SD+AD)*PRD) - gammavec = c((n/n-1),x) - gammavec2 = c(0,x) - if (x<0){ - gamma = max(gammavec2) - }else{ - gamma = min(gammavec) - } #x1 = ((PRU+PRD+PSU+PSD)*x-(PRU+PRD)*(1-wr))/(PSD+PSU+((PRD+PRU)*(wr))) #x1 = (x2^n)/(1+x2^n) - x1 = gamma*(rats1)+ratr1*wr - dPSU.dt = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU - dPSD.dt = ((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD - dPRU.dt = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU - dPRD.dt = ((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD - e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD - er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD - es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD + x1 = gamma*(rats1)+ratr1*h + dPSU.dt = (1-sigma)*n*gamma*((NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))/N)*((1-r)*delS+r*delP*(1-x1^n))+ (1-n*gamma)*((NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU + dPSD.dt = (1-sigma)*n*gamma*((ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))/N)*((1-r)*delS+r*delP*(1-x1^n))+ (1-n*gamma)*((ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD + dPRU.dt = ((NU+SU*(meanstrainsS3+omega1)+AU*(omega1+meanstrainsS5))/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU + dPRD.dt = ((ND+SD*(meanstrainsS2+omega1)+AD*(omega1+meanstrainsS6))/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD + #e = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD+((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD + #er = ((NU+SU+AU)/N)*((1-r)*delR+r*delP*(x1^n))-muRU*PRU+((ND+SD+AD)/N)*((1-r)*delR+r*delP*(x1^n))-muRD*PRD + #es = ((NU+SU+AU)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSU*PSU+((ND+SD+AD)/N)*((1-r)*delS+r*delP*(1-x1^n))-muSD*PSD #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) - dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d2+rho3)) - (del*att*SU) - dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2)-(SD/Tau)-(ratr*b*I*eps*SD*rho4) - (del*att*SD) + dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d2*(meanstrainsS1+omega1)+(eta))) - (del*att*SU) + dSD.dt = (ratr*b*I*eps*ND*rho2)+(b*I*eps*SU*d2*(meanstrainsS1+omega1))-(SD/Tau)-(ratr*b*I*eps*SD*eta) - (del*att*SD) #dAU.dt = rat*(AU+AD)*(1-((AU+AD)/L))+(AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d) - (del*att*AU) - dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(b*I*eps*omega1*AU*d3) - (del*att*AU) + dAU.dt = (AD/Tau)+(b*I*eps*SU*eta)-(b*I*eps*(omega1)*AU*d3) - (del*att*AU) #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*rho4)+(b*I*eps*omega1*AU*d3)-(AD/Tau) - (del*att*AD) + dAD.dt = (ratr*b*I*eps*SD*eta)+(b*I*eps*(omega1)*AU*d3)-(AD/Tau) - (del*att*AD) #dgamma.dt = (1/(PSU+PSD+((1-S)*(PRU+PRD))))*((1/rats)-(gamma))*(dPRD.dt+dPRU.dt) #dgamma.dt = ratr + rats*((gamma)-ratr)-((1/(PSU+PSD+((1-S)*(PRU+PRD))))*(dPRD.dt+dPRU.dt)*((gamma)-ratr)) #dgamma.dt = -(dPRD.dt+dPRU.dt)/(dPSD.dt+dPSU.dt) #dx.dt = (1/((PSU+PSD+PRU+PRD)+(e)))*((er)-(x*(e))+((es)*((x*(PSD+PSD+PRU+PRD)-(PRU+PRD))/(PSU+PSD)))) - dx.dt = (delP*(x1-x)+delR*(x-1))/((PSU+PSD)+delS) - dbeta.dt = ((L*sin((t/a)-kap*sin(t/a))*(kap*cos(t/a)-1))/a) - return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt,dx.dt,dbeta.dt))) + dx.dt = (delPx*(x1-x)+delR*(x-1))/((PSU+PSD)+delS*plus) + dbeta.dt = ((L*bavgnull*sin(((t/a)+ps)-kap*sin((t/a)+ps))*(kap*cos((t/a)+ps)-1))/a) + dPSUACCN2.dt = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(NU/N) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(NU/N)+(ND/Tau)*(PSDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PSUACCN2/NU)#*(1-rho3) + dPSDACCN2.dt = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(ND/N) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(ND/N) +(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 = ((1-r)*delR+r*delP*(x1^n))*(NU/N)+(ND/Tau)*(PRDACCN2/ND)+(-(b*I*eps*NU*(d1+rho1))-(muNU*NU*I*b*eps) - (del*att*NU))*(PRUACCN2/NU)#*(1-rho3) + dPRDACCN2.dt = (((1-r)*delR+r*delP*(x1^n))*(ND))+ (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 = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(SU*(meanstrainsS1+omega1)/N) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(SU*(meanstrainsS1+omega1)/N)+(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 = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(SD*(meanstrainsS2+omega1)/N) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(SD*(meanstrainsS2+omega1)/N)+(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 = ((1-r)*delR+r*delP*(x1^n))*(SU*(meanstrainsS3+omega1)) + (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 = ((1-r)*delR+r*delP*(x1^n))*(SD*(meanstrainsS4+omega1)) + (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 = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(AU*(meanstrainsS5+omega1)/N) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(AU*(meanstrainsS5+omega1))+(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 = (1-sigma)*n*gamma*((1-r)*delS+r*delP*(1-x1^n))*(AD*(meanstrainsS6+omega1)) + (1-n*gamma)*((1-r)*delS+r*delP*(1-x1^n))*(AD*(meanstrainsS6+omega1))+(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 = ((1-r)*delR+r*delP*(x1^n))*(AU*(meanstrainsS5+omega1))+(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 = ((1-r)*delR+r*delP*(x1^n))*(AD*(meanstrainsS6+omega1))+(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,dx.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))) }) } + #b = 10 nstrains = 150 @@ -270,17 +320,18 @@ rho1 = 0.5#prob of immunity gain from NU to SU rho2 = 0.5#prob of immunity gain from ND to SD rho3 = 1/nstrains #prob of immunity gain from SU to AU rho4 = rho3#=1/nstrains #prob of immunity gain from SD to AD -omega1 = 5/100 #prob of appearance of new antigen +omega1 = 1/100 #prob of appearance of new antigen omega2 = omega1 #prob of appearance of new antigen only in resistant strains - +bavgnull = 0.1 L = 1 +kap = -1 ws = 1 -wr = 1-0.05 +wr = 1-0.6 S = 1-wr #muSU = 1 -muSU=1/250#0.001 +muSU=1/150#0.001 #muSD = 20 muSD=0.99 #muRD = 1 @@ -288,28 +339,33 @@ muRD = muSU#0.001 #muRU = 1 muRU = muSU#0.001 Tau = 20#/365 -d1 = 0.87#*0.3 +d1 = 1#*0.3 d2 = d1#/6 d3 = d1#/10 #m = 3 K=10 -p2 = list(ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K) +eta=0.05 +p2 = list(eta=eta,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K, bavgnull=bavgnull) tmax=40000 t2 = seq(from=0,to=tmax,by=1) -N0 = c(10000,10000,1,1,333,333,333,333,333,333,0.34,1) +N0 = c(500,500,500,500,333,333,333,333,333,333,0.34,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200) out2 = ode(y=N0,times=t2,func=SEAI,parms=p2) plot(out2[,1],(out2[,2]),ylim=c(0,max(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")) -plot(out2[,1],(out2[,4]),type="l", ylim=c(0,max(out2[,4])),xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30") +plot(out2[,1],(out2[,4]),type="l", ylim=c(0,max(out2[,5])),xlim=c(0,9000),xlab = "Time", ylab="Population", main = "b=0.2239*365,d=0.30") lines(out2[,1],(out2[,5]),col="RED") legend("bottomright",legend=c("PRU","PRD"),bty="n",lwd = c(2,1),col=c("BLACK","RED")) +plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(0,9000),main="2 loci, low seasonality") +lines(emp2,emp1) + + plot(out2[,1],(out2[,6]),ylim=c(0,max(out2[,6])),type="l", xlab="Time", ylab="Population", main = "b=0.2239*365,d=0.30") lines(out2[,1],(out2[,7]),col="RED") legend("bottomright",legend=c("NU","ND"),bty="n",lwd=c(2,1),col=c("BLACK","RED")) @@ -318,26 +374,35 @@ plot(out2[,1],(out2[,8]),ylim=c(0,max(out2[,9])),type="l", xlab="Time",ylab="Pop lines(out2[,1],(out2[,9]),col="RED") legend("bottomright",legend=c("SU","SD"),bty="n",lwd=c(2,1),col=c("BLACK","RED")) -plot(out2[,1],out2[,10],ylim=c(0,max(out2[,10])),type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30") +plot(out2[,1],out2[,10],ylim=c(0,max(out2[,11])),type="l",xlab="Time",ylab="Asymptomatic Population (A)", main="b=0.2239*365,d=0.30") lines(out2[,1],out2[,11],col="RED") plot(out2[,1],out2[,12],type="l") #print(out2[,12][35000]) -plot(out2[,1],-out2[,13],type="l",xlim=c(0,2000)) +plot(out2[,1],abs(out2[,13]),type="l",xlim=c(0,2000)) + +L=0.01 + +N02=c(10000,10000,1,1,333,333,333,333,333,333,0.34,bavgnull+(bavgnull*L)) +p2 = list(ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K, bavgnull=bavgnull) +out7=ode(y=N0,times=t2,func=SEAI,parms=p2) + +plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(0,700)) +lines(out7[,1],(out7[,4]+out7[,5])/(out7[,2]+out7[,3]+out7[,4]+out7[,5])) -plot(out2[,1],(out2[,4]+out2[,5])/(out2[,2]+out2[,3]+out2[,4]+out2[,5]),type="l",xlim=c(30000,40000),ylim=c(0.75,0.85)) #out3eq = vector(mode = "numeric", length = 18) -out3eq2 = vector(mode = "numeric", length = 19) +out3eq2 = vector(mode = "numeric", length = 18) #out3eq3 = vector(mode = "numeric", length = 18) #out3eq4 = vector(mode = "numeric", length = 18) #out3eq5 = vector(mode="numeric", length=18) -#out3eq7 = vector(mode="numeric", length=18) +out3eq7 = vector(mode="numeric", length=18) #out4 = vector(mode="numeric", length=22) #bvalues = c(0,0.0001,0.0005,0.05,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.70,0.80,0.90,1) bvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4) +bavg = vector(mode="numeric",length=13) #bvalues = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2) kapvalues = c(-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1) #bvalues = c(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24) @@ -348,11 +413,14 @@ wrvalues = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1) dvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1) muSUvalues = c(0,0.10,0.15,0.20,0.25,0.27,0.30,0.32,0.35,0.37,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,1) delvalues = c(10,20,40,50,60,70,90,100,150,200,250,300,350,400,450,500,550,600,700,800,900,1000) -rho3values = c(0.0005,0.002,0.004,0.005,0.007,0.008,0.01,0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11) +#rho3values = c(0.0005,0.002,0.004,0.005,0.007,0.008,0.01,0.02,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11) +rho3values = c(0.0005,0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01) omegavalues = c(0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.024,0.026,0.028,0.030,0.032,0.034,0.036,0.038,0.040)#,0.042,0.044,0.046,0.048,0.050,0.052) -Lvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4) -for (i in 1:19){ - #del=delvalues[i] +#Lvalues = c(-3,-2.5,-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.5,2)#,3,4) +Lvalues = c(0.04,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9) +bavgvalues = c(-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.5,2) +for (i in 14:18){ + #del=delvales[i] #b=10^(bvalues[i]) #wr=wrvalues[i] #d1=dvalues[i] @@ -361,81 +429,245 @@ for (i in 1:19){ #muSU = muSUvalues[i]*0.1 #rho3 = rho3values[i] #rho4=rho3 - L=10^Lvalues[i] + bavgnull = 10^bavgvalues[i] + #bavg[i] = bavgnull*L + #L=Lvalues[i] #kap=kapvalues[i] #omega1 = omegavalues[i] - p2 = list(b=b,ws=ws,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att,kap=kap) + p2 = list(eta=eta,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att, K=K, bavgnull=bavgnull) + N0 = c(500,500,500,500,333,333,333,333,333,333,0.34,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200) + out3 = ode(y=N0,times=t2,func=SEAI,parms=p2) #out3eq[i] = out3[39999,3]+out3[39999,2] +out3[39999,4] + out3[39999,5] - out3eq2[i] = (out3[39999,4]+out3[39999,5])/(out3[39999,2] + out3[39999,3] + out3[39999,4] + out3[39999,5]) + #out3eq2[i] = mean((out3[39635:40000,4]+out3[39635:40000,5])/(out3[39635:40000,2] + out3[39635:40000,3] + out3[39635:40000,4] + out3[39635:40000,5])) #out3eq3[i] = (out3[3999,6]+out3[3999,7])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]) #out3eq4[i] = (out3[39999,8]+out3[39999,9])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11]) #out3eq5[i] =(out3[39999,10]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11]) #out4 [i] = (1/(K*(out3[399999,6]+out3[399999,7]+out3[399999,8]+out3[399999,9]+out3[399999,10]+out3[399999,11])))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K)) #out4[i] = (1/(K*(K)))*(-(((out3[3999,3]+out3[3999,2]+out3[3999,4] + out3[3999,5])/(out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]))-K)) #out3eq6[i] = (out3[3999,6]+out3[3999,7]+out3[3999,8]+out3[3999,9]+out3[3999,10]+out3[3999,11]) - #out3eq7[i] = (out3[39999,7]+out3[39999,9]+out3[39999,11])/(out3[39999,6]+out3[39999,7]+out3[39999,8]+out3[39999,9]+out3[39999,10]+out3[39999,11]) + out3eq7[i] = mean(1-exp(-(out3[39635:40000,2] + out3[39635:40000,3] + out3[39635:40000,4] + out3[39635:40000,5])/(out3[39635:40000,6]+out3[39635:40000,7]+out3[39635:40000,8]+out3[39635:40000,9]+out3[39635:40000,10]+out3[39635:40000,11]))) } #plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") #plot(bvalues,1-out4)#,ylab="Proportion Resistant") -plot(Lvalues,out3eq2) +plot(bavgvalues,out3eq7) #print(out3eq2) -print(out3eq2) +print(out3eq7) +print(bavg) + +plot(out3eq7,out3eq2) -mat6 = matrix(data=0,nrow=11,ncol=19) +#mat62 = matrix(data=0,nrow=11,ncol=18) +mat645 = matrix(data=0,nrow=10,ncol=18) #mat7 = matrix(data=0,nrow=11,ncol=17) -#mat8 = matrix(data=0,nrow=11,ncol=19) +#mat8 = matrix(data=0,nrow=11,ncol=18) #mat9 = matrix(data=0,nrow=11,ncol=18) #mat10 = matrix(data=0,nrow=11,ncol=18) #mat11 = matrix(data=0,nrow=11,ncol=18) -#mat12 = matrix(data=0,nrow=11,ncol=18) -for (i in 1:11){ - for (j in 1:19){ +#mat123 = matrix(data=0,nrow=11,ncol=18) +for (i in 1){ + for (j in 1:18){ #d = dvalues[i] # A = Avalues[j] - #rho3=rho3values[i] - #rho4=rho3 + rho3=rho3values[i] + rho4=rho3 + #omega1 = omegavalues[i] #wr = wrvalues[i] - kap=kapvalues[i]#muSU = muSUvalues[i]*0.01 + wr = 1-0.9 + bavgnull = 10^bavgvalues[j] + #kap=kapvalues[i]#muSU = muSUvalues[i]*0.01 #b = 10^(bvalues[j]) - L=10^Lvalues[j] + #L=Lvalues[j] #omega1 = omegavalues[j] - p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att) + + p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att,bavgnull=bavgnull) + N0 = c(500,500,500,500,333,333,333,333,333,333,0.34,bavgnull+bavgnull*L,200,200,200,200,200,200,200,200,200,200,200,200)#,100,100,70,70,150,150,150,150,500,500,500,500)#,500,500,500,500,150,150,150,150,200,200,200,200) out5 = ode(y=N0,times=t2,func=SEAI,parms=p2) # PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5] # mat[i,j] = (out5[3999,4] + out5[3999,5])/PTotal #mat2[i,j] = (out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5])/6000 #mat3[i,j] = out5[3999,6]/out5[3999,7] - mat6[i,j] = (out5[39999,4]+out5[39999,5])/(out5[39999,2] + out5[39999,3] + out5[39999,4] + out5[39999,5]) + #mat62[i,j] = (out5[39999,4]+out5[39999,5])/(out5[39999,2] + out5[39999,3] + out5[39999,4] + out5[39999,5]) + mat645[i,j] = mean((out5[39635:40000,4]+out5[39635:40000,5])/(out5[39635:40000,2] + out5[39635:40000,3] + out5[39635:40000,4] + out5[39635:40000,5])) #mat7[i,j] = out5[3999,2]+out5[3999,3]+out5[3999,4]+out5[3999,5] #mat8[i,j]=(out5[39999,10]+out5[39999,11])/(out5[39999,6]+out5[39999,7]+out5[39999,8]+out5[39999,9]+out5[39999,10]+out5[39999,11]) #mat9[i,j] = (out5[3999,8]+out5[3999,9])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) #mat10[i,j] = (out5[3999,6]+out5[3999,7])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) #mat11[i,j] = (out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - #mat12[i,j]=(out5[3999,7]+out5[3999,9]+out5[3999,11])/(out5[3999,6]+out5[3999,7]+out5[3999,8]+out5[3999,9]+out5[3999,10]+out5[3999,11]) - } + #mat122[i,j]=(out5[39999,7]+out5[39999,9]+out5[39999,11])/(out5[39999,6]+out5[39999,7]+out5[39999,8]+out5[39999,9]+out5[39999,10]+out5[39999,11]) + #mat123[i,j]=(out5[39999,7]+out5[39999,9]+out5[39999,11])/(out5[39999,6]+out5[39999,7]+out5[39999,8]+out5[39999,9]+out5[39999,10]+out5[39999,11]) + + } } -print(mat6) +print(mat645) #print(mat7) #print(mat8) #print(mat9) #print(mat10) #print(mat11) -#print(mat12) +#print(mat123) -FLUC<-function(t,L,kap,a){ - L = 10^-3 +FLUC<-function(t,L,kap,beta,a){ + L = 1 + beta = 10^-0.5 kap = 1 - a =1 - b = L*cos(kap*sin(t/a)-(t/a))+L + a =1/(2*pi) + b = -L*beta*cos(kap*sin(t/a)-(t/a))+beta*2 return(b) } integrate(FLUC,0,1) +mat63vec = c(mat63) +mat123vec=c(mat123) + +threerowmatrix2 = cbind(mat63vec,mat123vec,vec6t,vec7t) +hostparasite3.df <-as.data.frame(threerowmatrix2) + +colnames(hostparasite.df)<-c("freqResist","propD","wr","tr") +colnames(hostparasite2.df)<-c("freqResist","propD","wr","tr") +colnames(hostparasite3.df)<-c("freqResist","propD","wr","tr") + + +loc<-rep(c(1,2,3),c(nrow(hostparasite.df),nrow(hostparasite2.df),nrow(hostparasite3.df))) +allDat<-rbind(hostparasite.df,hostparasite2.df,hostparasite3.df) +allDat$loc<-loc +ggplot(data=allDat,aes(log10(tr), freqResist,col=as.factor(loc)))+ + geom_line()+ + facet_wrap(.~wr)+ + theme_cowplot() +ggplot(data=allDat,aes(log10(tr),propD,col=as.factor(loc)))+ + geom_line()+ + facet_wrap(.~wr)+ + theme_cowplot() + + +tr = vec7t +wr = vec6t +freqResist = mat63vec +ggplot(data=hostparasite3.df,aes(tr, freqResist))+geom_point()+geom_line(aes(color=as.factor(wr))) + +mat62vec = c(mat62) +print(mat62vec) + +mat122vec = c(mat122) +print(mat122vec) + +tworowmatrix = cbind(mat6vec,mat12vec) +print(tworowmatrix) + + +hostparasite.df <- as.data.frame(tworowmatrix) + +ggplot(data=hostparasite.df)+geom_point(mapping=aes(x=mat6vec,y=mat12vec))+geom_point(mapping=aes(color=vec6)) + +vec7 = matrix(data=0,nrow=11,ncol=18) + +vec6 = matrix(data=0,nrow=11,ncol=18) +for(i in 1:11){ + for(j in 1:18){ + vec6[i,j]=wrvalues[i] + vec7[i,j]=10^bavgvalues[j] + } +} +print(vec6) +print(vec7) +vec6t=c(vec6) +print(vec6t) +vec7t=c(vec7) +print(vec7t) + +threerowmatrix=cbind(mat62vec,mat122vec,vec6t,vec7t) +print(threerowmatrix) + +hostparasite2.df <- as.data.frame(threerowmatrix) + +ggplot(data=hostparasite2.df,mapping=aes(x=freqResist,y=propD))+geom_point(mapping=aes(color=as.factor(tr)))+geom_line(aes(color=as.factor(wr))) + +tr = vec7t +wr = vec6t +freqResist = mat62vec + +require(tidyverse) +ggplot(data=hostparasite2.df,aes(tr, propD))+geom_point()+geom_line(aes(color=as.factor(wr))) +ggplot(data=hostparasite2.df,aes(tr, freqResist))+geom_point()+geom_line(aes(color=as.factor(wr))) + +ggplot(hostparasite2.df[hostparasite2.df$tr>0,],aes(wr, propD,color=as.factor(tr)))+geom_line() + + + + + + + + + + + +bavgvalues = c(-2,-1.5, -1,-0.9,-0.8,-0.6,-0.4,-0.2,-0.1,0,0.2,0.4,0.6,0.8,1,1.2,1.4,1.5,2) +t3 = seq(from=0,to=20*365,by=1) +trajectory = matrix(data=0,nrow=19,ncol=20) +tstag = seq(from=365,to=20*365,by=365) +for(i in 1:19){ + for(j in tstag){ + #d1=dvalues[i] + #d2=d1 + #d3=d1 + d1=0.02 + d2=d1 + d3=d1 + bavgnull=10^bavgvalues[i] + #rho3 = rho3values[i] + #rho4=rho3 + p2 = list(ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d1=d1,d2=d2,d3=d3,del=del,muND=muND,muNU=muNU,rho1=rho1,rho2=rho2,rho3=rho3,rho4=rho4,omega1=omega1,omega2=omega2,L=L,att=att,bavgnull=bavgnull) + N0 = c(1,1,1000,1000,333,333,333,333,333,333,0.34,(bavgnull)+(bavgnull*L)) + out10 = ode(y=N0,times=t3,func=SEAI,parms=p2) + trajectory[i,(j/365)] = (out10[j,4]+out10[j,5])/(out10[j,2] + out10[j,3] + out10[j,4] + out10[j,5]) + } +} +plot(tstag,trajectory[1,],type="l",xlim = c(365,3000),xlab="Time since removal of drug (days)",ylab = "frequency of resistance") +lines(tstag,trajectory[2,]) +lines(tstag,trajectory[3,]) +lines(tstag,trajectory[4,]) +lines(tstag,trajectory[5,]) +lines(tstag,trajectory[6,]) +lines(tstag,trajectory[7,]) +lines(tstag,trajectory[8,]) +lines(tstag,trajectory[9,]) +lines(tstag,trajectory[10,]) +lines(tstag,trajectory[11,]) +lines(tstag,trajectory[12,]) +lines(tstag,trajectory[13,]) +lines(tstag,trajectory[14,]) +lines(tstag,trajectory[15,]) +lines(tstag,trajectory[16,]) +lines(tstag,trajectory[17,]) +lines(tstag,trajectory[18,]) +lines(tstag,trajectory[19,]) + + +rows<- c("FR1","FR2","FR3","FR4","FR5","FR6","FR7","FR8","FR9","FR10","FR11","FR12","FR13","FR14","FR15","FR16","FR17","FR18","FR19") +columns<-c("t1","t2","t3","t4","t5","t6","t7","t8","t9","t10","t11","t12","t13","t14","t15","t16","t17","t18","t19","t20") +trajectoryvec = c(t(trajectory)) +trajectory2 = cbind(timesvec,vec6t2,trajectoryvec) +trajectories.df <- as.data.frame(trajectory2)#,row.names = columns, col.names=rows) +rownames(trajectories.df) +times=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20) + +times = matrix(data=0,nrow=20,ncol=17) +for (i in 1:20){ + for(j in 1:17){ + times[i,j] = i + } +} +print(times) +timesvec = c(times) +print(timesvec) +ggplot(data=trajectories.df,mapping=aes(x=timesvec,y=trajectoryvec))+geom_line(aes(color=as.factor(vec6t2)))+xlim(1,10) +#print(trajectory[1,])