COVID-19 epidemics with Non-Pharmaceutical Interventions and zonal restraints

;; --------------------------------------------------------------------------
;; COVID-19 epidemic simulation with Non-Pharmaceutical Interventions,
;; including a novel zonal restraint NPI.
;;
;; A model by Sergio Rojas-Galeano and Lindsay Alvarez
;; v1.24 Copyright (c) June 2020 The authors
;; Correspondance email: srojas@udistrital.edu.co
;; Universidad Distrital Francisco Jose de Caldas, Bogota, Colombia
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License (GPLv3)
;; (see license at: https://www.gnu.org/licenses/gpl-3.0.txt)
;;
;; The model is made publicly available in the hope that it will be useful
;; to modelers, but WITHOUT ANY WARRANTY whatsoever (see license for details).
;; --------------------------------------------------------------------------

;; Define globals
globals [
  daytime?      ; is it day or night? ( 1 day + 1 night = ticks-day)
  sort-codes    ; a list of codes to identify each zone of a city
  shapes-list   ; a list of shapes to identify agents of the same zone
  colors-list   ; a list of colors to identify different zones
  day hour      ; the current day and hour since simulation started
  sicks-day     ; list of daily count of sicks
  n-cases       ; number of overall COVID-19 sick people
  n-confirmed   ; number of confirmed COVID-19 sick people
  n-tested      ; number of overall people tested by sentinel
  n-positives   ; number of people that tested COVID-19 positive
  doubling-time ; estimated days to double today's number of sicks
  R0 S0 R0_     ; reproduction number and initial Susceptibles
  p-contagion   ; probability of contagion on a single encounter
]

;; Define breeds for healthy and infected people
breed[ healthies healthy ]
breed[ sicks sick ]
breed[ deaths death ]
breed[ houses house ]
breed[ ambulances ambulance ]

;; Define attributes for patches
patches-own [
  zone                   ; the zone code of this patch
]

;; Define attributes of all the people
turtles-own [
  homebase               ; the home patch of this person
  scope                  ; the reach of distance away from home for this person
  speed                  ; how fast does he/she move? A speed of 0 is still
  unlocked?              ; is the person free from lockdown?
  tested?                ; was the person tested by sentinel test?
  contagions             ; how many people this person infected if he/she ever got sick
]

;; Define attributes for sick people only
sicks-own [
  days-to-cure           ; the actual duration of the illness for this person
  recovery               ; counter of days remaining to recovery
  confirmed?             ; is the person isolated at homebase?
  asymptomatic?          ; is the person asymptomatic?
  risky?                 ; is the person high-risk due to existing factor (obese, diabetes, older, etc)
  severe?                ; is the person in severe condition?
  deadly?                ; is the person in deadly (critical) condition?
  hospitalised?          ; has he/she been admitted to hospital?
  ICU-admitted?          ; has he/she been admitted to ICU?
]

;; Define attributes for healthy people only
healthies-own [
  immune?                ; has the person recovered and acquired immunity?
]

;; Define attributes for ambulances
ambulances-own [
  tests-stock            ; how many tests in stock
  tests-supply           ; numbers of tests to top-up (depends on amount of zone residents)
]

;; One iteration of simulation's steps

to go
  ifelse day < end-day [ tick ] [ stop ] ;screenshot stop ]
  lifestyle
  epidemic
  isolation
  quarantine
  distancing
  sentinel
  lockdown
  illness
  clock
  indicators
end 

;; Simulate everyday routine of people

to lifestyle
    ask (turtle-set healthies sicks with [not confirmed?] ) [
    if zone-enforcing? [
      if zone != [zone] of homebase [ face homebase fd 0.5 ]  ; if away from resident zone, head back
    ]
    forward speed                                             ; move ahead
    if distance homebase > scope [ face homebase ]            ; if too far from home, head back
    set heading heading + (random-float 3 - random-float 3)   ; change direction swiftly
  ]
end 

;; Spread the infection from sicks not yet isolated, to healthies

to epidemic
  ask sicks with [not confirmed?] [
    let counter 0
    ask healthies with [not immune?] in-radius .5 [
      ;; Disease spreads depending on contagion probability
      if random-float 1 < p-contagion [
        get-sick
        set counter counter + 1
      ]
    ]
    set contagions contagions + counter
  ]
