• SIR (Kermack and McKendrick 1927): an applicable epidemiological model.

• Many categories and extensions for SIR types of models. Its varieties can be found in social networks, stochastic interest rates, ageing population in ecology, knowledge diffusion in memetics etc.

• SIR stands for Susceptible (if previously unexposed to the pathogen), Infected (if currently colonized by the pathogen), and Recovered (if they have successfully cleared the infection).

• Consider transitions \(S \rightarrow I\) and \(I \rightarrow R\). For infections, it is generally observed that the amount of time spent in the infectious class (the “infectious period”) is distributed around some mean value.

• The percentage of total population: \(S + I + R =1\).

• Consider infectious period is Markovian. \(S \rightarrow I\), \(I \rightarrow R\) follow two seperated Poisson jump events.

• The force of infection \(\lambda\): is the per capita rate at which susceptible individuals contact the infection.

• Recovery rate \(\gamma\): \(\gamma I(t)\) is the infected proportion that recovers at time \(t\).

• In a time interval (\([t, t+ \Delta t)\)), the percentage of susceptible individuals is \[\Pr\{ \mbox{ susceptible individuals survive in }\Delta t \} =1-e^{-\lambda \Delta t}.\] By the relation between Poisson approximations and ODEs: \[ \frac{dS(t)}{dt} = -\lambda S (t).\\ \mbox{As } \left( \frac{ d (e^{-\lambda t})}{d t} = - \lambda e^{ - \lambda t} \right) = \left( \frac{dS(t)}{dt} = - \lambda S(t) \right) \mbox{ for } S(t) = e^{- \lambda t}.\]

• Endogenizing the parameter \(\lambda = \beta I(t)\), a fraction of infected proportion gives the force of infection. Then \[ \frac{dS(t)}{dt} = -\beta S (t) I(t).\]

• SIR model \[\frac{dS(t)}{dt} = - \beta S(t) I(t) \\ \frac{dI(t)}{dt} = \beta S(t)I(t) - \gamma I(t) \\ \frac{dR(t)}{dt} = \gamma I(t) \\ S(t)+I(t)+R(t)=1 \mbox{ for all }t.\]

solve.SIR = function(parms, times){
  SIR = function (times, init, parms) 
  {
       y = init; p = parms; 
       if (y[2]>1){y[2]=1}
       dy1 = - p["beta"] * y[1]* y[2]
       dy2 = p["beta"] *y[1]* y[2] - p["gamma"] * y[2]
       dy3 = p["gamma"] * (1- y[1] - y[3]) # y[1] + y[2] + y[3] = 1
       list(c(dy1, dy2, dy3))
  }
  state = c(0.99,0.01,0) # initial states
  return(ode(y=state,times=times, func=SIR, parms=parms))
}

library(deSolve)
parms=c(beta=0.05, gamma=0.01); times=seq(0,600,by=1)
init.mod=as.data.frame(solve.SIR(parms, times)); matplot(init.mod[,1],init.mod[,-1],type="l")

1 - Let the time increment between steps \(\Delta t\) be small and fixed.

