Purpose: construct an individual based model to simulate how fish move to track temperature-based growth opportunities across a habitat network, the resulting effects on growth, and the implications for long-term population dynamics.

Source dependencies

Load pre-calculated growth

Code 
load("data/wt.growth.array.RData")

6.1 Design notes

Configuration options

  • Stochasticity in movement is accounted for by one of three alternative mechanisms:
    • No stochastic movement/entirely deterministic: all movers do the same thing…they all can perfectly sense and compare growth opportunities
    • Probabilistic (softmax) rule: the softmax function assigns probabilities to each patch based on the growth difference, with higher probabilities for patches offering greater growth benefits.
    • Individual-level movement threshold: each fish has a unique threshold for movement, representing individual variation in boldness/dispersal propensity
  • Density-dependent food availaility: after fish have selected patches, calculate food availability based on patch density and baseline ration levels. Two alternative approaches:
    • Modified Fullerton approach: ration size decreases linearly with fish density, up to some (a priori defined) threshold above which there is no effect of there is no effect of density of ration size.
    • Hyperbolic reduction: describes a non-linear decrease in ration size with increasing fish density, where the rate of decrease slows as density increases.
  • Environment sensing: how omniscient are fish to temperature and density-dependent growth across the habitat patches?
    • abiotic_only: fish can only sense temperature-dependent GP, regardless of which habitat is occupied
    • density_current: fish sense density-dependent GP in current habitat, and compare to abiotic (temp-based) GP in the alternative habitat
    • density_all: fish can sense density-dependent GP in all habitats

Senstive parameters:

  • General
    • Habitat variables (temperature regimes and food availability in each patch)
    • Number of fish per strategy/total number of starting fish
  • Density-dependence
    • Fullerton upper density threshold
    • Hyperbolic half saturation density
  • Movement decision
    • Probabilistic/softmax tau (strength of threshold)
  • Movement cost
    • Allometric function b (metabolic exponent) and c (scaling constant)
  • Movement stochasticity
    • Probabilistic/softmax tau (sensitivity to growth difference)
    • Individual-level movement threshold, SD in individual boldness

Notes:

  • Currently, we just look-up growth from the precalculated growth array, but this is going to be less accurate for the IBM b/c precalculated growth rates are calculated over 0.25g mass increments, whereas growth in the IBM is continuous. Probably fine for now, but perhaps should borrow code from above/Fullerton IBM to increase accuracy. My guess is that this will dramatically increase run times.
  • If sense_environment == "density_all", the density used for the alternative patch is the current pre-movement density (i.e., what the fish can observe before anyone decides to move that day). This is the most tractable approximation — computing post-movement density would require solving a joint decision problem.
  • If using individual-level movement thresholds (i.e., boldness), residents can be given move_threshold = Inf to formally unify all three strategies under one mechanism, rather than branching by strategy label — keeps the code cleaner if we expand to more strategies later.

To do:

  • eventually, should probably collate all parameters into a single control script/page…that way we can see where all the sensitive/arbitrary parameters are
  • stochasticity in movement: it probably makes sense to combine stochastic movement processes…e.g, probabilistic and individual thresholds, and perhaps even random noise as well
  • Size-at-maturity could probably be sampled probabilistically, with noise. Preliminary simulation will just set a specific size threshold.
  • Same for timing of reproduction: rather than all mature fish reproducing on a single day, we could have them spread out over a given period, or make spawn timing size-dependent, with larger fish spawning first. Lots of empirical evidence for this. That would be very interesting because it would give the offspring of the larger/early spawners a competitive advantage over offspring of smaller/later spawners, likely magnifying the effect of size on reproductive success.

6.2 Parameters and initialization

Three movement strategies are simulated: fish that remain year-round in the warm patch (resident_warm), fish that remain year-round in the cold patch (resident_cold), and fish that choose the better-growing patch each day (optimal_mover). Ration sizes match the unequal-ration scenario above.

Configuration notes:

  • All fish start at same size, no variation
  • If relevant, variation in boldness parameter probably needs tuning

Misc. objects

Code 
waterTemps <- cbind("pid" = seq(1, length(seq(0.1, 25, 0.1))),"WT" = seq(0.1, 25, 0.1))
rations <- seq(0.001, 0.4, 0.001)
weights <- seq(0.25, 1500, 0.25)

Configure and initialize fish

Code 
# Configuration options:
move_stochastic <- "prob" # individual movement stochasticity: options = "none", "prob", "indiv"
food_densdepen <- "hyperbolic" # density dependence: options = "fullerton", "hyperbolic
sense_environment <- "density_current" # how omniscient are fish to temperature and density-dependent growth across the habitat patches?
                        #   - "abiotic_only" = Omniscient to abiotic (temp-based) GP only...this is how it currently is. 
                              #   - "density_current" = Sense density dependent GP in current habitat, and compare to abiotic GP in the alternative habitat (perhaps the most biologically realistic)
                              #   - "density_all" = Omniscient to density dependent GP in all habitats

# misc. parameters
n_fish_per_strat  <- 50      # fish per strategy
start_wt          <- 0.5     # initial weight (g)
wt_sd             <- 0.05    # SD of starting weight
MaxDensity4Growth <- 50      # growth is not depressed further above this density (currently this is arbitary)
movecost_c        <- 0.025   # movement cost allometric function: `c` = scaling constant; `
movecost_b        <- 0.6     # movement cost allometric function `a` = decay exponent (higher = steeper drop-off)
sigma_bold        <- 0.002   # SD of boldness distribution, for dispersal propensity (g/g/d units)
tau               <- 0.001   # sensitivity to growth difference (for probabilistic/softmax variation in movement)
T1_mort           <- 30      # temperature (°C) at which daily survival from thermal stress = 0.1
T9_mort           <- 25.8    # temperature (°C) at which daily survival from thermal stress = 0.9
K9_starv          <- 0.55    # relative condition (W_current/W_peak) at which starvation survival = 0.9
K1_starv          <- 0.45    # relative condition at which starvation survival = 0.1
K_warm <- K_cold  <- 500     # For hyperbolic reduction, strength of density dependence/half-saturation density: fish density at which ration is halved

# Fullerton has movement cost coded slightly differently, where movement cost is a function of the instantaneous growth rate, distance moved, and a constant...rather than simply a fixed cost as is the case here. But because we only have two patches that are spatially implicit, we can ignore this more complex parameterization...at least for now.

# ration_warm = 0.10, ration_cold = 0.07 (defined in "Unequal ration" section above)

# Lookup sequences — identical to the axes of the wt.growth array
wt_seq_ibm <- waterTemps[, "WT"]  # seq(0.1, 25, 0.1), 250 values       
ra_seq_ibm <- rations             # seq(0.001, 0.4, 0.001), 151 values

# set.seed(7843)
# initialize fish
fish_pop <- bind_rows(
  tibble(strategy = "resident_warm", patch = "warm",
         weight = pmax(min(weights), rnorm(n_fish_per_strat, start_wt, wt_sd))),
  tibble(strategy = "resident_cold", patch = "cold",
         weight = pmax(min(weights), rnorm(n_fish_per_strat, start_wt, wt_sd))),
  # optimal movers are randomly placed into either warm or cold patches
  tibble(strategy = "optimal_mover", patch = sample(c("warm", "cold"), size = n_fish_per_strat, replace = TRUE), #"warm",
         weight = pmax(min(weights), rnorm(n_fish_per_strat, start_wt, wt_sd)))) %>%
  mutate(pid = row_number(),
         # ration = if_else(patch == "warm", ration_warm, ration_cold),
         # Rations are calculated in the loop as functions of temperature and body size
         cmax_allometric = NA,
         pcmax_baseline = NA,
         pcmax_adjusted = NA,
         pcmax_adjusted_dd = NA,
         func_temp = NA,
         ration = NA,
         # Individual movement threshold — drawn once, fixed for life
         # Positive = reluctant to move (needs a clear advantage)
         # Negative = bold/prone to move (moves even when slightly disadvantaged)
         move_threshold = c(rep(Inf, n_fish_per_strat*2), # residents never move
                            rnorm(n_fish_per_strat, mean = 0, sd = sigma_bold)),
         peak_weight = weight,  # track historical max weight for condition factor
         age_days = 0L,         # age in days, incremented each iteration independent of d
         prob_surv = 1, survive = 1) %>%
  group_by(patch) %>% mutate(density = n()) %>% ungroup()
fish_pop
# A tibble: 150 × 16
   strategy     patch weight   pid cmax_allometric pcmax_baseline pcmax_adjusted
   <chr>        <chr>  <dbl> <int> <lgl>           <lgl>          <lgl>         
 1 resident_wa… warm   0.451     1 NA              NA             NA            
 2 resident_wa… warm   0.448     2 NA              NA             NA            
 3 resident_wa… warm   0.541     3 NA              NA             NA            
 4 resident_wa… warm   0.431     4 NA              NA             NA            
 5 resident_wa… warm   0.435     5 NA              NA             NA            
 6 resident_wa… warm   0.461     6 NA              NA             NA            
 7 resident_wa… warm   0.465     7 NA              NA             NA            
 8 resident_wa… warm   0.539     8 NA              NA             NA            
 9 resident_wa… warm   0.481     9 NA              NA             NA            
