Here, β represents the transmission rate, and N is the total population. 

Note: The transmission coefficient β is calculated using the formula:

 β = (average number of new infections / number of susceptible individuals) * (1 / duration of the time interval)

This formula gives you the average rate at which susceptible individuals become infected per unit of time. Let's illustrate this with an example:

Suppose you have data for a population of 10,000 individuals over a week. During that week, 300 new infections occurred. At the beginning of the week, there were 9,500 susceptible individuals.

Thus, average number of new infections = 300/7 = 42.86; & β = (42.86/9500)* (1/7) = 0.000476

So, in this example, the estimated transmission rate (β) is 0.000476 per susceptible individual per day.

Similarly, the flow of individuals from the exposed (E) to the infectious (I) compartment is described by:

where σ represents the rate at which exposed individuals become infectious. For example, if an individual becomes infectious after 4 days of exposure, σ will be 1/4 (or 0.25).

Once infectious, individuals either recover or succumb to the disease, captured by the equations:

Here, γ is the recovery rate, and μ is the mortality rate.

Note: If we divide the transmission rate (β) by the recovery rate (γ), it gives us R0 or the "basic reproduction number" - this number estimates the number of people who will be infected by the average infectious individual. In viral epidemiology;

• if R₀ > 1, then an outbreak of the virus is likely to become an epidemic

• If R₀ < 1, then an outbreak is likely to be contained

By solving these equations, we can simulate how a disease spreads through a population under different conditions. Let’s consider two conditions involving lockdown.

We first simulated the model under a "no lockdown" scenario (see code below), where the transmission rate (β) remained constant. The results were alarming but not unexpected. In a population of 1 million, infections peaked on day 112, with 126,516 individuals simultaneously infectious—a surge that would overwhelm any healthcare system. By the end of the simulation, 887,300 individuals recovered, 108,264 avoided infection, and 4,436 succumbed to the disease. This rapid rise and fall of infections demonstrated the virus’s natural progression when left unchecked.

library(deSolve)

library(ggplot2)

library(dplyr)

library(tidyr)

library(knitr)

# SEIR model function

SEIR <- function(time, current_state, params){

  with(as.list(c(current_state, params)), {

    N <- S + E + I + R

    dS <- -(beta * S * I) / N

    dE <- (beta * S * I) / N - sigma * E

    dI <- sigma * E - gamma * I - mu * I

    dR <- gamma * I

    dM <- mu * I

    return(list(c(dS, dE, dI, dR, dM)))

  })

}

# Parameters

params <- c(beta = 0.5, sigma = 0.25, gamma = 0.2, mu = 0.001)

# Initial state

initial_state <- c(S = 999999, E = 1, I = 0, R = 0, M = 0)

# Time points

times <- seq(0, 365, by = 1)

# Solve the model

model <- ode(initial_state, times, SEIR, params)

# Convert to data frame

model_df <- as.data.frame(model) %>%

  rename(Time = time) %>%

  pivot_longer(cols = -Time, names_to = "Compartment", values_to = "Population")

# Visualization

ggplot(model_df, aes(x = Time, y = Population, color = Compartment)) +

  geom_line(size = 1) +

  labs(title = "SEIR Model: Dynamics Over Time",

       x = "Time (Days)",

       y = "Population",

       color = "Compartment") +

  theme_minimal()

# Peak Infection

peak_infection <- model_df %>%

  filter(Compartment == "I") %>%

  slice(which.max(Population))

print(peak_infection)

Next, we introduced a lockdown into the model (see code) by reducing β during a 60-day period (day 90 to day 150). This adjustment reflected reduced human interactions during the lockdown. The impact was significant: infections peaked much later, on day 223, with a lower peak of 72,444 simultaneous infections. Additionally, the total number of deaths dropped to 4,195, and 156,886 individuals avoided infection entirely. The curve flattened, showcasing the immediate benefits of lockdowns in reducing strain on healthcare resources.

# Lockdown scenario

SEIR_lockdown <- function(time, current_state, params) {

  with(as.list(c(current_state, params)), {

    beta <- ifelse(time >= start_lockdown & time <= end_lockdown, 0.1, 0.5)

    N <- S + E + I + R

    dS <- -(beta * S * I) / N

    dE <- (beta * S * I) / N - sigma * E

    dI <- sigma * E - gamma * I - mu * I

    dR <- gamma * I

    dM <- mu * I

    return(list(c(dS, dE, dI, dR, dM)))

  })

}

# Lockdown Parameters

lockdown_params <- c(sigma = 0.25, gamma = 0.2, mu = 0.001, start_lockdown = 90, end_lockdown = 150)

lockdown_model <- ode(initial_state, times, SEIR_lockdown, lockdown_params)

# Convert to data frame

lockdown_df <- as.data.frame(lockdown_model) %>%

  rename(Time = time) %>%

  pivot_longer(cols = -Time, names_to = "Compartment", values_to = "Population")

# Comparison Visualization

combined_df <- bind_rows(

  model_df %>% mutate(Scenario = "No Lockdown"),

  lockdown_df %>% mutate(Scenario = "Lockdown")

)

ggplot(combined_df, aes(x = Time, y = Population, color = Compartment, linetype = Scenario)) +

  geom_line(size = 1) +

  labs(title = "SEIR Model: Impact of Lockdown",

       x = "Time (Days)",

       y = "Population",

       color = "Compartment",

       linetype = "Scenario") +

  theme_minimal()

# Summary Statistics

summary_stats <- combined_df %>%

  group_by(Scenario, Compartment) %>%

  summarize(

    Peak_Population = max(Population),

    Peak_Time = Time[which.max(Population)],

    .groups = "drop"

  )

kable(summary_stats, caption = "Summary Statistics of SEIR Model with and without Lockdown")