Agent-based model: AND logic gate for AMP secretion

;
; Phytophthora/AMP ABM (AND = Elicitin ∧ ROS)
; Two-phase defense: baseline containment + pathogen-specific amplification
;

patches-own [
  chemical              ; amount of AMP present
  phyto                 ; phyto mass blob
  elicitin              ; pathogen-derived PAMP produced with phyto
  ros                   ; reactive oxygen species (from AMP/phyto/glucans)
  glucan                ; β-1,3/1,6-glucan fragments released upon AMP attack
  toxin                 ; kill switch toxin

]

turtles-own [
  m                     ; mRNA level of the output gene
  y                     ; protein level of the output gene - defines on/off
  on?                   ; if on then secreting amp iif y >= y - threshold
  active?               ; activity sensor when no phyto present
]

globals[
  ; circuit parameters
  time-step
  a0 amax                ; transcription leak & max transcription gain
  beta dm dp             ; translation, mRNA decay, protein decay
  K1 K2 n1 n2            ; Hill parameters: K1/n1 elicitin, K2/n2 for ROS
  y-threshold            ; ON threshold on Y

  ;secretion in two steps
  q                      ; baseline secretion of AMP, alway on if secretion ON
  q_boost                ; extra secretion when AND gate is ON during attack

  ; dehydration (from 0 to 1) which simultaneously reduces diffusion,
  ; increases evaporation, and reduces transcriptional gain (proxy for stress)
  diffusion-rate
  evaporation-rate
  diff-min diff-max
  evap-min evap-max

  ; pathogen (phyto)
  C_eff                    ; AMP level considered effective
  phyto-initial
  phyto-initial-max        ; max value when size is = 1 biomass value in blob
  phyto-kill-rate          ; AMP driven kill rate
  phyto-decay-rate         ; natural decay of pathogen
  phyto-present?           ; true if any phyto remains

  target-protected-area    ; number of patches wanted >= C_eff
  time-to-protection       ; first tick when target is met or -1 if not yet

  phyto-radius-max       ; radius of the blob
  phyto-radius
  blob-roughness           ; 0=perfect disk, >0 adds noisy edge
  rough-min rough-max

  amp-max                  ; scale halo intensity
  amp-max-attack           ; scale halo intensity in phyto presence

  ; inputs
  elicitin-prod              ; production by phyto
  elicitin-decay
  ros-prod-amp            ; ROS made in response to AMP
  ros-prod-phyto          ; ROS from pathogen stress/presence
  ros-prod-glucan         ; ROS amplified by glucan fragments
  ros-decay
  glucan-release-rate     ; fragments released when AMP attacks phyto
  glucan-decay



  chemical-initial
  chemical-add

  reshuffle-period      ; ticks between non active/active random
  next-reshuffle-tick

  ;plotting
  amp-secreted-per-tick
  mean-m mean-y

  ; toxin for kill switch
  toxin-dose         ; uniform toxin level injected on all patches
  lethal-dose        ; threshold above which turtles die instantly

]

to-report hill[ s K n ]
  let denom ((K ^ n) + (s ^ n))
  if denom = 0 [ report 0 ]
  report ((s ^ n) / denom)
end 

to seed-phyto-blob
  let cx (min-pxcor + max-pxcor) / 2
  let cy (min-pycor + max-pycor) / 2
  ;; noisy disk include patch if distance < noisy radius
  ask patches [
    set phyto 0
    let d distancexy cx cy
    ;; noise multiplies the radius by a factor in [1 - blob-roughness, 1 + blob-roughness]
    let jitter (1 + blob-roughness * (random-float 2 - 1))
    if d <= phyto-radius * jitter [
      set phyto phyto-initial
    ]
  ]
end 

