Remove all objects from workspace.

Load add-on packages - deSolve - contains lsoda function - differential equation solver.

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
##X = 4916632  # Susceptible cattle
##Y = 1        # Infection cattle
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,150,by=0.5)

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, ylab="", ylim = c(0,1), 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("topright", legend = c("susceptible","Infected","immune"), col=c(2,3,4), lwd=2, cex = 1)

Intervention

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)

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(S~time,data=output,type="l",lwd=2,col=2,ylab="",ylim=c(0,1),main=paste("compliance=",i, "efficiency=",ii))
  lines(I~time,data=output,type="l",lwd=2,col=3,ylab="")
  lines(R~time,data=output,type="l",lwd=2,col=4,ylab="")
  legend("right",legend=c("S","I","R"),col=c(2,3,4),lwd=2,cex=0.7)
title ("FMD dynamic vaccination intervention", outer=TRUE)

}