Skip to content

Commit

Permalink
Updated with incorporation of resistant proportion and differing d va…
Browse files Browse the repository at this point in the history 
…lues
  • Loading branch information
John Chaillet committed Oct 26, 2021
1 parent a4d4fca commit fdb35f5
Showing 1 changed file with 44 additions and 21 deletions.
65 changes: 44 additions & 21 deletions ODEsImmunity.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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]))
Expand Down Expand Up @@ -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)
Expand All @@ -114,38 +120,55 @@ 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]
# 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[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)






0 comments on commit fdb35f5

Please sign in to comment.