end 

;; Make a single person ill

to get-sick
  ;; S->I (set the initial properties of a sick person)
  set breed sicks
  set shape item ([zone] of homebase) shapes-list
  set confirmed? false
  set tested? false
  set severe? false
  set deadly? false
  set hospitalised? false
  set ICU-admitted? false

  ;; Set stochastic attributes
  set days-to-cure abs int random-normal avg-duration 4         ; duration of illness for this sick
  set recovery 0                                                ; initial count of days until recovery
  set asymptomatic? (random 100) < %-asymptomatics              ; is this an asymptomatic patient?
  set risky? (random 100) < %-high-risk  and not asymptomatic?  ; is it high-risk (co-morbidities) given not asymptomatic?
  set color ifelse-value asymptomatic? [ yellow ] [ red ]

  ;; Update how many cases so far
  set n-cases n-cases + 1
end 

;; Ask confirmed sick people to stay home

to isolation
  ;; If quarantine on, enforce isolation of symptomatic
  ifelse case-isolation? [
    ask sicks with [ not asymptomatic? and not confirmed? ] [
      move-to homebase
      set confirmed? true
      set n-confirmed n-confirmed + 1         ; these sicks are now confirmed cases
    ]
  ]
  ;; Else lift isolation, except if was sent sentinel testing
  [
    ask sicks with [ not tested? and confirmed? ] [
      set confirmed? false
      set n-confirmed n-confirmed - 1         ; these sicks are set back to not confirmed cases
    ]
  ]

  ;; Decorate houses with color of isolated people
  ask houses [
    set color ifelse-value any? sicks-here with [ not asymptomatic? ] [ red ]
            [ ifelse-value any? sicks-here with [ asymptomatic? ] [ yellow ] [ white ] ]
  ]
end 

;; Ask housemates of isolated  people to stay home

to quarantine
  if home-quarantine? [
    ask sicks with [ confirmed? ] [
      let sick-home homebase
      ask other turtles with [ homebase = sick-home ] [ move-to homebase ]
    ]
  ]

  ;; Decorate houses with color of quarantined people
  ask houses [
    set color ifelse-value any? sicks-here with [ not asymptomatic? ] [ red ]
            [ ifelse-value any? sicks-here with [ asymptomatic? ] [ yellow ] [ white ] ]
  ]
end 

;; Avoid contacting nearby people (closer than than min-dist)

to distancing
  if social-distancing?
  [ let people (turtle-set healthies sicks)
    let willings (count people * %-willings / 100)
    ask n-of willings people [
      if any? other people in-radius 0.6 [  ; a radius of 0.6 is slightly higher than the 0.5 threshold to contagion
        let closest min-one-of other people [distance myself]
        face closest rt 180                 ; head in the opposite direction of closest person
      ]
    ]
  ]
end 

;; Ask ambulances to catch sicks by mass testing

to sentinel
  ifelse sentinel-testing? [
    ask ambulances [
      st
      if tests-stock > 0 [
        ;; Move ambulances around
        forward speed
        if zonal? [ if zone != [zone] of homebase [ face homebase fd 0.5 ] ]  ; keep ambulance in-zone
        set heading heading + (random-float 10 - random-float 10)             ; change ambulance's heading randomly

        ;; Pick nearbys and perform tests
        let target-people ((turtle-set healthies sicks) in-radius 5) with [not tested?]
        ask target-people  [
          if is-sick? self and not confirmed? [                               ; if sick and not isolated, target is positive
            move-to homebase set confirmed? true
            set n-positives n-positives + 1                                   ; update count positives tests so far
            set n-confirmed n-confirmed + 1                                   ; update count of confirmed cases so far
          ]
          set tested? true                                                    ; mark this person as tested
        ]
        set n-tested n-tested + (count target-people)                         ; update count of tests made so far
        set tests-stock tests-stock - (count target-people)                   ; update stock of tests available
      ]
    ]
  ]
  ;; Else, stop and hide ambulances
  [ ask ambulances [ ht ] ]
end 

;; Enforce everyone staying home, except for a few permits

