Plant Speciation

turtles-own [ tolerance fitness flower-time seedling? will-die? ]
patches-own [ metal barrier? ]

globals [
   day                             ;; monitor value
   year                            ;; current year
   old-year                        ;; previous year
end -of-days-counter             ;;
   transition-time?                ;;
   old-visualize-time-steps-state  ;; allows for switching between visualization modes during model run
   chance-death-per-year           ;; the probability that the plant will die this year
   chance-seed-dispersal           ;; the probability that a newborn plant grows in a patch adjacent to its parents', instead of in the same patch.
   average-number-offspring        ;; avg. number of seeds dispersed by the plant when it does disperse them
   percent-same-flowering-time
]

to setup
  clear-all
  ;; initialize globals
  set day 0
  set old-year 0
  set year 0
  set end-of-days-counter 0
  set percent-same-flowering-time 0
  set chance-death-per-year 10 ;; set to 10% by default
  set chance-seed-dispersal 50 ;; set to 50% by default
  set average-number-offspring 3
  set transition-time? false
  set old-visualize-time-steps-state visualize-time-steps

  ;color the frontier
  ask patches [
    set barrier? false
    setup-two-regions
    set pcolor calc-patch-color metal
  ]

  ;spawn the initial population -- carrying-capacity-per-patch allows more than plant at this patch
  ask patches with [pxcor = min-pxcor] [
    sprout plants-per-patch [
      if initial-tolerance = "all no tolerance" [set tolerance 0]
      if initial-tolerance = "all full tolerance" [set tolerance 100]
      if initial-tolerance = "random tolerances" [set tolerance random-float 100]
      set flower-time (365 / 2)
      set heading random 360
      fd random-float .5
      set fitness 1
      set shape "plant"
      set seedling? false
      set will-die? false
      set color calc-plant-color tolerance
    ]
  ]
  reset-ticks
end 

to setup-two-regions
  ;; set the metal value based on the width designated by the frontier-sharpness slider value
  set metal precision (100 / (1 + exp (frontier-sharpness * ((max-pxcor + min-pxcor) / 2 - pxcor)))) 0
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;; RUNTIME PROCEDURES ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to go
   check-labels
   ;; if visualization switched during model run, redraw plants
   if old-visualize-time-steps-state != visualize-time-steps [
     ask turtles [
       redraw-plants-as-full-sized-plants
     ]
   ]

  if visualize-time-steps = "years" [
    set day 0 ;; always at day 0 in years mode
    set year year + 1
    ask turtles [ redraw-plants-as-full-sized-plants  ]
    do-start-of-new-year-events
  ]

  if visualize-time-steps = "days" [
    visualize-bloom
    do-end-of-days-events
  ]
  tick
end 

to do-start-of-new-year-events
  if year > old-year [
    do-reproduction
    mark-turtles-to-kill
    kill-marked-turtles
    set old-year year
  ]
end 

to do-end-of-days-events
  ifelse day = 365 [  ;; at end of year
    if transition-time? = false  [
      do-reproduction
      mark-turtles-to-kill
      set transition-time? true
    ]
    ;; 10 ticks of transition time will elapse between years when visualize-time-steps is in "day" mode
    ;; this transition time is for animation events to show death and growth of new plants before the next
    ;; year begins. It is done before day 0 of the new year so that these visualizations do not interfere
    ;; with flower time speciation mechanisms in the model

    if transition-time? [
      set end-of-days-counter end-of-days-counter + 1
      ;; grow the new plants to be visible for the next season
      ask turtles with [seedling?] [visualize-seedling-growth]
      ;; shrink and fade the plants that are going to die
      ask turtles with [will-die?] [set size (1 - (end-of-days-counter / 10))]

      if end-of-days-counter > 9 [  ;; start a new year
        set year year + 1
        set end-of-days-counter 0
        set transition-time? false
        set day 0
        ask turtles [turn-seedlings-into-full-plants]
        kill-marked-turtles
      ]
    ]
  ]

  ;; the else event
  [set day day + 1]
