inbreeding and outbreeding depression model

; ============================================================
; INBREEDING & OUTBREEDING DEPRESSION
; An Agent-Based Model for Conservation Genetics
; THEORETICAL BACKGROUND:
;
; 1. INBREEDING DEPRESSION
;    Mechanism: deleterious recessive alleles become homozygous as relatedness
;    rises in small populations, exposing fitness costs.
;    Two submechanisms:
;      (a) Partial dominance: recessive deleterious alleles expressed in homozygotes
;      (b) Overdominance: heterozygote advantage lost under inbreeding
;    Key parameters: s (selection coefficient), h (dominance: h=0 fully recessive,
;    h=0.5 additive, h>0.5 dominant).
;    Real examples: Florida panther (heart defects, kinked tails),
;    Swedish adders (stillborn offspring), Mauritius kestrel (50% het lost),
;    Isle Royale wolves (vertebral malformations, 80%+ F-coefficient).
;    Literature: Charlesworth & Willis (2009) Nat Rev Genet 10:783-796
;                Frankham (1995) Ann Rev Genet 29:305-327
;
; 2. OUTBREEDING DEPRESSION
;    Mechanism: disruption of locally coadapted gene complexes when
;    genetically divergent populations interbreed.
;    Three submechanisms:
;      (a) Fixed chromosomal differences → F1 hybrid sterility
;      (b) Disruption of coadapted gene complexes → F1 or F2 fitness reduction
;      (c) Loss of local adaptation: offspring maladapted to parental environments
;    Expressed in F1 (immediate) and worsens in F2 (recombination breaks coadaptations).
;    Real examples: Tatra Mountain ibex (wrong calving season after hybridization),
;    Ipomopsis aggregata (spatial-scale-dependent outbreeding depression).
;    Literature: Edmands (2007) Mol Ecol 16:463-475
;                Frankham et al. (2011) Conserv Biol 25:55-63
;
; 3. GENETIC RESCUE
;    Introducing migrants into inbred populations masks deleterious recessives
;    via heterosis and restores genetic variation.
;    Risk: outbreeding depression if populations are too diverged (long isolation,
;    strong local adaptation, chromosomal rearrangements).
;    The "between a rock and a hard place" dilemma (Edmands 2007):
;    when is a species unique enough that outbreeding might cause problems?
;    Evidence: SUCCEEDS for Florida panther (Texas cougar intro), Illinois prairie
;    chicken, Swedish adders. FAILS risk: Tatra ibex, genetically distinct subspecies.
;    Guidelines (Frankham et al. 2011): Use gene flow if (1) inbreeding depression
;    evident, (2) source population is not chromosomally incompatible, (3) populations
;    not isolated >500 years or in markedly different environments.
;    Literature: Ralls et al. (2020) Biol Conserv 252
;                Hedrick & Garcia-Dorado (2016) Trends Ecol Evol 31:940-952
;
; 4. PURGING
;    Natural selection removes deleterious recessives when they are exposed in
;    homozygotes under inbreeding. Purging is most effective for LETHAL alleles
;    (high s, low h) but slow and incomplete for mildly deleterious ones.
;    Empirical verdict: purging cannot be relied on to prevent inbreeding depression;
;    it slows extinction slightly but does not eliminate the threat.
;    Literature: Frankham et al. (2001) Conserv Genet 2:279-285
;                Garcia-Dorado et al. (2016) Heredity 116:580
;
; 5. Ne/N RATIO
;    Effective population size (Ne) << census size (N). Average Ne/N ≈ 0.11
;    due to variance in reproductive success, unequal sex ratios, overlapping
;    generations. This means genetic erosion happens much faster than expected.
;    Critical thresholds: Ne ≈ 50 (short-term inbreeding), Ne ≈ 500 (long-term
;    evolutionary potential). The 50/500 rule (Franklin 1980; Soule 1980).
;    v2 computes Ne-A / Ne-B live from the sex-ratio component of this formula.
;
; GENOME STRUCTURE (diploid):
;   Positions 0 to (n-del-loci - 1):           Deleterious recessive loci
;     allele 0 = wildtype, allele 1 = deleterious
;   Positions n-del-loci to (n-del-loci + n-local-loci - 1): Local adaptation loci
;     Pop A local optimum = allele 0
;     Pop B local optimum = allele 1
;   Positions (n-del-loci + n-local-loci) to (total-1): Neutral loci
;     For heterozygosity/diversity tracking only
;
; FITNESS MODEL:
;   W = W_inbreeding * W_local_adaptation
;   W_inbreeding = PRODUCT over del loci of:
;       (1 - s)    if homozygous deleterious [aa]
;       (1 - h*s)  if heterozygous          [Aa]
;   W_local = 1 - local-adapt-strength * (n_mismatched_alleles / max_possible_mismatch)
;
; PARAMETERS (set via sliders):
;   n-populations       : 1 (isolated) or 2 (two-pop with possible gene flow)
;   initial-pop-A/B     : starting census sizes
;   carrying-cap-A/B    : ecological carrying capacities
;   n-del-loci          : number of deleterious loci (depth of genetic load)
;   del-allele-freq     : starting frequency of deleterious alleles
;   sel-coeff (s)       : fitness cost of homozygous deleterious genotype
;   dom-coeff (h)       : dominance; h=0.1 = nearly recessive (most realistic)
;   n-local-loci        : locally adapted loci (controls outbreeding risk)
;   local-adapt-strength: strength of local adaptation (0=none, 1=strong)
;   n-neutral-loci      : neutral loci for diversity tracking
;   mutation-rate       : rate of new deleterious mutations per locus per gamete
;   gene-flow-rate      : fraction of each pop that migrates to other pop per gen
;   birth-rate          : average offspring per female per generation
;   max-generations     : run length
;   n-rescue-migrants   : number of migrants in genetic rescue intervention
;
;   (new, code-defined, no widget required -- see setup for defaults)
;   spatial-mating?     : if true, mates/offspring are spatially local
;   mating-radius       : radius (patches) for local mate search when spatial-mating?
;   migrant-glow-duration : generations a migrant keeps its star shape
;
; HOW TO USE:
;   1. Choose a scenario preset button (recommended), or set sliders manually
;   2. Click SETUP
;   3. Click GO (runs continuously) or STEP (one generation at a time)
;   4. During a run, click GENETIC-RESCUE or FORCED-HYBRIDIZATION to test
;      conservation interventions interactively
;   5. Watch plots: population collapse, fitness decline, allele freq changes
;   6. Watch the view: color = ancestry blend, size = fitness, star = recent migrant
; ============================================================

