Instructions

  1. Complete the self-guided coding activity, then complete the exercises below.
  2. Print and staple your .html file and submit a hard copy for grading.

Part 1

Recall the discrete epidemic model from Section 1.4 of the text, which describes a flu epidemic in a population of \(N=1000\) individuals. Note that the \(R_n\) equation has been omitted. \[ \begin{aligned} S_{n+1} &= S_n - a S_nI_n,\\ I_{n+1} &= I_n + a S_n I_n - b I_n, \end{aligned} \] where \(n\) is measured in weeks.

Exercise 1

Suppose that 5 people are initially infected and after one week 9 are infected. Assume no one has recovered initially and that the average length of infection is 2 weeks. Define R variables \(S_0\), \(I_0\), \(b\), and \(a\) based on this information (let R do the hard work for you!). What is the value of the transmission coefficient?

N <- 1000
I0 <- 5
S0 <- N - I0
b <- 0.5
a <- (9 - I0 + (b*I0)) / (S0 * I0)

print(a)
## [1] 0.001306533

The value of the transmission coefficient is 0.001306533.

Exercise 2

Solve the model over 20 weeks and plot the numerical solutions on the same set of axes. Include appropriate axes labels and a legend. Describe the outcome of the epidemic, e.g. how many people remain healthy, when infection peaks, etc.

N <- 1000
I0 <- 5
S0 <- N - I0
b <- 0.5
a <- 0.001306533

update_values <- function(S, I, a, b) {
  Sn1 <- S - a * S * I
  In1 <- I + a * S * I - b * I
  return(c(Sn1, In1))
}

weeks <- 20
S <- numeric(weeks + 1)
I <- numeric(weeks + 1)
S[1] <- S0
I[1] <- I0

for (i in 1:weeks) {
  new_values <- update_values(S[i], I[i], a, b)
  S[i + 1] <- new_values[1]
  I[i + 1] <- new_values[2]
}

plot(1:(weeks + 1), S, type = "l", col = "blue", xlab = "Weeks", ylab = "Population",
     main = "Epidemic Model Simulation", ylim = c(0, 1000), lwd = 2, yaxt = "n")
lines(1:(weeks + 1), I, col = "red", lty = 2, lwd = 2)
legend("topright", legend = c("Susceptible", "Infected"), col = c("blue", "red"), lty = 1:2, lwd = 2)

axis(2, at = seq(0, 1000, by = 100))
axis(1, at = seq(0, weeks + 1, by = 1))

The outcome of the epidemic is that approximately 300 people remain healthy when infection peaks at which point the number of healthy individuals drops below the number of infected individuals as infection decreases. Once the number of infected individuals drops to 100, the number of healthy individuals is greater than the number of infected individuals as infection steadily declines.

Part 2

Now suppose that removed individuals can be come susceptible again with parameter \(f\). Including the equation for \(R\), the system becomes \[ \begin{aligned} S_{n+1} &= S_n - a S_n I_n + fR_n,\\ I_{n+1} &= I_n + a S_n I_n - b I_n,\\ R_{n+1} &= R_n + b I_n - fR_n. \end{aligned} \]

Based on previous work, we know there are two equilibrium points of the form \((S^*,I^*)\): \((N,0,0)\) and \(\displaystyle\left(\frac{b}{a},\frac{f (N-b/a)}{f+b}, N-\frac{b}{a} - \frac{f(N-b/a)}{f+b} \right)\).

Exercise 3

  1. Consider all the same parameters as in Part 1, and suppose \(f = 0.2\). Calculate the values of each of these equilibrium points. (Use R or type your answers below.)
N <- 1000
S0 <- N - I0
I0 <- 5
b <- 0.5
a <- 0.001306533
f <- 0.2

equilibrium1 <- c(N, 0, 0)

equilibrium2 <- c(b/a, (f * (N - b/a)) / (f + b), N - b/a - (f * (N - b/a)) / (f + b))

cat("Equilibrium Point 1:", equilibrium1, "\n")
## Equilibrium Point 1: 1000 0 0
cat("Equilibrium Point 2:", equilibrium2, "\n")
## Equilibrium Point 2: 382.6922 176.3737 440.9341

