A Trojan horse approach to medical intervention strategies

; Revisions:
; 2013-08-02
;   - accuracy now reports ( 100% - per-capita-mean-error )
;   - diffuse-individual moved into .nls file
;   - output-wrap moved into .nls file
;   - sigfigs moved into .nls file
; 2013-07-29
;   - line-breaking procedure to replace output-prints
;   - moved some generic code into includes file
;   - events are defined at beginning and applied behind the scenes (in included file)

; design choices:
; * the (full) model is agent-based
; * bacteria (wilds and cheaters) are mobile (Brownian)
; * their products (food, toxin, and bacteriocide) are stationary
; * interactions are local (per patch)
; * time steps forward in increments large enough for multiple events to occur
; * the full model fixes the fast (food, toxin, & bacteriocide) production rates
; * the System Dynamics Modeler contains a numeric approximation
; * the numeric model applies the QSSA to reduce the system, removing the products
; * the numeric model approximates the equations in Brown et al., 2009

__includes [
  "diffuse-individual.nls"
  "output-wrap.nls"
  "per-capita-tau-leaping.nls"
  "sigfigs.nls"
]


; the possible "breeds" of turtles
breed [wilds wild]
breed [cheaters cheater]


; turtles own event rates
patches-own [
  food
  toxin
  bacteriocide
]


globals [ 
  count-patches  ; assign to variable to improve speed; used in setup and plots
  beta-food
  delta-toxin
  gamma-bacteriocide
  fast-rate      ; = beta = delta = gamma (fast consumption rates for food, toxin, & bacteriocide)
  wild-virulence
]

to startup
  setup
end 

to setup
  clear-all
  set fast-rate 5 ; TO DO: how big does this need to be to satisfy QSSA?
  set wild-virulence 2 ; assume wild type twice as virulent as cheaters (as per Brown et al.)
  set beta-food fast-rate
  set delta-toxin fast-rate
  set gamma-bacteriocide fast-rate
  set count-patches count patches
  system-dynamics-setup ; setup system-dynamics-modeler for numeric comparison
  ; set breed shapes
  set-default-shape turtles "circle"
  ; create numbers of each breed to get the desired average density (per patch)
  create-wilds       ( W-num * count-patches ) [ set color green ]
  create-cheaters    ( C-num * count-patches ) [ set color red   ]
  ; shuffle them 
  ask turtles [ setxy random-xcor random-ycor ]
  ask patches [
    ; start with QSSA densities of food, toxin, & bacteriocide
    let count-turtles-here  count turtles-here
    let count-cheaters-here count cheaters-here
    if ( count-turtles-here > 0 ) [
      ; Poisson-distributed random variates
      set food random-poisson ( b-shared-benefit * count wilds-here / ( beta-food * count-turtles-here ) )
      set toxin random-poisson ( a-toxin-production * count-cheaters-here / ( delta-toxin * count-turtles-here ) )
      set bacteriocide random-poisson ( e-bacteriocinogen * count-cheaters-here / ( gamma-bacteriocide * count-turtles-here ) )
    ]
  ]
  reset-ticks
  ; hint
  output-wrap 26 ; set wrapping column
  output-wrap "An agent-based model of 'Social evolution in micro-organisms and a Trojan horse approach to medical intervention strategies' by Brown et al. (2009)."
  output-wrap "Press go to get started.\n"
  
  ; merge all events into one master list [2013-07-29]
  add-per-capita-event turtles                  ; N -> 2 N @ rate 1
    task [ 1 ] 
    task [ hatch 1 ]
  add-per-capita-event turtles                  ; N + N -> N @ rate 1
    task [ count other turtles-here ] 
    task [ die ]
  add-per-capita-event turtles                  ; N + F -> 2 N @ rate beta
    task [ beta-food * food ]
    task [ set food food - 1
           hatch 1           ]
  add-per-capita-event turtles                  ; N + T -> 0 @ rate delta
    task [ delta-toxin * toxin ]
    task [ set toxin toxin - 1    
           die                 ]
  add-per-capita-event wilds                    ; W -> 0 @ rate x
    task [ x-wild-cost ]
    task [ die ]
  add-per-capita-event wilds                    ; W -> W + F @ rate b
    task [ b-shared-benefit ]
    task [ set food food + 1 ]
  add-per-capita-event wilds                    ; W + B -> 0 @ rate gamma
    task [ gamma-bacteriocide * bacteriocide ]
    task [ set bacteriocide bacteriocide - 1    
           die                               ]
  add-per-capita-event cheaters                 ; C -> 0 @ rate q
    task [ q-cheater-cost ]
    task [ die ]
  add-per-capita-event cheaters                 ; C -> C + T @ rate a
    task [ a-toxin-production ]
    task [ set toxin toxin + 1 ]
  add-per-capita-event cheaters                 ; C -> C + B @ rate e
    task [ e-bacteriocinogen ]
    task [ set bacteriocide bacteriocide + 1 ]
  add-per-capita-event cheaters                 ; C + B -> C @ rate gamma
    task [ gamma-bacteriocide * bacteriocide ]
    task [ set bacteriocide bacteriocide - 1 ]
end 

to go
  if ticks = 0 [ 
    reset-timer 
    output-wrap "You may adjust all parameters while the simulation runs."
    output-wrap "Compare the numeric, non-spatial (faint) vs. agent-based (bold) dynamics in the graphs."
    output-wrap ( word "It is assumed that the wild type is " wild-virulence " times as virulent as the cheater.\n" )
  ]
  system-dynamics-go ; numeric model
  set dt per-capita-tau-leap inaccuracy dt
  diffuse-individual turtles diffusion-const dt
  do-plots
  if (dt = 0) [  ; simulation over
    output-wrap (word "Finished in " timer " seconds.")
    stop 
  ]
end 

to do-plots
  set-current-plot "dynamics"
  system-dynamics-do-plot
  update-plots ; update other plots ("Edit" the plots to see code)
end 

to-report diffusion-const
  report ifelse-value well-mixed [ 0 ][ 10 ^ log-diffusion ]
end 

to-report inaccuracy
; complement of accuracy, depends on err-tolerance
  report 10 ^ ( err-tolerance )
end