globals [
  ;; Population sizes
  pop-A-size
  pop-B-size
  total-pop-size

  ;; Fitness metrics
  mean-fitness            ;; grand mean across all individuals
  mean-fitness-A          ;; population A only
  mean-fitness-B          ;; population B only
  mean-inbred-comp        ;; mean W_inbreeding component
  mean-local-comp         ;; mean W_local_adaptation component
  genetic-load            ;; = 1 - mean-fitness

  ;; Genetic diversity metrics
  mean-F                  ;; mean inbreeding coefficient (genome homozygosity proxy)
  mean-het                ;; mean genome-wide heterozygosity = 1 - mean-F

  ;; Allele frequencies
  del-freq-A              ;; deleterious allele frequency in pop A
  del-freq-B              ;; deleterious allele frequency in pop B
  local-opt-freq-A        ;; locally optimal allele freq in pop A (should stay near 1.0)
  local-opt-freq-B        ;; locally optimal allele freq in pop B (should stay near 1.0)

  ;; Between-population differentiation
  fst                     ;; approximate Fst (Nei 1973 analog from heterozygosity)

  ;; [v2] Effective population size (sex-ratio component, Wright 1938)
  Ne-A
  Ne-B

  ;; State
  generation
  any-extinct?            ;; has any population gone extinct?

  ;; [v2] new tunables (plain code globals -- optional widgets, see header)
  spatial-mating?          ;; if true: local mate choice + offspring near mother
  mating-radius            ;; patch radius for local mate search
  migrant-glow-duration    ;; generations a migrant keeps its star shape
]

turtles-own [
  ;; Diploid genome: two haploid chromosomes, each a list of 0/1 alleles
  ;; Total length = n-del-loci + n-local-loci + n-neutral-loci
  chrom1    ;; maternal chromosome
  chrom2    ;; paternal chromosome

  ;; Fitness components (decomposed for analysis)
  fitness         ;; overall individual fitness [0, 1]
  inbred-comp     ;; W_inbreeding: fitness cost from deleterious recessive alleles
  local-comp      ;; W_local: fitness cost from local adaptation mismatch

  ;; Genetic statistics
  F-ind           ;; individual inbreeding coeff proxy: fraction of homozygous loci
  het-ind         ;; = 1 - F-ind
  del-load        ;; total deleterious alleles carried (both chromosomes, all del loci)

  ;; Identity
  pop-id          ;; 0 = population A, 1 = population B (current ecological context)
  sex             ;; 0 = female, 1 = male

  ;; [v2] visual / lineage tracking
  ancestry-A      ;; 0..1, continuous genetic ancestry from pop A founders,
                  ;; inherited as mean of both parents -- independent of pop-id,
                  ;; so a migrant's descendants visibly blend over generations
  migrant-glow    ;; >0 while this individual should show as a migrant (star shape)
]

;; ============================================================
;; REPORTER: total genome length
;; ============================================================

to-report n-total-loci
  report n-del-loci + n-local-loci + n-neutral-loci
end 

;; ============================================================
;; SETUP
;; ============================================================

to setup
  clear-all
  set generation 0
  set any-extinct? false

  ;; [v2] defaults for new (optional-widget) parameters
  set spatial-mating? false      ;; OFF by default -- preserves original dynamics
  set mating-radius 3
  set migrant-glow-duration 6

  ask patches [ set pcolor black ]
  if n-populations = 2 [
    ;; [v2] faint background tint per side, purely cosmetic
    ask patches with [pxcor < 0] [ set pcolor 101 ]
    ask patches with [pxcor > 0] [ set pcolor 11 ]
    ask patches with [pxcor = 0] [ set pcolor 5 ]
  ]

  create-turtles initial-pop-A [
    init-turtle 0
    setxy (min-pxcor + random-float (abs min-pxcor)) random-ycor
  ]

  if n-populations = 2 [
    create-turtles initial-pop-B [
      init-turtle 1
      setxy random-float max-pxcor random-ycor
    ]
  ]

  reset-ticks        ;; ← must come BEFORE update-monitors
  update-monitors    ;; ← now ticks exists and can be read safely
