Appendix_S1.Rmd

---
title: |
  | Appendix S1
  |
  | Optimal energy allocation trade-off driven by size-dependent physiological and demographic responses to warming
  |
  | Ecology
author: Viktor Thunell$^1$$^*$, Anna Gårdmark$^1$, Magnus Huss$^1$ & Yngvild Vindenes$^2$

output:
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

1. Department of Aquatic Resources, Swedish University of Agricultural Sciences, Box 7018, 750 07 Uppsala, Sweden

2. Centre for Ecological and Evolutionary Synthesis (CEES), Department of Biosciences, University of Oslo, 0316 Oslo, Norway

\* Author to whom correspondence should be addressed. Current address:
Viktor Thunell, 1. Department of Aquatic Resources (SLU Aqua), Swedish University of Agricultural Sciences, Box 7018, 750 07 Uppsala, Sweden. Email: viktor.thunell@slu.se


# Table of contents

Section S1. Size and temperature dependent maintenance, consumption, allocation and their effects on the energy budget.

Section S2. Fitting the DEB-model at mean temperature to data on size at age and reproduction in Windermere pike.

Section S3. Effects of mass-independent energy allocation on the DEB and modeled size at age and fecundity.

Section S4. Effect of number of size classes on long-term population growth rate ($\lambda$) and stable population size structure ($w$).

Section S5. Sensitivity analysis of long-term population growth rate ($\lambda$) to $\kappa_0$: Methods and analytical decomposition into contributions from demographic functions.

Section S6. Growth trajectories and fecundity at optimal energy allocation ($\kappa_0^*$).

Section S7. Fitness landscape and stable size structure for survival probability contrasts.

\newpage

## Section S1. Size and temperature dependent maintenance, consumption, allocation and their effects on the energy budget.

### Section S1.1. Size and temperature dependent maintenance and consumption and the energy budget

The effects of mass and temperature on maintenance and consumption in combination with mass dependent allocation determine somatic growth and reproductive reserve in the Dynamic Energy Budget Integral Projection Model (DEB-IPM, Fig. S1.1). In the DEB-part, the effect of mass on maintenance and consumption is modeled using a power function (Fig. S1.1a), the interacting effects of temperature and mass (equation 3, main text) on maintenance rate is exponential (Fig. S1.1b) and the effect of temperature on consumption rate is unimodal (Fig. S1.1b, solid line). We assume allocation to reproductive reserve to increase with mass according to data (nearly linearly, see S2) and not with consumption rate (which is sub-linear). This means that allocation to growth ($\kappa(m)$ with a mass intercept $\kappa_0$) decreases with mass (Fig. S1.1c). The resulting DEB, which includes assimilation efficiency ($\alpha =0.4$), thus describes the proportion of consumed energy that each compartment use up (Fig. S1.1d).

```{r, echo = FALSE, warning = FALSE}
test_Pars <- c(T = 287, kappa = deb.optim$par[1], Y = 1)

cons_E <- (eps1*x^eps2)*rC_T_Pad(test_Pars['T'])

E_bud <- bind_cols(Size=x,
                   Assimilation = alpha*cons_E,
                   Maintenance =  rho1*x^rho2*rM_T_A(test_Pars['T']),
                   Growth =    kap_fun(x,test_Pars["kappa"])*alpha*cons_E - rho1*x^rho2*rM_T_A(test_Pars['T']),
                   Reserve = (1-kap_fun(x,test_Pars["kappa"]))*alpha*cons_E )

E_bud_long <- pivot_longer(E_bud, cols = c(2:ncol(E_bud)), names_to = "comp", values_to = "gram")
E_bud_long$comp <- factor(E_bud_long$comp, levels=c( "Reserve", "Growth", "Maintenance", "Assimilation") )

Eb <- E_bud_long %>%
  ggplot() +
  geom_area(aes(Size, gram, fill = comp)) +
  coord_cartesian(xlim = c(min(x),18000)) +
  scale_fill_manual(values = c( "#252525", "#cccccc", "#969696", "#636363" ), name="Compartment") +
  xlab("Mass [g]" ) +
  ylab("Energy [g]" ) +
  ggtitle('(d)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold")) 

### Plot Size functions
Srate <- bind_cols(Size=x,
                   Consumption=eps1*x^eps2,
                   Maintenance=rho1*x^rho2)
Srate_long <- pivot_longer(Srate, cols=c(2:3), names_to = "Rate", values_to = "Effect")

Sr <- Srate_long %>%
  ggplot() +
  geom_line(aes(Size, Effect, linetype = Rate), size=0.8) +
  scale_x_continuous(expand = c(0,0), name = "Mass [g]", limits = c(0,20000), breaks = seq(0,20000,5000)) +
  scale_linetype_manual(values = c("solid","dashed"),
                        labels = c("Consumption","Maintenenance")) +
  ylab("Effect [g/day]") +
  ggtitle('(a)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"), 
        legend.position="none")

### Plot Temp functions
Trate <- bind_cols(Temp=273:300,
                   Consumption=rC_T_Pad(273:300),
                   Maintenance1=rM_T_A(273:300))
Trate_long <- pivot_longer(Trate, cols=c(2:3), names_to = "Rate", values_to = "Effect")

Tr <- Trate_long %>%
ggplot() +
  geom_line(aes(Temp, Effect, linetype = Rate), size=0.8) +
  scale_x_continuous(expand = c(0,0), name = "Temperature [K]", limits = c(274,300), breaks = seq(275,300,5))+
  scale_linetype_manual(values = c("solid","dashed"), labels = c("Consumption","Maintenenance")) +
  ylab("") +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold")) 

Kf_long <- rbind( tibble(Size=x, val = 1-kap_fun(x,kappa=deb.optim$par[1]), k="Rep. reserve"),
  tibble(Size=x, val=kap_fun(x,kappa=deb.optim$par[1]), k="Growth"))
Kf <-
ggplot(Kf_long) +
  geom_line( aes(Size, val, color=k), size=0.8, ) +
  ylab(expression(paste("Allocation,", italic(" k(m)"),))) +
  scale_color_manual(values = c( "#cccccc" ,"#252525" ), name="") +
  scale_x_continuous(expand = c(0,0), name = "Mass [g]", limits = c(0,20000), breaks = seq(0,20000,5000)) +
  ggtitle('(c)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none",
        legend.key.width = unit(1, 'cm')) 
```

```{r, echo = FALSE, warning = FALSE, fig.height = 4, fig.width = 6}
(Sr + Tr) / (Kf + Eb)
```

Figure S1.1. Size- and temperature dependence of the processes that constitute the dynamic energy budget part of the model. The mass (a) and temperature (b) scaling of maintenance (a and b, dotted line for a 1 gram individual and dashed line for a 100 gram) and consumption (a and b, solid line) in combination with mass dependent energy allocation of reserve and growth (c, via $\kappa(m)$) determines the size dependence of the energy budget (d).

\newpage