to setup
  clear-all
  set-default-shape turtles "rod2"

  set time-step 1
  set reshuffle-period 30
  set next-reshuffle-tick reshuffle-period

  set a0 0.8
  set amax 1.2
  set beta 0.5
  set dm 0.10
  set dp 0.02
  set K1 0.5
  set K2 0.8
  set n1 2.0
  set n2 2.0
  set y-threshold 2
  set q 0.2
  set q_boost 0.8
  set C_eff 0.3
  set phyto-initial-max 1.0
  set phyto-radius-max 15
  set phyto-radius   (phyto-size * phyto-radius-max)
  set phyto-initial  (phyto-size * phyto-initial-max)
  set phyto-kill-rate 0.4
  set phyto-decay-rate 0
  set target-protected-area 200
  set chemical-initial 0.0
  set chemical-add 0.0

  ; dehydration set in slider derives diffusion and evaporation
  set diff-min 10      ; lowest diffusion when very dry
  set diff-max 95      ; highest diffusion when very wet
  set evap-min 5       ; lowest evaporation when very wet
  set evap-max 80      ; highest evaporation when very dry

  ; from current dehydration slider
  set diffusion-rate   (diff-max - (diff-max - diff-min) * dehydration)
  set evaporation-rate (evap-min + (evap-max - evap-min) * dehydration)

  set amp-max 0.5
  set amp-max-attack 1.5

  set elicitin-prod    0.02
  set elicitin-decay 0.0005
  set ros-prod-amp    0.05
  set ros-prod-phyto   0.02
  set ros-prod-glucan 0.06
  set ros-decay 0.10
  set glucan-release-rate  0.04
  set glucan-decay 0.01


  set toxin-dose      1.0
  set lethal-dose     0.5


  ask patches [
    set chemical chemical-initial
    set phyto 0
    set elicitin 0
    set ros 0
    set glucan 0
     set toxin 0
    set pcolor black
  ]

  set time-to-protection -1

  set rough-min 0.00    ; more compact when wet
  set rough-max 1    ; drying means rougher
  set blob-roughness (rough-min + (rough-max - rough-min) * dehydration)

  seed-phyto-blob

  ; bacteria is randomly scattered
  create-turtles population [
    setxy random-xcor random-ycor
    set size 2
    set m 0
    set y random-float 1
    set on? false
    set active? false
    ; initial color while pathogen exists = OFF > orange
    set color orange
  ]


  set phyto-present? any? patches with [ phyto > 0.0001 ]
  ask patches [
    ifelse phyto > 0 and phyto-initial > 0
      [ let b min list 1 (phyto / phyto-initial)
        set pcolor rgb round (120 * b) round (30 * b) round (25 * b) ]
    [ set pcolor black ]
  ]


  ; each turtle has a 50% chance of being active at start when no phyto present
  if not phyto-present? [
    ask turtles [ set active? (random-float 1 < 0.5) ]
  ]

  if plot? [
    clear-all-plots
  ]

  reset-ticks
end 

;;;;;;;;;;;;;;;;;;;;;
;;; Go procedures ;;;
;;;;;;;;;;;;;;;;;;;;;

to go  ; forever button

  set phyto-present? any? patches with [ phyto > 0.0001 ]

  set amp-secreted-per-tick 0

  if not phyto-present? [
    if ticks >= next-reshuffle-tick [
      ask turtles [ set active? (random-float 1 < 0.5) ]
      set next-reshuffle-tick ticks + reshuffle-period
    ]
  ]

  ask turtles [ cell-step ]

  ; update each tick so can see changes without resetting
  set diffusion-rate   (diff-max - (diff-max - diff-min) * dehydration)
  set evaporation-rate (evap-min + (evap-max - evap-min) * dehydration)

  ; global background AMP injection
  if chemical-add > 0 [
    ask patches [ set chemical chemical + chemical-add * time-step ]
  ]

  ; dehydration reduces diffusion, increases evaporation
  ; use the dehydration-derived rates directly
  diffuse chemical (diffusion-rate / 100)

  ask patches [
    set chemical (chemical * (100 - evaporation-rate) / 100)
  ]


  ;   ; AMP-driven kill + natural decay
  ;  ask patches [
  ;    let kill phyto-kill-rate * chemical
  ;    set phyto max list 0 ( phyto - (phyto-decay-rate * phyto + kill) * time-step )
  ;  ]
  ;  ; this method not strong enough too linear, need stringer kill effect for visible results on blob


  ask patches [
    let Hkill ((chemical ^ 3) / ((0.8 ^ 3) + (chemical ^ 3) + 1e-9))  ;; MIC≈0.8, n=3
    let frac-kill 0.25 * Hkill        ;; max 25% of phyto per tick near/above MIC
    let dphyto (phyto-decay-rate * phyto + frac-kill * phyto)
    set phyto max list 0 (phyto - dphyto * time-step)
  ]

  ; input dynamics on patches
  ask patches [
    ;; elicitin from pathogen present?
    set elicitin max list 0 (elicitin
      + elicitin-prod * phyto * time-step
      - elicitin-decay * elicitin * time-step)

    ; glucan fragments released when AMP attacks phyto
    set glucan max list 0 (glucan
      + glucan-release-rate * chemical * phyto * time-step
      - glucan-decay * glucan * time-step)

    ; ROS from AMP, pathogen stress, and glucan-augmented signaling
    let amp_ros_term (ifelse-value (phyto > 1e-6) [ ros-prod-amp * chemical ] [ 0 ])
    set ros max list 0 (ros
      + (amp_ros_term + ros-prod-phyto * phyto + ros-prod-glucan * glucan) * time-step
      - ros-decay * ros * time-step)
  ]

  ; any phyto remains?
  set phyto-present? any? patches with [ phyto > 0.0001 ]


  ; plots
  let turtle-count count turtles
  let frac_on (ifelse-value (turtle-count = 0) [0]
    [ (count turtles with [ on? ]) / turtle-count ])
  let frac_active (ifelse-value (turtle-count = 0) [0]
    [ (count turtles with [ active? ]) / turtle-count ])

  let total_amp   sum [ chemical ] of patches
  let total_phyto sum [ phyto ] of patches
  let total_elic  sum [ elicitin ] of patches
  let total_ros   sum [ ros ] of patches
  let total_glu   sum [ glucan ] of patches
  let protected_area count patches with [ chemical >= C_eff ]

  set mean-m (ifelse-value (turtle-count = 0) [0] [ mean [ m ] of turtles ])
  set mean-y (ifelse-value (turtle-count = 0) [0] [ mean [ y ] of turtles ])

  ; capture first time protection target is reached
  if (time-to-protection = -1 and protected_area >= target-protected-area) [
    set time-to-protection ticks
  ]