end 

;; make babies

to do-reproduction
  ask turtles [
    ;; construct list of potential mates based on pollen-radius slider
    let potential-mates []
    let nearby-turtles turtles
    if (pollen-radius < (max-pxcor - min-pxcor)) [set nearby-turtles turtles in-radius pollen-radius]
    let num-candidates 10
    ifelse count nearby-turtles < 10 [ set potential-mates (sort nearby-turtles) ]
                                     [ set potential-mates (sort (n-of 10 nearby-turtles))
                                       set potential-mates (fput self potential-mates)]

    ;; pick mate randomly weighted by compatibility
    let compatibilities map [compatibility self ?] potential-mates
    let mate (pick-weighted (potential-mates) (compatibilities))

    ;; spawn children
    hatch (random-poisson average-number-offspring) [
      set seedling? true
      set will-die? false
      ifelse visualize-time-steps = "days" [ set size 0.0] [set size 1.0]

      ;combine parents' genes and give average value of to the child
      if genetics-model = "avg. genotype" [
        set tolerance (tolerance + [tolerance] of mate) / 2
        set flower-time (flower-time + [flower-time] of mate) / 2
      ]

      ;mutate tolerance gene
      if (random-float 100) < chance-tolerance-mutation [
        set tolerance (tolerance + (random-normal 0 tolerance-mutation-stdev))
      ;; set tolerance (tolerance + random tolerance-mutation - random tolerance-mutation)
      ]

      ;mutate flowering time gene
      if (random-float 100) < chance-flower-time-mutation [
        set flower-time (flower-time + (random-normal 0 flower-time-mutation-stdev))
      ]

      ;; keeps values from going above a min and max
      if tolerance < 0 [ set tolerance  0 ]
      if tolerance > 100 [ set tolerance 100 ]
      if flower-time < 0 [ set flower-time 0 ]
      if flower-time > 365 [ set flower-time 365 ]

      ;change color to reflect metal tolerance
      set color calc-plant-color tolerance
      migrate-this-plant
    ]

   if plant-type = "annual" [set will-die? true] ;end of the generation
  ]
end 

;kill ill-adapted turtles and fix solve overpopulation

to mark-turtles-to-kill
  ask turtles [
    let t tolerance / 100
    let m metal / 100
    ;; Fitness is a linear function dependent on tolerance whose slope and y-intercept
    ;; vary linearly with respect to an increases in metal amount.
    ;; This linear function would have the following slopes in various metal levels:
    ;; (i.e. negative slope in clean ground -> high tolerance is bad
    ;;       positive slope in dirty ground -> high tolerance is good
    ;;       zero     slope in between      -> no benefit or disadvantage to any level of tolerance )

    ;; This is a model of a "tradeoff", where specializing in one variation of trait is advantageous
    ;; in one environmental extreme, but specializing in another variation of the trait is advantageous in a different
    ;; environmental extreme. Intermediate "hybridization" or averaging between both variations is disadvantageous
    ;; in both environments or at the very least it is not advantageous in either extreme environment.
    ;; Such tradeoff models can lead to speciation when other traits permit a population to reproductively
    ;; fragment and isolate itself into non-interbreeding sub populations.

    ;; This makes for a hyperbolic paraboloid or saddle shaped function that is dependent on metal amount and
    ;; tolerance.  A general form of this fitness function would be the following:
    ;; set fitness ((1 +  (A * t * m + B * t * m - C * t * m) - ( A * t + B * m) ) )
    ;; where fitness is 1 at clean ground and no tolerance
    ;; A is the penalty (0 to 1) for having tolerance in clean ground, therefore fitness is (1 - A)
    ;; B is the penalty (0 to 1) for having the highest level of metal in the ground and no tolerance, therefore fitness is (1 - B)
    ;; C is the penalty (0 to 1) for having the highest tolerance in the highest level of metal, therefore fitness is (1 - C)
    ;; As long as C is less than both B and A, then you will have a fitness function that can be used in this section
    ;; The fitness function has been hard coded here to use A = .4 and B = .4 and C = 0

       set fitness (1 - m) * (1 - .4 * t) + m * (1 - .4 * (1 - t))

    ;; survival probability based on fitness
    if (random-float 1) > (fitness) [set will-die? true]
    ;; survival probability based on fixed number of additional deaths-a-year
    if random-float 100 < chance-death-per-year [set will-die? true]
  ]

  ;; In overpopulated patches, kill the least fit plant until we are down to the carrying capacity
  ask patches [
    let overpopulation ((count turtles-here) - plants-per-patch)
    if overpopulation > 0 [
      ask (min-n-of overpopulation turtles-here [fitness]) [
        set will-die? true
      ]
    ]
  ]