### Section S1.2. Effects of temperature on energy allocation to growth and reproduction

The effects of mass and temperature on maintenance and consumption rate (Fig. S1.1) make the rate of energy allocated to growth and reproductive reserve unimodal over temperature (Fig. S1.2). The optimum temperature for energy spent on reproduction is a function of consumption, but not of mass, and is therefore 291 Kelvin (Fig. S1.2, dotted line, Fig. S1.1b, solid line). The optimum temperature for somatic growth however, decreases with mass (Fig. S1.2, solid line) due to the disparate effects of size on maintenance and consumption rate combined with a unimodal effect of temperature on consumption (Lindmark et al. 2022).

```{r, echo = FALSE, warning = FALSE, fig.height = 3, fig.width = 6}
Temp <- seq(280, 300, 0.2)
test_Pars <- c(T = 287, 
             kappa = deb.optim$par[1], 
             Y = 1)

T_rates2 <- NULL
for(i in Temp) {
  maint <- rho1*(x^rho2)*rM_T_A(i)
  intake  <- eps1*(x^eps2)*rC_T_Pad(i) 
  T_rates2 <- 
    rbind(T_rates2,cbind(Temp = i, mass = x, maint, intake, 
                         growth_E = kap_fun(x,test_Pars["kappa"])*alpha*test_Pars["Y"]*intake - maint,
                         repro_E = alpha* test_Pars["Y"]*(1-kap_fun(x,test_Pars["kappa"]))*intake))
  }
T_rates <- T_rates2 
maxgT<- as.data.frame(T_rates) %>% 
        group_by(mass) %>%
        slice_max(0.35*(intake-maint))

T_rates_long <- pivot_longer(as_tibble(T_rates), 5:6, names_to = "E_type", values_to = "g_day")
T_rates_long["g_g_day"] <- T_rates_long$g_day/T_rates_long$mass


T_rates_long %>%
  filter(as.factor(mass) %in% c(T_rates_long$mass[c(5,50,500,5000)])) %>%
  ggplot(.) +
    geom_line(size=0.7 , aes(Temp, g_day, color = as.factor(mass), 
                         linetype = as.factor(E_type))) +
    scale_color_manual(values = c( "#252525", "#636363", "#969696", "#cccccc" ), name = "Mass [g]", labels = round(c(T_rates_long$mass[c(5,50,500,5000)]))) +
    scale_linetype_manual(values = c(  "solid", "dotted"), name = "Compartment", labels = c( "Growth", "Rep. reserve")) +
    ylab("Allocation rate [g/day]") +
    xlab("Temperature [K]") +
    ylim(0,NA)+
    theme_light(base_size=8)
  
```

Figure S1.2. Energy allocation rate to reproductive reserve (dotted lines) and somatic growth (solid lines) for four different sizes.

\newpage

## Section S2. Fitting the DEB-model at mean temperature to data on size at age and reproduction in Windermere pike

Four parameters in the DEB ($\epsilon_{1,2}$ and $\kappa_{0,m}$) are free and we estimated their values by numerical optimization of the DEB-function to fit observed size at age and fecundity of Windermere pike. The size at age data on Windermere pike is based on back-calculated length (Winfield et al. 2013b). The data on fecundity represents individuals caught between 1963-2003 and for growth 1944-2003 and therefore encompass a wide temperature range (Winfield et al. 2013, 2013a, 2013b). We wish to reduce such variation to get closer to average fecundity and growth trajectories at mean temperature (287 K) as we add temperature effects in the DEB via consumption and maintenance. Therefore we removed data on individuals that experienced high or low growing season temperatures (April-September, single seasons for fecundity and mean life time temperature for somatic growth). We removed data with temperatures from the first and fourth quartile of the growing season temperature range leaving 16191 of 34215 size at age measurements (3573 of 7939 individual growth trajectories) and 1139 of 3696 individuals measured for fecundity. The difference between parameter estimates using the selected data and all data was $0$ for $\epsilon_1$, $0.38$ for $\epsilon_2$, $-0.01$ for $\kappa_0$ and $0$ for $\kappa_m$.

We convert length to weight based on our estimated length to weight relationship (eq. 11, main text). Fig. S2 shows how our estimates compare with data for the mean (growth trajectories, Fig. S2a) and variance (Fig. S2c) of growth and of reproduction (S2b, fecundity measured in October-December) at mean temperature (287 K).

For the optimization of free parameters, we used package optimParalell() (a parallel core wrapper for "optim()") in R using the "L-BFGS-B" algorithm (Byrd et. al. 1995). For starting values in the optimization, we use values that provides visually good fits of the model to the data and are reasonable for fish: $\epsilon_1=1$, $\epsilon_2=0.5$, $\kappa_0=0.8$ and $\kappa_m=2$. We calculated confidence intervals from 1000 bootstrap estimates of the modeled data.

\newpage

```{r, echo = FALSE, warning = FALSE, fig.height = 4, fig.width = 6}
a <- data.frame() # 10 individual growth traj, 20+1 years of growth
for(i in 1:10){ 
  a[i,1] <- sample(x = x, size = 1 , replace = TRUE, prob = age1size(test_Pars))
    for(j in 2:20){ # 20 years of growth
      a[i,j] <- sample(x = x, size = 1 , replace = TRUE, prob = growthfun(a[i,j-1],ratefun(a[i,j-1], test_Pars), test_Pars))
  }
}

b <- data.frame() # 10 individual growth traj, 20+1 years of growth
for(i in 1:10){ 
  b[i,1] <- sample(x = x, size = 1 , replace = TRUE, prob = age1size(Pars = c(T = 291,                                           kappa = deb.optim$par[1], Y = 1) ))
    for(j in 2:20){ # 20 years of growth
      b[i,j] <- sample(x = x, size = 1 , replace = TRUE, prob = growthfun(b[i,j-1], ratefun(b[i,j-1], c(T = 291,kappa = deb.optim$par[1], Y = 1)), Pars = c(T = 291,kappa = deb.optim$par[1], Y = 1)))
  }
}

ranTraja <- bind_cols(age=1:20,TK=287,t(a))
ranTrajb <- bind_cols(age=1:20,TK=291,t(b))
ranTraj <- bind_rows(ranTraja,ranTrajb)
ranTraj_long <- pivot_longer(ranTraj, cols=3:12, names_to = "ind", values_to = "mass") 

ranT <- 
ggplot() +
  geom_line(data=data.g, aes(Age, ltow(Length), group = as.factor(Ind), color = "WData")) +
  geom_line(data=ranTraj_long, size=0.4,  aes(age, mass, group = interaction(as.factor(ind), as.factor(TK)), color = as.factor(TK), color="DEB-IPM")) +
  scale_color_manual(name = '', labels = c("DEB-IPM 287 K","DEB-IPM 291 K","Wind. data"), values = c('#91bfdb','#fc8d59',"black")) +
  ylab(expression(paste(italic("g"),"(",italic("m"["s+1"]),italic(",T"),")"))) +
  xlab("Age [years]") +
  ggtitle('(c)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.key.width = unit(1, 'cm')) 
```