end 

to init-turtle [pid]
  set pop-id pid
  set sex random 2
  set chrom1 new-genome pid
  set chrom2 new-genome pid
  set ancestry-A ifelse-value (pid = 0) [1.0] [0.0]  ;; [v2] founders are "pure"
  set migrant-glow 0
  calc-fitness
  set shape "circle"
  refresh-appearance
end 

;; ============================================================
;; GENOME BUILDING
;; ============================================================

to-report new-genome [pid]
  ;; Builds one haploid chromosome for an individual from population pid
  ;; pid determines which allele is "locally optimal" at local adaptation loci
  let g []

  ;; --- Deleterious loci [indices 0 .. n-del-loci - 1] ---
  ;; Start at approximate mutation-selection equilibrium
  ;; frequency of deleterious allele q ≈ del-allele-freq (set by user)
  repeat n-del-loci [
    set g lput (ifelse-value (random-float 1 < del-allele-freq) [1] [0]) g
  ]

  ;; --- Local adaptation loci [indices n-del-loci .. n-del-loci + n-local-loci - 1] ---
  ;; Pop A optimal allele = 0 (pid = 0)
  ;; Pop B optimal allele = 1 (pid = 1)
  ;; Each population starts with ~90% locally adapted allele (well-adapted)
  let opt-allele pid
  repeat n-local-loci [
    set g lput (ifelse-value (random-float 1 < 0.90) [opt-allele] [1 - opt-allele]) g
  ]

  ;; --- Neutral loci [remaining positions] ---
  ;; Used purely for diversity metrics (Fst, heterozygosity, genetic drift)
  repeat n-neutral-loci [
    set g lput (random 2) g
  ]

  report g
end 

;; ============================================================
;; FITNESS CALCULATION (must be called in turtle context)
;; ============================================================

to calc-fitness
  ;; -------------------------------------------------------
  ;; Component 1: INBREEDING DEPRESSION
  ;; Dominance model (Charlesworth & Willis 2009):
  ;;   Homozygous deleterious [aa]:  W *= (1 - s)
  ;;   Heterozygous [Aa]:            W *= (1 - h*s)
  ;;   Wildtype [AA]:                no effect
  ;; With h << 0.5, deleterious alleles are mostly recessive,
  ;; meaning they are hidden in heterozygotes but exposed in
  ;; inbred homozygotes — the core mechanism of inbreeding depression.
  ;; -------------------------------------------------------
  let w-del 1.0
  let cnt-del 0
  let i 0

  repeat n-del-loci [
    let a1 item i chrom1
    let a2 item i chrom2
    if a1 = 1 [ set cnt-del cnt-del + 1 ]
    if a2 = 1 [ set cnt-del cnt-del + 1 ]
    ifelse (a1 = 1 and a2 = 1) [
      ;; Homozygous deleterious: full selection
      set w-del w-del * (1 - sel-coeff)
    ] [
      if (a1 = 1 or a2 = 1) [
        ;; Heterozygous: partial expression (h = dominance coefficient)
        set w-del w-del * (1 - dom-coeff * sel-coeff)
      ]
    ]
    set i i + 1
  ]
  set inbred-comp w-del
  set del-load cnt-del

  ;; -------------------------------------------------------
  ;; Component 2: LOCAL ADAPTATION / OUTBREEDING DEPRESSION
  ;; Each individual's fitness is reduced proportionally to
  ;; how many of their local adaptation alleles DON'T match
  ;; their current population's environmental optimum.
  ;; This captures the disruption of coadapted gene complexes
  ;; (Edmands 2007, Lynch 1991): when divergent populations
  ;; hybridize, offspring carry allele combinations that were
  ;; never tested by selection together, reducing fitness.
  ;; F2 outbreeding depression (recombination effects) is
  ;; partially captured since alleles segregate independently.
  ;; -------------------------------------------------------
  let w-local 1.0
  if local-adapt-strength > 0 and n-local-loci > 0 [
    let mismatches 0
    set i 0
    repeat n-local-loci [
      let li n-del-loci + i
      let a1 item li chrom1
      let a2 item li chrom2
      ;; Pop A optimum = 0, Pop B optimum = 1
      ;; Mismatched alleles reduce fitness
      if a1 != pop-id [ set mismatches mismatches + 1 ]
      if a2 != pop-id [ set mismatches mismatches + 1 ]
      set i i + 1
    ]
    let max-mm (2 * n-local-loci)
    set w-local max list 0 (1 - (local-adapt-strength * mismatches / max-mm))
  ]
  set local-comp w-local

  ;; -------------------------------------------------------
  ;; Overall fitness: multiplicative model
  ;; -------------------------------------------------------
  let total-fit w-del * w-local
  set fitness max list 0 (min list 1 total-fit)

  ;; -------------------------------------------------------
  ;; Inbreeding coefficient proxy
  ;; True F uses pedigree information; here we use genome-wide
  ;; homozygosity as a proxy (reflects identity-by-descent over time)
  ;; F_hat = proportion of loci where both alleles are identical
  ;; -------------------------------------------------------
  let hom-count 0
  set i 0
  repeat n-total-loci [
    if (item i chrom1 = item i chrom2) [ set hom-count hom-count + 1 ]
    set i i + 1
  ]
  set F-ind hom-count / n-total-loci
  set het-ind 1 - F-ind