10 resident_wa… warm   0.508    10 NA              NA             NA            
# ℹ 140 more rows
# ℹ 9 more variables: pcmax_adjusted_dd <lgl>, func_temp <lgl>, ration <lgl>,
#   move_threshold <dbl>, peak_weight <dbl>, age_days <int>, prob_surv <dbl>,
#   survive <dbl>, density <int>

6.3 Simulation loop

For each day of year, optimal movers compare the bioenergetics growth rate at both patches — the vectorised equivalent of calling fncGrowthPossible() for each fish — and switch only when the movement-cost-penalised gain in the alternative patch exceeds staying (but this may differ among individuals depending on how movement stochasticity is configured). Warm and cold patch residents stay in their pre-defined patch. fncGrowthFish() then looks up the actual growth rate for every fish at their current patch, weight, and ration, and weights are updated accordingly.

Code 
n_days  <- nrow(habitat_df) # number of days in the simulation
n_fish  <- nrow(fish_pop)   # number of fish in the simulation
fish_pop_init <- fish_pop   # snapshot of full population before any mortality

# Pre-allocate storage (rows = fish pid, columns = days)
# Rows are indexed by fish pid (1:n_fish), which is stable even as fish die.
# Dead fish slots remain NA for all days after death.
ggd_matrix         <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # standardized growth rates experienced
gd_matrix          <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # absolute growth rates
temp_matrix        <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # water temps experienced 
patch_matrix       <- matrix(NA_character_, nrow = n_fish, ncol = n_days)  # patch occupied 
mass_matrix        <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # weight/mass trajectories
surv_matrix        <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # survived (1) or died (0) each day
pcmax_base_matrix  <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # baseline P_Cmax threshold 
pcmax_adj_matrix   <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # P_Cmax adjusted for temperature dependence
pcmax_adjdd_matrix <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # P_Cmax adjusted for tempderature and density dependence
func_temp_matrix   <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # Temperature dependence function
ration_matrix      <- matrix(NA_integer_,   nrow = n_fish, ncol = n_days)  # ration size/consumption in g/d/d
condition_matrix   <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # relative condition (W_current / W_peak)
age_matrix         <- matrix(NA_real_,      nrow = n_fish, ncol = n_days)  # age
switches           <- integer(n_fish)                                      # count patch switches per fish (indexed by pid)


