From fdb35f584252dc79a450226176304786462b9d1e Mon Sep 17 00:00:00 2001 From: John Chaillet Date: Tue, 26 Oct 2021 12:30:32 -0400 Subject: [PATCH] Updated with incorporation of resistant proportion and differing d values --- ODEsImmunity.R | 65 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/ODEsImmunity.R b/ODEsImmunity.R index ec2303e..b901009 100644 --- a/ODEsImmunity.R +++ b/ODEsImmunity.R @@ -17,32 +17,33 @@ SEAI<-function(t,y,p){ N = y[5]+y[6]+y[7]+y[8]+y[9]+y[10] m = (PSU+PSD+PRU+PRD)/N I = 1-exp(-m) - rat = 0.8 + ratr = (PRU+PRD)/(PSU+PSD+PRU+PRD) + #rats = (PSU+PSD)/(PSU+PSD+PRU+PRD) eps = (((PSD+PSU+((PRD+PRU)*(wr))))/(PSD+PSU+PRD+PRU)) dPSU.dt = -(b*(PSU/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSD/N)*ws*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSU)*(PSU) dPSD.dt = -(b*(PSD/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PSU/N)*ws*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muSD)*(PSD) dPRU.dt = -(b*(PRU/N)*wr*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRD/N)*(NU+SU+AU)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRU)*(PRU) dPRD.dt = -(b*(PRD/N)*wr*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (b*(PRU/N)*(ND+SD+AD)*(1/K)*(((PSU+PSD+PRU+PRD)/N)-K)) - (muRD)*(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*(d+rho1))-(muNU*NU) - (del*att*NU) - dND.dt = (b*I*eps*NU*d)-(ND/Tau)-(b*I*eps*ND*rho2)-(muND*ND) - (del*att*ND) - dSU.dt = (b*I*eps*NU*rho1)+(SD/Tau)-(b*I*eps*SU*(d+rho3)) - (del*att*SU) - dSD.dt = (b*I*eps*ND*rho2)+(b*I*eps*SU*d)-(SD/Tau)-(b*I*eps*SD*rho4) - (del*att*SD) + 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*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) #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*d) - (del*att*AU) + dAU.dt = (AD/Tau)+(b*I*eps*SU*rho3)-(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 = (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) return(list(c(dPSU.dt,dPSD.dt,dPRU.dt,dPRD.dt,dNU.dt,dND.dt,dSU.dt,dSD.dt,dAU.dt,dAD.dt))) }) } -b = 9 +b = 0.9 nstrains = 20 att = 0.001#.000000000001 -del = 50 -muND = 0.1 +del = 60 +muND = 0.05 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 @@ -67,16 +68,18 @@ muRD = 0.05 #muRU = 1 muRU = 0.05 Tau = 20#/365 -d = 0.3 +d1 = 0.3 +d2 = 0.01 +d3 = 0.01 #m = 3 K=10 -p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,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) t2 = seq(from=0,to=5000,by=1) -N0 = c(500,500,500,500,333,333,333,333,333,333) +N0 = c(1,1,1000,1000,333,333,333,333,333,333) out2 = ode(y=N0,times=t2,func=SEAI,parms=p2) #plot(out2[,1],(out2[,2]+out2[,3]+out2[,4]+out2[,5])) @@ -106,6 +109,9 @@ plot(out2[,1],(out2[,8]+out2[,9])/(out2[,10]+out2[,11])) out3eq = vector(mode = "numeric", length = 22) out3eq2 = vector(mode = "numeric", length = 22) out3eq3 = vector(mode = "numeric", length = 22) +out3eq4 = vector(mode = "numeric", length = 22) +out3eq5 = vector(mode="numeric", length=22) +out3eq6 = 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(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) 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) @@ -114,25 +120,33 @@ for (i in 1:22){ b=bvalues[i]*1 #wr=wrvalues[i] #d=dvalues[i] - p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,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) out3 = ode(y=N0,times=t2,func=SEAI,parms=p2) out3eq[i] = out3[3999,2] + out3[3999,3]+out3[3999,4] + out3[3999,5] out3eq2[i] = (out3[3999,4]+out3[3999,5])/(out3[3999,2] + out3[3999,3] + out3[3999,4] + out3[3999,5]) - out3eq3[i] = out3[3999,10]+out3[3999,11] + out3eq3[i] = out3[3999,6]+out3[3999,7] + out3eq4[i] = out3[3999,8]+out3[3999,9] + out3eq5[i] = out3[3999,10]+out3[3999,11] + 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]) } -plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") -#plot(bvalues,out3eq3) +#plot(bvalues,out3eq2,xlab="b/10",ylab="Proportion resistant") +plot(bvalues,out3eq2)#,ylab="Proportion Resistant") +#print(out3eq2) print(out3eq2) -#print(out3eq3) mat6 = matrix(data=0,nrow=22,ncol=22) +#mat7 = matrix(data=0,nrow=22,ncol=22) +#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) for (i in 1:22){ for (j in 1:22){ #d = dvalues[i] # A = Avalues[j] wr = wrvalues[i] - b = bvalues[j]*10 + b = bvalues[j]*1 p2 = list(b=b,ws=ws,wr=wr,muSU=muSU,muSD=muSD,muRD=muRD,muRU=muRU,Tau=Tau,d=d,m=m) out5 = ode(y=N0,times=t2,func=SEAI,parms=p2) # PTotal = out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5] @@ -140,12 +154,21 @@ for (i in 1:22){ #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[3999,4]+out5[3999,5])/(out5[3999,2] + out5[3999,3] + out5[3999,4] + out5[3999,5]) - + #mat7[i,j] = out5[3999,2]+out5[3999,3]+out5[3999,4]+out5[3999,5] + #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]) } } print(mat6) +#print(mat7) +#print(mat8) +#print(mat9) +#print(mat10) +#print(mat11) + -