end 

;; Turtle reporter: count alleles matching a given optimal allele at local loci

to-report count-local-opt-alleles [opt]
  let cnt 0
  let i 0
  repeat n-local-loci [
    let li n-del-loci + i
    if item li chrom1 = opt [ set cnt cnt + 1 ]
    if item li chrom2 = opt [ set cnt cnt + 1 ]
    set i i + 1
  ]
  report cnt
end 

;; ============================================================
;; [v2] APPEARANCE
;; Centralizes all color/size logic that used to be duplicated
;; across init-turtle, do-gene-flow, reproduce-one-pop, and the
;; intervention buttons. Call this any time fitness or ancestry
;; changes and the turtle should be redrawn.
;; ============================================================

to refresh-appearance
  ;; Color: continuous blue (pop-A ancestry) <-> red (pop-B ancestry) gradient.
  ;; This is independent of current pop-id/location, so migrants keep their
  ;; true ancestry color even after moving, and hybrid offspring visibly
  ;; blend toward purple rather than snapping between two fixed colors.
  let r 255 * (1 - ancestry-A)
  let b 255 * ancestry-A
  set color (list r 40 b)

  ;; Size: scales with fitness, so an inbreeding vortex is visible as
  ;; individuals shrinking, not just the population shrinking in count.
  set size 0.5 + 1.3 * fitness
end 

;; ============================================================
;; MAIN SIMULATION LOOP
;; ============================================================

to go
  if ticks >= max-generations [ stop ]
  if count turtles = 0 [ set any-extinct? true  stop ]

  ;; 1. GENE FLOW between populations (if enabled)
  if n-populations = 2 and gene-flow-rate > 0 [
    do-gene-flow
  ]

  ;; 2. REPRODUCTION (both populations)
  reproduce-all-pops

  ;; [v2] fade migrant "star" markers back to circles over time
  ask turtles with [migrant-glow > 0] [
    set migrant-glow migrant-glow - 1
    if migrant-glow = 0 [ set shape "circle" ]
  ]

  ;; 3. MONITORING
  update-monitors
  tick
end 

to step
  ;; Single-generation advance (useful for close observation)
  go
end 

;; ============================================================
;; GENE FLOW
;; ============================================================

to do-gene-flow
  ;; Symmetric migration: a fraction gene-flow-rate of each pop migrates
  ;; This represents natural dispersal or conservation-managed translocation

  let setA turtles with [pop-id = 0]
  let setB turtles with [pop-id = 1]
  if not any? setA or not any? setB [ stop ]

  let n-AtoB max list 0 (round (count setA * gene-flow-rate))
  let n-BtoA max list 0 (round (count setB * gene-flow-rate))

  ;; Ensure at least 1 individual stays in source population
  set n-AtoB min list n-AtoB (count setA - 1)
  set n-BtoA min list n-BtoA (count setB - 1)

  ;; Migrate A → B
  if n-AtoB > 0 [
    ask n-of n-AtoB setA [
      set pop-id 1
      setxy random-float max-pxcor ycor
      set migrant-glow migrant-glow-duration  ;; [v2]
      set shape "star"                        ;; [v2]
      calc-fitness  ;; fitness recalculated in new environmental context
      refresh-appearance  ;; [v2]
    ]
  ]

  ;; Migrate B → A
  if n-BtoA > 0 [
    ask n-of n-BtoA setB [
      set pop-id 0
      setxy (min-pxcor + random-float (abs min-pxcor)) ycor
      set migrant-glow migrant-glow-duration  ;; [v2]
      set shape "star"                        ;; [v2]
      calc-fitness
      refresh-appearance  ;; [v2]
    ]
  ]
end 

;; ============================================================
;; REPRODUCTION
;; ============================================================

to reproduce-all-pops
  ;; Important: capture agentsets BEFORE any turtles are killed
  ;; These are evaluated after gene-flow has occurred
  let setA turtles with [pop-id = 0]
  let setB turtles with [pop-id = 1]

  if any? setA [ reproduce-one-pop setA 0 carrying-cap-A ]
  if n-populations = 2 and any? setB [
    reproduce-one-pop setB 1 carrying-cap-B
  ]
end 