for (d in 1:n_days) {
  t <- habitat_df$dayofsim[d] # day of simulation
  fish_pop$age_days <- fish_pop$age_days + 1L  # increment age for all living fish

  # 1. GET/FIND PATCH TEMPERATURE AND RATION -----------------------------------------------------
  
  # Temperature indices for this day — shared by all fish
  T_warm      <- get_patchtemp(t, "warm")             # get warm patch temp on day t
  T_cold      <- get_patchtemp(t, "cold")             # get cold patch temp on day t
  wt_idx_warm <- which.min(abs(wt_seq_ibm - T_warm))  # water temp index for warm patch on day t
  wt_idx_cold <- which.min(abs(wt_seq_ibm - T_cold))  # water temp index for cold patch on day t
  
  # Ration indices for this day - shared by all fish
  # R_warm      <- get_patchration(t, "warm")
  # R_cold      <- get_patchration(t, "cold")
  # ra_idx_warm <- which.min(abs(ra_seq_ibm - R_warm))  # ration index for warm patch on day t
  # ra_idx_cold <- which.min(abs(ra_seq_ibm - R_cold))  # ration index for cold patch on day t
  
  # Calculate habitat, temperature, and body size dependent rations following convo. with Jonny
  fish_pop <- fish_pop %>%
    mutate(func_temp = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)),
           pcmax_baseline = ifelse(patch == "warm", get_patchpcmax(t, "warm"), get_patchpcmax(t, "cold")),
           pcmax_adjusted = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline),
           cmax_allometric = fncAllomCmax(weight)
           # ration = cmax_allometric * pcmax_adjusted
           )

  
  # 2. CHOOSE PATCH: OPTIMAL MOVERS ---------------------------------------------------------------

  # Mirrors fncGrowthPossible() logic, vectorized across all movers at once.
  # Patch choice senses baseline temperature and ration, but not density
  mover_rows <- which(fish_pop$strategy == "optimal_mover")  # get row indices for optimal movers
  
  if (length(mover_rows) > 0) {
    # Mass indices for each mover (clamped to array bounds)
    ma_idx_vec <- pmax(1L, pmin(4500L, round(fish_pop$weight[mover_rows]))) # get row index for current mass
    
    ## 2.A.
    in_warm       <- fish_pop$patch[mover_rows] == "warm"  # is the fish in warm patch? T/F
    move_cost_vec <- fncMoveCost_allometric(fish_pop$weight[mover_rows], c = movecost_c, b = movecost_b) # calculate size-dependent cost of movement 

    # Pre-compute each mover's hypothetical ration in BOTH patches for patch choice sensing.
    # This initial step is done without accounting for the effect of density on p_cmax/ration
    # pcmax scalars are day-level constants; cmax_allometric is already in fish_pop from step 1.
    pcmax_warm_adj  <- min(fncTempDepend(T_warm), get_patchpcmax(t, "warm"))
    pcmax_cold_adj  <- min(fncTempDepend(T_cold), get_patchpcmax(t, "cold"))
    ra_warm_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj
    ra_cold_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj
    ra_idx_warm_vec <- pmax(1L, pmin(400L, map_int(ra_warm_vec, ~which.min(abs(ra_seq_ibm - .x)))))
    ra_idx_cold_vec <- pmax(1L, pmin(400L, map_int(ra_cold_vec, ~which.min(abs(ra_seq_ibm - .x)))))

    if (sense_environment == "abiotic_only") {
      ### 2.A.1. Individuals sense only temperature-based growth potential across all habitats
      g_warm_vec <- wt.growth[cbind(wt_idx_warm, ra_idx_warm_vec, ma_idx_vec)]
      g_cold_vec <- wt.growth[cbind(wt_idx_cold, ra_idx_cold_vec, ma_idx_vec)]
      g_stay     <- ifelse(in_warm, g_warm_vec, g_cold_vec) # growth if staying in current patch
      g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec) # growth if moving to alternate patch

    } else if (sense_environment == "density_current") {
      ### 2.A.2. Individuals sense density-dependent GP in current habitat, abiotic GP in alternative habitat
      # Density-adjusted individual rations for each patch (hyperbolic: Cmax * P_cmax * (K / (K + N)))
      # Adjust P_Cmax by density to effects of body size on the density penalty
      n_warm             <- sum(fish_pop$patch == "warm")
      n_cold             <- sum(fish_pop$patch == "cold")
      ra_warm_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj * (K_warm / (K_warm + n_warm - 1))
      ra_cold_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj * (K_cold / (K_warm + n_cold - 1))
      ra_idx_warm_dd_vec <- pmax(1L, pmin(400L, map_int(ra_warm_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))
      ra_idx_cold_dd_vec <- pmax(1L, pmin(400L, map_int(ra_cold_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))

      # Warm fish sense dd-adjusted warm / abiotic cold; cold fish sense abiotic warm / dd-adjusted cold
      g_warm_vec <- wt.growth[cbind(wt_idx_warm, ifelse(in_warm, ra_idx_warm_dd_vec, ra_idx_warm_vec), ma_idx_vec)]
      g_cold_vec <- wt.growth[cbind(wt_idx_cold, ifelse(in_warm, ra_idx_cold_vec, ra_idx_cold_dd_vec), ma_idx_vec)]
      g_stay     <- ifelse(in_warm, g_warm_vec, g_cold_vec)
      g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec)

    } else if (sense_environment == "density_all") {
      ### 2.A.3. Individuals sense density-dependent GP across all habitats
      # Density-adjusted individual rations for each patch (hyperbolic: Cmax * P_cmax * (K / (K + N)))
      # Adjust P_Cmax by density to effects of body size on the density penalty
      n_warm             <- sum(fish_pop$patch == "warm")
      n_cold             <- sum(fish_pop$patch == "cold")
      ra_warm_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_warm_adj * (K_warm / (K_warm + n_warm - 1))
      ra_cold_dd_vec     <- fish_pop$cmax_allometric[mover_rows] * pcmax_cold_adj * (K_cold / (K_warm + n_cold - 1))
      ra_idx_warm_dd_vec <- pmax(1L, pmin(400L, map_int(ra_warm_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))
      ra_idx_cold_dd_vec <- pmax(1L, pmin(400L, map_int(ra_cold_dd_vec, ~which.min(abs(ra_seq_ibm - .x)))))

      # Both patches sensed with density-adjusted individual rations
      g_warm_vec <- wt.growth[cbind(wt_idx_warm, ra_idx_warm_dd_vec, ma_idx_vec)]
      g_cold_vec <- wt.growth[cbind(wt_idx_cold, ra_idx_cold_dd_vec, ma_idx_vec)]
      g_stay     <- ifelse(in_warm, g_warm_vec, g_cold_vec)
      g_move_net <- ifelse(in_warm, g_cold_vec - move_cost_vec, g_warm_vec - move_cost_vec)
    }
    
    prev_patch <- fish_pop$patch[mover_rows]
    
    ## 2.B. Impose stochasticity in movement
    if (move_stochastic == "none") {
      ### 2.B.1. Deterministic: All movers do the same thing
      # Update patch based on whether moving is beneficial. 
      # 1. If fish in is warm patch and growth potential of moving exceeds growth potential of staying, move to cold. 
      # 2. If fish in is cold patch and growth potential of moving exceeds growth potential of staying, move to warm. 
      fish_pop$patch[mover_rows] <- case_when(
        in_warm  & g_move_net > g_stay ~ "cold",
        !in_warm & g_move_net > g_stay ~ "warm", 
        TRUE                           ~ fish_pop$patch[mover_rows]
        )
    } else if (move_stochastic == "prob") {
      ### 2.B.2. Probabilistic (softmax) rule: Instead of a hard threshold, fish switch with a probability that increases with the growth advantage. One parameter: τ (sensitivity; small = nearly deterministic, large = nearly random).
      # Use cost-adjusted net growth rates: staying is free, moving incurs move_cost_vec.
      # g_stay     = growth in current patch (no cost)
      # g_move_net = growth in alternate patch - move_cost_vec
      # Re-map to warm/cold perspective for fncMoveSoftmax:
      gwarm_eff <- ifelse(in_warm, g_stay, g_move_net) # effective growth if occupying warm patch
      gcold_eff <- ifelse(in_warm, g_move_net, g_stay) # effective growth if occupying cold patch
      p_warm <- fncMoveSoftmax(gwarm = gwarm_eff, gcold = gcold_eff, tau = tau) # calculate probability of selecting warm habitat
      choose_warm <- runif(length(mover_rows)) < p_warm # choose warm if p_warm exceeds a random number between 0 and 1
      fish_pop$patch[mover_rows] <- if_else(choose_warm, "warm", "cold") # update patch
    } else if (move_stochastic == "indiv") {
      ### 2.B.3. Individual-level movement threshold: each fish has a unique threshold for movement, representing individual variation in boldness/dispersal propensity
      fish_pop$patch[mover_rows] <- case_when(
      in_warm  & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "cold",
      !in_warm & (g_move_net - g_stay) > fish_pop$move_threshold[mover_rows] ~ "warm",
      TRUE                                                                   ~ prev_patch) #fish_pop$patch[mover_rows])
    }
    
    # Tally switches (indexed by pid so counts remain correct as fish die)
    switched <- fish_pop$patch[mover_rows] != prev_patch
    switches[fish_pop$pid[mover_rows]] <- switches[fish_pop$pid[mover_rows]] + as.integer(switched)

  } # end patch choice

  # Recompute func_temp and pcmax_adjusted based on post-choice patch.
  # Required because some fish may have switched patches during step 2;
  # step 1 values reflect the pre-choice patch and would be stale for switchers.
  fish_pop <- fish_pop %>%
    mutate(
      func_temp      = ifelse(patch == "warm", fncTempDepend(T_warm), fncTempDepend(T_cold)),
      pcmax_baseline = ifelse(patch == "warm", get_patchpcmax(t, "warm"), get_patchpcmax(t, "cold")),
      pcmax_adjusted = ifelse(func_temp < pcmax_baseline, func_temp, pcmax_baseline)
    )

  # 3. UPDATE HABITAT QUALITY (DENSITY-DEPENDENT RATION) -----------------------------------------

  # Update density after fish have moved
  fish_pop <- fish_pop %>% group_by(patch) %>% mutate(density = n()) %>% ungroup()
  
  if (food_densdepen == "fullerton") {
    ### 3a. Fullerton approach
    # get density effect on ration
    fdens <- fdens_raw <- fish_pop$density # get fish "density" per patch (here, simply abundance b/c space is not specified)
    fdens[fdens > MaxDensity4Growth] <- MaxDensity4Growth # high densities cap out at MaxDensity4Growth
    density.effect <- fncRescale((1 - c(fdens, 0.01, MaxDensity4Growth)), c((0.5),1)) # calculate density effect scalar
    density.effect <- density.effect[-c(length(density.effect), (length(density.effect) - 1))] # remove the last 2 temporary values
    # update ration in patches to account for density dependence in food availability
    fish_pop <- fish_pop %>% mutate(ration = ration * density.effect)
  } else if (food_densdepen == "hyperbolic") {
    ### 3b. Hyperbolic reduction approach — scale each fish's individual P_Cmax by density, then calculate ration as C_max * P_Cmax 
    fish_pop <- fish_pop %>%
      mutate(
        k                 = if_else(patch == "warm", K_warm, K_cold),
        pcmax_adjusted_dd = pcmax_adjusted * (k / (k + density - 1)),
        ration            = cmax_allometric * pcmax_adjusted_dd
      ) %>%
      select(-k)
  }
  
  
  # 4. GROW FISH (WISCONSIN BIOENERGETICS) -------------------------------------------------------
  
  # Growth lookup for all fish via fncGrowthFish. Growth in g/g/d
  growth_df <- fncGrowthFish(NA, as.data.frame(fish_pop), t)
  
  # Store experienced rates before weight update.
  # All indexing uses fish pid (rows 1:n_fish) so slots remain stable as fish die.
  temp_matrix[fish_pop$pid, d]  <- growth_df$WT.actual
  ggd_matrix[fish_pop$pid, d]   <- growth_df$growth                     # growth in g/g/d
  gd_matrix[fish_pop$pid, d]    <- growth_df$growth * fish_pop$weight   # growth in g/d
  
  # Update weight: W_{t+1} = W_t + (growth_rate * W_t)
  fish_pop$weight <- fish_pop$weight + (growth_df$growth * fish_pop$weight)
  
  # Update peak weight: ratchet up only, never down
  fish_pop$peak_weight <- pmax(fish_pop$peak_weight, fish_pop$weight)
  
  
  # 5. SURVIVAL BASED ON FISH WEIGHT AND GROWTH --------------------------------------------------
  
  # Weight and growth dependent -- growth effect depends on range of growth rates experienced on any given day
  # surv.list <- fncSurvive(growth_df)
  # growth_df <- growth_df %>%
  #   mutate(prob_surv = surv.list[[1]],
  #          survive   = surv.list[[2]])
  
  # Weight and growth dependent -- growth effect is rescaled to a specific range
  p_sg    <- fncSurviveSimp(growth_df)[[1]] # fncSurviveSimp(growth_df)[[1]]
  
  # Temperature-dependent survival: multiplied with size/growth probability so both act simultaneously
  p_temp  <- fncSurviveTemp(growth_df$WT.actual, T1 = T1_mort, T9 = T9_mort)

  # Condition-based (starvation) survival: relative condition = current / peak weight
  condition <- fish_pop$weight / fish_pop$peak_weight
  p_starv   <- fncSurviveStarve(condition, K9 = K9_starv, K1 = K1_starv)
  
  # Age-based survival (senescence): survival probability declines past 5 years old, fish older than 8 have ~0 chance of survival.
  p_age <- fncSurviveAge(fish_pop$age_days / 365)
  
  # Calculate combined daily survival rate
  prb.srv   <- pmin(p_sg * p_temp * p_starv * p_age, 1)
  
  # Minimum weight floor: fish that drop below hatch weight die immediately (backstop)
  prb.srv[fish_pop$weight < start_wt] <- 0
  
  survivors <- rbinom(n = nrow(growth_df), size = 1, prob = prb.srv)
  
  growth_df <- growth_df %>%
    mutate(prob_surv = prb.srv,
           survive   = survivors)
  
  # Weight dependent only --- the issue with this is that basically nobody dies. See GrowthScapes working notes in OneNote
  # surv.list <- fncSurviveSimp(growth_df)
  # growth_df <- growth_df %>% 
  #   mutate(prob_surv = surv.list[[1]], 
  #          survive   = surv.list[[2]])
  
  
  # 6. STORE RESULTS AND REMOVE NON-SURVIVORS ----------------------------------------------------
  
  # Primary: Store updated weight, patch, and survival outcome for this day
  mass_matrix[fish_pop$pid, d]  <- fish_pop$weight
  patch_matrix[fish_pop$pid, d] <- fish_pop$patch
  surv_matrix[fish_pop$pid, d]  <- growth_df$survive
  age_matrix[fish_pop$pid, d]  <- fish_pop$age_days/365 # store age in years
  
  # Secondary: store bioenergetics inputs
  pcmax_base_matrix[fish_pop$pid, d]  <- fish_pop$pcmax_baseline
  pcmax_adj_matrix[fish_pop$pid, d]   <- fish_pop$pcmax_adjusted
  pcmax_adjdd_matrix[fish_pop$pid, d] <- fish_pop$pcmax_adjusted_dd 
  func_temp_matrix[fish_pop$pid, d]   <- fish_pop$func_temp
  ration_matrix[fish_pop$pid, d]      <- fish_pop$ration
  condition_matrix[fish_pop$pid, d]   <- condition
  
  # Remove non-survivors — only living fish carry forward to next iteration
  fish_pop <- fish_pop[growth_df$survive == 1, ]
  
  # exit loop if all fish have died
  if (nrow(fish_pop) == 0) break
}