; update plots
  if plot? [
; 1. Secretion
  set-current-plot "AMP Secretion"
  set-current-plot-pen "secreted"
  plot amp-secreted-per-tick

  ; Population states
  set-current-plot "Population states (on/off)"
  set-current-plot-pen "frac-on"      plot frac_on
  set-current-plot-pen "frac-active"  plot (ifelse-value phyto-present? [0] [frac_active])

  ; Concentrations
  set-current-plot "Field concentrations"
  set-current-plot-pen "AMP"      plot total_amp
  set-current-plot-pen "Phyto"    plot total_phyto
  set-current-plot-pen "Elicitin" plot total_elic
  set-current-plot-pen "ROS"      plot total_ros
  set-current-plot-pen "Glucan"   plot total_glu

  ]

  ;; recolor patches: show Phytophthora when present; otherwise show AMP old code does not work
  ;  ask patches [
  ;    ifelse (phyto > 0)
  ;    [ set pcolor scale-color blue phyto 0 phyto-initial ]
  ;    [ set pcolor scale-color green chemical 0.1 5 ]
  ;  ]
;
;  ask patches [
;    ;;new colouring that works
;    let g round (255 * chemical / amp-max)          ;; AMP intensity (green)
;    let b round (255 * (ifelse-value (phyto-initial > 0) [phyto / phyto-initial] [0]))    ;; Phyto intensity (blue)
;    set g max list 0 (min list 255 g)
;    set b max list 0 (min list 255 b)
;
;    let R 255 - (min list 255 (g + b))   ;; both AMP and Phyto dim red
;    let G2 255 - b                        ;; Phyto dims green
;    let B2 255 - g                        ;; AMP dims blue
;
;    set pcolor rgb R G2 B2
  ;  ]

  ask patches [
    let ampscale (ifelse-value phyto-present?
      [ amp-max-attack ]
      [ amp-max ])
    let ampN  min list 1 ((ifelse-value (ampScale > 0) [ chemical / ampScale ] [ 0 ]) ^ 0.6)
    let phyN  min list 1 (ifelse-value (phyto-initial > 0) [ phyto / phyto-initial ]     [ 0 ])
    let g round (255 * ampN)          ;; AMP contribution (green channel)
    let pR round (120 * phyN)         ;; phyto maroon R
    let pG round ( 30 * phyN)         ;; phyto maroon G
    let pB round ( 25 * phyN)         ;; phyto maroon B
    let R pR
    let G2 pG + g
    let B pB

    set pcolor rgb (max list 0 (min list 255 R))
    (max list 0 (min list 255 G))
    (max list 0 (min list 255 B))
  ]
  tick