to reproduce-one-pop [pop-set pid cap]
  ;; --- Check viability: need both sexes ---
  let females pop-set with [sex = 0]
  let males   pop-set with [sex = 1]

  if not any? females or not any? males [
    ;; Population goes extinct (no mates)
    ask pop-set [ die ]
    set any-extinct? true
    stop
  ]

  ;; --- STAGE 1: Generate offspring list ---
  ;; Each female produces birth-rate potential offspring
  ;; Parents are still alive at this point
  let n-fem count females
  let offspring-data []

  ;; [v2 perf fix] build mating lists + fitness sums ONCE per generation
  ;; instead of re-sorting/re-summing on every single offspring draw
  let female-list [self] of females
  let male-list [self] of males
  let fem-total-fit sum [fitness] of females
  let male-total-fit sum [fitness] of males

  repeat round (n-fem * birth-rate) [
    let mom roulette-wheel-select-list female-list fem-total-fit

    ;; [v2] optional spatially-local mate choice
    let dad nobody
    ifelse spatial-mating? and mom != nobody [
      let nearby-males males with [distance mom <= mating-radius]
      ifelse any? nearby-males [
        set dad roulette-wheel-select-list ([self] of nearby-males) (sum [fitness] of nearby-males)
      ] [
        set dad roulette-wheel-select-list male-list male-total-fit
      ]
    ] [
      set dad roulette-wheel-select-list male-list male-total-fit
    ]

    if mom != nobody and dad != nobody [
      ;; Form gametes via independent segregation at all loci
      let g1 apply-mutation (make-gamete ([chrom1] of mom) ([chrom2] of mom))
      let g2 apply-mutation (make-gamete ([chrom1] of dad) ([chrom2] of dad))
      ;; [v2] inherited ancestry = mean of both parents' ancestry
      let anc (([ancestry-A] of mom) + ([ancestry-A] of dad)) / 2
      ;; [v2] remember mother's location for spatial placement
      let mx [xcor] of mom
      let my [ycor] of mom
      set offspring-data lput (list g1 g2 pid anc mx my) offspring-data
    ]
  ]

  ;; --- STAGE 2: Kill parent generation ---
  ;; Non-overlapping generations: parents do not survive reproduction
  ask pop-set [ die ]

  ;; --- STAGE 3: Create offspring and apply viability selection ---
  ;; Each offspring survives with probability = its own fitness
  ;; This is the mechanism through which inbreeding depression and
  ;; outbreeding depression actually reduce population fitness
  foreach offspring-data [ entry ->
    create-turtles 1 [
      set chrom1 item 0 entry
      set chrom2 item 1 entry
      set pop-id  item 2 entry
      set ancestry-A item 3 entry   ;; [v2]
      let mom-x item 4 entry        ;; [v2]
      let mom-y item 5 entry        ;; [v2]
      set sex random 2
      set migrant-glow 0            ;; [v2] offspring are not migrants
      set shape "circle"
      calc-fitness

      ;; Viability selection: survival probability = its own fitness
      ;; Low-fitness inbred or maladapted individuals die here
      ifelse random-float 1 < fitness [
        ;; Survived viability selection
        refresh-appearance  ;; [v2] color by ancestry, size by fitness
        ifelse spatial-mating? [
          ;; [v2] place near mother (jittered), clamped to the correct side
          ;; and within world bounds -- builds persistent spatial clusters
          let jx mom-x + (random-float 2 - 1)
          let jy min list max-pycor (max list min-pycor (mom-y + (random-float 2 - 1)))
          ifelse pop-id = 0 [
            set jx min list -0.1 (max list min-pxcor jx)
          ] [
            set jx max list 0.1 (min list max-pxcor jx)
          ]
          setxy jx jy
        ] [
          ifelse pop-id = 0 [
            setxy (min-pxcor + random-float (abs min-pxcor)) random-ycor
          ] [
            setxy random-float max-pxcor random-ycor
          ]
        ]
      ] [
        die  ;; Failed viability selection (inbreeding or outbreeding depression)
      ]
    ]
  ]

 ;; --- STAGE 4: Density regulation (random culling) ---
;; Random rather than fitness-sorted, so genetic load
;; can accumulate — this allows the inbreeding vortex to emerge
let new-pop turtles with [pop-id = pid]
let n-current count new-pop
if n-current > cap [
  ask n-of (n-current - cap) new-pop [ die ]
]
end 

;; ============================================================
;; GAMETE FORMATION & MUTATION
;; ============================================================

to-report make-gamete [c1 c2]
  ;; Random segregation (independent assortment) at each locus
  ;; Each allele is drawn independently from either parent chromosome
  ;; Note: this effectively treats each locus as on a separate chromosome.
  ;; For studying linked selection or recombination effects, this could
  ;; be extended to model physical chromosomes with recombination.
  let gamete []
  repeat n-total-loci [
    let k length gamete
    set gamete lput (ifelse-value (random 2 = 0) [item k c1] [item k c2]) gamete
  ]
  report gamete
end 

to-report apply-mutation [gamete]
  ;; Forward mutation: wildtype (0) → deleterious (1) at deleterious loci
  ;; Rate: mutation-rate per locus per gamete
  ;; Back mutation (1 → 0) is neglected (rate typically 100x lower)
  ;; This maintains mutation-selection balance and allows studying
  ;; mutational meltdown in very small populations
  let g gamete
  let i 0
  repeat n-del-loci [
    if random-float 1 < mutation-rate [
      set g replace-item i g 1
    ]
    set i i + 1
  ]
  report g