6.4 Summarize results

Collate results into long tibbles and summarized tibbles

Code 
# assign DOY column names
colnames(mass_matrix)     <- as.character(habitat_df$dayofsim)
colnames(patch_matrix)    <- as.character(habitat_df$dayofsim)
colnames(temp_matrix)     <- as.character(habitat_df$dayofsim)
colnames(ggd_matrix)      <- as.character(habitat_df$dayofsim)
colnames(surv_matrix)     <- as.character(habitat_df$dayofsim)
colnames(age_matrix)      <- as.character(habitat_df$dayofsim)

colnames(pcmax_base_matrix)    <- as.character(habitat_df$dayofsim)
colnames(pcmax_adj_matrix)     <- as.character(habitat_df$dayofsim)
colnames(pcmax_adjdd_matrix)   <- as.character(habitat_df$dayofsim)
colnames(func_temp_matrix)     <- as.character(habitat_df$dayofsim)
colnames(ration_matrix)        <- as.character(habitat_df$dayofsim)
colnames(condition_matrix)     <- as.character(habitat_df$dayofsim)


# Convert weight matrix to long format, join with patch matrix
# Use fish_pop_init to label rows — fish_pop at this point only contains survivors
mass_long <- as_tibble(mass_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "weight") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(
    as_tibble(patch_matrix) %>%
      bind_cols(fish_pop_init %>% select(pid)) %>%
      pivot_longer(-pid, names_to = "dayofsim", values_to = "patch") %>%
      mutate(dayofsim = as.integer(dayofsim)),
    by = c("pid", "dayofsim")
  ) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# Convert temp matrix to long format
temp_long <- as_tibble(temp_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "temp") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert standardized growth rate matrix to long format
ggd_long <- as_tibble(ggd_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "ggd") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert survival matrix to long format (1 = survived, 0 = died, NA = already dead)
surv_long <- as_tibble(surv_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "survived") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert baseline P_Cmax matrix to long format
pcmax_base_long <- as_tibble(pcmax_base_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "pcmax_base") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert temp-adjusted P_Cmax matrix to long format
pcmax_adj_long <- as_tibble(pcmax_adj_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "pcmax_adj") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert temp and density adjusted matrix to long format
pcmax_adjdd_long <- as_tibble(pcmax_adjdd_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "pcmax_adjdd") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert survival matrix to long format (1 = survived, 0 = died, NA = already dead)
func_temp_long <- as_tibble(func_temp_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "func_temp") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# convert survival matrix to long format (1 = survived, 0 = died, NA = already dead)
ration_long <- as_tibble(ration_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "ration") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

condition_long <- as_tibble(condition_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "condition") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

age_long <- as_tibble(age_matrix) %>%
  bind_cols(fish_pop_init %>% select(pid, strategy)) %>%
  pivot_longer(-c(pid, strategy), names_to = "dayofsim", values_to = "age") %>%
  mutate(dayofsim = as.integer(dayofsim)) %>%
  left_join(habitat_df %>% select(dayofsim, date), by = "dayofsim")

# bind all individual data
ibm_long <- mass_long %>%
  left_join(temp_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(ggd_long %>% select(-date),  by = c("pid", "strategy", "dayofsim")) %>%
  left_join(surv_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(pcmax_base_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(pcmax_adj_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(pcmax_adjdd_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(func_temp_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(ration_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(condition_long %>% select(-date), by = c("pid", "strategy", "dayofsim")) %>%
  left_join(age_long %>% select(-date), by = c("pid", "strategy", "dayofsim"))

# summarize by strategy
ibm_summary <- ibm_long %>%
  group_by(strategy, dayofsim, date) %>%
  summarise(
    mean_weight      = mean(weight, na.rm = TRUE),
    sd_weight        = sd(weight,   na.rm = TRUE),
    mean_temp        = mean(temp,   na.rm = TRUE),
    sd_temp          = sd(temp,     na.rm = TRUE),
    mean_ggd         = mean(ggd,    na.rm = TRUE),
    sd_ggd           = sd(ggd,      na.rm = TRUE),
    prop_warm        = mean(patch == "warm", na.rm = TRUE),
    n_alive          = sum(survived, na.rm = TRUE),
    mean_pcmax_bas   = mean(pcmax_base, na.rm = TRUE),
    sd_pcmax_base    = sd(pcmax_base, na.rm = TRUE),
    mean_pcmax_adj   = mean(pcmax_adj, na.rm = TRUE),
    sd_pcmax_adj     = sd(pcmax_adj, na.rm = TRUE),
    mean_pcmax_adjdd = mean(pcmax_adjdd, na.rm = TRUE),
    sd_pcmax_adjdd   = sd(pcmax_adjdd, na.rm = TRUE),
    mean_func_temp   = mean(func_temp, na.rm = TRUE),
    sd_func_temp     = sd(func_temp, na.rm = TRUE),
    mean_ration      = mean(ration, na.rm = TRUE),
    sd_ration        = sd(ration, na.rm = TRUE),
    mean_condition   = mean(condition, na.rm = TRUE),
    sd_condition     = sd(condition, na.rm = TRUE),
    # mean_age         = mean(age, na.rm = TRUE),
    # sd_age           = sd(age, na.rm = TRUE),
    .groups          = "drop"
  )

ibm_summary
# A tibble: 11,691 × 23
   strategy dayofsim date       mean_weight sd_weight mean_temp sd_temp mean_ggd
   <chr>       <int> <date>           <dbl>     <dbl>     <dbl>   <dbl>    <dbl>
 1 optimal…        1 2024-05-01       0.534    0.0533      13.0    3.28   0.0728
 2 optimal…        2 2024-05-02       0.597    0.0366      13.3    3.34   0.0740
 3 optimal…        3 2024-05-03       0.641    0.0391      13.5    3.34   0.0726
 4 optimal…        4 2024-05-04       0.685    0.0408      13.4    3.40   0.0707
 5 optimal…        5 2024-05-05       0.733    0.0441      13.6    3.50   0.0689
 6 optimal…        6 2024-05-06       0.783    0.0481      13.8    3.50   0.0674
 7 optimal…        7 2024-05-07       0.834    0.0519      13.9    3.55   0.0659
 8 optimal…        8 2024-05-08       0.890    0.0560      14.2    3.60   0.0646
 9 optimal…        9 2024-05-09       0.946    0.0598      14.3    3.65   0.0631
10 optimal…       10 2024-05-10       1.00     0.0633      14.5    3.65   0.0618
# ℹ 11,681 more rows
# ℹ 15 more variables: sd_ggd <dbl>, prop_warm <dbl>, n_alive <int>,
#   mean_pcmax_bas <dbl>, sd_pcmax_base <dbl>, mean_pcmax_adj <dbl>,
#   sd_pcmax_adj <dbl>, mean_pcmax_adjdd <dbl>, sd_pcmax_adjdd <dbl>,
#   mean_func_temp <dbl>, sd_func_temp <dbl>, mean_ration <dbl>,
#   sd_ration <dbl>, mean_condition <dbl>, sd_condition <dbl>

6.5 Check bioenergetics functionality

Ration size by weight: points should generally track allometric scaling of Cmax (halved to account for baseline/treshold P_Cmax, dashed line). Noise is due to variation in temperature and density. This looks fine…points fall away from 0.5xC_max during very cold periods when temperature restricts physiological capacity. The three separate “cliffs” show that you get the same behavior at the end of the year when temps are limited but fish have attained different mean body sizes depending on their life-history strategy

Code 
constants <- fncReadConstants() 
xs <- seq(from = min(ibm_long$weight, na.rm = TRUE), to = max(ibm_long$weight, na.rm = TRUE), length.out = 100)
ys <- constants$Consumption$CA * (xs ^ constants$Consumption$CB)

p1 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = weight, y = ration, color = strategy)) +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "black") +
  theme_bw() + theme(legend.position = "top") +
  ylab("Ration size (g/g/d)") + xlab("Fish mass (g)")

p2 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = weight, y = ration, color = temp)) +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "black") +
  theme_bw() + theme(legend.position = "top") +
  ylab("Ration size (g/g/d)") + xlab("Fish mass (g)")

egg::ggarrange(p2, p1, nrow = 1)

f(T) should perfectly track the temperature dependence function b/c it is only based on temperature and constants that are common to all individuals

Code 
ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1), aes(x = temp, y = func_temp, color = weight)) +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + 
  ylab("f(Temperature)") + xlab("Temperature (deg. C)")

