Cellular Automata

The necessary library

library(ggplot2)
library(grid)
library(dplyr)
library(magrittr)

Geral explanation

This is a cellular automata model inspired by the SIR (Susceptible-Infected-Recovered) epidemiological model. In our model, we have susceptible, infected, vaccinated and dead individuals in order to understand the effect of vaccination on the population in different epidemiological scenarios. In this simulation, we view cells as houses in a large neighborhood. The dynamics are determined as follows: susceptible cells can either get vaccinated or get infected. Vaccinated cells remain vaccinated. Infected cells can become susceptible again or die. Dead cells remain dead.

About the parameters

grid_size <- 100       # Size of the grid (100x100) equal to 10000 cells
n_steps <- 20          # Number of simulation steps
set.seed(13)           # Defining the random number sequence

About the probabilities

p_vaccinated is the probability of a cell be vaccinated
p_removed is the probability of a cell be removed
p_infected is the probability of a cell be infected
p_spread is the probability of spreading the infection to neighbors

Scenarios

Scenario I (Pro-vaccine Population)

In this scenario, the vaccination is widely accepted, the disease has a moderate mortality rate, the chance of initial infection is low and transmission to neighbors is moderate.

pro_vaccine <- list(p_vaccinated = 0.8,   
                    p_removed = 0.3,    
                    p_infected = 0.05,    
                    p_spread = 0.5)

Scenario II (Vaccine Rejection)

In this scenario, the vaccination is widely rejected, the disease has a moderate mortality rate, the chance of initial infection is low and transmission to neighbors is moderate.

vaccine_rejection <- list(p_vaccinated = 0.02, 
                          p_removed = 0.3,    
                          p_infected = 0.05,    
                          p_spread = 0.5)

Scenario III (Contagious Disease)

contagious <- list(p_vaccinated = 0.2,   
                   p_removed = 0.3,    
                   p_infected = 0.05,    
                   p_spread = 0.9)

Scenario IV (Epidemiological Collapse)

In this scenario, the vaccination is widely rejected, the disease has a high mortality rate, the chance of initial infection is low and transmission to neighbors is high.

collapse <- list(p_vaccinated = 0.1,   
                 p_removed = 0.9,    
                 p_infected = 0.05,    
                 p_spread = 0.95)

All the scenarios

scenarios <- list(pro_vaccine, vaccine_rejection, contagious, collapse)

Function to plot the city grid

plot_city <- function(city, step) {
df <- expand.grid(x = 1:grid_size, y = 1:grid_size)
df$state <- as.vector(city)
df$state <- factor(df$state, levels = c(0, 1, 2, 3), labels = c("Vaccinated", "Susceptible", "Infected", "Removed")) # Ensuring all states are included in the factor levels

ggplot(df, aes(x = x, y = y, fill = state)) +
    geom_tile() + 
    scale_fill_manual(
      values = c("white", "#DAF7A6", "#DE3163", "black"),  # Defining colors for states
      drop = FALSE  # Preventing ggplot from dropping missing states
    ) +
    labs(
      title = paste("Scenario 1 (Pro-vaccine population) - Step", step),
      fill = "State"
    ) +
    theme_minimal() +
    theme(
      axis.title = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      legend.position = "bottom",
      legend.text = element_text(size = 12, face = "italic"), 
      legend.title = element_text(size = 14, face = "bold"),  
      plot.title = element_text(size = 16, face = "bold")
    ) +
    scale_x_continuous(expand = c(0, 0)) +  # No space around grid
    scale_y_continuous(expand = c(0, 0))    # No space around grid
}

Function to spread infection to neighbors on an absorbent edge

spread_to_neighbors <- function(i, j, grid_size, city, p_spread, n_city) {
  neighbors <- list(
    c(i-1, j), c(i+1, j),      # up and down
    c(i, j-1), c(i, j+1),      # left and right
    c(i-1, j-1), c(i-1, j+1),  # upper diagonals
    c(i+1, j-1), c(i+1, j+1)   # lower diagonals
  )
  for (n in neighbors) {
    ni <- n[1]
    nj <- n[2]
     if (ni >= 1 && ni <= grid_size && nj >= 1 && nj <= grid_size) { # check if the neighbor is inside the grid (absorbing edge)
      if (city[ni, nj] == 1) {  # only susceptible may be infected
        if (runif(1) < p_spread) {
          n_city[ni, nj] <- 2  # neighbor become infected
        }
      }
    }
  }
  return(n_city)
}

Function to save data

data_saved <- function(city, step, grid_size) {
  df <- expand.grid(x = 1:grid_size, y = 1:grid_size)
  df$state <- as.vector(city)
  df$step <- step 
  return(df)
  }