to lockdown
  if hour mod 24 = 0 [                                              ; lockdown starts or ends at 00:00 of next day
    ifelse total-lockdown?[
      let people (turtle-set healthies sicks) with [not unlocked?]  ; for people not in an unlocked zone
      ask people [ move-to homebase set speed 0 ]                   ; sent everyone home, not confirmed
      foreach sort-codes [
        z -> let residents people with [ zone = z ]
        let permits (count residents * %-permits / 100)
        ask n-of permits residents [ set speed step-size ]
      ]
    ]
    ;; Else unlock everyone
    [ ask (turtle-set healthies sicks) [ set speed step-size set unlocked? false ] ]
  ]
end 

;; Lift total lockdown by unlocking zones gradually

to unlock-zone
  let unlocked-zone random zones
  ask (turtle-set healthies sicks) with [zone = unlocked-zone] [ set speed step-size set unlocked? true ]
end 

;; Simulate illness evolution in each person

to illness
  if ticks mod ticks-day = 0 [                        ; illness unfolds in a daily basis

    ;; I->R (sicks recovering after days-to cure)
    ask sicks [
      set recovery recovery + 1                       ; countdown of days left to recover
      if recovery >= days-to-cure [ get-well ]        ; if recovery time over, the person is cured
    ]

    ;; I->E (sicks getting extinct, that is, dead)
    ask sicks with [not asymptomatic?] [
      ;; Define probability of death
      let p 0.005                                     ; P(Extinct ∣ Symptomatic)
      if deadly? and not ICU-admitted? [ set p 1.0  ] ; P(Extinct ∣ Deadly, no ICU)
      if deadly? and     ICU-admitted? [ set p 0.5  ] ; P(Extinct ∣ Deadly, ICU)
      if severe? and not hospitalised? [ set p 0.2  ] ; P(Extinct ∣ Severe, Hospitalised)
      if severe? and     hospitalised? [ set p 0.05 ] ; P(Extinct ∣ Severe, no Hospitalised)

      ;; Death may occur any day within days-to-cure, with triangular probability distribution peaking at half the days-to-cure
      let day-chance cdf-triangle recovery days-to-cure ; CDF of triangle distribution

      ;; Simulate death
      let r random-float 1.0
      if r < (p * day-chance) [ pass-away ]
    ]

    ;; I+Severe
    ask sicks with [not asymptomatic? and not (severe? or deadly?)] [
      let p 0.0
      if not risky? [ set p 0.05 ]                    ; P(Severe ∣ Symptomatic, no Risky)
      if risky?     [ set p 0.25 ]                    ; P(Severe ∣ Symptomatic, Risky)

      ;; Severity may arise any day within days-to-cure, following an uniform distribution
      let day-chance recovery / days-to-cure          ; CDF of uniform distribution

      ;; Simulate severity
      let r random-float 1.0
      if r < (p * day-chance) [ set severe? true ]
    ]

    ;; I+Deadly
    ask sicks with [severe?] [
      let p 0.0
      if not risky? [ set p 0.05 ]                    ; P(Deadly ∣ Severe, no Risky)
      if risky?     [ set p 0.25 ]                    ; P(Deadly ∣ Severe, Risky)

      ;; Deadliness may arise any day within days-to-cure, following an uniform distribution
      let day-chance recovery / days-to-cure          ; CDF of uniform distribution

      ;; Simulate deadliness
      let r random-float 1.0
      if r < (p * day-chance) [ set deadly? true set severe? false ]
    ]

    ;; I+Hospitalised
    let occupied count sicks with [severe? and hospitalised?]
    let required count sicks with [severe? and not hospitalised?]
    let incoming min list (hospital-beds - occupied) required                   ; how many can be hospitalised?
    ask n-of incoming sicks with [severe? and not hospitalised?] [              ; hospitalised as many as possible
      set hospitalised? true
      set tested? true                                                          ; so they can't move even when isolation off
      if not confirmed? [ set confirmed? true set n-confirmed n-confirmed + 1 ] ; if wasn't confirmed yet, confirm and update count
    ]

    ;; I+ICU-admitted
    set occupied count sicks with [deadly? and ICU-admitted?]
    set required count sicks with [deadly? and not ICU-admitted?]
    set incoming min list (ICU-beds - occupied) required                        ; how many can be admitted to ICU?
    ask n-of incoming sicks with [deadly? and not ICU-admitted?] [              ; ICU-admit as many as possible
      set ICU-admitted? true
      set tested? true                                                          ; so they can't move even when isolation off
      if not confirmed? [ set confirmed? true set n-confirmed n-confirmed + 1 ] ; if wasn't confirmed yet, confirm and update count
    ]
  ]