```{r, echo = FALSE, warning = FALSE, fig.height = 5, fig.width = 7}
load("optim_pred.RData")

gci <- as_tibble(pred$gci,c(20000,20000))
colnames(gci) = c("low","high")

fci <- as_tibble(pred$fci,c(20000,20000))
colnames(fci) = c("low","high")

groz <- 
ggplot() +
  geom_point(data = growthTemp_select, aes(Age, mass), alpha=0.2) +
  geom_point(data = data.g, aes(Age, body_mass), alpha=0.2, shape=4) +
  geom_line(data = pred$g, aes(Age, Mass, color = "#91bfdb"), size=0.85) +
  geom_ribbon(data = gci, aes(x = 1:(nrow(gci)), ymin = low, ymax = high), fill = "grey", alpha=0.5) +
  ylab(expression(paste(italic("g"),"(",italic("m"["s+1"]),italic(",T"),")"))) +
  xlab("Age [years]") +
scale_color_manual(values = c('#91bfdb'), name="Temp [K]") +
  ggtitle('(a)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none",
        legend.key.width = unit(1, 'cm'))   

repz <- 
ggplot() +
  geom_point(data=fecTemp_select, aes(mass, Egg_number*Egg_weight/e_m), alpha=0.2) +
  geom_point(data=data.f, aes(ltow(Length), Egg_mass/e_m), alpha=0.2, shape=4) +
  geom_line(data = pred$f, aes(mass, fecund), color = "#91bfdb", size=0.85) +
  geom_ribbon(data = fci, aes(x = pred$f$mass, ymin = low, ymax = high), fill = "grey", alpha=0.5) +
    ylab(expression(paste(italic("f"),"(",italic("m"["s"]),italic(",T"),")"))) +
  xlab(expression(paste("Mass ",italic("m")[s]," [g]"))) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none",
        legend.key.width = unit(1, 'cm'))

(groz + repz) / ranT 
```

Figure S2. Comparison of the temperature dependent DEB-IPM and data on Windermere pike. Circles in a and b are data points used for parameter  optimization (crosses are removed data, see text) and lines in c) correspond to individual growth trajectories (all data). The predicted mean size at age (a, shaded line area are confidence intervals from 1000 bootstrap estimates) and size specific fecundity (b, eggs produced in year s+1 by individuals of size  $m_s$ in year s) in the DEB-IPM represent $\kappa_0=0.89$ (i.e. 89 % of daily assimilated energy is allocated to growth and maintenance). Indvidual growth trajectories in c) (black, corresponds to point data in a) are based on individual back calculated growth from wing bones from Windermere pike and 10 random growth trajectories produced by the DEB-IPM at 287 K (blue) and 291 K (orange).

\newpage

## Section S3. Effects of mass-independent energy allocation

The standard DEB-model assumes that energy allocation is independent of body mass, i.e. $\kappa$ is constant. Below we show the basis for our assumption that allocation to reproductive reserve increases with body mass (at the cost of somatic growth). 

We evaluated the energy budget over size and the modeled growth trajectories (mean size at age) and size specific fecundity for the mass dependent and energy allocation and  the standard DEB-model. We visually compared the growth trajectories and size dependent fecundity based on parameter estimates from the optimization procedure described in Section S2 assuming either a mass-dependent (eq. 5 main text, Fig. S2a, b) or independent (constant $\kappa$, Fig. S3.1 and 2) allocation. 
The main difference between the two energy allocation assumptions in terms of how they compare to data on Windermere Northern pike is that the mass independent energy allocation overestimates fecundity for smaller individuals and underestimates fecundity for larger sizes. Based on this difference we choose the mass dependent allocation for the model in the main text. 


```{r, echo = FALSE, warning = FALSE, fig.height = 2, fig.width = 6}
kap_fun <- function(m, kappa, ha=deb.optim$par[2]){kappa} #for body size indep. kappa

test_Pars_mi =  c(T = 287, kappa = 0.8254219, Y = 1)

eps1=1.2356264 
eps2=0.4394062

cons_E <- (eps1*x^eps2)*rC_T_Pad(test_Pars_mi['T'])

E_bud <- bind_cols(Size=x,
                   Assimilation = alpha*cons_E,
                   Maintenance =  rho1*x^rho2*rM_T_A(test_Pars_mi['T']),
                   Growth =    kap_fun(x,test_Pars_mi["kappa"])*alpha*cons_E - rho1*x^rho2*rM_T_A(test_Pars_mi['T']),
                   Reserve = (1-kap_fun(x,test_Pars_mi["kappa"]))*alpha*cons_E )

E_bud_long <- pivot_longer(E_bud, cols = c(2:ncol(E_bud)), names_to = "comp", values_to = "gram")
E_bud_long$comp <- factor(E_bud_long$comp, levels=c( "Reserve", "Growth", "Maintenance", "Assimilation") )

Eb_const<-E_bud_long %>%
  ggplot() +
  geom_area(aes(Size, gram, fill = comp)) +
  coord_cartesian(xlim = c(min(x),max(x))) +
  scale_fill_manual(values = c( "#252525", "#cccccc", "#969696", "#636363" ), name="Compartment") +
  xlab("Mass [g]" ) +
  ylab("Energy [g]" ) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.key.width = unit(1, 'cm')) 

Kf_long <- rbind( tibble(Size=x, val = 1-kap_fun(x,kappa=test_Pars_mi["kappa"]), k="Rep. reserve"),
  tibble(Size=x, val=kap_fun(x,kappa=deb.optim$par[1]), k="Growth"))

Kf_const <-
ggplot(Kf_long) +
  geom_line( aes(Size, val, color=k), size=0.8) +
  ylab(expression(paste("Allocation,", italic(" k(m)"),))) +
  scale_color_manual(values = c( "#cccccc" ,"#252525" ), name="") +
  scale_x_continuous(expand = c(0,0), name = "Mass [g]", limits = c(0,20000),breaks = seq(0,20000,5000)) + 
  ggtitle('(a)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none") 

Kf_const + Eb_const
```

Figure S3.1. Effects of a mass independent energy allocation (a, constant $\kappa$) on the energy budget (b, compare to Fig. S1.1c and d).

