Profile Picture

Simulations of Infected Individuals by COVID-19 in Portugal using the SEIR Model

Simulating some COVID-19 scenarios in Portugal with the SEIR Model

While on reddit, I came across this paper called “Phase-adjusted estimation of the number of Coronavirus Disease 2019 cases in Wuhan, China”, where the authors applied a model to rum some simulations and reason about if the public health interventions indeed suffice or not.

In the paper, the model used is the SEIR Model - that I had never heard about until today (well, I started this by 9PM on 2020-03-13, and now it is 2AM 2020-03-14)… Until I gave it a try, some googling and youtubing until I understand it.

In this post, I will apply this model with the same parameters as described in the paper to the Portuguese numbers, and check what are the results.

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

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.

Becoming Exposed

Every day, individuals from a population can become exposed. The only individuals that can become exposed are the susceptible individuals, and there are factors that influence this exposition:

  • For a susceptible individual to become exposed, the individual has to have contact with an infectious individual
  • Even when in contact with an infectious individual, there is a given transmission probability of the pathogen
  • As the more contacts a susceptible individual has during the day, the greater is the probability for contacting with an infectious individual
  • As the more contacts an infectious individual has during the day, the greater is the probability for contacting with a susceptible individual
  • Has more susceptible individuals became exposed, there are less individuals to become exposed out of the entire population

Becoming Infected

Every day, there will be exposed individuals that become infectious individuals. This amount of new infectious individuals is a process controlled by the parameter infectious_rate, which indicates the rate at which exposed individuals became infectious individuals.

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!

Modelling Transmission Rate

In this paper, the transmission_rate is modeled with two factors (as also shown in this paper, in the section R0 IN SIMPLE MODELS:

  • \(R_0\): this value represent the amount of individuals a given infectious person can infect (an average) by the begin of the epidemic (that0s why it is subscript with 0, representing time = 0)
  • Recovery Rate: To take into consideration the rate of infectious individuals that become recovery individuals

Therefore, we have:

\[ Transmission Rate = R_0 \times Recovery Rate \]

Simulating the change over time of Susceptible, Exposed, Infectious and Removed Cases

This simulation occurs for a given amount_of_days, where the \(R_0\) can vary over time (days_with_new_r0) - for example, the Portuguese Government decreed for schools to be closed, as well as disco bars, while taking other measures in restaurant and malls, which results in a lower \(R_0\) from that point on.

seir_model_huwen <- function(
  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/5.2){
  
  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){
    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
  }
  return(list(
    susceptible_individuals = susceptible_individuals,
    exposed_individuals = exposed_individuals,
    infectious_individuals = infectious_individuals,
    removed_individuals = removed_individuals
  ))
}

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)

portuguese_simulation_covid19_data_generation <- function(amount_of_days = 30, days_with_new_r0 = list()){
  

  simulation_huwen_seir_data = seir_model_huwen(
    amount_of_days = amount_of_days,
    days_with_new_r0 = days_with_new_r0
  )
  
  simulation_huwen_seir_df <- data.frame(
    days = 1:length(simulation_huwen_seir_data$susceptible_individuals) + as.Date("2020-03-05"),
    susceptible_individuals = simulation_huwen_seir_data$susceptible_individuals,
    exposed_individuals = simulation_huwen_seir_data$exposed_individuals,
    infectious_individuals = simulation_huwen_seir_data$infectious_individuals,
    removed_individuals = simulation_huwen_seir_data$removed_individuals
  )
    
  return(simulation_huwen_seir_df)
  
}
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
)

amount_of_days = 30 * 4 # six months

simulation_without_isolation = portuguese_simulation_covid19_data_generation(amount_of_days = amount_of_days)
simulation_with_moderate_isolation = portuguese_simulation_covid19_data_generation(amount_of_days = amount_of_days, days_with_new_r0 = moderate_decrese_transmission_rate)
simulation_with_aggressive_isolation = portuguese_simulation_covid19_data_generation(amount_of_days = amount_of_days, days_with_new_r0 = agressive_decrease_transmission_rate)

#portugal has around 4.2 intensive care beds for each 100.000 persons
portuguese_intensive_care_beds_entire_population = 10374822*4.2/100000

comparison_simulation_isolation_measures = data.frame(
  days = simulation_without_isolation$days,
  infections_without_isolation = simulation_without_isolation$infectious_individuals,
  infectious_with_moderate_isolation = simulation_with_moderate_isolation$infectious_individuals,
  infectious_with_agressive_isolation = simulation_with_aggressive_isolation$infectious_individuals,
  portugal_hospital_beds = rep(portuguese_intensive_care_beds_entire_population,length(simulation_without_isolation$susceptible_individuals))
)

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!

Comments