end 

;; Compute Cumulative Distribution Function (CDF) of triangular distribution at a given number of days elapsed (x) within illness period (width)

to-report cdf-triangle [x width]
  ;; Height of the inner triangle to obtain CDF, i.e. the value of the PDF on x
  let height 2 * (width - abs (width - (2 * x))) / (width ^ 2)

  ;; If x is before peak, CDF is the area of the inner triangle with base x
  ifelse x <= (width / 2) [
    let base x
    report height * base / 2
  ]
  ;; Else if x is after peak, CDF is to total area (1) minus area of inner triangle with base (width - x)
  [
    let base (width - x)
    report 1 - (height * base / 2)
  ]
end 

;; Become healthy again

to get-well
  let old-shape shape
  set breed healthies
  set immune? true
  set shape old-shape
  set color white
  set speed step-size
end 

;; Become dead

to pass-away
  if not confirmed? [ set confirmed? true set n-confirmed n-confirmed + 1 ]  ; some cases are confirmed only at death
  set breed deaths
  set shape "x"
  set speed 0
  set color black
end 

;; Update counters of days and hours

to clock
  set day int (ticks / ticks-day)           ; track of number of days elapsed since beginning
  set hour int ((ticks / ticks-day) * 24)   ; track of number of hours elapsed since beginning
  sunset
end 

;; Record epidemic indicators

to indicators
  ;; Daily statistics are sicks per day and doubling-time
  if ticks mod ticks-day = 0 [
    set sicks-day lput (count sicks) sicks-day

    ;; Compute doubling time
    ;; see (Bakir, 2016, "Compound Interest Doubling Time Rule: Extensions and Examples from Antiquities.")
    if any? sicks[
      ifelse (last sicks-day) > (item (day - 1) sicks-day)
      [ set doubling-time ln 2 / ln ((last sicks-day) / (item (day - 1) sicks-day + 1E-10)) ]
      [ set doubling-time "Inf" ]
    ]
  ]

  if ticks mod (ticks-day / 24) = 0 [
    ;; Estimate R0 based on the final size of susceptibles in population
    ;; see (Dietz, 1993, "The estimation of the basic reproduction number for infectious diseases.")
    if S0 != 0 [
      let Sn count healthies with [ not immune? ] / pop-size
      set R0 (ln (S0 + 1E-10) - ln (Sn + 1E-10)) / (S0 - Sn + 1E-10)
    ]
    ;; Estimate R0_ with an individual-based approach, as the average of contagions per sick
    if any? turtles with [contagions > 0] [
      set R0_ mean [contagions] of turtles with [contagions > 0]
    ]
  ]
end 

;; Inoculate infection on one agent

to infect
  if any? healthies with [not immune?][
    ask one-of healthies with [not immune?] [ get-sick ]
  ]

  ;; Record the initial count of susceptibles for R0 computation
  set S0 (count healthies with [not immune?]) / pop-size
end 

;; Contagion by an outbreak

to outbreak
  ;; Infect a percentage of population
  repeat pop-size * (%-spread / 100) [ infect ]

  ;; Record the initial count of susceptibles for R0 computation
  set S0 (count healthies with [not immune?]) / pop-size
end 

;; Light up or down the city depending on daytime

to sunset
  if see-sunset? [
    if (ticks mod (ticks-day / 2)) = 0 [
      ifelse daytime?
      [ ask patches [ set pcolor pcolor + 4 ] ]
      [ ask patches [ set pcolor pcolor - 4 ] ]
      set daytime? not daytime?
    ]
  ]
end 

;; Setup procedures

to setup
  clear-all
  reset-ticks
  setup-globals
  setup-city
  setup-people
  setup-ambulances
  set daytime? false
end 

;; Initialize global variables