2 - Transition probability amongst \(3\) variables \[\Pr (\mbox{Jump from S to I } = 1 ) \approx \beta S I , \, \, \Pr (\mbox{Jump from I to R } = 1 ) \approx \gamma I\]. We can summarize these into a matrix \[ [\{ \# \mbox{Jump from S to I} \}, \{ \# \mbox{Jumps from I to R} \}] \left[\begin{array}{ccc} -1 & 1 & 0\\ 0 & -1 & 1 \\ \end{array}\right] \\ = [\{ \# \mbox{Jumps leave S} \}, \{ \# \mbox{Jumps enter I and leave I} \}, \{ \# \mbox{Jumps enter R} \}].\]

3 - For \(\Delta t\), the transition random numbers are generated by \[\mbox{Transition from S to I } \sim \mathcal{Poi}\left( \beta S I \Delta t \right),\, \, \mbox{Transition from I to R } \sim \mathcal{Poi}\left( \gamma I \Delta t \right) \]

To summarize the previous scheme: \[S(t + \Delta t) = S(t) - \# \mbox{Jumps leave S in }\Delta t \\ I(t + \Delta t) = I(t) + \# \mbox{Jumps enter I in }\Delta t - \# \mbox{Jumps leave I in }\Delta t\\ R(t + \Delta t ) = R(t) + \# \mbox{Jumps enter R in} \Delta t\]

• Jump unit is \(1\). But \(S, I, R \leq 1\). Rescale \(S, I, R\) by \(1000\).

set.seed(2019); StateTrans = matrix(c(-1, 1, 0,
                                  0, -1, 1),nrow=2,ncol=3, byrow = TRUE)
mcState = list(); i = 1; vState = c(0.99,0.01,0)*1000; clock = 0; beta=0.05; gamma=0.01;
while(vState[2]>0 & vState[1]>0){
  mcState[[i]] = c(vState,clock)
  y1 = vState[1]; y2 = vState[2]; y3 = vState[3]
  vecLambda = c(beta*y1*y2, gamma*y2*1000); deltaTime = 1/sum(vecLambda)
  Jumps = rpois(2, vecLambda*deltaTime)
  vState = vState + Jumps%*%StateTrans
  vState[vState<0] = 0;  i = i+1;  clock = clock + deltaTime
}
y1mc = sapply(mcState, "[[", 1)/1000; y2mc = sapply(mcState, "[[", 2)/1000
y3mc = sapply(mcState, "[[", 3)/1000; vtime = sapply(mcState, "[[", 4)*1000

par(mfrow=c(1,3)); plot(vtime, y1mc, type="s"); lines(init.mod$time, init.mod$`1`, col="red")
plot(vtime, y2mc, type="s"); lines(init.mod$time, init.mod$`2`, col="red"); 
plot(vtime, y3mc, type="s"); lines(init.mod$time, init.mod$`3`, col="red")

• The previous Poisson approximation is a low-cost computation.

• In reality, the data is random and discrete. Similar to the stochastic Poisson approximation.

• But quite likely the underlying mechanism (law of motion or endogenized force) follows a deterministic structure.

• Estimate the stochastic model. Use the estimated parameters to simulate both deterministic and stochastic counterparts.

• Extract 10 observations from the stochastic simulation.

poiData = data.frame(vtime[seq(2,2000,by=200)], y1mc[seq(2,2000,by=200)], 
             y2mc[seq(2,2000,by=200)], y3mc[seq(2,2000,by=200)])
names(poiData)=c("time", "1", "2","3")
poiData

##          time     1     2     3
## 1    1.680672 0.989 0.011 0.000
## 2   78.591937 0.852 0.115 0.033
## 3  104.893173 0.715 0.210 0.075
## 4  122.950383 0.565 0.322 0.113
## 5  139.054743 0.440 0.389 0.171
## 6  155.644680 0.306 0.445 0.249
## 7  174.979037 0.192 0.477 0.331
## 8  199.842631 0.099 0.458 0.443
## 9  237.479727 0.043 0.332 0.625
## 10 306.152218 0.021 0.186 0.793

matplot(init.mod[,1],init.mod[,-1],type="l", ylab="Probability", xlab="Time")
points(poiData$time,poiData$`1`, col="black"); points(poiData$time,poiData$`2`, col="red")
points(poiData$time,poiData$`3`, col="green")

library(FME) 
obj.func = function(theta){
  times=seq(0,600,by=1)
  out = solve.SIR(parms=theta, times=times)
  return(modCost(obs=poiData, model=out))
}
initParm = c(beta=0.08, gamma=0.05);
fit=modFit(f=obj.func, p=initParm); fit$par

##        beta       gamma 
## 0.044318713 0.009685511

• Spatiotemporal dynamics happens for infection processes.

• Spatial redistribution of \(S\), \(I\), \(R\) of a fraction \(\epsilon\) of all individuals to other neighboring patches.

epsilon = 0.25; ny = nx = 10; xy = expand.grid(x = 1:nx, y = 1:ny)
dst  =  as.matrix(dist(xy)) #make distance matrix

(First law of geography: everything is related to everything else, but near things are more related than distant things.)

image(dst)

redist = matrix(0, nrow = ny*nx, ncol = ny*nx)
#populate the matrix so each of the 8 neighbors gets their share
redist[dst<1.5] = epsilon/8; diag(redist) = 1-epsilon #the remaining fraction
image(redist)

redist[1:15, 1:5];

##          [,1]    [,2]    [,3]    [,4]    [,5]
##  [1,] 0.75000 0.03125 0.00000 0.00000 0.00000
##  [2,] 0.03125 0.75000 0.03125 0.00000 0.00000
##  [3,] 0.00000 0.03125 0.75000 0.03125 0.00000
##  [4,] 0.00000 0.00000 0.03125 0.75000 0.03125
##  [5,] 0.00000 0.00000 0.00000 0.03125 0.75000
##  [6,] 0.00000 0.00000 0.00000 0.00000 0.03125
##  [7,] 0.00000 0.00000 0.00000 0.00000 0.00000
##  [8,] 0.00000 0.00000 0.00000 0.00000 0.00000
##  [9,] 0.00000 0.00000 0.00000 0.00000 0.00000
## [10,] 0.00000 0.00000 0.00000 0.00000 0.00000
## [11,] 0.03125 0.03125 0.00000 0.00000 0.00000
## [12,] 0.03125 0.03125 0.03125 0.00000 0.00000
## [13,] 0.00000 0.03125 0.03125 0.03125 0.00000
## [14,] 0.00000 0.00000 0.03125 0.03125 0.03125
## [15,] 0.00000 0.00000 0.00000 0.03125 0.03125

S = I = matrix(NA, nrow = ny*nx, ncol = 520)
S[,1] = 1; I[,1] = 0; I[100,1] = 1

local.dyn = function(t, S, I) {
  beta = 0.05*(1+0.5*cos(pi*t/20))
  I = S*(1-exp(-beta*I))
  S = 0.99*S- I + 5
  list(S = S,I = I)
}
for (t in 2:520) {
   #local growth:
   tmp = local.dyn(t,S=S[,t-1],I=I[,t-1])
   #spatial movement
   S[,t] = redist%*%tmp$S
   I[,t] = redist%*%tmp$I
}
x = xy[,1]; y = xy[,2]; Icubed = I^(1/4)/(max(I[,10:520]^(1/4)))

symbols(x,y,fg=2,circles=Icubed[,10],inches=FALSE,bg=2,xlab="",ylab="")

symbols(x,y,fg=2,circles=Icubed[,20],inches=FALSE,bg=2,xlab="",ylab="")

• Consider the following SIR model with births. Individuals are recruited directly into the susceptible class at birth. \[\frac{dS(t)}{dt} = \mu (1 - S(t)) - \beta S(t) I(t) \\ \frac{dI(t)}{dt} = \beta S(t)I(t) - (\mu + \gamma) I(t) \\ \frac{dR(t)}{dt} = \gamma I(t) - \mu R(t) \\ S(t)+I(t)+R(t)=1 \mbox{ for all }t.\] The new parameter \(\mu\) is a birth rate.

• Simulate it. Poissonize it. Estimate it. Couple it.

• END