Abstract_Berlin

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_base

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)

How prevalence decrease in 1 year

Vac_efficency = 0.2
compliance = seq(0.1,1,0.1)
data.store = matrix(NA,nrow = length(compliance), ncol=4)
for (i in 1:length(compliance)) {
  Vac_compliance = compliance[i]
  alpha_value = Vac_compliance*Vac_efficency
  parameter_list = c(beta=beta_value,gamma=gamma_value,alpha=alpha_value)
  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)
    data.store[i,]<- output[366,]
}
data.store

##       [,1]         [,2]         [,3]      [,4]
##  [1,]  365 7.871990e-04 6.799673e-15 0.9992128
##  [2,]  365 6.195207e-07 3.824713e-16 0.9999994
##  [3,]  365 6.135463e-10 2.102460e-16 1.0000000
##  [4,]  365 9.376707e-13 1.769369e-16 1.0000000
##  [5,]  365 5.358876e-15 1.543050e-16 1.0000000
##  [6,]  365 6.635720e-16 1.328439e-16 1.0000000
##  [7,]  365 3.004891e-16 1.201958e-16 1.0000000
##  [8,]  365 1.846611e-16 1.107966e-16 1.0000000
##  [9,]  365 1.306580e-16 1.045264e-16 1.0000000
## [10,]  365 9.938608e-17 9.938608e-17 1.0000000

plot (compliance,data.store[,3], xlab="compliance", ylab="prevalence")
lines (compliance,data.store[,3],type="l", lwd=2, col="purple")

Abstract_Berlin

Naree Ketusing

February 17, 2016

Modelling Final Project

Narges