Modeling the effect of outside travellers on epidemic spread in an agent-based simulation

This model is based on the model/simulation found found here.

The agent

We will augment the agent with an additional state–not in network. This will let us have ‘vacationer’ or ‘traveller’ agents that leave the network and no longer impact anyone.

  1. Unexposed
  2. Asymptomatic but infected/contagious
  3. Symptomatic and contagious
  4. Symptomatic and not contagious
  5. Post-COVID Immune
  6. Naturally immune (will not contract)
  7. Death
  8. Left network
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Loading required package: statnet.common
## 
## Attaching package: 'statnet.common'
## The following object is masked from 'package:base':
## 
##     order
## Loading required package: network
## network: Classes for Relational Data
## Version 1.16.0 created on 2019-11-30.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##                     Mark S. Handcock, University of California -- Los Angeles
##                     David R. Hunter, Penn State University
##                     Martina Morris, University of Washington
##                     Skye Bender-deMoll, University of Washington
##  For citation information, type citation("network").
##  Type help("network-package") to get started.
## sna: Tools for Social Network Analysis
## Version 2.5 created on 2019-12-09.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:sna':
## 
##     betweenness, bonpow, closeness, components, degree, dyad.census,
##     evcent, hierarchy, is.connected, neighborhood, triad.census
## The following objects are masked from 'package:network':
## 
##     %c%, %s%, add.edges, add.vertices, delete.edges, delete.vertices,
##     get.edge.attribute, get.edges, get.vertex.attribute, is.bipartite,
##     is.directed, list.edge.attributes, list.vertex.attributes,
##     set.edge.attribute, set.vertex.attribute
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
## 
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
## 
##     nasa

Create transition matrix.

makeAgent <- function(psychstate,biostate,age=30,traveler=0,is_infected = 0)
{
  
  return (list(psychstate=psychstate,
               biostate=biostate,
               age=age,
               nextbiostate=NA,
               biostatecountdown=NA,traveler=traveler,is_infected = is_infected))
}

# * 1. Unexposed
# * 2. Asymptomatic but infected/contagious
# * 3. Symptomatic and contagious
# * 4. Symptomatic and not contagious
# * 5. Post-COVID Immune
# * 6. Naturally immune (will not contract)
# * 7. Death
# * 8. Left network


bioTransition <- matrix(0,STATES,STATES)
bioMin <- matrix(1,STATES)      #state time minimum
bioMax <- matrix(1,STATES)      #state time maximum



bioMin[2] <- 3            #infected but asymptomatic for 3 to 10 days
bioMax[2] <- 10         
bioTransition[2,3] <- .1  #transition to infected with symptoms
bioTransition[2,5] <- .90  #transition to no longer contagious/cured


bioMin[3] <-    3             #symptoms + contagion
bioMax[3] <- 8                #symptoms + contagion max
bioTransition[3,4] <-  .95    #transitioon to no longer contagious
bioTransition[3,7] <-  .05    #transitioon to death state 


bioMin[4] <- 1          #symptoms bot no longer contagiious
bioMax[4] <- 7
bioTransition[4,5] <- 1  #Transition to 'immune' cured state.




setAgentState<- function(agent, biostate)
{
  agent$biostate <- biostate
  if(sum(bioTransition[biostate,])>0) # this state transitions to something else.
  {
    ##which state do we go to?
     agent$biostatecountdown <- sample(x=seq(bioMin[biostate],bioMax[biostate]),1) #how long will we state in this state?
     agent$nextbiostate <- sample(1:STATES, prob=bioTransition[agent$biostate,],size=1)
     
  } else{
   agent$biostatecountdown <- NA
   agent$nextbiostate <- NA   ##just so we can tell if the agent is finished.
  }
  return(agent) 
}

transitionAgent<- function(agent)
{
   return(setAgentState(agent,agent$nextbiostate))
}