The values of the first equilibrium are when S = 1000, I = 0, and R = 0, and the values of the second equilibrium are when S = 382.6923, I = 176.3736, and R = 440.9341.

  1. Plot the solution to the DDS over 52 weeks for each of the following sets of initial conditions. You should have exactly four graphs, each with three curves (S,I, and R). Use the same \(y\)-axis limits for all graphs, from \(0\) to \(N\). Include an appropriate legend and axes labels.
  • \(S_0 = 995,\,I_0 = 5,\,R_0=0\)
  • \(S_0 = 350,\,I_0 = 150,\,R_0 = 500\)
  • \(S_0 = 150,\,I_0 = 350,\,R_0 = 500\)
  • \(S_0 = 400,\,I_0 = 150,\,R_0 = 500\)
update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.2

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

Based on your results, how might you classify the stability of the two equilibrium points for this new model?

The second equilibrium point for this new model is stable because the second set of equilibrium values are always achieved regardless of the initial conditions whereas the first equilibrium point is unstable.

Exercise 4

Suppose there’s a magical potion that decreases the original transmission rate by 85%, i.e. \(a\) is now 15% of its original value. Repeat the previous exercise for the following initial conditions.

  • \(S_0 = 995,\,I_0 = 5,\,R_0=0\)
  • \(S_0 = 350,\,I_0 = 150,\,R_0 = 500\)
  • \(S_0 = 150,\,I_0 = 350,\,R_0 = 500\)
  • \(S_0 = 400,\,I_0 = 150,\,R_0 = 500\)
update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533 * 0.15
b <- 0.5
f <- 0.2

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

Based on your results, how might you classify the stability of the two equilibrium points with this modified transmission coefficient? How do you interpret these results?

With the modified transmission coefficient, the first equilibrium point is stable and the second equilibrium point is unstable as the values of the first equilibrium point are always achieved regardless of the initial conditions. It is clear from the results that modifying the transmission rate changes which equilibrium point the model approaches.

Exercise 5

Using a numerical approach, explore how the parameter \(f\) affects the long-term behavior of this DDS. Begin with the original values for \(N\), \(a\), and \(b\), and use a combination of initial conditions and \(f\) values to explore the model. Comment on and interpret your results.

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.1

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.1", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.2

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.2", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.3

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.3", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.4

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.4", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.5

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.5", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.6

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.6", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.7

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where. = 0.7", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.8

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.8", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 0.9

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 0.9", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

update_values <- function(S, I, R, a, b, f) {
  Sn1 <- S - a * S * I + f * R
  In1 <- I + a * S * I - b * I
  Rn1 <- R + b * I - f * R
  return(c(Sn1, In1, Rn1))
}

a <- 0.001306533
b <- 0.5
f <- 1

weeks <- 52
time_points <- 0:weeks

initial_conditions <- matrix(c(995, 5, 0, 
                               350, 150, 500, 
                               150, 350, 500, 
                               400, 150, 500), ncol = 3, byrow = TRUE)

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))

for (i in 1:4) {
  S <- numeric(weeks + 1)
  I <- numeric(weeks + 1)
  R <- numeric(weeks + 1)

  S[1] <- initial_conditions[i, 1]
  I[1] <- initial_conditions[i, 2]
  R[1] <- initial_conditions[i, 3]

  for (n in 1:weeks) {
    new_values <- update_values(S[n], I[n], R[n], a, b, f)
    S[n + 1] <- new_values[1]
    I[n + 1] <- new_values[2]
    R[n + 1] <- new_values[3]
  }

  plot(time_points, S, type = "l", col = "blue", ylim = c(0, N), ylab = "Population", xlab = "Weeks", main = paste("Scenario where f = 1", i))
  lines(time_points, I, col = "red")
  lines(time_points, R, col = "green")

  legend("topright", legend = c("Susceptible", "Infected", "Recovered"), col = c("blue", "red", "green"), lty = 1:2, lwd = 2, cex = 0.8)
}

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

It is clear from the above results that increasing the value of the parameter f results in the model achieving the second equilibrium point sooner except for when the parameter f = 0 when the model only achieves the first equilibrium point.