Initializating the grid with all cells susceptible, with the central cell infected

city <- matrix(1, nrow = grid_size, ncol = grid_size)
city[grid_size/2,grid_size/2]=2

Viewing the initial condition

print(plot_city(city, 0))

Choosing which scenario we will simulate

Note that our example will be scenario I, pro-vaccine.

chosen_scenario <- scenarios[[1]] # The chosen scenario
print(chosen_scenario)

$p_vaccinated
[1] 0.8

$p_removed
[1] 0.3

$p_infected
[1] 0.05

$p_spread
[1] 0.5

saved_data_list <- list()

Simulation loop

for (step in 1:n_steps) {                                 # Do until the final "time"
  n_city <- city                                          # Create a copy of the city to update states
  for (i in 1:grid_size) {                                # Do for all lines
    for (j in 1:grid_size) {                              # Do for all columns
      if (city[i, j] == 0) {                              # If the cell be equal to 0 (vaccinated)
        n_city[i, j] <-  0                                # Stay 0 (vaccinated)
      } else if (city[i, j] == 3) {                       # If the cell be equal to 3 (dead)
        n_city[i, j] <-  3                                # Stay 3 (dead)
      } else if (city[i, j] == 1) {                       # If the cell be equal to 1 (susceptible)
        if (runif(1) < chosen_scenario$p_infected) {      # And the random number is within the probability of being infected 
          n_city[i, j] <-  2                              # The state changes to 2 (infected)
        } else if 
        (runif(1) < chosen_scenario$p_vaccinated) {       # If the random number is within the probability of being vaccinated
          n_city[i, j] <-  0                              # The state changes to 0 (vaccinated)
        }
    } else if (city[i, j] == 2) {                         # If the cell be equal to 2 (infected)
        if (runif(1) < chosen_scenario$p_removed) {       # If the random number is within the probability of being removed
          n_city[i, j] <- 3                               # The state changes to 3 (removed)
        } else {
          n_city <- spread_to_neighbors(i, 
                                        j, 
                                        grid_size, 
                                        city, 
                                        chosen_scenario$p_spread,
                                        n_city)           # spread the infection to susceptible neighbors
        }
      }
    }
  }
  city <- n_city                                          # Updating the matrix   
  saved_data_list[[step]] <- data_saved(city, 
                                        step, 
                                        grid_size) 
  print(plot_city(city, step))                            # Plot the current state of the city
  #Sys.sleep(0.5)                                         # Pause for a short time to visualize the simulation
 }

Viewing the results

results_scenario_one <- do.call(rbind, saved_data_list)
summary_one <- results_scenario_one %>%
  group_by(step) %>%
  summarise(
    vaccinated = sum(state == 0),
    susceptible = sum(state == 1),
    infected = sum(state == 2),
    removed = sum(state == 3)
  )

summary_one <- summary_one %>%
  bind_rows(data.frame(step = 0,
                       vaccinated = 0,
                       susceptible = 9999,
                       infected = 1,
                       removed = 0))

ggplot(summary_one) +
  geom_line(aes(x = step, y = vaccinated, color = "Vaccinated"), size = 1) +
  geom_line(aes(x = step, y = susceptible, color = "Susceptible"), size = 1) +
  geom_line(aes(x = step, y = infected, color = "Infected"), size = 1) +
  geom_line(aes(x = step, y = removed, color = "Removed"), size = 1) +
  labs(title = "Scenario 1 - (Pro-vaccine Population)",
       x = "Time",
       y = "Number of sites",
       color = "State") +  
  scale_color_manual(values = c("Vaccinated" = "#e5e8e8", 
                                "Susceptible" = "#DAF7A6", 
                                "Infected" =  "#DE3163", 
                                "Removed" = "black")) +  
  theme_minimal()+
  theme(
    panel.grid = element_blank(),
    axis.title = element_text(size = 12, face = "bold"),  
    axis.text = element_text(size = 10),                  
    axis.ticks = element_line(size = 0.5),
    legend.position = "bottom",
    legend.text = element_text(size = 12, face = "italic"), 
    legend.title = element_text(size = 14, face = "bold"),  
    plot.title = element_text(size = 16, face = "bold")
  )

In this scenario, we observed that the high rate of early vaccination reduced the number of susceptible sites and consequently the number of infected and dead sites.

This is the end of the scenario I simulation, but you can try the other scenarios.

See you!