updateAgent<- function(agent)
{
  if(!is.na(agent$biostatecountdown))
  {
   agent$biostatecountdown <- agent$biostatecountdown -1
    if(agent$biostatecountdown <=0)  ##new state
    {
       agent <- transitionAgent(agent)

    }
  }
   return(agent)
}

setAgentStatusInfected<- function(agent,status=1)
{
   agent$is_infected <- status
   return(agent)
}

Social model

This is not quite right, because we’d really expect interactions to be clumpy–you interact with the same people every day, but this model is like the opinion dynamics model in that anyone can interact with anyone. Let’s make a network using preferential attachment to represent that social world. We can make a couple of them and add them together to represent how we each have different types of relationships and interactions–you may be the president of your drama club and so are a central member, but you are also in a pottery class where you don’t talk to very many people (although your teacher does).

Adding a new agent to the network

Here, we write code that will grow the network by adding new agents.

Background

A key step in contingency planning for a possible pandemic in a community is barring entry of travelers in the community. The potential for infectious diseases to spread in a community rapidly through well-connected infected travelers in the community is very high. Past researches looked into the effect of the international travel ban on reducing the spread of disease in countries. Partial restriction of travel helps in slowing the spread if the global number of cases of infection is low (Hollingsworth et al., 2006). According to Epstein et al. (2007), international air travel restrictions may provide a small but important delay in the spread of a pandemic, especially if other disease control measures are implemented during a timeframe. Using the mathematical model, Cooper et al. (2006) concluded that restrictions on air travel would achieve very little. This is probably because a virus may transmit from one person to another very quickly and affects many people. Once a major outbreak is underway, banning flights from affected cities would be effective at significantly delaying spread only if almost all travel between cities could be stopped almost as soon as an outbreak was detected in each city.

One of the main reasons against oversea travel ban is because it has a big negative economic impact in a country if we stop adults traveling oversea (Lam et al., 2011). But for small and moderate length communities whose economic activities circulated within the community, a travel restriction may be helpful during a pandemic. If the health care systems for these communities are not prepared to deal with a pandemic, it is better to enforce some sort of travel ban so that the health care system does not overwhelm with patients. This will help to keep a check on the number of patients. In the traditional SIR model (Smith & Moore, 2004), it is not defined how adding travelers to the scenario will affect the spread of disease. So, we created an agent-based social model for a community where we will be adding travelers to the local populations to see how a disease spreads in the whole community. We will be looking into the worst-case scenario in the model, where travelers will enter the community but nobody leaves within a timeframe from the community. We will vary different variables in the simulations to see the disease progress in the community within a specific timeframe.

References

  • Cooper, B. S., Pitman, R. J., Edmunds, W. J., & Gay, N. J. (2006). Delaying the international spread of pandemic influenza. PLoS Medicine, 3(6).
  • Epstein, J. M., Goedecke, D. M., Yu, F., Morris, R. J., Wagener, D. K., & Bobashev, G. V. (2007). Controlling pandemic flu: The value of international air travel restrictions. PloS One, 2(5).
  • Hollingsworth, T. D., Ferguson, N. M., & Anderson, R. M. (2006). Will travel restrictions control the international spread of pandemic influenza? Nature Medicine, 12(5), 497–499.
  • Lam, E. H., Cowling, B. J., Cook, A. R., Wong, J. Y., Lau, M. S., & Nishiura, H. (2011). The feasibility of age-specific travel restrictions during influenza pandemics. Theoretical Biology and Medical Modelling, 8(1), 44.
  • Smith, D., & Moore, L. (2004). The SIR Model for Spread of Disease: The Differential Equation Model. Loci.(Originally Convergence.) Https://Www. Maa. Org/Press/Periodicals/Loci/Joma/the-Sir-Model-for-Spread-of-Disease-the-Differential-Equation-Model.

Travel ban lift and allowing traveler to the community

Social connection of the travelers