end 

;; ============================================================
;; FITNESS-PROPORTIONATE (ROULETTE WHEEL) SELECTION
;; [v2] Now takes a pre-built list + pre-summed fitness so callers
;; can compute these ONCE per generation instead of once per draw.
;; ============================================================

to-report roulette-wheel-select-list [agent-list total-fit]
  ;; Returns one agent selected with probability proportional to fitness
  ;; This implements sexual selection: fitter individuals mate more
  ;; Combined with viability selection, this gives total selection pressure
  if empty? agent-list [ report nobody ]

  ;; Fallback: if all fitness = 0, pick randomly (rare edge case)
  if total-fit <= 0 [ report one-of agent-list ]

  let threshold random-float total-fit
  let cumul 0.0
  let winner nobody

  foreach agent-list [ a ->
    if winner = nobody [
      set cumul cumul + ([fitness] of a)
      if cumul >= threshold [ set winner a ]
    ]
  ]

  ;; Failsafe (numerical precision edge case)
  if winner = nobody [ set winner last agent-list ]
  report winner
end 

;; ============================================================
;; MONITOR UPDATES
;; ============================================================

to update-monitors
  set generation ticks
  set pop-A-size count turtles with [pop-id = 0]
  set pop-B-size count turtles with [pop-id = 1]
  set total-pop-size count turtles

  if total-pop-size = 0 [
    set mean-fitness 0
    set mean-F 1
    set mean-het 0
    set genetic-load 1
    set mean-fitness-A 0  set del-freq-A 0  set local-opt-freq-A 0  ;; [v2 fix]
    set mean-fitness-B 0  set del-freq-B 0  set local-opt-freq-B 0  ;; [v2 fix]
    set Ne-A 0  set Ne-B 0                                          ;; [v2]
    stop
  ]

  ;; Grand-mean fitness stats
  set mean-fitness mean [fitness] of turtles
  set mean-F mean [F-ind] of turtles
  set mean-het mean [het-ind] of turtles
  set genetic-load 1 - mean-fitness
  set mean-inbred-comp mean [inbred-comp] of turtles
  set mean-local-comp mean [local-comp] of turtles

  ;; Population A specific stats
  ifelse pop-A-size > 0 [
    set mean-fitness-A mean [fitness] of turtles with [pop-id = 0]
    ;; Deleterious allele frequency: (total del alleles) / (total alleles at del loci)
    set del-freq-A (sum [del-load] of turtles with [pop-id = 0]) /
                   (2 * n-del-loci * pop-A-size)
    ;; Local adaptation: fraction of alleles that match pop A's environment (allele 0)
    if n-local-loci > 0 [
      set local-opt-freq-A (sum [count-local-opt-alleles 0] of turtles with [pop-id = 0]) /
                           (2 * n-local-loci * pop-A-size)
    ]
    ;; [v2] Ne from sex-ratio component (Wright 1938): Ne = 4*Nf*Nm / (Nf+Nm)
    let nfA count turtles with [pop-id = 0 and sex = 0]
    let nmA count turtles with [pop-id = 0 and sex = 1]
    set Ne-A ifelse-value (nfA + nmA > 0) [ (4 * nfA * nmA) / (nfA + nmA) ] [0]
  ] [
    set mean-fitness-A 0  set del-freq-A 0  set local-opt-freq-A 0  set Ne-A 0  ;; [v2 fix]
  ]

  ;; Population B specific stats
  ifelse pop-B-size > 0 [
    set mean-fitness-B mean [fitness] of turtles with [pop-id = 1]
    set del-freq-B (sum [del-load] of turtles with [pop-id = 1]) /
                   (2 * n-del-loci * pop-B-size)
    if n-local-loci > 0 [
      set local-opt-freq-B (sum [count-local-opt-alleles 1] of turtles with [pop-id = 1]) /
                           (2 * n-local-loci * pop-B-size)
    ]
    let nfB count turtles with [pop-id = 1 and sex = 0]
    let nmB count turtles with [pop-id = 1 and sex = 1]
    set Ne-B ifelse-value (nfB + nmB > 0) [ (4 * nfB * nmB) / (nfB + nmB) ] [0]
  ] [
    set mean-fitness-B 0  set del-freq-B 0  set local-opt-freq-B 0  set Ne-B 0  ;; [v2 fix]
  ]

  ;; Approximate Fst between populations (Nei 1973 analog)
  ;; Fst = (H_total - H_within) / H_total
  ;; Positive Fst = populations are more different from each other than random
  ;; Fst → 0: panmictic / rescued    Fst → 1: completely diverged
  if n-populations = 2 and pop-A-size > 0 and pop-B-size > 0 [
    let h-A mean [het-ind] of turtles with [pop-id = 0]
    let h-B mean [het-ind] of turtles with [pop-id = 1]
    let h-within (h-A * pop-A-size + h-B * pop-B-size) / total-pop-size
    let h-total mean-het
    set fst ifelse-value (h-total > 0.001) [(h-total - h-within) / h-total] [0]
  ]