```{r, echo = FALSE, warning = FALSE, fig.height = 3, fig.width = 6}
Temp <- seq(285,293,2)
t <- 25

traj = data.frame() 
 for(i in 1:length(Temp)){
   traj[1,i] <- x[which.max(age1size(Pars = c(T = Temp[i], kappa = test_Pars_mi[["kappa"]], Y = 1) ))]
    for(j in 2:t){ 
      traj[j,i] <- x[which.max(growthfun(traj[j-1,i], ratefun(traj[j-1,i], Pars = c(T = Temp[i], kappa = test_Pars_mi[["kappa"]], Y = 1)),
                                         Pars = c(T = Temp[i], kappa = test_Pars_mi[["kappa"]], Y = 1)))]
   }
 }
colnames(traj) <- Temp[1:5]
traj<-cbind(Age=1:nrow(traj),traj)
trajG_mi <- pivot_longer(traj, cols=c(2:6),names_to = "Temp", values_to = "mass" )
trajG_mi$Temp <- as.double(trajG_mi$Temp)

groT_const <- ggplot(trajG_mi) +
  geom_point(data=data.g, aes(Age, body_mass), alpha=0.1) +
  geom_line(data = trajG_mi, size=0.85, aes(Age, mass, color = as.factor(Temp))) +
  ylab(expression(paste("Mass ", italic("\u03BC")["g"]," [g]"))) +
  xlab("Age [years]") +
  xlim(0,20) +
  ylim(0,20000) +
  scale_color_manual(values = c('#4575b4','#91bfdb','#fee090','#fc8d59','#d73027'), name="Temp [K]") +
ggtitle('(a)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none",
        legend.key.width = unit(1, 'cm')) 

#### Plot size dependent fecundity, f(m_s,T) via repfun() ####
# Note that we need to double (*2) the fecundity when plotting using repfun(). Repfun account only for female eggs.
rep_mi=data.frame()
for(i in 2:ncol(traj)){ # the temps, from 2nd last column in traj
  for(j in 1:nrow(traj)){ # the sizes in traj, i.e. Age traj
    if(traj$Age[j]==1){
      rep_tplus1 = cbind(traj[j,i], repfun(e_m, ratefun(e_m, Pars = c(T = Temp[i-1], kappa = test_Pars_mi[["kappa"]], Y = 1)), 
                                            c(T = Temp[i-1], kappa = test_Pars_mi[["kappa"]], Y = 1))*2, Temp[i-1], test_Pars_mi[["kappa"]], "Y"=1)} 
    else{ rep_tplus1 = cbind(traj[j,i], repfun(traj[j-1,i], ratefun(traj[j-1,i], Pars = c(T = Temp[i-1], kappa = test_Pars_mi[["kappa"]], Y = 1)),
                                                c(T = Temp[i-1], kappa = test_Pars_mi[["kappa"]], Y = 1))*2, Temp[i-1], test_Pars_mi[["kappa"]], "Y"=1)}
    colnames(rep_tplus1) = c("mass","fec","Temp", "Kappa", "Y")
    rep_mi = rbind(rep_mi, rep_tplus1)
  }
}
repT_const <-
  ggplot() +
  geom_point(data=data.f, aes(body_mass, Egg_mass/e_m), alpha=0.1) +
  geom_line(data=as_tibble(rep_mi), size=0.85, aes(mass, fec, color = as.factor(Temp)))+
  scale_color_manual(values = c('#4575b4','#91bfdb','#fee090','#fc8d59','#d73027'), name="Temp [K]")+
  ylab(expression(paste(italic("f"),"(",italic("m"["s"]),italic(",T"),")"))) +
  xlab(expression(paste("Mass ",italic("m")[s]," [g]"))) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.key.width = unit(1, 'cm')) 

#reset rates to main model
eps1=deb.optim$par[3]
eps2=deb.optim$par[4]
kap_fun <- function(m, kappa, ha=deb.optim$par[2]){kappa*exp(-m/(ha*max(x)))} #for body size dep. kappa

groT_const + repT_const
```

Figure S3.2. Effects of a mass independent energy allocation (constant $\kappa$) on growth trajectories (a, mean size at age) and size specific fecundity (b) compared to data on Windermere pike (points in a and b, compare S2a and b).

\newpage

## Section S4. Effect of number of size classes on long-term population growth rate ($\lambda$) and stable population size structure ($w$)

When analyzing the IPM, the number of size classes in the population (or matrix mesh size) should not affect the estimation of $\lambda$ and stable size structure ($w$). We therefore estimate $\lambda$ for DEB-IPMs with $0<n<1400$ (Fig. S3) to find the minimum number size classes that does not affect $\lambda$ to an accuracy of three decimal points nor the qualitative shape of the $w$. This resulted in $n=1000$. 

```{r, echo = FALSE, warning = FALSE, fig.height = 4, fig.width = 6}
#  ns <- rev(seq(200,1400,200))
#  ts <- c(285,287,289)
#  nlam=data.frame()
#  nwlam=NULL
#  nvlam=list()
#  for(j in ts) {
#    for(i in ns){
#        n <- i
#        x <- seq(mmin,mmax, length=n)
#        dx <- x[2] - x[1] # step size (grams) in the continuos size
#        lam <- wvlambda.projection(K.matrix(Pars <- c(T = j,  kappa = deb.optim$par[1], Y = 1) ))
#        nlam <- rbind(nlam,c(j,i,lam[[1]]))
#        nwlam <- rbind(nwlam,cbind(j,i,x,lam[[2]][2:nrow(lam[[2]]),]))
#        #nvlam <- c(nvlam,c(cbind(i,x,lam[[3]])))
#    }
#  }
# 
# colnames(nlam) <- c("T","n","lambda")
# colnames(nwlam) <- c("T","n","x","w")

#save(nlam, file = "nlam.RData")
#save(nwlam, file = "nwlam.RData")
load("nlam.RData")
load("nwlam.RData")
nw.labs <- c("285 K", "287 K", "289 K")
names(nw.labs) <- c("285", "287", "289")
nW <-
  as_tibble(nwlam) %>%
  filter(n %in% c(seq(200,1400, by=400))) %>%
ggplot(., aes(x,w, color=as.factor(n))) +
  geom_line(size=0.5) +
  facet_grid(.~T, scales="free", labeller = labeller(T = nw.labs)) +
  xlab(expression(paste("Mass ",italic("m"["s"])," [g]"))) +
  xlim(0,1500) +
  scale_color_discrete(name = expression(paste("# size classes" ))) + 
  ylab("Stable structure w") +
  theme_light(base_size=8) +
  theme(strip.text = element_text(color = "black"),
        strip.background = element_blank())

nL <-
as_tibble(nlam) %>%
  ggplot(., aes(n,lambda, color = as.factor(T) )) +
  geom_line(size = 0.5) +
  scale_color_manual(values = c('#4575b4','#91bfdb','#fee090','#fc8d59','#d73027'), name="Temp [K]")+
  ylab(expression(paste(italic(lambda)))) +
  xlim(200,NA) +
  xlab(expression(paste("# size classes"))) + 
  theme_light(base_size=8)

nL / nW

n <- 1000
x <- seq(mmin,mmax, length=n)
dx <- x[2] - x[1] # step size (grams) in the continuous size
```

Figure S4. Effect of number of size classes ($n$) on long-term population growth rate ($\lambda$, top panel) and stable size structure ($w$, bottom panel). We use $n=1000$ in the main analysis of the DEB-IPM.

\newpage