--- title: "*Final Model - Mathematical Models Applied to Ecology*" author: "Amanda B Loureiro & Emanuel B S A Cambraia & Juliana R M Ciccheto" date: "26 Feb 2025" format: html editor: visual cap-location: bottom title-block-banner: true code-tools: true theme: custom.scss toc: true toc-title: "Sections" editor_options: chunk_output_type: console --- # Cellular Automata ### The necessary library ```{r, warning=FALSE, message=FALSE} #| Warning: false #| Message: false #| Include: false library(ggplot2) library(grid) library(dplyr) library(magrittr) ``` ### Geral explanation This is a **cellular automata model** inspired by the **SIR (Susceptible-Infected-Recovered) epidemiological model**. In our model, we have **susceptible, infected, vaccinated and dead** individuals in order to understand the **effect of vaccination** on the population in different epidemiological scenarios. In this simulation, we view cells as houses in a large neighborhood. The dynamics are determined as follows: susceptible cells can either get vaccinated or get infected. Vaccinated cells remain vaccinated. Infected cells can become susceptible again or die. Dead cells remain dead. ### About the parameters ```{r} grid_size <- 100 # Size of the grid (100x100) equal to 10000 cells n_steps <- 20 # Number of simulation steps set.seed(13) # Defining the random number sequence ``` ::: callout-tip ### About the probabilities **p_vaccinated** is the probability of a cell be vaccinated **p_removed** is the probability of a cell be removed **p_infected** is the probability of a cell be infected **p_spread** is the probability of spreading the infection to neighbors ::: # Scenarios ## Scenario I (Pro-vaccine Population) In this scenario, the **vaccination** is **widely accepted**, the **disease** has a **moderate mortality rate**, the chance of **initial infection** is **low** and **transmission** to neighbors is **moderate**. ```{r} pro_vaccine <- list(p_vaccinated = 0.8, p_removed = 0.3, p_infected = 0.05, p_spread = 0.5) ``` ## Scenario II (Vaccine Rejection) In this scenario, the **vaccination** is **widely rejected**, the **disease** has a **moderate mortality rate**, the chance of **initial infection** is **low** and **transmission** to neighbors is **moderate**. ```{r} vaccine_rejection <- list(p_vaccinated = 0.02, p_removed = 0.3, p_infected = 0.05, p_spread = 0.5) ``` ## Scenario III (Contagious Disease) In this scenario, the **vaccination** is **widely rejected**, the **disease** has a **moderate mortality rate**, the chance of **initial infection** is **low** and **transmission** to neighbors is **high**. ```{r} contagious <- list(p_vaccinated = 0.2, p_removed = 0.3, p_infected = 0.05, p_spread = 0.9) ``` ## Scenario IV (Epidemiological Collapse) In this scenario, the **vaccination** is **widely rejected**, the **disease** has a **high mortality rate**, the chance of **initial infection** is **low** and **transmission** to neighbors is **high**. ```{r} collapse <- list(p_vaccinated = 0.1, p_removed = 0.9, p_infected = 0.05, p_spread = 0.95) ``` ## All the scenarios ```{r} scenarios <- list(pro_vaccine, vaccine_rejection, contagious, collapse) ``` ### Function to plot the city grid ```{r} #| Warning: false #| Message: false plot_city <- function(city, step) { df <- expand.grid(x = 1:grid_size, y = 1:grid_size) df$state <- as.vector(city) df$state <- factor(df$state, levels = c(0, 1, 2, 3), labels = c("Vaccinated", "Susceptible", "Infected", "Removed")) # Ensuring all states are included in the factor levels ggplot(df, aes(x = x, y = y, fill = state)) + geom_tile() + scale_fill_manual( values = c("white", "#DAF7A6", "#DE3163", "black"), # Defining colors for states drop = FALSE # Preventing ggplot from dropping missing states ) + labs( title = paste("Scenario 1 (Pro-vaccine population) - Step", step), fill = "State" ) + theme_minimal() + theme( axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "bottom", legend.text = element_text(size = 12, face = "italic"), legend.title = element_text(size = 14, face = "bold"), plot.title = element_text(size = 16, face = "bold") ) + scale_x_continuous(expand = c(0, 0)) + # No space around grid scale_y_continuous(expand = c(0, 0)) # No space around grid } ``` ### Function to spread infection to neighbors on an absorbent edge ```{r} #| Warning: false spread_to_neighbors <- function(i, j, grid_size, city, p_spread, n_city) { neighbors <- list( c(i-1, j), c(i+1, j), # up and down c(i, j-1), c(i, j+1), # left and right c(i-1, j-1), c(i-1, j+1), # upper diagonals c(i+1, j-1), c(i+1, j+1) # lower diagonals ) for (n in neighbors) { ni <- n[1] nj <- n[2] if (ni >= 1 && ni <= grid_size && nj >= 1 && nj <= grid_size) { # check if the neighbor is inside the grid (absorbing edge) if (city[ni, nj] == 1) { # only susceptible may be infected if (runif(1) < p_spread) { n_city[ni, nj] <- 2 # neighbor become infected } } } } return(n_city) } ``` ### Function to save data ```{r} #| Warning: false data_saved <- function(city, step, grid_size) { df <- expand.grid(x = 1:grid_size, y = 1:grid_size) df$state <- as.vector(city) df$step <- step return(df) } ``` ### Initializating the grid with all cells susceptible, with the central cell infected ```{r} #| Warning: false city <- matrix(1, nrow = grid_size, ncol = grid_size) city[grid_size/2,grid_size/2]=2 ``` ### Viewing the initial condition ```{r} #| Warning: false #| Message: false print(plot_city(city, 0)) ``` ### Choosing which scenario we will simulate Note that our example will be scenario I, pro-vaccine. ```{r} #| Warning: false chosen_scenario <- scenarios[[1]] # The chosen scenario print(chosen_scenario) saved_data_list <- list() ``` # Simulation loop ```{r} #| Warning: false for (step in 1:n_steps) { # Do until the final "time" n_city <- city # Create a copy of the city to update states for (i in 1:grid_size) { # Do for all lines for (j in 1:grid_size) { # Do for all columns if (city[i, j] == 0) { # If the cell be equal to 0 (vaccinated) n_city[i, j] <- 0 # Stay 0 (vaccinated) } else if (city[i, j] == 3) { # If the cell be equal to 3 (dead) n_city[i, j] <- 3 # Stay 3 (dead) } else if (city[i, j] == 1) { # If the cell be equal to 1 (susceptible) if (runif(1) < chosen_scenario$p_infected) { # And the random number is within the probability of being infected n_city[i, j] <- 2 # The state changes to 2 (infected) } else if (runif(1) < chosen_scenario$p_vaccinated) { # If the random number is within the probability of being vaccinated n_city[i, j] <- 0 # The state changes to 0 (vaccinated) } } else if (city[i, j] == 2) { # If the cell be equal to 2 (infected) if (runif(1) < chosen_scenario$p_removed) { # If the random number is within the probability of being removed n_city[i, j] <- 3 # The state changes to 3 (removed) } else { n_city <- spread_to_neighbors(i, j, grid_size, city, chosen_scenario$p_spread, n_city) # spread the infection to susceptible neighbors } } } } city <- n_city # Updating the matrix saved_data_list[[step]] <- data_saved(city, step, grid_size) print(plot_city(city, step)) # Plot the current state of the city #Sys.sleep(0.5) # Pause for a short time to visualize the simulation } ``` ### Viewing the results ```{r, warning=FALSE, message=FALSE} #| Code-fold: true results_scenario_one <- do.call(rbind, saved_data_list) summary_one <- results_scenario_one %>% group_by(step) %>% summarise( vaccinated = sum(state == 0), susceptible = sum(state == 1), infected = sum(state == 2), removed = sum(state == 3) ) summary_one <- summary_one %>% bind_rows(data.frame(step = 0, vaccinated = 0, susceptible = 9999, infected = 1, removed = 0)) ggplot(summary_one) + geom_line(aes(x = step, y = vaccinated, color = "Vaccinated"), size = 1) + geom_line(aes(x = step, y = susceptible, color = "Susceptible"), size = 1) + geom_line(aes(x = step, y = infected, color = "Infected"), size = 1) + geom_line(aes(x = step, y = removed, color = "Removed"), size = 1) + labs(title = "Scenario 1 - (Pro-vaccine Population)", x = "Time", y = "Number of sites", color = "State") + scale_color_manual(values = c("Vaccinated" = "#e5e8e8", "Susceptible" = "#DAF7A6", "Infected" = "#DE3163", "Removed" = "black")) + theme_minimal()+ theme( panel.grid = element_blank(), axis.title = element_text(size = 12, face = "bold"), axis.text = element_text(size = 10), axis.ticks = element_line(size = 0.5), legend.position = "bottom", legend.text = element_text(size = 12, face = "italic"), legend.title = element_text(size = 14, face = "bold"), plot.title = element_text(size = 16, face = "bold") ) ``` In this scenario, we observed that the high rate of early vaccination reduced the number of susceptible sites and consequently the number of infected and dead sites. ### This is the end of the scenario I simulation, but you can try the other scenarios. ![](_110911141_gettyimages-151032915.jpg.webp) See you!