travel_ban_lift <- function( enable_outside_travel_restriction_day,second_wave_day,
                             connection_number = 1, will_affect_traveler = 0, numDays = 50,
                              people_leave_status = 1,
                             number_of_people_leave = 4,
                             naturalImmunity = .01,numInteractions_value = 14,
                              contagionProb = .05,####5% chance
                             sampleFromNetwork = .98,traveler_percentage_mean = 1,
                             traveler_percentage_ceiling = 2, traveler_increase = 0,
                             boundary_of_traveler_change = 2, traveler_interval_day = 5, 
                             enable_second_wave = 1,
                             enable_outside_travel_restriction =0,need_to_update_graph_cord = 0 ){ 



 
  socialnetwork <- original_network_of_agents
  numAgents <- numAgents_initial
  people_left <-0
  

  
  plotNetwork <-T
  
  # If i create the social network here, works ok
  #socialnetwork <-makeNetwork(numAgents,numsets=1,power=.5,steps=2)
  #print(socialnetwork)
  
  #pool <<-pool_initial
  
  .GlobalEnv$pool <<- pool_initial
  
  numInteractions <-  rep(numInteractions_value,numDays) ##interaction during no travel bans
  numInteractions[enable_outside_travel_restriction_day:second_wave_day] <- as.integer(numInteractions_value*(1/3))  ##interaction during travel restriction
  
  
  local_population_infection_number_local_variable <<- local_population_infection_number
  traveler_infection_number_local_variable <- 0

  if(plotNetwork)
  {
    cc <-mygplot_updated(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
    
    xrange <-   range(cc[,1])
    xval <- xrange[2] + (xrange[2]-xrange[1])/3
    yval <- mean(cc[,2])
    sd <- (xrange[2]-xrange[1]) /10
  }  

  
  

  disthistory <- matrix(NA,ncol=STATES,nrow=numDays)
  
 # pool <<-pool_initial

  
  for(day in 1:numDays)
  {
      
      
       ###################agent leaving Start
    
      #innetwork <- sapply(.GlobalEnv$pool,function(y){!(y$biostate==7||y$biostate==8)})
      if(!enable_outside_travel_restriction){
        if(people_leave_status){
          ####Hard coding day for people leave
           #.GlobalEnv$pool <- list.remove(.GlobalEnv$pool,  c(1))
          
          if(day == (second_wave_day+7)){
                  
                  number_of_people_leave <- as.integer((1/100)*numAgents) ##1% will leave the system
                  
                      
                    u <- 1:numAgents
                
                    for (val in u) {
                      
                      if((.GlobalEnv$pool[[val]]$biostate != 8) ){
                
                        if((.GlobalEnv$pool[[val]]$biostate != 7)){
                          x<-u
                          
                        }
                      }
                     
                    }    
                    
                    a <- sample(x, number_of_people_leave)
                    
                            
                    for (i in seq_along(a)) {
                      
                      #print(a[i])
                      x <- a[i]
                      .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          
                      socialnetwork[x,] <- 0
                      socialnetwork[,x] <-0
                    }                
                           
                  
                  
                 
          }          
          if(day == (second_wave_day+14)){
                  
                  number_of_people_leave <- as.integer((1/100)*numAgents) ##1% will leave the system
                  
                      
                    u <- 1:numAgents
                
                    for (val in u) {
                      
                      if((.GlobalEnv$pool[[val]]$biostate != 8) ){
                
                        if((.GlobalEnv$pool[[val]]$biostate != 7)){
                          x<-u
                          
                        }
                      }
                     
                    }    
                    
                    a <- sample(x, number_of_people_leave)
                    
                            
                    for (i in seq_along(a)) {
                      
                      #print(a[i])
                      x <- a[i]
                      .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          
                      socialnetwork[x,] <- 0
                      socialnetwork[,x] <-0
                    }                
                           
                  
                  
                 
          }              

          # if(day == (second_wave_day+5)){
          #   
          #   d <- 1:numAgents
          #   number_traveler_found
          #   repeat{
          #     statements...
          #     if(condition){
          #       break
          #     }
          #   }            
          #   
          #   
          #   
          #       x <- 18
          #       .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          # 
          #       socialnetwork[x,] <- 0
          #       socialnetwork[,x]<-0
          # }

        }

      }

      ###################agent leaving ENd
    

          innetwork <- sapply(.GlobalEnv$pool,function(y){!(y$biostate==7||y$biostate==8)})
  
  ##who are you going to talk to today.
      sneezers <- rep((1:numAgents)[innetwork],each=numInteractions[day])
      sneezedons <- rep(NA,length(sneezers))
      
      for(i in 1:length(sneezers))
        {
        if(runif(1)<(1-sampleFromNetwork))
        {
    
          sneezedons[i] <- sample((1:numAgents)[innetwork],1)
        }else{
          
          probs <- socialnetwork[sneezers[i],]* innetwork
          sneezedons[i] <- sample(1:numAgents,prob=probs,1)
          
        }
       }
      
        
    for(i in 1:length(sneezers))
    {
      
      agent1 <- .GlobalEnv$pool[[ sneezers[i] ]]
      agent2 <- .GlobalEnv$pool[[ sneezedons[i] ]]
    
      
      ##this constitutes the rules of infection.
      if((agent1$biostate==2 || agent1$biostate==3 ) & agent2$biostate==1 & runif(1)<contagionProb)
      {
           .GlobalEnv$pool[[ sneezedons[i] ]] <- setAgentState(agent2,2)##infect!
      }
        
    }
      ##increment each agent 1-day.
      for(i in 1:numAgents)
      {
         .GlobalEnv$pool[[i]] <- updateAgent(.GlobalEnv$pool[[i]])
         
         ###keep track number of local population infection
         
           
           if(.GlobalEnv$pool[[i]]["traveler"] == 0 ){
             if(.GlobalEnv$pool[[i]]["biostate"] == 2){
               if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
                 local_population_infection_number_local_variable <-local_population_infection_number_local_variable+1
                .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
               }
               
             }
            
           
         }         
           if(.GlobalEnv$pool[[i]]["traveler"] == 1 ){
             if(.GlobalEnv$pool[[i]]["biostate"] == 2){
               if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
                 traveler_infection_number_local_variable <-traveler_infection_number_local_variable+1
                .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
               }
               
             }
            
           
         }            
         
         
      }
    
      states <- sapply(.GlobalEnv$pool,FUN=function(x){x$biostate})
      distrib <- table(factor(states,levels=1:STATES))
      disthistory[day,] <- distrib
     # print(disthistory[day,])
      
      
     
  
      
       #####################################
      ############################
      ##################################
      ################
      ################################   
      ####################Adding New agents as travelers in the network
      
      #traveler_interval_day <-sample(c(3,2,5),1)
      
     # traveler_increase <- sample(c(0,1),1)
     
      
    
      if(traveler_increase == 0){
        if(day == second_wave_day){
          
          if(enable_second_wave == 1){
            traveler_percentage_mean <- as.integer(traveler_percentage_ceiling/1.5)
            
             
          }
          
          
        }
      }
      ##travel restriction on
  
      
      
      ###if second wave already happned due to premature removal of travel ban. traveler will come in future days
      if((enable_second_wave == 1) && (day == second_wave_day)){
        
        enable_outside_travel_restriction <-0
        
      }
      else if((enable_second_wave == 1) && (day > second_wave_day)){
        
        enable_outside_travel_restriction <-0 ###if we do not want to close down the travel after second wave
      }
      else if(enable_outside_travel_restriction_day == 0){
        
         enable_outside_travel_restriction <-0
         
      }else if(enable_outside_travel_restriction_day == day){
         enable_outside_travel_restriction <-1
        
      }else if(enable_outside_travel_restriction == 1){
        
        enable_outside_travel_restriction <-1
        
      }
      else{
        
        enable_outside_travel_restriction <-0
        
      }
      
      

      
  
      ####
      if(((!(day %% traveler_interval_day)) ||
           (day==second_wave_day)) && 
           (enable_outside_travel_restriction==0)){
        
            if(traveler_percentage_mean < 1){

              traveler_percentage_mean =  as.integer(traveler_percentage_ceiling/2)
            }
            else if(traveler_percentage_mean>traveler_percentage_ceiling){
              
              traveler_percentage_mean = traveler_percentage_ceiling
            }
     
        
          if(traveler_percentage_mean > 0){
            number_of_traveler <-as.integer((traveler_percentage_mean/100) * numAgents_initial)
            current_number_traveler_remaining <- number_of_traveler
            if(traveler_increase < 1){

              traveler_percentage_mean = traveler_percentage_mean - boundary_of_traveler_change

            }else{

              traveler_percentage_mean = traveler_percentage_mean + boundary_of_traveler_change

            

            }
            ###enabling travel restruition after second wave
            if(enable_second_wave == 1){
              if((day==second_wave_day) && (enable_outside_travel_restriction_day != 0)){
                
                enable_outside_travel_restriction<-1
              }
            }
            
            ###for evenly distributing travel to network I run a do while loop
            #print("number_of_traveler")
            #print(number_of_traveler)
            #print("-----------------------start-------------------")
            repeat{
              current_number_traveler <- sample(c(1,2,3,4,5),1)
                
                #
  
              
              if(current_number_traveler>current_number_traveler_remaining){
                
                current_number_traveler <- current_number_traveler_remaining
              }
              
             # print(current_number_traveler)
  
              
              if(current_number_traveler > 0)
              {
                starting_numagents <- numAgents + 1
                ending_numagent <- numAgents+current_number_traveler
                
                ####connection_number indicates new agent connected how many number of previous agents, currently it is a number between 1-3
                
                socialnetwork <- addAgentToNetwork(numAgents = current_number_traveler,
                                                   prev_network = socialnetwork,
                                                   connection_number = connection_number,
                                                   self_connected=1,.GlobalEnv$pool)
               # print(socialnetwork)
                
                newxy <- cbind(rnorm(current_number_traveler,mean=xval,sd=sd),
                               rnorm(current_number_traveler, mean=yval,sd=sd))
                
                cc <- rbind(cc,newxy)
                for(i in starting_numagents:ending_numagent)
                  {
                    ###printing to see how are travelers
                    # print("adding")
                     .GlobalEnv$pool[[i]] <<- makeAgent(psychstate=1,
                                            biostate=sample(c(1,6),
                                            p=c(1-naturalImmunity, naturalImmunity),1),
                                            traveler = 1)
                   # print("added")
                    ################to make randomize infect -> traveler
                     
                     will_be_infected <- sample(c(0,1), 1)
                     
                    
                     
                    # print(paste("infecting",i))
                     #####Tell shane about if setAgentState is 2, then number of local remains ok, but if state is 1 then it changes
                     if(will_affect_traveler == 1){
                       #print("poit a")
                       tmp <-  setAgentState(.GlobalEnv$pool[[i]],2) ##infect this person
                       #print("point b")
                       .GlobalEnv$pool[[i]] <<-tmp
                       
                     }else{
          
                      # print("point b1")
                      # tmp <- setAgentState(.GlobalEnv$pool[[i]],1) ##do not infect this person
                      # print("point b2")
                      # .GlobalEnv$pool[[i]] <<- tmp
                      # print("point b3")
                     }
                    # print("infected")
                     
                     # 
                     
                     
                 
                     
                     
                }
                numAgents <- ending_numagent
                need_to_update_graph_cord <- 1
              }
              current_number_traveler_remaining <- current_number_traveler_remaining-current_number_traveler
              if(current_number_traveler_remaining < 1){
                 #print("-----------------------end-------------------")
                break
              }
              
            }
        }
      }
      
     if(plotNetwork)
     {
        ## print(1)
       
       ##show this to SHane
        # cc <-mygplot(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
        # mygplot(coord=NULL,socialnetwork,states,main=paste("Day",day))
       
       if(need_to_update_graph_cord){
        # print("UPDATE NEEDED")
        
         mygplot_updated(cc,socialnetwork,states,main=paste("Day",day, " ", "Traveler Added", number_of_traveler))
        
         
       }else{
         
         mygplot_updated(cc,socialnetwork,states,main=paste("Day",day))
       }
        need_to_update_graph_cord <- 0
       
        
         
         
         
  
     }    
      
      
  ###vary week for travel restriction
     
    
      
         
  }

  total_local <-0
  affected_locals <- 0
  affected_travelers <- 0
  local_death <- 0
  total_traveler <-0
  traveler_death <- 0
  ok_travler <- 0
  ok_local<- 0
  #print(disthistory)
     
   for(i in 1:numAgents)
   {
     
     if(pool[[i]]["traveler"] == 0){
       
       #print(i)
       
       total_local <-total_local+1
       if(pool[[i]]["biostate"] == 2 || pool[[i]]["biostate"] == 3 || pool[[i]]["biostate"] == 4)
         {
           affected_locals <-affected_locals + 1
       }
       else if(pool[[i]]["biostate"] == 7){
         
         local_death <-local_death+1
       }
       else{
         
         ok_local <- ok_local+1
       }
       
     }
     else{
       total_traveler <- total_traveler+1
       if(pool[[i]]["biostate"] == 2 || pool[[i]]["biostate"] == 3 || pool[[i]]["biostate"] == 4)
         {
           affected_travelers <-affected_travelers + 1
       }
       else if(pool[[i]]["biostate"] == 7){
         
         traveler_death <-traveler_death+1
       }  
       else{
         
         ok_travler <- ok_travler+1
       }       
       
       
     }
   }
  
    #.GlobalEnv$pool <- list.remove(.GlobalEnv$pool, 1)
    #print(.GlobalEnv$pool)
  print("TOTAL")
  print(numAgents)
  print("Total Local")
  print(total_local)
  # print("Local Recovered/Not affected")
  # print(ok_local)
  # print("Local Affected")
  # print(affected_locals)
  print("Local Death")
  print(local_death)
  print("-----")
  print("Total Traveler")
  print(total_traveler)  
  
  # print("Traveler Recovered/Not affected")
  # print(ok_travler)
  # print("Traveler Affected")
  # print(affected_travelers)
  print("Traveler Death")
  print(traveler_death)
  
  
  print("Local infected at somepoint")
  print(local_population_infection_number_local_variable)
  print("Traveler infected at somepoint")
  print(traveler_infection_number_local_variable) 
  print("Total infected at somepoint")
  print(traveler_infection_number_local_variable + local_population_infection_number_local_variable) 
  
 # print(traveler_percentage_mean)
  

  


  
  return(disthistory)
  
  


 
}

We will run the model for a community of 700 local agents, for 100 days

The simulations will be divided into major parts. Interactions among agents will be one-third of normal interactions during travel ban

In the first part, we will see how travlers affect the local populations when no from the local populaton if infeceted.

Travel ban

In this simulation, a ‘travel ban’ is implemented from Day 7 and removed early (In this case from day 49). Each traveler will have only three connections with local agents. For this simulation, we will look at the spread of the disease through the network over time.

## [1] "TOTAL"
## [1] 791
## [1] "Total Local"
## [1] 700
## [1] "Local Death"
## [1] 4
## [1] "-----"
## [1] "Total Traveler"
## [1] 91
## [1] "Traveler Death"
## [1] 0
## [1] "Local infected at somepoint"
## [1] 536
## [1] "Traveler infected at somepoint"
## [1] 84
## [1] "Total infected at somepoint"
## [1] 620