## Section S5. Sensitivity analysis of long-term population growth rate ($\lambda$) to $\kappa_0$: Methods and analytical decomposition into contributions from demographic functions

### Methodological approach

Our sensitivity analysis aims to explain how $\kappa_0$ is optimized at different temperatures, i.e. how the 'configuration' of demographic functions at optimal values of $\kappa$ may shift with temperatures. To this end we decompose the sensitivity of $\lambda$ into the respective size dependent sensitivity contribution from the demographic functions for  size at age 1, somatic growth and fecundity at three temperatures (287, 289, 291 K). We use the chain rule of derivation based on numerical approximations of the derivatives for mass and temperature specific demographic functions with respect to $\kappa_0$. We use numerical approximation as the integrals in the DEB (that approximates the energy budget and consecutive growth and fecundity) are also numerical approximations. However, we also present analytical expressions for the sensitivity.

### Notation when deriving sensitivity

To describe this approach, we here use a simplified matrix notation instead of integrals (the underlying functions and analysis are the same as the main text DEBIPM). The expressions describing the derivation of the sensitivity below is thus a discrete version of the IPM with two main stages (eggs, and 'non-eggs', i.e. individuals with a given size, from age 1 and above), resulting in the following four main parts of the projection matrix $\mathbf{K}$:

```{=Latex}
\begin{align}\label{Kmat1}
\mathbf{K}&=\left[\begin{array}{c|c}
\text{Egg to egg}&\text{Reproduction}\\
\hline
\text{Egg to 'non-egg'}&  \text{Size transitions for 'non-eggs'}
\end{array}\right].
\end{align}
```

As all eggs transition to 1-year olds (or die) next year, the upper left element of the matrix is zero. The other parts of the matrix are defined by the DEBIPM, the discretized version assumes the model consists of $n+1$ stages, where the first stage is the egg stage and the last stages from 2 to $n+1$ are size stages. The term $g_{ij}$ is the probability of transitioning from size stage $j$ to size stage $i$ and is determined by the underlying DEB-function $g(m_{s+1};m_s,T)$. The  projection matrix is then given by

```{=Latex}
\begin{align}\label{Kmat2}
\mathbf{K}&=\left[\begin{array}{c|ccc}
0&f_2a_2&\cdots & f_{n+1}a_{n+1} \\
\hline
a_1g_{21}&g_{22}a_2&\cdots & g_{2,n+1}a_{n+1} \\
\vdots &\vdots &\ddots & \vdots \\
a_1g_{n+1,1}&g_{n+1,2}a_2&\cdots & g_{n+1,n+1}a_{n+1} \\
\end{array}\right].
\end{align}
```

Here, $g_{i1}$ is the probability that a surviving egg grows to stage $i$ at age 1 (determined by $o(m_{s+1},T)$ in the DEB-models) and is thus calculated differently from  $g_{ij}$ for $j>1$. However, they both represent the same transition process (growth transitions) and have the same symbol here $g_{ij}$. Similarly, $a_j$ is survival probability, given by the parameter $o_{surv}$ for $j=1$ and calculated by the $a(m_s, T)$ function for $j>1$.

The $\kappa_0$ parameter affects growth $g_{ij}$ and fecundity $f_j$ in stage $j$, while survival $a_j$ is not dependent on this parameter. 

### Sensitivity analysis

#### Decomposition of sensitivity of $\lambda$ to $\kappa_0$

Using this projection matrix $\mathbf{K}$ we can decompose the sensitivity $\frac{d \lambda}{d\kappa_0}$ (which should be zero at the optimum),  into contributions from the different vital rates across different stages. Applying the chain rule (and the fact that $\frac{d \lambda}{dK_{ij}}=v_iw_j$ where $v_i$ is the reproductive value or left eigenvector of the projection matrix), we get

```{=Latex}
\begin{align}\label{sensitkappa}
\frac{d \lambda}{d\kappa_0}&=\sum_{i=1}^{n+1}\sum_{j=1}^{n+1}\frac{d \lambda}{dK_{ij}}\frac{d K_{ij}}{d\kappa_0}\\
&=\sum_{i=1}^{n+1}\sum_{j=1}^{n+1}v_iw_j\frac{d K_{ij}}{d\kappa_0}\\
&=\sum_{j=2}^{n+1}v_1w_ja_j\frac{df_j}{d\kappa_0}+\sum_{i=2}^{n+1}v_iw_1a_1\frac{dg_{i1}}{d\kappa_0}+\sum_{i=2}^{n+1}\sum_{j=2}^{n+1}v_iw_ja_j\frac{dg_{ij}}{d\kappa_0}.
\end{align}
```

In the last equation, the first term corresponds to contributions from fecundity, the second to contributions from first year growth, and the third to contributions from growth of individuals age 1 and older.  The contributions from survival are zero because survival depends only indirectly on $\kappa_0$.

This decomposition requires the three sensitivities  $\frac{df_j}{d\kappa_0}$ ($j>1$), $\frac{dg_1j}{d\kappa_0}$ ($j>1$), and $\frac{dg_{ij}}{d\kappa_0}$ ($i,j>1$). These are calculated in the next section.  These derivations are on the scale of demographic functions that are defined on continuous scale, which requires that we switch to notation using continuous functions.

### Analytical expression for $\frac{dg_{ij}}{d\kappa_0}$

The growth function $g(m_{s+1};m_s,T)$ is assumed to be a lognormal distribution, where in our case the mean $\mu=\mu(m)$ depends on $\kappa_0$ (mean defined on the absolute scale) and the standard deviation $\sigma$ depends on the mean (and thereby also on $\kappa_0$), $\sigma^2(\mu)=\beta^2 e^{-2\nu \mu}$.  Note that in the above derivations $g_{ij}$ is the discretized version of this function, where for a given bin size $dy$ we get  $g_{ij}=g(m_{s+1};m_s,T)dy/\int g(m_{s+1};m_s,T)dy$ (and $\sum_i g_{ij}=1$). 

To simplify calculations, we can define the following function $h(\mu)$ and its derivative:

```{=Latex}
\begin{align}\label{lognormal}
h(\mu)&=1+\frac{\sigma^2}{\mu^2}=1+\frac{\beta^2 e^{-2\nu \mu}}{\mu^2}=1+\mu^{-2}\beta^2 e^{-2\nu \mu}\\
h'(\mu)&=-2\mu^{-3}\beta^2e^{-2\nu\mu}+e^{-2\nu\mu}(-2\nu)\beta^2\mu^2=-2\beta^2e^{-2\nu\mu}\mu^{-2}(\mu^{-1}+1).
\end{align}
```

For this model we can write the lognormal distribution as a function of $\mu=\mu(m)$:

```{=Latex}
\begin{align}\label{lognormal2}
g(m_{s+1},\mu) &=\frac{1}{m_{s+1}}\left[2\pi\ln\left(1+\frac{\sigma^2}{\mu^2}\right)\right]^{-\frac{1}{2}}\exp\left(-\frac{\left[\ln m_{s+1}-\ln\left(\mu\left[1+\frac{\sigma^2}{\mu^2}\right]^{-\frac{1}{2}}\right)\right]^2}{2\ln\left(1+\frac{\sigma^2}{\mu^2}\right)}\right)\\
&=\frac{1}{m_{s+1}}\left[2\pi\ln h(\mu) \right]^{-\frac{1}{2}}\exp\left(-\frac{\left[\ln m_{s+1}-\ln\left(\mu\left[h(\mu)\right]^{-\frac{1}{2}}\right)\right]^2}{2\ln h(\mu) }\right)\\
&=\frac{1}{m_{s+1}}g_a(\mu)\exp\left(-\frac{g_b(\mu)}{g_c(\mu)}\right),
\end{align}
```

where $g_a(\mu)=\left[2\pi\ln h(\mu) \right]^{-\frac{1}{2}}$, $g_b(\mu)=\left[\ln m_{s+1}-\ln\left(\mu\left[h(\mu)\right]^{-\frac{1}{2}}\right)\right]^2$, and $g_c(\mu)= 2\ln h(\mu)$.  The derivatives of these are

```{=Latex}
\begin{align*}
g_a'(\mu)&=-\frac{1}{2}\left[2\pi\ln h(\mu) \right]^{-\frac{3}{2}}2\pi[h(\mu)]^{-1}h'(\mu)=-\frac{\pi}{h(\mu)}\left[2\pi\ln h(\mu) \right]^{-\frac{3}{2}}h'(\mu),
\end{align*}
```

```{=Latex}
\begin{align*}
g_b'(\mu)&= 2\left[\ln m_{s+1}-\ln\left(\mu\left[h(\mu)\right]^{-\frac{1}{2}}\right)\right]\frac{-1}{\mu\left[h(\mu)\right]^{-\frac{1}{2}}}\left[\left[h(\mu)\right]^{-\frac{1}{2}}-\frac{\mu}{2}\left[h(\mu)\right]^{-\frac{3}{2}}h'(\mu)\right]\\
&= 2\left[\ln m_{s+1}-\ln\left(\mu\left[h(\mu)\right]^{-\frac{1}{2}}\right)\right]\left[ \frac{h'(\mu)}{2h(\mu)}-\frac{1}{\mu}\right],
\end{align*}
```

and

```{=Latex}
\begin{align*}
g_c'(\mu)&= \frac{2h'(\mu)}{h(\mu)}.
\end{align*}
```

Using the chain rule, the derivative of $g(m_{s+1},\mu)$ with respect to $\kappa_0$  is given by

```{=Latex}
\begin{align}\label{dG.dkappa2}
\frac{d g(m_{s+1},\mu )}{d \kappa_0}&=\frac{d g(m_{s+1},\mu)}{d \mu}\frac{d \mu}{d \kappa_0}
\end{align}
```

The first term is given by

```{=Latex}
\begin{align}\label{dg.dmug2}
\frac{d g(m_{s+1},\mu;\sigma)}{d\mu} &=\frac{1}{m_{s+1}}g'_a(\mu)\exp\left(-\frac{g_b(\mu)}{g_c(\mu)}\right)+\frac{1}{m_{s+1}}g_a(\mu)\exp\left(-\frac{g_b(\mu)}{g_c(\mu)}\right)\left(\frac{-g_b'(\mu)g_c(\mu)+g_c'(\mu)g_b(\mu)}{g_b(\mu)^2}\right). 
\end{align}
```

This can be found by plugging in the expressions defined above. The second term of eq. \eqref{dG.dkappa2} is the derivative of the mean with respect to $\kappa_0$, and is given by

```{=Latex}
\begin{align}\label{dmu.dkappa2}
\frac{d \mu}{d\kappa_0} &=\frac{d}{d\kappa_0}\int_0^{t_l}[\kappa_0e^{-\frac{m}{20000\kappa_m}}\alpha C(m,T)-M(m,T)]dt\\
&=\int_0^{t_l}\alpha C(m,T)e^{-\frac{m}{20000\kappa_m}}dt.
\end{align}
```


### Analytical expression for $\frac{dg_{1j}}{d\kappa_0}$

The distribution of size at age 1 follows a lognormal distribution similar to \eqref{lognormal2}, but the mean and variance are defined differently for first year growth. Here the variance is a constant, and the integral defining the mean is over the interval 0 to $t_l/2$.  Thus, the derivative of $\mu$ with respect to $\kappa_0$ is given by

```{=Latex}
\begin{align}\label{dmu.dkappa2}
\frac{d \mu(m_s)}{d\kappa_0} &=\frac{d}{d\kappa_0}\int_0^{t_l/2}[\kappa_0e^{-\frac{m}{20000\kappa_m}}\alpha C(m,T)-M(m,T)]dt\\
&=\int_0^{t_l/2}\alpha C(m,T)e^{-\frac{m}{20000\kappa_m}}dt.
\end{align}
```

In this case we get $h(\mu)=1+\frac{\sigma^2}{\mu^2}$ and $h'(\mu)=-2\sigma^2\mu^{-3}$, which can be plugged into all of the above expressions calculated for the growth sensitivity  (with correct values of $\mu$ and $\sigma$).  

### Analytical expression for $\frac{df_j}{d\kappa_0}$

Fecundity is given by (assuming enough energy available for reproduction)

```{=Latex}
\begin{align}\label{fecundity}
f(m_s,T)&=\frac{1}{2e_m}\int_0^{t_l}\frac{dR}{dt}dt=\frac{1}{2e_m}\int_0^{t_l}(1-\kappa_0e^{-\frac{m}{20000\kappa_m}}\alpha C(m,T))dt.
\end{align}
```

The sensitivity with respect to $\kappa_0$ is thus given by

```{=Latex}
\begin{align}\label{fecundity}
\frac{df(m_s,T)}{d\kappa_0}&= -\frac{1}{2e_m}\int_0^{t_l}e^{-\frac{m}{20000\kappa_m}}\alpha C(m,T))dt.
\end{align}
```

\newpage

## Section S6. Growth trajectories and fecundity at optimal energy allocation ($\kappa_0^*$)

Size dependent growth and fecundity for optimal energy allocation to somatic growth $(\kappa_0^*)$ differs from those presented in Fig. 3 in the main text (i.e. with $\kappa_0=0.89$ which was used to fit the model). Compared to $\kappa_0=0.89$, optimal allocation in the Windermere pike DEBIPM (main text Fig. 2) results in decreased size specific fecundity (Fig. S5b) and reduced growth  (mean size at age, Fig. S5A) for temperatures above 287 K, and increased growth below or equal to 287 K.  compared to $\kappa_0=0.89$. 