end 

;; ============================================================
;; CONSERVATION INTERVENTION BUTTONS
;; ============================================================

to genetic-rescue
  ;; GENETIC RESCUE: Introduce n-rescue-migrants from pop B into pop A
  ;; Simulates conservation translocation / assisted gene flow.
  ;; Expected outcome:
  ;;   - If local-adapt-strength LOW:  rapid fitness recovery (rescue succeeds)
  ;;   - If local-adapt-strength HIGH: initial fitness drop (outbreeding depression)
  ;; This is the core management question modeled here.
  ;;
  ;; [v2] Migrants now MOVE from B to A (matching do-gene-flow) rather than
  ;; being cloned, so pop B genuinely pays a cost for each rescue -- a real
  ;; conservation trade-off instead of an infinite free lunch.
  if n-populations < 2 [
    user-message "Set n-populations = 2 to use genetic rescue!"
    stop
  ]
  if pop-B-size = 0 [
    user-message "Population B is extinct — no source migrants available!"
    stop
  ]
  if pop-A-size = 0 [
    user-message "Population A is extinct — too late for rescue."
    stop
  ]

  let n-m min list n-rescue-migrants pop-B-size
  set n-m min list n-m (pop-B-size - 1)  ;; [v2] keep pop B non-empty when possible

  ask n-of n-m turtles with [pop-id = 1] [
    set pop-id 0
    set sex random 2
    set migrant-glow migrant-glow-duration  ;; [v2]
    set shape "star"                        ;; [v2]
    calc-fitness  ;; recalculate in pop A context (outbreeding depression visible here)
    setxy (min-pxcor + random-float (abs min-pxcor)) random-ycor
    refresh-appearance  ;; [v2]
  ]
  update-monitors
end 

to forced-hybridization
  ;; FORCED HYBRIDIZATION: Flood pop A with 30 migrants from pop B
  ;; Tests extreme gene flow scenario for studying outbreeding depression.
  ;; If local-adapt-strength is high, expect immediate fitness collapse in pop A
  ;;
  ;; [v2] Migrants now MOVE from B to A, same reasoning as genetic-rescue above.
  if n-populations < 2 [
    user-message "Set n-populations = 2!"
    stop
  ]
  if pop-B-size = 0 [ stop ]

  let n-m min list 30 pop-B-size
  set n-m min list n-m (pop-B-size - 1)  ;; [v2] keep pop B non-empty when possible

  ask n-of n-m turtles with [pop-id = 1] [
    set pop-id 0
    set sex random 2
    set migrant-glow migrant-glow-duration  ;; [v2]
    set shape "star"                        ;; [v2]
    calc-fitness
    setxy (min-pxcor + random-float (abs min-pxcor)) random-ycor
    refresh-appearance  ;; [v2]
  ]
  update-monitors
end 

to bottleneck-pop-A
  ;; BOTTLENECK: Suddenly reduce pop A to 10 individuals
  ;; Simulates a catastrophe, habitat loss, or founder event
  ;; Watch F rapidly increase after this event
  if pop-A-size <= 10 [ stop ]
  let n-to-keep 10
  ask min-n-of (pop-A-size - n-to-keep) turtles with [pop-id = 0] [fitness] [
    die
  ]
  update-monitors
end 

;; ============================================================
;; SCENARIO PRESETS
;; ============================================================

to preset-inbreeding-vortex
  ;; Parameters tuned to produce visible inbreeding vortex
  ;; within 200 generations with random density regulation:
  ;;   - Small population (Ne < 20 effective breeders)
  ;;   - High initial genetic load (del-allele-freq 0.20)
  ;;   - Strong selection on homozygotes (s = 0.35)
  ;;   - Very recessive alleles (h = 0.05) so load stays hidden
  ;;     until inbreeding exposes it — maximum surprise effect
  ;;   - Low birth rate so population cannot easily recover
  set n-populations 1
  set initial-pop-A 30
  set initial-pop-B 100
  set carrying-cap-A 30
  set carrying-cap-B 200
  set n-del-loci 30
  set del-allele-freq 0.20
  set sel-coeff 0.35
  set dom-coeff 0.05
  set n-local-loci 0
  set local-adapt-strength 0.0
  set n-neutral-loci 10
  set mutation-rate 0.0002
  set gene-flow-rate 0.0
  set birth-rate 3
  set max-generations 200
  set n-rescue-migrants 5
  setup
end 

to preset-outbreeding-depression
  ;; SCENARIO 2: Two populations with strong local adaptation
  ;; Run 50+ generations to let pops diverge, then click forced-hybridization
  ;; Watch: immediate fitness drop in pop A (outbreeding depression)
  ;; Key insight: how long must populations be isolated before hybridization harms them?
  ;; Try: lower local-adapt-strength and compare hybridization outcome
  set n-populations 2
  set initial-pop-A 100
  set initial-pop-B 100
  set carrying-cap-A 100
  set carrying-cap-B 100
  set n-del-loci 5
  set del-allele-freq 0.03
  set sel-coeff 0.05
  set dom-coeff 0.05
  set n-local-loci 30
  set local-adapt-strength 0.90
  set n-neutral-loci 10
  set mutation-rate 0.000001
  set gene-flow-rate 0.0
  set birth-rate 3
  set max-generations 250
  set n-rescue-migrants 30
  setup