P_Cmax by temp: points should perfect track the temperature dependence function, but max out at the threshold P_Cmax set in each habitat.

Code 
ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1), aes(x = temp, y = pcmax_adj, color = weight)) +
  geom_abline(slope = 0, intercept = 0.5, color = "grey40", linetype = "dashed") +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + 
  ylab("Temperature adjusted P_Cmax") + xlab("Temperature (deg. C)")

Now show P_cmax adjusted by both temperature and density: points should generally track temperature-dependence function, capped at baseline/threshold P_Cmax. Noise is due to variation in body size and density. Smaller fish feeding at higher densities typically experience more food limitation…this appears to occur more commonly in the warm habitat during late spring.

Code 
p1 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = temp, y = pcmax_adjdd, color = weight)) +
  geom_abline(slope = 0, intercept = 0.5, color = "grey40", linetype = "dashed") +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + ylim(0.07, 0.5) + theme(legend.position = "top") +
  ylab("Temperature and density adjusted P_Cmax") + xlab("Temperature (deg. C)")

p2 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = temp, y = pcmax_adjdd, color = density)) +
  geom_abline(slope = 0, intercept = 0.5, color = "grey40", linetype = "dashed") +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + ylim(0.07, 0.5) + theme(legend.position = "top") +
  ylab("Temperature and density adjusted P_Cmax") + xlab("Temperature (deg. C)")

p3 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = temp, y = pcmax_adjdd, color = patch)) +
  geom_abline(slope = 0, intercept = 0.5, color = "grey40", linetype = "dashed") +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + ylim(0.07, 0.5) + theme(legend.position = "top") +
  ylab("Temperature and density adjusted P_Cmax") + xlab("Temperature (deg. C)")

p4 <- ggplot() +
  geom_point(data = ibm_long %>% filter(survived == 1) %>% group_by(dayofsim, patch) %>% mutate(density = n()) %>% ungroup(), 
             aes(x = temp, y = pcmax_adjdd, color = dayofsim)) +
  geom_abline(slope = 0, intercept = 0.5, color = "grey40", linetype = "dashed") +
  # geom_line(aes(x = xs, y = ys), linetype = "dashed", color = "grey40") +
  # geom_line(aes(x = xs, y = ys*0.5), linetype = "dashed", color = "red") +
  theme_bw() + ylim(0.07, 0.5) + theme(legend.position = "top") +
  ylab("Temperature and density adjusted P_Cmax") + xlab("Temperature (deg. C)")

egg::ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2)

6.6 Visualize results

6.6.1 Size trajectories

Weight trajectories by strategy (ribbon = ±1 SD across individuals).

I wonder how realistic overwinter mass losses are…