```{r, echo = FALSE, warning = FALSE, fig.height = 3, fig.width = 6}
OptPars <- 
  as.data.frame(maxl_1[,1:3]) %>%
  filter(T %in% c(285,287,289,291,293))

trajO <- data.frame()
for(i in 1:nrow(OptPars)){
  trajO[1,i] <- x[which.max(age1size(Pars = unlist(OptPars[i,])))]
  for(j in 2:t){ # 25 years of growth
    trajO[j,i] <- x[which.max(growthfun(trajO[j-1,i], ratefun(trajO[j-1,i], unlist(OptPars[i,])), Pars = unlist(OptPars[i,])))]
  }
}
colnames(trajO) <- OptPars[1:5,1]
trajO <- cbind(Age=1:nrow(trajO),trajO)
trajGOpt <- pivot_longer(trajO, cols=c(2:6),names_to = "Temp", values_to = "mass" )
trajGOpt["Allocation"] <- "Optimal"
trajG["Allocation"] <- "0.89"
trajBoth <- rbind(trajGOpt,trajG)

groTopt <- 
  trajBoth %>%
  ggplot(.) +
  geom_line(size=0.7,  aes(Age, mass, color = as.factor(Temp), linetype = as.factor(Allocation))) + #alpha = strat, 
  scale_color_manual(values = c('#4575b4','#91bfdb','#fee090','#fc8d59','#d73027'), name="Temp [K]")+
  scale_linetype_manual(name="Allocation", values = c("solid", "dashed")) +
  ylab(expression(paste("Mass ", italic("\u03BC")["g"]," [g]"))) +
  xlab("Age [years]") +
  xlim(0,20)+
  ylim(0,max(trajBoth$mass)) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.position="none",
        legend.key.width = unit(1, 'cm')) 

# Plot optimal fecundity from DEB and in Lake Windermere. Note that we need to double (*2) the fecundity when plotting using repfun(). Repfun account only for female eggs.
repO=data.frame()
for(i in 2:ncol(trajO)){ # the temps, from 2nd to 6th column in trajO
  for(j in 1:nrow(trajO)){ # the sizes in trajo, i.e. Age trajO
    if(trajO$Age[j]==1){
      rep_tplus1 = cbind(trajO[j,i], repfun(e_m, ratefun(e_m, unlist(OptPars[i-1,])), unlist(OptPars[i-1,]))*2, OptPars[i-1,] )} #use the OptPars corresponding to temps in columns in trajO
    else{ rep_tplus1 = cbind(trajO[j,i], repfun(trajO[j-1,i], ratefun(trajO[j-1,i], unlist(OptPars[i-1,])), unlist(OptPars[i-1,]))*2, OptPars[i-1,] )} #j-1 because last years size predicts fecundity this year
    colnames(rep_tplus1) <- c("mass","fec","Temp", "Kappa", "Y")
    repO <- rbind(repO, rep_tplus1)
  }
}

repO["Allocation"] <- "Optimal"
rep["Allocation"] <- "0.89"
repBoth <- rbind(repO,rep)

repTopt <- 
ggplot() +
  geom_line(data=as_tibble(repBoth), size=0.85, aes(mass, fec, color = as.factor(Temp), linetype= as.factor(Allocation))) +
  scale_color_manual(values = c('#4575b4','#91bfdb','#fee090','#fc8d59','#d73027'), name="Temp [K]") +
  scale_linetype_manual(name="Allocation", values = c("solid", "dashed")) +
  ylab(expression(paste(italic("f"),"(",italic("m"["s"]),italic(",T"),")"))) +
  xlab(expression(paste("Mass ",italic("m")[s]," [g]"))) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(plot.title = element_text(size = 9, face = "bold"),
        legend.key.width = unit(1, 'cm')) 

groTopt + repTopt
```

Figure S6. Mean size at age (a) and size specific fecundity (b) predicted from the optimal allocation strategy $(\kappa_0^*)$, Fig. 2 main text) at five temperatures and contrasted by growth and fecundity for $\kappa_0=0.89$ (i.e. $\kappa_0$ used for fitting the model).

\newpage

## Section S7. Fitness landscape and stable size structure for survival probability contrasts