end 

to cell-step  ; turtle procedure
  ; immediate kill if local toxin >= lethal dose
  if [toxin] of patch-here >= lethal-dose [
    die
    stop
  ]

  ; if no phyto, go into radom ON/OFF = random active?
  if not phyto-present? [
    set on? active?

    if secretion? [
      ask patch-here [
        let q_idle (0.05 * q)                 ; 5% leak when OFF
        let add (ifelse-value [on?] of myself
          [ q * time-step ]               ; ON
          [ q_idle * time-step ])         ; OFF
        set chemical chemical + add
        set amp-secreted-per-tick amp-secreted-per-tick + add
      ]
    ]

    rt (random-float 6 - 3)
    if not can-move? 0.01 [ set heading towardsxy 0 0 ]
    fd 0.01
    ;; simple idle coloring
    set color (ifelse-value on? [cyan] [orange])
    stop
  ]

  ; sense inputs on current patch
  let s_elicitin [elicitin] of patch-here
  let s_ros      [ros]      of patch-here

  ; AND gate & dehydration-attenuated transcription gain
  let H1 hill s_elicitin K1 n1
  let H2 hill s_ros      K2 n2
  let gE (0.6 + (1.0 - 0.6) * (1 - dehydration))  ; moisture effect on expression gain: drier -> smaller gE, dehydration = 0 then gE=1.0

  ; ODE tau-leap with tiny uniform noise = adding noise, stochastic gene expression
  let dm_dt ((a0 * gE) + (amax * gE *  (H1 * H2))) - ((dm + 0.01) * m)
  ; F(H1,H2) = H1*H2 -> AND, if we want NAND -> change to  F(H1,H2) = 1 - (H1*H2)
  let dY_dt ((beta * m)) - ((dp + 0.01) * y)

  set m max (list 0 (m + dm_dt * time-step + (random-float 0.2 - 0.1)))
  set y max (list 0 (y + dY_dt * time-step + (random-float 0.2 - 0.1)))

 ; set on? (y >= y-threshold)   ;; ON/OFF decision OLD

  ; new
  ; probabilistic ON: quicker near the blob, slower farther away
  ifelse( y < y-threshold) [
    set on? false]
  [
    ;; p = exp(-d/6) with 0.1 floor; near center p≈1, far away p→0.1
    let d distancexy ((min-pxcor + max-pxcor) / 2) ((min-pycor + max-pycor) / 2)
    let p max list 0.1 (exp (- d / 6))
    if random-float 1 < p [ set on? true ]
  ]


;  ;; two-phase secretion: baseline containment always; boost only if AND is ON and invasion persists OLD
;  if secretion? [
;    ask patch-here [
;      let add (q * time-step)
;      if [on?] of myself and [phyto-present?] of myself [
;        set add add  + ( q_boost * time-step)
;      ]
;      set chemical chemical + add
;      set amp-secreted-per-tick amp-secreted-per-tick + add
;    ]
;  ]

  ;new
  if secretion? [
    if phyto-present? [
      ; distance-weighted secretion only during invasion
      let d distancexy ((min-pxcor + max-pxcor) / 2) ((min-pycor + max-pycor) / 2)
      let scale max list 0.15 (exp (- d / 6))   ;; stronger falloff; 0.05 floor
      ask patch-here [
        ; OFF = no baseline during invasion; ON = baseline + boost, both scaled
        let add (ifelse-value [on?] of myself
          [ (q * scale + q_boost * scale) * time-step ]
          [ 0 ])
        set chemical chemical + add
        set amp-secreted-per-tick amp-secreted-per-tick + add
      ]
    ]
  ]


  ; gentle wobble and inwards push
  rt (random-float 6 - 3) ;; uniform wobble
  if not can-move? 0.01 [ set heading towardsxy 0 0 ]
  fd 0.01

  ifelse on? [ set color orange ] [ set color cyan ]
  ; COLOURS:
  ; phyto present: ON: yellow, OFF: red
  ; phyto cleared: mission complete: cyan
  ifelse phyto-present? [
    ifelse on?
    [ set color yellow]   ; ON & pathogen present
    [ set color orange ]      ; OFF & pathogen present
  ] [
    set color cyan          ; pathogen cleared  mission complete
  ]
end 

to kill-switch
  ; inject lethal toxin
  ask patches [ set toxin toxin-dose ]
end