Appendix
Naree Ketusing
November 30, 2015
Remove all objects from workspace.
Load add-on packages - deSolve - contains lsoda function - differential equation solver.
Modelling Final Project
Construct SIR Model
FMD_model = function (timepoint,state,parameters)
{
S = state [1] # Susceptible bovine
I = state [2] # Infected bovine
R = state [3] # Immune bovine
with (as.list (parameters),
{ dS = (-beta*S*I) + (gamma*I) - (alpha*S)
dI = (+beta*S*I) - (gamma*I)
dR = (alpha*S)
results = c(dS,dI,dR)
list (results)}
)
}Parameter
Vac_efficency = 0
Vac_compliance = 0
beta_value = 0.115 # transmission rate
gamma_value = 0.1 # recovery_rate per day
alpha_value = Vac_efficency * Vac_compliance
parameter_list = c(beta=beta_value,gamma=gamma_value,alpha=alpha_value)Initial value
X = 4080805 # Susceptible cattle
Y = 835827 # Infection cattle (17% of total pop)
##X = 4916631 # Susceptible cattle
##Y = 1 # Infection cattle
## Total pop = 4916632
Z = 0 # Immune cattle
N = X+Y+Z
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)simulation duration
timepoint = seq (0,365,by=1)Run model
output = lsoda(initial_value,timepoint,FMD_model,parameter_list)Plot FMD Dynamic (Base_model)
par (mfrow=c(1,3), oma = c(0,0,1,0))
plot (S~time, data=output, type="l", lwd = 2, col = 2, ylab = "Susceptible bovine")
plot (I~time, data=output, type="l", lwd = 2, col = 3, ylab = "Infected bovine")
plot (R~time, data=output, type="l", lwd = 2, col = 4, ylab = "Immune bovine")
title ("FMD dynamic_Base model", outer=TRUE)
Combine plot
plot (S~time, data=output, type="l", lwd=2, col=2,xlab="days", ylab="prevalence", ylim = c(0,0.9), main ="FMD dynamic_Base model")
lines (I~time, data=output, type="l", lwd=2, col=3)
lines (R~time, data=output, type="l", lwd=2, col=4)
legend("right", legend = c("susceptible","Infected","immune"), col=c(2,3,4), lwd=2, cex = 1)
Intervention
par(mfrow = c(1,1), oma = c( 0, 0, 1, 0 ))
{Vac.eff = 0.2
vac.com = 1
timepoint = seq (0,500,by=1)
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = lsoda(initial_value,timepoint,FMD_model,parameter_list)
plot(I~time,data=output,type="l",lwd=2,col=510,ylab="",ylim=c(0.1,0.2),main=paste("compliance=",100, "efficiency=",20))
legend("right",legend=c("Infected"),col=c(510),lwd=2,cex=0.7)
title ("FMD dynamic vaccination intervention", outer=TRUE)
summary(output, at=500)
}
## S I R
## Min. 8.300000e-01 1.305000e-01 0
## 1st Qu. 8.647000e-01 1.305000e-01 0
## Median 8.688000e-01 1.312000e-01 0
## Mean 8.649000e-01 1.351000e-01 0
## 3rd Qu. 8.695000e-01 1.353000e-01 0
## Max. 8.695000e-01 1.700000e-01 0
## N 5.010000e+02 5.010000e+02 501
## sd 8.195713e-03 8.195713e-03 0{Vac.eff = 0.2
vac.com = 1
timepoint = seq (0,2000,by=1)
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = lsoda(initial_value,timepoint,FMD_model,parameter_list)
plot(I~time,data=output,type="l",lwd=2,col=510,ylab="",ylim=c(0.1,0.2),main=paste("compliance=",100, "efficiency=",20))
legend("right",legend=c("Infected"),col=c(510),lwd=2,cex=0.7)
title ("FMD dynamic vaccination intervention", outer=TRUE)
summary(output, at=500)
}
## S I R
## Min. 8.300000e-01 1.304000e-01 0
## 1st Qu. 8.695000e-01 1.304000e-01 0
## Median 8.696000e-01 1.304000e-01 0
## Mean 8.684000e-01 1.316000e-01 0
## 3rd Qu. 8.696000e-01 1.305000e-01 0
## Max. 8.696000e-01 1.700000e-01 0
## N 2.001000e+03 2.001000e+03 2001
## sd 4.563583e-03 4.563583e-03 0Intervention
Vac.eff = c(0.2,0.4,0.6,0.8,0.9,1)
vac.com = c(0.5,0.6,0.7,0.8,0.9,1)
timepoint = seq (0,365,by=1)
par(mfrow = c(2,2), oma = c( 0, 0, 1, 0 ))
for (i in vac.com)
for (ii in Vac.eff) {
Vac_compliance = i
Vac_efficency = ii
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = lsoda(initial_value,timepoint,FMD_model,parameter_list)
plot(I~time,data=output,type="l",lwd=2,col=510,xlab="",ylab="Prevalence",ylim=c(0,0.2),main=paste("compliance=",i, "efficiency=",ii))
title ("FMD dynamic vaccination intervention", outer=TRUE)
}
Intervention (Vaccine 20% efficacy)
Run model
timepoint = seq (0,120,by=1)
par(mfrow = c(1,1), oma = c( 0, 0, 1, 0 ))
for (i in 1:6){
if (i==1) Vac_compliance = 0.5
if (i==2) Vac_compliance = 0.6
if (i==3) Vac_compliance = 0.7
if (i==4) Vac_compliance = 0.8
if (i==5) Vac_compliance = 0.9
if (i==6) Vac_compliance = 1
Vac_efficency = 0.2
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = as.data.frame(lsoda(initial_value,timepoint,FMD_model,parameter_list))
if (i==1){
plot(output$time,output$I,type="l",lwd=2,col=1,xlab="Time(days)",ylab="Prevalence",ylim=c(0,0.2),
main="vaccination efficacy 20%")
}else{
lines(output$time,output$I,lwd=2,col=i)
}
}
legend("topright",title="% compliance", legend=c("50","60","70","80","90","100"),bty="n",lwd=2,col=c(1,2,3,4,5,6))
Intervention (Vaccine 60% efficacy)
Run model
timepoint = seq (0,120,by=1)
par(mfrow = c(1,1), oma = c( 0, 0, 1, 0 ))
for (i in 1:6){
if (i==1) Vac_compliance = 0.5
if (i==2) Vac_compliance = 0.6
if (i==3) Vac_compliance = 0.7
if (i==4) Vac_compliance = 0.8
if (i==5) Vac_compliance = 0.9
if (i==6) Vac_compliance = 1
Vac_efficency = 0.6
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = as.data.frame(lsoda(initial_value,timepoint,FMD_model,parameter_list))
if (i==1){
plot(output$time,output$I,type="l",lwd=2,col=1,xlab="Time(days)",ylab="Prevalence",ylim=c(0,0.2),
main="vaccination efficacy 60%")
}else{
lines(output$time,output$I,lwd=2,col=i)
}
}
legend("topright",title="% compliance", legend=c("50","60","70","80","90","100"),bty="n",lwd=2,col=c(1,2,3,4,5,6))
Intervention (Vaccine 80% efficacy)
Run model
timepoint = seq (0,120,by=1)
par(mfrow = c(1,1), oma = c( 0, 0, 1, 0 ))
for (i in 1:6){
if (i==1) Vac_compliance = 0.5
if (i==2) Vac_compliance = 0.6
if (i==3) Vac_compliance = 0.7
if (i==4) Vac_compliance = 0.8
if (i==5) Vac_compliance = 0.9
if (i==6) Vac_compliance = 1
Vac_efficency = 0.8
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = as.data.frame(lsoda(initial_value,timepoint,FMD_model,parameter_list))
if (i==1){
plot(output$time,output$I,type="l",lwd=2,col=1,xlab="Time(days)",ylab="Prevalence",ylim=c(0,0.2),
main="vaccination efficacy 80%")
}else{
lines(output$time,output$I,lwd=2,col=i)
}
}
legend("topright",title="% compliance", legend=c("50","60","70","80","90","100"),bty="n",lwd=2,col=c(1,2,3,4,5,6))
Intervention (Vaccine 90% efficacy)
Run model
timepoint = seq (0,120,by=1)
par(mfrow = c(1,1), oma = c( 0, 0, 1, 0 ))
for (i in 1:6){
if (i==1) Vac_compliance = 0.5
if (i==2) Vac_compliance = 0.6
if (i==3) Vac_compliance = 0.7
if (i==4) Vac_compliance = 0.8
if (i==5) Vac_compliance = 0.9
if (i==6) Vac_compliance = 1
Vac_efficency = 0.9
alpha_value = Vac_efficency*Vac_compliance
initial_value = c(S = X*(1-alpha_value)/N, I = Y/N, R= (Z+X*alpha_value)/N)
output = as.data.frame(lsoda(initial_value,timepoint,FMD_model,parameter_list))
if (i==1){
plot(output$time,output$I,type="l",lwd=2,col=1,xlab="Time(days)",ylab="Prevalence",ylim=c(0,0.2),
main="vaccination efficacy 90%")
}else{
lines(output$time,output$I,lwd=2,col=i)
}
}
legend("topright",title="% compliance", legend=c("50","60","70","80","90","100"),bty="n",lwd=2,col=c(1,2,3,4,5,6))