SEIR Model

In the context Epidemics dynamics, there is a modeling approach called SEIR.

This techniques describes a process that has 4 stages:

SEIRS Model

• Susceptible Individuals - These are the individuals who have not been infected, but can be exposed to the infection
• Exposed Individuals - These are the individuals who have acquired the infection, but are still not transmitting the disease
• Infectious Individual - These are the individuals who have the pathogen and are active transmitting the disease
• Resolved Individual - These are the individuals that have either recovered or have died

In this stages, the following interactions can happen:

• Becoming Exposed : A Susceptible Individual can become an Exposed Individual
• Becoming Infected : An Exposed Individual can become an Infectious Individual
• Becoming Resolved: An Infectious Individual will be a Resolved Individual

Let’s assume that Recovered Individuals become immune to the pathogen, and therefore they cannot be a Susceptible Individual again. Indeed, the latter don’t even matter because we are going to focus on the number of infected - and for this, the only interactions that matter are Becoming Exposed and Becoming Infected.

Implementing a SEIR Model

Here you can see what I was able to implement from the paper’s differential equations. If you are reading this and the paper, and you have found something that doesn’t seems quite right, please get in touch with me!

seir_model_huwen <- function(
  start_date = NULL, # Start date of the simulation
  amount_of_days=30,
  days_with_new_r0=list(),
  initial_susceptible_individuals = 10374822, #https://en.wikipedia.org/wiki/Portuguese_people
  initial_infectious_individuals = 13,
  initial_exposed_individuals_factor = 20,
  initial_removed_individuals = 0,
  recovery_rate = 1/18,
  initial_r0 = 2.6,
  infection_rate = 1/4.2,
  intensive_care_beds = NULL){
  
  susceptible_individuals = rep(0, amount_of_days)
  exposed_individuals = rep(0, amount_of_days)
  infectious_individuals = rep(0, amount_of_days)
  removed_individuals = rep(0, amount_of_days)
  
  
  susceptible_individuals[1] = initial_susceptible_individuals 
  infectious_individuals[1] = initial_infectious_individuals
  exposed_individuals[1] = infectious_individuals[1]*initial_exposed_individuals_factor
  removed_individuals[1] = initial_removed_individuals
  
  current_transmission_rate = initial_r0 * recovery_rate
  
  for(step in 1:(amount_of_days-1)){
    if(step %in% names(days_with_new_r0)){
      current_transmission_rate = days_with_new_r0[[as.character(step)]]*recovery_rate
    }
    all_population = susceptible_individuals[step] + exposed_individuals[step] + infectious_individuals[step] + removed_individuals[step]
    
    new_exposed_individuals = current_transmission_rate*(susceptible_individuals[step]/all_population)*infectious_individuals[step]

    new_infectious_individuals = infection_rate*exposed_individuals[step]
    new_removed_individuals = recovery_rate*infectious_individuals[step]
    
    susceptible_individuals[step + 1] = susceptible_individuals[step] - new_exposed_individuals
    exposed_individuals[step + 1] = exposed_individuals[step] + new_exposed_individuals - new_infectious_individuals
    infectious_individuals[step + 1] = infectious_individuals[step] + new_infectious_individuals - new_removed_individuals
    removed_individuals[step + 1] = removed_individuals[step] + new_removed_individuals
  }
  
  is_valid_date = !is.na(lubridate::parse_date_time(start_date,"ymd"))
  if(!is.null(start_date) && is_valid_date){
    days = 0:(amount_of_days-1) + as.Date(start_date)
  }else{
    days = 1:amount_of_days
  }

  simulation_data = data.frame(
    days = days,
    susceptible_individuals = susceptible_individuals,
    exposed_individuals = exposed_individuals,
    infectious_individuals = infectious_individuals,
    removed_individuals = removed_individuals
  )
    
  if(!is.null(intensive_care_beds)){
    simulation_data$intensive_care_beds = rep(intensive_care_beds,amount_of_days)
  }
  
  return(simulation_data)
}

Simulating for 3 isolation measures

Now let’s run the simulation for 3 different isolation measures scenarios:

• No isolation measures (this occurred for the first weeks in Portugal, but around 13 March measures started to be taken)
• Moderate isolation measures: Let’s assume that at the beginning no isolation measures were taken. Then by day 20 of the epidemic, after the implementation of isolation measures, the \(R_0\) was able to drop to 2. The same thing happened after yet another set of isolation measures, dropping the \(R_0\) to 1.5. In day 50 it was able to reach an \(R_0\) of 1, and at day 65 the achieved \(R_0\) was 0.7.
• Aggressive isolation measures: Let’s assume that at the beginning no isolation measures were taken. Then by day 15, isolation measures were implemented that allowed the \(R_0\) to drop to 1.5. Then at day 30, more measures allowed the \(R_0\) to drop to 1. The last isolation measures implementation is this case would occur in day 45, dropping the \(R_0\) to 0.7.

