Agent-based Simulation of An SIR epidemic model for COVID-19 spread with fuzzy parameter: the case of Indonesia

globals [
  S           ; Susceptible proportion
  I           ; Infected proportion
  R           ; Recovered proportion
  beta        ; Infection rate
  mu_c        ; COVID death rate
  gamma       ; Recovery rate
  R0          ; Basic reproduction number

  ; Fuzzy parameters
  ;omega_min   ; Minimum virus load (default: 10)
  ;omega_0     ; Mid virus load (default: 100)
  ;omega       ; Current virus load (input via slider)

  ; Control parameters
  ;tau         ; Vaccine effectiveness (0-1)
  ;piy          ; Health protocol compliance (0-1)
  ;theta       ; Treatment effectiveness (0-1)
  ;vartheta    ; Medical factor (0-1, default: 0.9)

  ; Basic rates
  ;mu          ; Natural birth/death rate (default: 0.00625)
  ;mu0_c       ; Base COVID death rate (default: 0.00022114)
  ;gamma0      ; Base recovery rate (default: 0.001042)
]

breed [susceptibles susceptible]
breed [infecteds infected]
breed [recovereds recovered]

to setup
  clear-all
  ; Initialize proportions from paper data (scaled to proportions)
  set S 268757171 / 269600000  ; ≈0.996
  set I 457735 / 269600000     ; ≈0.0017
  set R 385094 / 269600000     ; ≈0.0014

  ; Initialize plots
  reset-ticks
end 

to go
  ; Update fuzzy parameters
  compute-beta
  compute-mu_c
  compute-gamma

  ; Calculate derivatives
  let dS (mu - (beta * S * I) - (mu + tau + piy) * S)
  let dI ((beta * S * I) - (mu + mu_c + theta + gamma) * I)
  let dR ((theta + gamma) * I + (piy + tau) * S - mu * R)

  ; Update proportions using Euler method
  set S S + dS * 0.1  ; Smaller time step for stability
  set I I + dI * 0.1
  set R R + dR * 0.1

  ; Calculate R0
  set R0 (beta * mu * (1 - tau) * (1 - piy)) /
         ((piy + tau + mu) * (theta + mu_c + gamma + mu))

  ; Ensure proportions stay reasonable
  set S max list 0 S
  set I max list 0 I
  set R max list 0 R

  tick
  update-plots
end 

to compute-beta
  ; Calculate infection rate based on virus load
  ifelse omega <= omega_min [
    set beta 0
  ] [
    ifelse omega < omega_0 [
      set beta ((omega - omega_min) * (1 - tau) * (1 - piy)) /
              (omega_0 - omega_min)
    ] [
      set beta (1 - tau) * (1 - piy)
    ]
  ]
end 

to compute-mu_c
  ; Calculate COVID death rate
  ifelse omega < omega_0 [
    set mu_c (((1 - vartheta) - mu0_c) * (1 - theta) * (omega / omega_0)) + mu0_c
  ] [
    set mu_c (1 - vartheta) * (1 - theta) + theta * mu0_c
  ]
end 

to compute-gamma
  ; Calculate recovery rate
  ifelse omega < omega_0 [
    set gamma ((gamma0 - 1) * (1 - theta) * (omega / omega_0)) + 1
  ] [
    set gamma gamma0 * (1 - theta) + theta
  ]
end 

to updates-plots
  set-current-plot "Population Proportions"
  plotxy ticks S
  set-current-plot-pen "Infected"
  plotxy ticks I
  set-current-plot-pen "Recovered"
  plotxy ticks R
end