Remove all objects

Load package

## 
## Attaching package: 'deSolve'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot

Base model

Construct model dS/dt = -??SI - ??S + ØV + ??R dI/dt = +??SI - µI - ??I dR/dt = +??I - ??R dV/dt = +??S +¶i - ¶oV - ØV

Base_model = function (t,state,parms)
{
  S = state [1]
  I = state [2]
  R = state [3]
  V = state [4]
  with (as.list (parms),
        { dS = (-beta*S*I)-(alpha*S)+(slashed_O*V)+(omega*R)
          dI = (beta*S*I)-(mue*I)-(gamma*I)
          dR = (gamma*I)-(omega*R)
          dV = (alpha*S)+(pie_in)-(pie_out*V)-(slashed_O*V)
          results = c(dS,dI,dR,dV)
          list (results)})
}

Parameters

vac.eff = 0
vac.com = 0
beta = 0.115 # transmission rate (animals/d)
gamma = 0.1 # recovery rate (animals/d)
mue = 0.0041 # death rate (paper p'kachen)
alpha = vac.eff*vac.com
slashed_O = 0.0038 # rate of animals returned to their susceptible after the protected period from vaccination (animals/d)
omega = 0.0038 # rate of animals returned to their susceptible after the protected period from natural infection (animals/d)
pie_in = 0.5 # move in rate (animals/day)
pie_out = 0.5 # move out rate (animals/day)


parms = c(beta=beta, gamma=gamma, mue=mue, alpha=alpha,slashed_O=slashed_O, omega=omega, pie_in=pie_in, pie_out=pie_out)

Initial value

x = 4080805 #suscep cattle
y = 835827 # infected cattle (Prev. 17%)
z = 0 # Vaccinated cattle
w = 0 # recovery
N = x+y+w+z

initial_value = c(S=x*(1-alpha)/N, I=y/N, R= w/N, V=(x*alpha)/N)

Duration

timepoint =  seq(0,365,1)

Run model

#lsoda (y,times,func,parms)
output = lsoda(initial_value,timepoint,Base_model,parms)

Plot

plot(output, main=c("S","I","R","V"))

plot(S~time,output, col="green")
lines(I~time,output, col="red")
lines(R~time,output, col="purple")
lines(V~time,output, col="brown")

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)