In the main analyses, we contrast the size- and temperature dependent survival used in our model using two alternative survival probability scenarios (Fig. 4c main text). Here we present the corresponding fitness landscape ($\lambda$), $\kappa_0^*$ (Fig. S7a,b) and stable structure $w$  over temperature (Fig. S7c,d) for the temperature independent survival contrast $a(m_s)$ (Fig. S7a,c) and constant survival $a=0.68$ (Fig. S7b,d). In the first scenario, we use temperature independent survival $a(m_s)$ that corresponds to survival at mean temperature (287 K) in Windermere where temperature does not affect the survival asymptote for mass ($a(m_s)$, (Fig. 4d). In the second scenario we use constant (both size- and temperature-independent) survival probability that equals this asymptote at mean temperature (a=0.68, Fig. 4e).

```{r, echo = FALSE, warning = FALSE, fig.height = 8, fig.width = 8}
mycolors <- colorRampPalette(colors=c("#feebe2","#ce1256"))(9)

maxl_2 <- as.data.frame(Res_lam_2) %>% 
  group_by(T) %>%
  slice_max(Lambda)
maxl_2["Sur_f"] <- "NotempVin"

FigS6a <- 
  as_tibble(Res_lam_2) %>%
  filter(kappa > 0.65) %>%
  filter(T < 292.25 & T > 283.5) %>%
  ggplot(., aes(T,kappa)) +
  geom_raster(aes(fill=Lambda)) +
  geom_line(data=maxl_2,aes(T,kappa),size=0.8) +
  scale_fill_gradientn( 
    limits = c(min(as.data.frame(Res_lam_2[Res_lam_2$kappa>0.65 & Res_lam_2$T<292.25, ][,'Lambda'])),max(as.data.frame(Res_lam_2[,'Lambda']))),
    colours = c("black","white","#f1eef6","#ce1256"),
    values = rescale(c(min(as.data.frame(Res_lam_2[Res_lam_2$kappa>0.65 & Res_lam_2$T<292.25, ][,'Lambda'])),0.9999,1, max(as.data.frame(Res_lam_2[,'Lambda']))),
                     to=c(0,1)),
    breaks = c(min(as.data.frame(Res_lam_2[Res_lam_2$kappa>0.65 & Res_lam_2$T<292.25, ][,'Lambda'])),1,max(as.data.frame(Res_lam_2[,'Lambda']))),
    labels=c(round(min(as.data.frame(Res_lam_2[Res_lam_2$kappa>0.65 & Res_lam_2$T<292.25, ][,'Lambda'])), 2),
             1, round(max(as.data.frame(Res_lam_2[,'Lambda'])),2)),
    name = expression(lambda)) +
  scale_x_continuous(expand = c(0,0), name = "Temperature [K]", breaks = seq(283,291,2)) +
  scale_y_continuous(expand = c(0,0), name=expression(paste(kappa[0]," (Growth allocation)"))) +
  coord_cartesian(xlim = c(283.75, 292)) +
  ggtitle('(a)') +
  theme_light(base_size=8) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())
  
Res_w_2 <- read.delim(text = getURL("https://raw.githubusercontent.com/VThunell/Allocation-Warming-DEB-IPM/main/data/Res_w_2_2022-10-07.txt"), sep=",")
colnames(Res_w_2)[4:ncol(Res_w_2)] <- sub("X", "", colnames(Res_w_2)[4:ncol(Res_w_2)])

Res_w_2_long <- pivot_longer(as.data.frame(Res_w_2), cols = c(4:ncol(Res_w_2)),
                           names_to = "Size", values_to ="biom")

FigS6c <- 
  as_tibble(Res_w_2_long) %>%
  filter(kappa %in% round(deb.optim$par[1],2)) %>%
  mutate(biom_log=log(biom)) %>% 
  ggplot(., aes(as.double(Size),T, z=biom_log, colour=biom_log))+
  geom_contour_filled(breaks = c(min(log(Res_w_2_long$biom)), 
 log(1e-10),log(1e-9),log(1e-8),log(5e-8),log(1e-7),log(5e-7),log(1e-6),
                                 max(log(Res_w_2_long$biom)))) + 
  scale_fill_manual(values = mycolors, name="Relative densities", 
                     labels=c(paste("<",1e-10), 
                              1e-9,1e-8,5e-8,1e-7,5e-7,1e-6,
                              paste("<",1))) + 
  scale_y_continuous(expand = c(0,0), name = "Temperature [K]", breaks = seq(283,291,2)) +
  scale_x_continuous(expand = c(0,0), name = expression(paste("Mass ",italic("m"["s"])," [g]")), limits = c(0,20000))+
  coord_cartesian(ylim = c(283.75, 292)) +
  ggtitle('(c)') +
  theme_light(base_size=8) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

maxl_3 <- as.data.frame(Res_lam_3) %>% 
  group_by(T) %>%
  slice_max(Lambda)
maxl_3["Sur_f"] <- "flat"

FigS6b <- 
  as_tibble(Res_lam_3) %>%
  filter(kappa > 0.65) %>%
  filter(T < 292.25 & T > 283.5) %>%
  ggplot(., aes(T,kappa)) +
  geom_raster(aes(fill=Lambda)) +
  geom_line(data=maxl_3,aes(T,kappa),size=0.8) +
  scale_fill_gradientn(
    limits = c(min(as.data.frame(Res_lam_3[Res_lam_3$kappa>0.65 & Res_lam_3$T<293.5, ][,'Lambda'])),max(as.data.frame(Res_lam_3[,'Lambda']))),
    colours = c("#f1eef6","#ce1256"),
    values = rescale(c(min(as.data.frame(Res_lam_3[Res_lam_3$kappa>0.65 & Res_lam_3$T<293.5, ][,'Lambda'])), max(as.data.frame(Res_lam_3[,'Lambda']))),
                     to=c(0,1)),
    breaks = c(min(as.data.frame(Res_lam_3[Res_lam_3$kappa>0.65 & Res_lam_3$T<293.5, ][,'Lambda'])),max(as.data.frame(Res_lam_3[,'Lambda']))),
    labels=c(round(min(as.data.frame(Res_lam_3[Res_lam_3$kappa>0.65 & Res_lam_3$T<293.5, ][,'Lambda'])), 2), round(max(as.data.frame(Res_lam_3[,'Lambda'])),2)),
    name = expression(lambda)) +
  scale_x_continuous(expand = c(0,0), name = "Temperature [K]", breaks = seq(283,291,2)) +
  scale_y_continuous(expand = c(0,0), name=expression(paste(kappa[0]," (Growth allocation)"))) +
  coord_cartesian(xlim = c(283.75, 292)) +
  ggtitle('(b)') +
  theme_light(base_size=8) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())
  

# a=0.68, i.e flat survival scenario
Res_w_3 <- read.delim(text = getURL("https://raw.githubusercontent.com/VThunell/Allocation-Warming-DEB-IPM/main/data/Res_w_3_2022-10-07.txt"), sep=",")
colnames(Res_w_3)[4:ncol(Res_w_3)] <- sub("X", "", colnames(Res_w_3)[4:ncol(Res_w_3)])

Res_w_3_long <- pivot_longer(as.data.frame(Res_w_3), cols = c(4:ncol(Res_w_3)),
                           names_to = "Size", values_to ="biom") 

FigS6d <- 
  as_tibble(Res_w_3_long) %>%
  filter(kappa %in% round(deb.optim$par[1],2)) %>%
  mutate(biom_log=log(biom)) %>% 
  ggplot(., aes(as.double(Size),T, z=biom_log, colour=biom_log))+
  geom_contour_filled(breaks = c(min(log(Res_w_3_long$biom)), 
 log(1e-10),log(1e-9),log(1e-8),log(5e-8),log(1e-7),log(5e-7),log(1e-6),
                                 max(log(Res_w_3_long$biom)))) + 
  scale_fill_manual(values = mycolors, name="Relative densities", 
                     labels=c(paste("<",1e-10), 
                              1e-9,1e-8,5e-8,1e-7,5e-7,1e-6,
                              paste("<",1))) +  
  scale_y_continuous(expand = c(0,0), name = "Temperature [K]", breaks = seq(283,291,2)) +
  scale_x_continuous(expand = c(0,0), name = expression(paste("Mass ",italic("m"["s"])," [g]")), limits = c(0,20000))+
  coord_cartesian(ylim = c(283.75, 292)) +
  ggtitle('(d)') +
  theme_light(base_size=8) +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

(FigS6a/FigS6c) | (FigS6b/FigS6d)
```

Figure S7. Fitness landscape (a,b) with solid line representing the optimal allocation strategy ($\kappa_0^*$), and the stable structure $w$ over temperature (c,d) for survival contrast $a(m_s)$(a,c) and $a=0.68$ (b,d).

## References

Byrd, R. H., Lu, P., Nocedal, J. and Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing, 16, 1190–1208. doi:10.1137/0916069.

Lindmark, M., J. Ohlberger, and A. Gårdmark. 2022. Optimum growth temperature declines with body size within fish species. Glob Chang Biol 28:2259-2271.

Winfield, I. J., and J. M. Fletcher. 2013. Windermere lake temperature data 1944-2002. . NERC Environmental Information Data Centre. https://doi.org/10.5285/9520664c-eb4d-4700-b064-5d215d23e595.

Winfield, I. J., J. M. Fletcher, and J. B. James. 2013a. Pike fecundity data 1963-2002., NERC Environmental Information Data Centre. https://doi.org/10.5285/b8886915-14cb-44df-86fa-7ab718acf49a.

Winfield, I. J., J. M. Fletcher, and J. B. James. 2013b. Pike growth data 1944-1995. NERC Environmental Information Data Centre. https://doi.org/10.5285/637d60d6-1571-49af-93f7-24c1279d884d.