Code 
p1 <- ibm_summary %>%
  ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_ribbon(aes(x = date, color = strategy, fill = strategy, ymin = mean_weight - sd_weight, ymax = mean_weight + sd_weight),
              alpha = 0.15, color = NA) +
  geom_line(aes(x = date, color = strategy, fill = strategy, y = mean_weight)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_fill_manual(values  = c(resident_warm = "red", 
                                resident_cold = "blue",
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Mean weight ± SD (g)")

p2 <- mass_long %>%
  ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = weight, color = strategy, group = pid), alpha = 0.5) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Individual weight (g)")

ggpubr::ggarrange(p1, p2, nrow = 1, common.legend = TRUE)

Now zoom into the beginning of the simulation, to look growth trajectories of early recruits

Code 
p1 <- ibm_summary %>%
  filter(date <= min(date) + days(100)) %>%
  ggplot() +
  geom_ribbon(aes(x = date, color = strategy, fill = strategy, ymin = mean_weight - sd_weight, ymax = mean_weight + sd_weight),
              alpha = 0.15, color = NA) +
  geom_line(aes(x = date, color = strategy, fill = strategy, y = mean_weight)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_fill_manual(values  = c(resident_warm = "red", 
                                resident_cold = "blue",
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Mean weight ± SD (g)") +
  geom_abline(slope = 0, intercept = 0, linetype = "dashed")

p2 <- mass_long %>%
  filter(date <= min(date) + days(100)) %>%
  ggplot() +
  geom_line(aes(x = date, y = weight, color = strategy, group = pid), alpha = 0.5) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Individual weight (g)") +
  geom_abline(slope = 0, intercept = 0, linetype = "dashed")

ggpubr::ggarrange(p1, p2, nrow = 1, common.legend = TRUE)

Now just look at warm residents

Code 
p1 <- ibm_summary %>%
  filter(strategy == "resident_warm", !is.na(mean_weight)) %>%
  ggplot() +
  # geom_rect(data = summer_shading,
  #           aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
  #           fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_ribbon(aes(x = date, color = strategy, fill = strategy, ymin = mean_weight - sd_weight, ymax = mean_weight + sd_weight),
              alpha = 0.15, color = NA) +
  geom_line(aes(x = date, color = strategy, fill = strategy, y = mean_weight)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_fill_manual(values  = c(resident_warm = "red", 
                                resident_cold = "blue",
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Mean weight ± SD (g)")

p2 <- mass_long %>%
  filter(strategy == "resident_warm", !is.na(weight)) %>%
  ggplot() +
  # geom_rect(data = summer_shading,
  #           aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
  #           fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = weight, color = strategy, group = pid), alpha = 0.5) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Individual weight (g)")

ggpubr::ggarrange(p1, p2, nrow = 1, common.legend = TRUE)

6.6.2 Survival

6.6.2.1 Cumulative survival

Plot daily counts of survivors and overall survival by strategy (plots are the same, survival is standardized by number of initial fish per strategy)

Code 
ibm_summary %>% 
  ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = n_alive / n_fish_per_strat, color = strategy)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    x     = "Date",
    y     = "Proportion alive",
    color = "Strategy",
    title = "Cumulative survival by strategy"
  ) +
  theme_bw() + theme(panel.grid = element_blank())

Code 
# ggpubr::ggarrange(p1, p2, common.legend = TRUE)

Summarize survival outcomes by strategy:

Code 
ibm_summary %>% 
  filter(dayofsim == n_days) %>%
  mutate(survive_rate = n_alive / n_fish_per_strat) %>%
  select(strategy, n_alive, survive_rate)
# A tibble: 3 × 3
  strategy      n_alive survive_rate
  <chr>           <int>        <dbl>
1 optimal_mover       0            0
2 resident_cold       0            0
3 resident_warm       0            0

Interpretation below is based on past model runs with different parameterization/configuration. But what is consistent across runs is that (1) summer is a challenging time for survival in the warm habitat, (2) winter/spring is challenges survival in the cold habitat, and (3) movers always have the highest survival.

sense_environment appears to have a large effect on survival outcomes, but doesn’t change overall patterns described above.

  • resident_warm has the lowest overall survival, consistent with it being the poorest-growing strategy.
    • mid-year acceleration in resident_warm mortality coincides with slower/negative growth in mid-summer, which is exactly what the function should be responding to.
    • similarly, late-year acceleration in mortality coincides with negative growth at this time
  • resident_cold has moderate overall survival. Survival declines the early part of the year when growth is negative and fish are small. Very little mortality occurs during the summer, consistent with theory on the role of perennially cold habitat providing summer refugia. By the end of the year, most fish have escaped the survival bottleneck associated with small body size
  • optimal_mover has the highest overall survival (~96%), reflecting its growth advantage keeping individual survival probabilities near 1.0 most of the year.

6.6.2.2 Context

When are fish dying and how large are they/how fast are they growing relative to the survivor mean?

Code 
n_days <- ncol(mass_matrix)

# Per-fish fate summary (used for panel C)
fish_fate <- fish_pop_init %>%
  mutate(
    survived   = if_else(!is.na(mass_matrix[pid, n_days]), "Survived", "Died"),
    days_alive = rowSums(!is.na(mass_matrix))
  )

# fish_fate %>% count(strategy, survived)
Code 
strategy_order <- c("resident_warm", "resident_cold", "optimal_mover")

# Panels A & B: compare dying fish to surviving fish on the same day
# -----------------------------------------------------------------------
# Dying fish: weight and growth rate on their day of death (survived == 0)
death_records <- ibm_long %>%
  filter(survived == 0) %>%
  select(pid, strategy, date, dayofsim, weight, ggd) %>%
  mutate(fate = "Died")

# Death days across all strategies
death_days <- death_records %>% distinct(dayofsim)

# Surviving fish (those alive at end): weight and growth on those same death days
survivor_pids <- fish_fate %>% filter(survived == "Survived") %>% pull(pid)

survivor_records <- ibm_long %>%
  filter(pid %in% survivor_pids) %>%
  inner_join(death_days, by = "dayofsim") %>%
  select(pid, strategy, dayofsim, weight, ggd) %>%
  mutate(fate = "Survived")

# Mean of ALL survivors (across strategies) on each death day
survivor_means <- survivor_records %>%
  group_by(dayofsim) %>%
  summarise(mean_wt_surv  = mean(weight, na.rm = TRUE),
            mean_ggd_surv = mean(ggd,    na.rm = TRUE),
            .groups = "drop")

# Difference: dying fish value minus all-survivor mean on same day
df_diff <- death_records %>%
  left_join(survivor_means, by = "dayofsim") %>%
  mutate(
    wt_diff  = weight - mean_wt_surv,
    ggd_diff = ggd    - mean_ggd_surv,
    strategy = factor(strategy, levels = strategy_order)
  )

# -----------------------------------------------------------------------
p_wt <- ggplot(df_diff, aes(x = strategy, y = wt_diff)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_boxplot(width = 0.4, outlier.shape = NA, fill = "grey") +
  geom_jitter(width = 0.15, size = 2, alpha = 0.7) +
  labs(x = NULL, y = "Weight difference (g)",
       title = "Dying fish wt. minus mean survivor wt. (same day)") +
  theme_bw() + theme(panel.grid = element_blank())

p_ggd <- ggplot(df_diff, aes(x = strategy, y = ggd_diff)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_boxplot(width = 0.4, outlier.shape = NA, fill = "grey") +
  geom_jitter(width = 0.15, size = 2, alpha = 0.7) +
  labs(x = NULL, y = "Growth rate difference (g/g/d)",
       title = "Dying fish growth minus mean survivor growth (same day)") +
  theme_bw() + theme(panel.grid = element_blank())

p_days <- ggplot() +
  geom_rect(data = summer_shading,
            aes(ymin = xmin, ymax = xmax, xmin = ymin, xmax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_jitter(data = death_records,
              aes(x = strategy, y = date, color = strategy), 
              width = 0.15, size = 2.5, alpha = 0.8) +
  scale_x_discrete(limits = strategy_order) +
  labs(x = NULL, y = "Day of death (day of year)",
       title = "When did non-survivors die?") +
  theme_bw() + theme(legend.position = "none", panel.grid = element_blank()) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen"))

egg::ggarrange(p_wt, p_ggd, p_days, ncol = 1)

Yes, overall, the survival module is behaving correctly: it’s largely selecting against small, slow-growing fish in a biologically interpretable way.

resident_warm fish that die in summer are smaller and actively losing mass relative to the population-wide survivor mean on those same days. The tight clustering near zero for resident_cold and optimal_mover reflects that deaths in those strategies tend to occur early in the year when fish are small and the all-strategy survivor mean hasn’t yet diverged much.

Timing of deaths (panel C) also makes biological sense: + resident_cold deaths cluster early when fish are small and baseline survival is lowest, and fish are also experiencing negative growth due to low temps + optimal_mover deaths also cluster early for the same reasons + resident_warm deaths cluster mid-summer as the warm patch becomes suboptimal. Likely not as much early season mortality b/c growth potential is just a bit higher

Now, let’s investigate the context of optimal_mover mortalities. Where and when to they occur and what are the thermal and growth conditions that they experienced?

Code 
mover_deaths <- ibm_long %>%
  filter(strategy == "optimal_mover", survived == 0) %>%
  select(pid, date, dayofsim, patch, weight, ggd)

dayofsim_max <- max(mover_deaths$dayofsim) + 10

temps <- habitat_df %>%
  filter(dayofsim <= dayofsim_max) %>%
  select(dayofsim, date, temp_warm, temp_cold) %>%
  pivot_longer(c(temp_warm, temp_cold), names_to = "patch", values_to = "temp") %>%
  mutate(patch = if_else(patch == "temp_warm", "warm", "cold"))

p_temp <- ggplot(temps, aes(x = date, y = temp, color = patch)) +
  geom_line(linewidth = 0.8) +
  geom_rug(
    data = mover_deaths, aes(x = date, color = patch),
    sides = "b", linewidth = 1.2, length = unit(0.06, "npc"), inherit.aes = FALSE
  ) +
  scale_color_manual(values = c(warm = "tomato", cold = "steelblue")) +
  labs(x = NULL, y = "Temperature (°C)", color = "Patch",
       title = "Optimal mover mortality: patch context and timing") +
  theme_bw() + theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

p_death <- ggplot(mover_deaths, aes(x = date, y = ggd, color = patch, size = weight)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_point(alpha = 0.85) +
  scale_color_manual(values = c(warm = "tomato", cold = "steelblue")) +
  scale_size_continuous(range = c(3, 8), name = "Weight (g)") +
  labs(x = "Date", y = "Growth rate at death (g/g/d)", color = "Patch") +
  theme_bw()

egg::ggarrange(p_temp, p_death, ncol = 1, heights = c(1, 1.2))

The early cluster of deaths is dominated by warm-patch fish with negative growth rates and small body sizes — these fish are actively losing mass and dying. Temperatures in both patches are essentially identical and near-zero during this period, so the warm patch offers no thermal advantage, but movers are nonetheless ending up there and dying. The one cold-patch early death follows the same pattern.

The two late deaths (days ~123–133) are different in character: both have positive growth rates and are much larger, suggesting they’re stochastic rather than condition-driven mortality.

The key diagnostic question this raises: why are movers preferentially occupying the warm patch in early winter when temperatures are identical and growth is negative in both? That points to investigating the patch-choice logic at low temperatures — the growth potential lookup may not be differentiating between patches when temperatures are near-flat, leading to essentially arbitrary placement. This makes sense, especially when combined with the softmax/probabilistic stochasticity in patch choice…fish won’t switch patches if the growth advantage is tiny…so they stay in the warm habitat where densities are ~high and

6.6.2.3 Drivers of survival

Now, look at how the mechanisms of survival manifest over time, for each strategy.

Code 
# Build survival probability components for all fish across all alive days
all_surv_df <- tibble(
  pid       = rep(fish_pop_init$pid, n_days),
  dayofsim  = rep(1:n_days, each = nrow(fish_pop_init)),
  condition = as.vector(condition_matrix),
  temp      = as.vector(temp_matrix),
  weight    = as.vector(mass_matrix),
  growth    = as.vector(ggd_matrix),
  age       = as.vector(age_matrix)
) |>
  filter(!is.na(condition)) |>
  left_join(fish_pop_init |> select(pid, strategy), by = "pid") |>
  mutate(
    p_temp  = fncSurviveTemp(temp, T1 = T1_mort, T9 = T9_mort),
    p_starv = fncSurviveStarve(condition, K9 = K9_starv, K1 = K1_starv),
    p_sg    = fncSurviveFixed(pick(weight, growth))[[1]],
    p_age   = fncSurviveAge(age)
  ) |>
  left_join(habitat_df |> select(dayofsim, date), by = "dayofsim") |>
  pivot_longer(c(p_temp, p_starv, p_sg, p_age),
               names_to = "source", values_to = "p_survive") |>
  mutate(
    source   = recode(source,
      p_temp  = "Temperature (p_temp)",
      p_starv = "Starvation (p_starv)",
      p_sg    = "Size/growth (p_sg)",
      p_age   = "Age (p_age)"),
    strategy = factor(strategy,
      levels = c("resident_warm", "resident_cold", "optimal_mover"))
  )

mycols <- brewer.pal(4, "Set2")
surv_source_colors <- c(
  "Temperature (p_temp)" = mycols[1],
  "Starvation (p_starv)" = mycols[2],
  "Size/growth (p_sg)"   = mycols[3],
  "Age (p_age)"          = mycols[4]
)

# Individual fish: one line per fish per mortality source, faceted by strategy
# all_surv_df |>
#   ggplot(aes(x = date, y = p_survive, color = source,
#              group = interaction(pid, source))) +
#   geom_rect(data = summer_shading,
#             aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
#             fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
#   geom_line(alpha = 0.2, linewidth = 0.3) +
#   scale_color_manual(values = surv_source_colors) +
#   #scale_y_continuous(limits = c(0, 1)) +
#   facet_wrap(~strategy, nrow = 3, scales = "free_y") +
#   theme_bw() + theme(panel.grid = element_blank()) +
#   xlab("Date") + ylab("Daily survival probability") +
#   labs(title = "Mortality sources by strategy — individual fish",
#        color = "Mortality source")

# Group means: daily mean per mortality source, faceted by strategy
all_surv_df |>
  group_by(strategy, date, source) |>
  summarise(mean_p = mean(p_survive), .groups = "drop") |>
  ggplot(aes(x = date, y = mean_p, color = source)) +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(linewidth = 0.9) +
  scale_color_manual(values = surv_source_colors) +
  #scale_y_continuous(limits = c(0, 1)) +
  facet_wrap(~strategy, nrow = 3, scales = "free_y") +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Mean daily survival probability") +
  labs(title = "Mortality sources by strategy — group means",
       color = "Mortality source")

  • Warm residents: Starvation (p_starv) is almost entirely responsible for warm resident mortality during summer. It remains at 1.0 through spring (fish growing well, condition = 1.0), then collapses rapidly in late June as the warm patch temperatures drive up metabolic costs faster than fish consume food — condition drops below K9=0.50 and p_starv crashes to ~0.35 within a matter of days.
    • Temperature (p_temp) plays almost no role. Because the warm patch peaks at 25°C — just below the T9 threshold of 25.8°C — p_temp barely budges from 1.0 throughout the survival window.
    • Relative to to p_starv, p_sg stays flat at ~1.0 throughout — by the time summer arrives, warm residents have grown large enough that size/growth survival is essentially at its ceiling.
  • Cold residents: size dependent survival, p_sg, drives mortality among new recruits. p_starv drives mortality during subsequent winters. p_temp plays almost no role as temperatures in the cold habitat never approach lethal limits
  • Optimate movers: similar patterns as cold residents, but smaller effects of winter starvation b/c these fish occupy the warm habitat where slightly warmer temps all fish to feed at slightly rates, better meeting metabolic demands

Condition/starvation is a one of the primary drivers or survival. Plot relative condition over time, for each strategy.

Code 
ibm_summary |>
  # filter(strategy != "resident_warm") |>   # warm residents extinct early — plot separately
  ggplot(aes(x = date, color = strategy, fill = strategy)) +
  geom_rect(data = summer_shading,
            aes(xmin=xmin, xmax=xmax, ymin=-Inf, ymax=Inf),
            fill="grey", alpha=0.3, inherit.aes=FALSE) +
  geom_ribbon(aes(ymin = mean_condition - sd_condition,
                  ymax = mean_condition + sd_condition),
              alpha = 0.15, color = NA) +
  geom_line(aes(y = mean_condition)) +
  geom_hline(yintercept = c(0.50, 0.30), linetype="dashed", alpha=0.5) +
  annotate("text", x = as.Date("2024-06-01"), y = 0.51, hjust=0,
           label="K9 = 0.50", size=3, alpha=0.7) +
  annotate("text", x = as.Date("2024-06-01"), y = 0.31, hjust=0,
           label="K1 = 0.30", size=3, alpha=0.7) +
  scale_color_manual(values=c(resident_warm = "red", resident_cold="blue", optimal_mover="darkgreen")) +
  scale_fill_manual(values=c(resident_warm = "red", resident_cold="blue",  optimal_mover="darkgreen")) +
  scale_y_continuous(limits=c(0,1)) +
  theme_bw() + theme(panel.grid=element_blank()) +
  xlab("Date") + ylab("Relative condition (W / W_peak)") +
  labs(title="Mean condition ± SD by strategy over time",
       subtitle="Dashed lines = starvation anchor points (K9, K1); shaded = June–Aug",
       color="Strategy", fill="Strategy")

Optimal movers and cold residents maintain relatively high condition, and we only see mass loss during the winter when very cold temperatures limit physiological capacity. Neither strategy every approaches the starvation anchor points, meaning that the starvation mechanism is correctly not firing during ~normal seasonal weight loss. In contrast, condition for warm residents declines steeply during the first summer, when food limitations (threshold P_cmax) are insufficient to meet temperature-driven physiological demands.

6.6.3 Age and weight-at-age

Plot histograms of maximum age by strategy

Code 
max_age_df <- ibm_long |>
  dplyr::filter(survived == 1) |>
  dplyr::group_by(strategy, pid) |>
  dplyr::summarize(max_age = max(age), .groups = "drop")

strategy_colors <- c(resident_warm = "red", resident_cold = "blue", optimal_mover = "darkgreen")

ggplot(max_age_df, aes(x = max_age, fill = strategy)) +
  geom_histogram(binwidth = 0.25, color = "white", linewidth = 0.2) +
  scale_fill_manual(values = strategy_colors) +
  facet_wrap(~ strategy, ncol = 1, scales = "free_y") +
  theme_bw() + theme(panel.grid = element_blank(), legend.position = "none") +
  xlab("Maximum age (years)") + ylab("Number of fish") +
  labs(title = "Maximum age by strategy")

Plot weight at age

Code 
ibm_long %>% filter(age %in% c(0:10)) %>%
  ggplot(aes(x = age, y = weight, color = strategy, fill = strategy, group = pid)) +
  geom_line() +
  scale_color_manual(values=c(resident_warm = "red", resident_cold="blue", optimal_mover="darkgreen")) +
  scale_fill_manual(values=c(resident_warm = "red", resident_cold="blue",  optimal_mover="darkgreen")) +
  labs(title = "Individual weight-at-age by strategy", x = "Age (years)", y = "Weight (g)", color="Strategy", fill="Strategy") +
  theme_bw() + theme(panel.grid=element_blank())

6.6.4 Abundance and biomass

Plot abundance through time

Code 
ibm_summary %>% 
  ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = n_alive, color = strategy)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    x     = "Date",
    y     = "Number alive",
    color = "Strategy",
    title = "Population size by strategy"
  ) +
  theme_bw() + theme(panel.grid = element_blank())

Plot biomass through time as abundance * mean weight

Code 
ibm_summary %>% 
  ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = (n_alive*mean_weight)/1000, color = strategy)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_y_continuous(limits = c(0, NA)) +
  labs(
    x     = "Date",
    y     = "Biomass (kg)",
    color = "Strategy",
    title = "Biomass by strategy"
  ) +
  theme_bw() + theme(panel.grid = element_blank())

6.6.5 Movement

6.6.5.1 Population means

Patch occupancy for optimal movers — proportion in warm patch through the year.

  • Stochasticity in movement
    • If move_stochastic = "prob", movers only really start to pile into the warm habitat during the shoulder seasons when benefits of warm habitat far exceed the cold habitat.
    • If move_stochastic = "Indiv", …
  • Environment sensing
    • If sense_environment = "density_current" or sense_environment = "density_all", movers bounce rapidly between the two habitats as they chase the apparently better habitat, but then experience reduced growth due to strong density dependence. This appears to be driven by fact that most movers are doing the same thing.
Code 
ibm_summary %>%
  filter(strategy == "optimal_mover") %>%
  ggplot(aes(x = date, y = prop_warm)) +
  geom_line() +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
  theme_bw() +
  xlab("Date") + ylab("Proportion of optimal movers in warm patch")

Code 
ibm_summary %>%
  filter(strategy == "optimal_mover", date < as.Date("2024-12-31")) %>%
  ggplot(aes(x = date, y = prop_warm)) +
  geom_line() +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "gray50") +
  theme_bw() +
  xlab("Date") + ylab("Proportion of optimal movers in warm patch")

6.6.5.2 Individual movement patterns

Build a long-form patch history for optimal movers:

Code 
patch_long <- as_tibble(patch_matrix) |>
  bind_cols(fish_pop_init |> select(pid, strategy)) |>
  pivot_longer(-c(pid, strategy), names_to = "day", values_to = "patch") |>
  mutate(dayofsim = as.integer(sub("V", "", day))) |>
  left_join(habitat_df |> select(dayofsim, date), by = "dayofsim") |>
  filter(!is.na(patch)) |>
  select(-day)

mover_patch <- patch_long |> filter(strategy == "optimal_mover")

Plot 1 — Raster of all optimal movers. Each row is one fish, each column is a day, fill = patch. Fish are ordered by their overall proportion of time in the warm patch. Synchronised individuals will produce horizontal bands switching at the same time; independent movers will produce a noisier mosaic.

At seasonal and annual timescales, fish show similar movement patterns–tracking the primary gradient in growth potential and it shifts seasonally between warm and cold patches. But we also see individual variation at finer timescales, driven by stochasticity from the softmax patch selection function. Another observation is that the allometric scaling of the cost of movement appears to be working as expected. Early in the simulation, small fish switch infrequently (only when the benefit of switching is very large). But later on, larger fish switch more frequently, navigating smaller magnitude differences in density-dependent growth potential.

Code 
# Order fish by proportion of time spent in warm patch
pid_order <- mover_patch |>
  group_by(pid) |>
  summarise(prop_warm = mean(patch == "warm"), .groups = "drop") |>
  arrange(prop_warm)

mover_patch |>
  mutate(pid = factor(pid, levels = pid_order$pid)) |>
  ggplot(aes(x = date, y = pid, fill = patch)) +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
            fill = "grey", alpha = 0.4, inherit.aes = FALSE) +
  geom_tile() +
  scale_fill_manual(values = c(warm = "#d95f02", cold = "#1f78b4")) +
  theme_bw() +
  theme(panel.grid  = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  xlab("Date") + ylab("Individual fish (ordered by % time in warm)") +
  labs(title = "Patch occupancy — all optimal movers",
       subtitle = "Shaded = June–Aug; fish ordered bottom-to-top by increasing warm-patch use",
       fill = "Patch")

Code 
# mover_patch |>
#   filter(date < as.Date("2024-06-01")) %>%
#   mutate(pid = factor(pid, levels = pid_order$pid)) |>
#   ggplot(aes(x = date, y = pid, fill = patch)) +
#   # geom_rect(data = summer_shading,
#   #           aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
#   #           fill = "grey", alpha = 0.4, inherit.aes = FALSE) +
#   geom_tile() +
#   scale_fill_manual(values = c(warm = "#d95f02", cold = "#1f78b4")) +
#   theme_bw() +
#   theme(panel.grid  = element_blank(),
#         axis.text.y = element_blank(),
#         axis.ticks.y = element_blank()) +
#   xlab("Date") + ylab("Individual fish (ordered by % time in warm)") +
#   labs(title = "Patch occupancy — zoomed to first month of simulation",
#        subtitle = "Shaded = June–Aug; fish ordered bottom-to-top by increasing warm-patch use",
#        fill = "Patch")

Plot 2 — Individual traces for 10 randomly sampled movers. Patch is encoded as a binary y-position (warm = 1, cold = 0) with a small per-fish y-offset to separate overlapping lines. A fish that moves independently of its peers will switch at a different time from the others.

Code 
set.seed(4921)
sample_pids <- sample(unique(mover_patch$pid), 10)

mover_patch |>
  filter(pid %in% sample_pids) |>
  mutate(patch_y   = if_else(patch == "warm", 1, 0),
         pid_label = factor(pid)) |>
  ggplot(aes(x = date, y = patch_y, color = pid_label, group = pid_label)) +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_step(alpha = 0.8, linewidth = 0.6) +
  scale_y_continuous(breaks = c(0, 1), labels = c("Cold", "Warm")) +
  theme_bw() + theme(panel.grid = element_blank()) +
  xlab("Date") + ylab("Patch") +
  labs(title = "Patch occupancy — 10 randomly sampled optimal movers",
       subtitle = "Each line is one fish; divergence between lines = independent decision-making",
       color = "Fish ID")

Plot 3 — Switch event timing. For each fish, record every day a patch switch occurs and plot as a rug. Vertical clustering = synchronised switching; diffuse scatter = independent behaviour.

Code 
switch_events <- mover_patch |>
  arrange(pid, date) |>
  group_by(pid) |>
  mutate(switched = patch != lag(patch)) |>
  filter(switched) |>
  ungroup()

switch_events |>
  mutate(to_patch = patch) |>
  ggplot(aes(x = date, y = factor(pid))) +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = -Inf, ymax = Inf),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_point(aes(color = to_patch), size = 1.2, alpha = 0.7) +
  scale_color_manual(values = c(warm = "#d95f02", cold = "#1f78b4")) +
  theme_bw() +
  theme(panel.grid   = element_blank(),
        axis.text.y  = element_blank(),
        axis.ticks.y = element_blank()) +
  xlab("Date") + ylab("Individual fish") +
  labs(title = "Patch switch events — optimal movers",
       subtitle = "Each dot = a switch; colour = destination patch. Vertical alignment = synchronised switching",
       color = "Destination")

Fish are generally acting similarly at seasonal and annual timescales, but on any given day fish are largely acting independently. To introduce additional variation in movement strategies, consider the following:

  • Inrease tau, the sensitivity to growth differences specified by the softmax function (steepness). Increasing tau would make all fish more randomly noisy in their decisions (similar to if we add simple random variation) but also may have the unintended consequence of making fish more reluctant to move (i.e., fish need a larger growth benefit to initiate a movement). But tau also applies identically to every fish…so any reduced synchrony in movement decisions would be entirely stochastic, rather than biologically meaningful individual variation
  • To generate structured among-individual differences in movement, we need to give each fish a persistent, individual-specific response to the same environment. move_stochastic = "indiv" with sigma_bold is the mechanism for this. Each fish draws a fixed movement threshold from a distribution, so some fish are consistently “bold” (move even when the gain is marginal) and others are consistently “reluctant” (only move when the growth advantage is large). That creates genuine individual variation in movement decisions that persists through time — not just day-to-day noise.

6.6.6 Experienced temperature

Temperature experienced by optimal movers

Code 
p1 <- ibm_summary %>%
  filter(strategy == "optimal_mover") %>%
  ggplot(aes(x = date, color = strategy, fill = strategy)) +
  geom_ribbon(aes(ymin = mean_temp - sd_temp, ymax = mean_temp + sd_temp),
              alpha = 0.15, color = NA) +
  geom_line(aes(y = mean_temp)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_fill_manual(values  = c(resident_warm = "red", 
                                resident_cold = "blue",
                                optimal_mover = "darkgreen")) +
  theme_bw() + #ylim(0,10000) +
  xlab("Date") + ylab("Mean experienced temperature (deg. C)")


p2 <- temp_long %>%
  filter(strategy == "optimal_mover") %>%
  ggplot(aes(x = date, y = temp, color = strategy, group = pid)) +
  geom_line(alpha = 0.5) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  theme_bw() + #ylim(0,10000) +
  xlab("Date") + ylab("Individual experienced temperature (deg. C)")

ggpubr::ggarrange(p1, p2, nrow = 1, common.legend = TRUE)

6.6.7 Experienced growth

Experienced standardized growth rates over time, by strategy.

Code 
p1 <- ibm_summary %>%
  # filter(strategy == "optimal_mover") %>%
  ggplot(aes(x = date, color = strategy, fill = strategy)) +
  geom_ribbon(aes(ymin = mean_ggd - sd_ggd, ymax = mean_ggd + sd_ggd),
              alpha = 0.15, color = NA) +
  geom_line(aes(y = mean_ggd)) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  scale_fill_manual(values  = c(resident_warm = "red", 
                                resident_cold = "blue",
                                optimal_mover = "darkgreen")) +
  theme_bw() + #ylim(0,10000) +
  geom_abline(slope = 0, intercept = 0, linetype = "dashed") +
  xlab("Date") + ylab("Mean growth rate (g/g/d)")

p2 <- ggd_long %>%
  ggplot(aes(x = date, y = ggd, color = strategy, group = pid)) +
  geom_line(alpha = 0.5) +
  scale_color_manual(values = c(resident_warm = "red", 
                                resident_cold = "blue", 
                                optimal_mover = "darkgreen")) +
  theme_bw() + #ylim(0,10000) +
  xlab("Date") + ylab("Individual growth rate (g/g/d)")

ggpubr::ggarrange(p1, p2, nrow = 1, common.legend = TRUE)

6.6.8 Boldness

Compile end of simulation results

Code 
fish_results <- fish_pop %>%
  mutate(
    final_weight   = mass_matrix[, n_days],
    n_switches     = switches
    # boldness       = move_threshold   # positive = shy, negative = bold
  )

fish_results %>%
  group_by(strategy) %>%
  summarise(mean_wt = round(mean(final_weight), 1),
            mean_switches = round(mean(n_switches), 1))

Plot only if move_stochastic = "indiv"

Code 
if (move_stochastic == "indiv") {
  
movers <- fish_results %>% filter(strategy == "optimal_mover")

# Panel A: boldness vs final weight
p_wt <- ggplot(movers, aes(x = boldness, y = final_weight)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "steelblue", linewidth = 0.8) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  annotate("text", x =  0.006, y = min(movers$final_weight) + 5,
           label = "shy →",  hjust = 0, size = 3.5, color = "grey40") +
  annotate("text", x = -0.006, y = min(movers$final_weight) + 5,
           label = "← bold", hjust = 1, size = 3.5, color = "grey40") +
  labs(x = "Boldness (move_threshold)", y = "Final weight (g)",
       title = "A: Boldness vs. end-of-season weight")

# Panel B: boldness vs number of patch switches
p_sw <- ggplot(movers, aes(x = boldness, y = n_switches)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "tomato", linewidth = 0.8) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  labs(x = "Boldness (move_threshold)", y = "No. patch switches",
       title = "B: Boldness vs. switching frequency")

# Panel C: switching frequency vs final weight  
p_sw_wt <- ggplot(movers, aes(x = n_switches, y = final_weight)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = TRUE, color = "mediumpurple", linewidth = 0.8) +
  labs(x = "No. patch switches", y = "Final weight (g)",
       title = "C: Switching frequency vs. end-of-season weight")

cowplot::plot_grid(p_wt, p_sw, p_sw_wt, ncol = 1)

}