end 

to preset-rescue-succeeds
  ;; SCENARIO 3: Genetic rescue in inbred pop with LOW local adaptation risk
  ;; Run 30-40 gens until inbreeding sets in, then click genetic-rescue
  ;; Watch: fitness rebounds, F drops, del-freq drops as heterosis masks del alleles
  ;; This is the Florida panther / Swedish adder scenario
  set n-populations 2
  set initial-pop-A 20
  set initial-pop-B 200
  set carrying-cap-A 25
  set carrying-cap-B 200
  set n-del-loci 25
  set del-allele-freq 0.15
  set sel-coeff 0.22
  set dom-coeff 0.05
  set n-local-loci 5
  set local-adapt-strength 0.10
  set n-neutral-loci 10
  set mutation-rate 0.0001
  set gene-flow-rate 0.0
  set birth-rate 4
  set max-generations 200
  set n-rescue-migrants 5
  setup
end 

to preset-between-rock-and-hard-place
  ;; SCENARIO 4: The dilemma! Pop A is inbred AND locally adapted
  ;; After ~30 gens, click genetic-rescue. Does rescue help or harm?
  ;; Vary local-adapt-strength slider: where is the tipping point?
  ;; This is Edmands (2007) "Between a rock and a hard place" — the key
  ;; conservation question: when is a species unique enough that outbreeding
  ;; might cause problems vs. when will inbreeding certainly cause extinction?
  set n-populations 2
  set initial-pop-A 25
  set initial-pop-B 150
  set carrying-cap-A 25
  set carrying-cap-B 150
  set n-del-loci 20
  set del-allele-freq 0.12
  set sel-coeff 0.20
  set dom-coeff 0.10
  set n-local-loci 15
  set local-adapt-strength 0.50   ;; <--- ADJUST THIS to explore the dilemma
  set n-neutral-loci 10
  set mutation-rate 0.0001
  set gene-flow-rate 0.0
  set birth-rate 4
  set max-generations 200
  set n-rescue-migrants 5
  setup
end 

to preset-metapopulation
  ;; SCENARIO 5: Two small populations with ongoing gene flow
  ;; Demonstrates optimal gene flow rate: enough to prevent inbreeding
  ;; depression, not so much as to erase local adaptation
  ;; Try: vary gene-flow-rate from 0 to 0.2 and compare outcomes
  ;; Watch Fst: it should decline as gene flow increases
  set n-populations 2
  set initial-pop-A 40
  set initial-pop-B 40
  set carrying-cap-A 50
  set carrying-cap-B 50
  set n-del-loci 20
  set del-allele-freq 0.10
  set sel-coeff 0.15
  set dom-coeff 0.10
  set n-local-loci 10
  set local-adapt-strength 0.40
  set n-neutral-loci 10
  set mutation-rate 0.0001
  set gene-flow-rate 0.05   ;; <--- ADJUST: 0=isolated, 0.2=high connectivity
  set birth-rate 4
  set max-generations 200
  set n-rescue-migrants 5
  setup
end 

to preset-purging-study
  ;; SCENARIO 6: Can purging save an inbred population?
  ;; High s = lethal alleles purge fast; low s = mild alleles purge slowly
  ;; Demonstrates: purging is real but insufficient (Frankham et al. 2001)
  ;; Compare with scenario 1 by changing sel-coeff and dom-coeff
  ;; Strong purging conditions: high s (0.3-0.5), low h (0.02-0.05)
  ;; Weak purging: low s (0.05), moderate h (0.2-0.3)
  set n-populations 1
  set initial-pop-A 25
  set initial-pop-B 100
  set carrying-cap-A 25
  set carrying-cap-B 200
  set n-del-loci 30
  set del-allele-freq 0.12
  set sel-coeff 0.40      ;; <--- High s promotes purging of lethal alleles
  set dom-coeff 0.02      ;; <--- Very low h = nearly fully recessive = better purging
  set n-local-loci 5
  set local-adapt-strength 0.0
  set n-neutral-loci 10
  set mutation-rate 0.00005
  set gene-flow-rate 0.0
  set birth-rate 5
  set max-generations 200
  set n-rescue-migrants 5
  setup
end 

to preset-spatial-family-clusters
  ;; [v2] SCENARIO 7: Same inbreeding vortex setup as scenario 1, but with
  ;; spatial-mating? turned on. Mates are drawn locally and offspring settle
  ;; near their mother, so you should see visible spatial clusters of related
  ;; (similarly-colored, since ancestry is uniform here) individuals emerge,
  ;; and local pockets of relatedness driving faster localized fitness decline
  ;; than the well-mixed version -- isolation-by-distance within a single
  ;; "population" patch, not just between pop A and pop B.
  preset-inbreeding-vortex
  set spatial-mating? true
  set mating-radius 3
end