Statsproject

set.seed(1345)

Tmax <- 100
N <- 1000

p <- 1.5 / 1000
q <- 1 - p

S <- numeric(Tmax + 1)
I <- numeric(Tmax + 1)
A <- numeric(Tmax + 1)

immune <- rep(0, N)

S[1] <- 998
I[1] <- 2
A[1] <- 2
immune[1:2] <- 3

for (t in 1:Tmax) {

  immune <- immune - 1
  immune[immune < 0] <- 0

  susceptible_ids <- which(immune == 0)
  S[t] <- length(susceptible_ids)

  prob_infect <- 1 - q^I[t]

  new_inf <- rbinom(1, S[t], prob_infect)
  I[t + 1] <- new_inf

  if (new_inf > 0) {
    infected_new_ids <- sample(susceptible_ids, new_inf)
    immune[infected_new_ids] <- 3
  }

  A[t + 1] <- A[t] + new_inf
}

barplot(I,
        xlab = "Time",
        ylab = "New infections",
        main = "Simulation with 3-Step Immunity")

When immunity lasts only 3 steps, people return to the susceptible group fast. This leads to repeated waves of infection because people can get infected again soon after recovering.

set.seed(1345)

Tmax <- 100
N <- 1000

p <- 1.5 / 1000
q <- 1 - p

S <- numeric(Tmax + 1)
I <- numeric(Tmax + 1)
A <- numeric(Tmax + 1)

immune <- rep(0, N)

S[1] <- 998
I[1] <- 2
A[1] <- 2
immune[1:2] <- 8

for (t in 1:Tmax) {

  immune <- immune - 1
  immune[immune < 0] <- 0

  susceptible_ids <- which(immune == 0)
  S[t] <- length(susceptible_ids)

  prob_infect <- 1 - q^I[t]

  new_inf <- rbinom(1, S[t], prob_infect)
  I[t + 1] <- new_inf

  if (new_inf > 0) {
    infected_new_ids <- sample(susceptible_ids, new_inf)
    immune[infected_new_ids] <- 8
  }

  A[t + 1] <- A[t] + new_inf
}

barplot(I,
        xlab = "Time",
        ylab = "New infections",
        main = "Simulation with 8-Step Immunity")

When immunity lasts 8 steps, people stay immune for a longer time. Fewer people return to the susceptible group, so the waves are smaller, and the spread slows down compared to the 3-step model.