to setup-globals
  ;; Codes to identify zones
  set sort-codes n-values zones [ i -> i ]

  ;; Residents of each zone are identified with different shapes and color
  set shapes-list [ "dog" "cat" "cow skull" "fish 2" "chess knight" "ant 2" "wolf 7" "bird"
    "butterfly 2" "rabbit" "bug" "spider" "ghost" "monster" "person" "fish" "sheep" "bee 2"]
  set colors-list [ cyan pink violet green brown turquoise gray sky lime violet blue
    brown gray sky pink turquoise ]

  ;; Initalise daily sicks, positives, total-tested, total and confirmed cases
  set sicks-day [0]
  set n-positives 0 set n-tested  0 set n-cases 0 set n-confirmed 0


  ;; Probablity of contagion depending on mask wearing combination
  if not mask-sicks? and not mask-healthies? [ set p-contagion 0.9 ]  ; P(Contagion ∣ sick, healthy)
  if not mask-sicks? and     mask-healthies? [ set p-contagion 0.5 ]  ; P(Contagion ∣ sick, healthy+mask )
  if     mask-sicks? and not mask-healthies? [ set p-contagion 0.08 ] ; P(Contagion ∣ sick+mask, healthy)
  if     mask-sicks? and     mask-healthies? [ set p-contagion 0.03 ] ; P(Contagion ∣ sick+mask, healthy+mask)

  ;; Initial value for doubling time
  set doubling-time 0.0
end 

;; Split the city into zones and setup facilities

to setup-city
  set-patch-size 1.9

  ;; Layout a grid of patch square zones by assingning auxiliary codes
  let step world-width / sqrt zones
  ask patches[
    let units int ((pxcor + max-pxcor) / step)
    let tens  int ((pycor + max-pycor) / step)
    set zone 10 * tens + units
  ]

  ;; Now map those auxiliary codes to consecutive sort-codes and corresponding colors
  let aux-codes sort remove-duplicates [zone] of patches
  (foreach aux-codes sort-codes [ [old new] -> ask patches with [zone = old] [set zone new] ])
  ask patches [ set pcolor (item zone colors-list) + 1 ]
  ;foreach sort-codes [ c -> ask one-of patches with [zone = c] [set plabel zone] ] ; display zone code on view area

  ;; If housing, create houses
  if housing? [
    create-houses pop-size / 4 [
      setxy random-xcor random-ycor
      set color white set shape "house" set size 9
    ]
  ]
end 

;; Create initial population

to setup-people
  create-healthies pop-size [
    ;; Assign as homebase either a house if any, otherwise his initial position
    ifelse housing?
    [  set homebase one-of houses
       move-to homebase
    ][ setxy random-xcor random-ycor
       set homebase patch-here
    ]

    ;; Assign color and shape depending on homebase's zone
    set color item zone colors-list - 2
    set shape item zone shapes-list
    set size 9
    set heading random 360
    set scope one-of (list 25 50 100 200) ; choose radomly movement range
    set immune? false
    set unlocked? false
    set tested? false
    set contagions 0
    set speed step-size
    forward 3           ; peek out of house to be seen at the beginning
  ]
end 

;; Create ambulances

to setup-ambulances
  foreach range zones [ z ->
    create-ambulances ambulances-zone [
      move-to one-of patches with [ zone = z ]
      set homebase patch-here
      set shape "ambulance 2" set size 12 set color white
      set speed 4 * step-size        ; ambulances move faster than people
      let residents (turtle-set healthies sicks) with [ zone = z ]

      ;; The amount of tests supply depends on residents population
      set tests-supply round (count residents * (%-tests / 100) / ambulances-zone )
      set tests-stock tests-supply   ; initialise stock of tests
    ]
  ]
end 

;; Replenish the stock of available tests in ambulances

to supply-ambulances
  ask ambulances [ set tests-stock tests-stock + tests-supply ]
end 

;; Replenish an infinite stock of tests in ambulances

to supply-infinity
  ask ambulances [ set tests-stock 1000000 ]
end 

;; Save a screenshot of the interface (view+plots, etc) in a file

to screenshot
  display
  ;export-interface (word "COVID19-NPI-ABM-" random 1000 random 1000 ".png")
end