end 

to migrate-this-plant ;; turtle procedure
  ;; Some plants grow in patch adjacent to the parent plant. This represents migration from seed dispersal
  if (random-float 100) < chance-seed-dispersal [
    move-to one-of neighbors
    rt random 360  fd random-float 0.45  ;;spread out plants that are on the same patch
  ]
end 

to kill-marked-turtles
    ask turtles with [will-die?] [die]
end 

to turn-seedlings-into-full-plants
         if seedling? [set seedling? false]
         redraw-plants-as-full-sized-plants
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;; visualization procedures ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to check-labels
  ask patches [
  ifelse show-labels-as = "metal in soil"
    [set plabel metal]
    [set plabel ""]
  ]
  ask turtles [
    if show-labels-as = "metal tolerance" [set label precision tolerance 0]
    if show-labels-as = "flower time"  [set label precision flower-time 0 ]
    if show-labels-as = "none" [set label ""]
  ]
end 

to redraw-plants-as-full-sized-plants ;; turtle procedure
  set shape "plant"
  set size 1
end 

to visualize-seedling-growth ;; turtle procedure
  if seedling? [set size (end-of-days-counter / 10) ]
end 

to visualize-bloom
  ask turtles [
    ;; If day of the year is greater than this plants flower time and less than the flower time plus the length of flower time, then it is in the
    ;; the flowering time window and a flower should be located here
    ifelse day >= (flower-time  ) and day <= (flower-time + (flower-duration)  )
       [ set shape "flower"]
       [ set shape "plant" ]
  ]
end 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;; reporters  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; make a random choice from  weighted by 
;; where highly weighted choices are more likely to be selected than others

to-report pick-weighted [ options weights ]
  let wsum 0
  foreach weights [
    set wsum (wsum + ?)
  ]
  let wret wsum * (random-float 1)
  let ret 0
  set wsum 0
  foreach weights [
    set wsum (wsum + ?)
    if wsum > wret [ report (item ret options) ]
    set ret (ret + 1)
  ]
end 

;; compatibility is a decreasing function of the difference in flowering times
;; since plants that have more of an overlap between flower-times are more likely
;; they will pollinate one another

to-report compatibility [ t1 t2 ]
  let diff abs ([flower-time] of t1 - [flower-time] of t2)
  ifelse diff < flower-duration [ report (flower-duration - diff) ] [ report 0 ]
end 

;; calculate the patch color based on the presence of metal in the soil

to-report calc-patch-color [ m ]
  report rgb 0 (255 * (1 - (m / 100)) / 2) (255 * (m / 100) / 2)
end 

to-report calc-plant-color [ t ]
  let black-pcolor rgb 0 0 0
  ifelse barrier?
    [report black-pcolor]
    [report rgb 0 (255 * (1 - (t / 100)) ) (255 * (t / 100) )]
end 


; Copyright 2012 Uri Wilensky.
; See Info tab for full copyright and license.