Let’s run the simulation for a total of 4 months, and compare the number of infections over time - while taking into consideration that Portugal has around 500 intensive care beds (and some of those are currently being c by patients with other conditions than COVID-19)

moderate_decrese_transmission_rate = list(
  "20" = 2,
  "35" = 1.5,
  "50" = 1,
  "65" = 0.7
)

agressive_decrease_transmission_rate = list(
  "15" = 1.5,
  "30" = 1,
  "45" = 0,7
)

portuguese_intensive_care_beds_entire_population = 10374822*4.2/100000
initial_infectious_individuals = 13
amount_of_days = 30 * 9 # six months
start_date = "2020-03-02"

portuguese_simulation_covid19_without_isolation = seir_model_huwen(
      amount_of_days = amount_of_days,
      start_date = start_date,
      initial_infectious_individuals = initial_infectious_individuals,
      initial_susceptible_individuals = 10374822, #https://en.wikipedia.org/wiki/Portuguese_people
      initial_r0 = 3.1,
      intensive_care_beds = portuguese_intensive_care_beds_entire_population
    )
portuguese_simulation_covid19_with_moderate_isolation = seir_model_huwen(
      amount_of_days = amount_of_days,
      start_date = start_date,
      days_with_new_r0 = moderate_decrese_transmission_rate,
      initial_infectious_individuals = initial_infectious_individuals,
      initial_susceptible_individuals = 10374822, #https://en.wikipedia.org/wiki/Portuguese_people
      initial_r0 = 3.1,
      intensive_care_beds = portuguese_intensive_care_beds_entire_population
    )
portuguese_simulation_covid19_with_aggressive_isolation = seir_model_huwen(
      amount_of_days = amount_of_days,
      start_date = start_date,
      days_with_new_r0 = agressive_decrease_transmission_rate,
      initial_infectious_individuals = initial_infectious_individuals,
      initial_susceptible_individuals = 10374822, #https://en.wikipedia.org/wiki/Portuguese_people
      initial_r0 = 3.1,
      intensive_care_beds = portuguese_intensive_care_beds_entire_population
    )

#portugal has around 4.2 intensive care beds for each 100.000 persons


comparison_simulation_isolation_measures = data.frame(
  days = portuguese_simulation_covid19_without_isolation$days,
  infections_without_isolation = portuguese_simulation_covid19_without_isolation$infectious_individuals,
  infectious_with_moderate_isolation = portuguese_simulation_covid19_with_moderate_isolation$infectious_individuals,
  infectious_with_agressive_isolation = portuguese_simulation_covid19_with_aggressive_isolation$infectious_individuals,
  portugal_hospital_beds = portuguese_simulation_covid19_with_aggressive_isolation$intensive_care_beds
)

portuguese_simulation_plot <- plotly::subplot(
  plotly::plot_ly(comparison_simulation_isolation_measures, x = ~days) %>%
    plotly::add_lines(y = ~infections_without_isolation, name="Without Isolation", line = list(color = 'black')) %>%
    plotly::add_lines(y = ~portugal_hospital_beds, name="Portugal Intensive Care Beds", line = list(color = 'red')),
  plotly::plot_ly(comparison_simulation_isolation_measures, x = ~days) %>%
    plotly::add_lines(y = ~infectious_with_moderate_isolation, name="Moderate Isolation", line = list(color = 'orange')) %>%
    plotly::add_lines(y = ~infectious_with_agressive_isolation, name="Aggressive Isolation", line = list(color = 'blue')) %>%
    plotly::add_lines(y = ~portugal_hospital_beds, name="Portugal Intensive Care Beds", line = list(color = 'red')),
  nrows=2,
  shareX=T
)

Conclusion

We ran three simulations, with the main parameters of the model being obtained from this paper (which takes into consideration the situation in China).

It is important to understand that these are just simulations taking into consideration those parameters, and the current number of confirmed infectious people in Portugal by 2020-03-13.

The forecast potential of these simulation is, at the best of my knowledge, pretty fuzzy due to the lack of data currently to understand what would be more fitted parameters for Portugal.

Under these conditions, we can see that dropping the \(R_0\) by acting now - starting to take isolation measures now - can have a huge difference in the months to come.

Please, let’s be considerative these next weeks so that everyone of us can see each other again, and hi-5 people and hug people, as fast as possible.

We are in this together people.

Stay Safe!