Identification of temperatures and diets maximizing growth at every week of development for the tropical banded cricket (Gryllodes sigillatus)

Introduction

Background

In the context of climate change, the sustainability of animal production like beef, lamb and pork is increasingly criticized as too polluting and needing high input resources. A new category of agricultural product that compares favorably to traditionally farmed animals are insects. They require little resources to produce an end product with high nutritional value while having a small environmental footprint. However, industrial production of insects is still in its infancy. Production methods need to be ameliorated to lower production costs (Halloran et al., 2018; Anderson et al., 2020).

Crickets are in this category of insects with promising future as food and feed that has not perfected its production process. The current technique to produce crickets sets a stable temperature throughout the growing period. That temperature usually has been confirmed to maximize mass at adulthood. Other industrially produced animals, such as chickens, are not raised under stable temperature, but a temperature that maximizes growth at every life stage. This results in a decrease in temperature with time leading to overall reduced temperature during production. Similarly for diets, the earlier stages are fed higher percentage of protein and it decreases with time. For cricket farming, the diet is set for their whole development with a stable diet formulation. There could be room to improve the farming of crickets if we could identify how performance, here growth, is influenced by temperature and diets at every week of development (NRC, 1994; May, 1998).

Above is the tropical banded cricket (Gryllodes sigillatus).

Objective

The objective is to identify if specific temperatures and diet formulations maximize growth for every week of development and if so which. If such parameters are found, create a new profile to continuously maximize growth during development.

Method

To achieve this objective, the tropical banded cricket (Gryllodes sigillatus) was reared at 6 temperatures (26, 29, 32, 35, 38˚C) and fed 7 diet formulation of protein to carbohydrate (4:1, 3:1, 2:1, 1:1, 1:2, 1:4 and 1:6; P:C), separately. They were reared from hatch in individual containers with feed, water and shelter. The animals were assessed once a week to check for weight, development (instar number) and survival status. The experiment was stopped at 6 weeks, emulating a production cycle.

Software used

To analyze these data, the open source software R was used. R is a coding language that uses different packages of functions. To analyse the data, I used a series of packages essential to my analysis.

Data Import, Hygiene and Tidying

Import

The first step in the analysis is to import the data. There are two data frames needed to run the upcoming code.

The first one has the information on the individual crickets. Their individual number (cricket), the incubator they are reared in, the temperature of said incubator and the diet number fed to the cricket. crick contains 390 rows representing the 390 crickets that were initially used in the experiment. Each treatment had 30 crickets initially except 32-1:1 and 1:2 with 45 crickets.

Rows: 390
Columns: 4
$ cricket     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,…
$ incubator   <dbl> 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32…
$ temperature <dbl> 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26…
$ diet        <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4…

The second data set contains the the information taken every week on the animals. The columns represent each individual observation, individual cricket number, the day of the experiment, the week of the experiment, the sex/maturity of the animal, the development stage (instar), the empty container weight and the weight of the container with the cricket (cricket_container) inside. It contains 2382 rows representing the measurements of the 390 crickets every week for 6 weeks, minus possible mortalities.

Rows: 2,463
Columns: 8
$ observation       <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
$ cricket           <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1…
$ day               <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ week              <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ sex               <chr> "juvenile", "juvenile", "juvenile", "juvenile", "juv…
$ instar            <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ container         <dbl> 0.53552, 0.53544, 0.53549, 0.53544, 0.53549, 0.53546…
$ cricket_container <dbl> 0.53641, 0.53640, 0.53628, 0.53610, 0.53649, 0.53651…

Hygiene Check

I then want to check if the different values used in the data sets are the values I expect.

First for crick

Cricket	Range	Incubator	Range	Temperature	Range	Diet	Range
390	[1, 390]	11	[31, 41]	6	[26, 41]	7	[1, 7]

Second for sig

Observation	Range	Cricket	Range	Day	Range	Week	Range	Instar	Range
2463	[1, 2507]	390	[1, 390]	7	[0, 42]	7	[0, 6]	8	[1, 8]

The number of values expected and range for each column is right for the cricket description and needs no corrections.

The weekly mass measurements on the crickets changed with time. To check the data, we will need to separate it into the different days of measurements

Day	Container	Range	Cricket+Container	Range
0	390	[0.53543, 53558]	390	[0.5361, 53639]
7	373	[0.53302, 0.53661]	373	[0.53663, 0.54771]
14	360	[0.53575, 6.38875]	360	[0.53826, 639855]
21	351	[6.3824, 638722]	351	[6.3907, 6.64018]
28	332	[6.38245, 6.38649]	332	[6.39203, 6.78781]
35	330	[6.38478, 6.38914]	330	[6.4025, 6.90416]
42	327	[6.382, 6.38967]	327	[6.41362, 6.95444]

The number of observation only decreases showing moralities during the experiment which is normal. The range of values for the mass of the container and the container and cricket are not as expected for every week. Some of those values are too high. Having run the experiment, I know that no values measured were above 7g. By isolating any values over 7g, we may be able to identify problematic values.

Observation	Cricket	Day	Week	Sex	Instar	Container	Cricket+Container
86	86	0	0	juvenile	1	5.35530e-01	5.36390e+04
2045	319	14	2	juvenile	4	6.38424e+00	6.39855e+05
33	33	0	0	juvenile	1	5.35580e+04	5.36490e-01
891	96	21	3	female	6	6.38722e+05	6.47354e+00
2163	346	21	3	female	6	6.38351e+05	6.47673e+00

Above are the values identifier out of the expected range. There are two in the Container column and 3 in the Cricket+Container column. The error made is rapidly obvious when comparing to other values. When the data was entered, the comma was forgotten. To fix this we need to find those values and then divide them by 10000 to get the mass in grams.

Code

sig <- sig %>% 
  mutate(container = if_else(container > 7, container / 100000, container)) # Divides the values that are too high by 100000 to replace the comma at the right place.
sig <- sig |>
  mutate(cricket_container = if_else(cricket_container > 7, cricket_container / 100000, cricket_container))# Divides the values that are too high by 100000 to replace the comma at the right place.

# Let's check if it corrected our issue.
a <- sig |>
  select(day, container, cricket_container) |>
  get_summary_stats_by_day()
# All values are now in the expected range and we can continue with our analysis.

colnames(a) <- c("Day", "Container", "Range", "Cricket+Container", "Range")

kable(a)

Day	Container	Range	Cricket+Container	Range
0	390	[0.53543, 0.53616]	390	[0.5361, 0.5371]
7	373	[0.53302, 0.53661]	373	[0.53663, 0.54771]
14	360	[0.53575, 6.38875]	360	[0.53826, 6.46507]
21	351	[6.3824, 6.38976]	351	[6.3907, 6.64018]
28	332	[6.38245, 6.38649]	332	[6.39203, 6.78781]
35	330	[6.38478, 6.38914]	330	[6.4025, 6.90416]
42	327	[6.382, 6.38967]	327	[6.41362, 6.95444]

The number of values for each columns and the range match expected values and then require no corrections.

An important check to make is to insure a cricket id number does not come up twice in a week. It would indicate a mistake in the data entry or source data.

Code

grouped_data <- sig %>%
  group_by(week) %>% # For every week of the experiment
  filter(duplicated(cricket)) # Filter for duplicated values in the 'cricket' column

num_duplicates <- nrow(grouped_data)

cat("There are", num_duplicates, "duplicated values.\n")

There are 0 duplicated values.

To be able to continue data cleaning, we need to join the two data sets to be able to have the information on the treatment of every cricket. Additionally, protein to carbohydrate ratios are not easy to deal with outside of factors. To be able to use the different diets as factor and numerical, we will give all the diets a percentage of protein here named “dietperc”.

Code

# Joining the two data sets by cricket to create the data set that we will use in this script.
sig <- full_join(sig, crick, by = "cricket")

# Protein to carbohydrate ratios are not easy to deal with outside of factors. To be able to use the different diets as factor and numerical, we will give all the diets a percentage of protein here named "dietperc".
sig <- sig |>
mutate(dietperc = case_when(
  diet == "1" ~ 62.4,
  diet == "2" ~ 58.8,
  diet == "3" ~ 51.9,
  diet == "4" ~ 38.3,
  diet == "5" ~ 25.9,
  diet == "6" ~ 16.6,
  diet == "7" ~ 11.8
))

Values of Interest

At this point, our data frame only have the mass of the empty container and the mass of the container and cricket combined. We won’t be able to go far with just that. We now need to generate some values from the base measurements. To achieve this, we will need to do some calculations.

Calculations

The mass of the individual cricket is the main value to calculate. For every observation, we take the values of container subtract it to cricket_container and then multiply it by 1000 to get mass in milligram. From mass, we will be calculating two new values.

The first is daily growth for every week of the experiment. It is the mass added to a cricket between week x and week x+1. It is then divided by 7 for the number of days in a week. This way, we get a number that would be comparable with studies measuring crickets at different intervals. It is calculated as: growth = (mass - lag(mass)) / 7) .

The second value is proportional growth for every week of the experiment. It is the mass added over a week as a percentage of the mass from the previous week. It is calculated as: propgrowth = (mass - lag(mass)) / lag(mass) *100) .

Verification of Distributions

To make sure the three new values calculated have a proper distribution, let’s graph them by separating by the two trials (temperature and diet) then temperature or diet and day.

Figure 1. Distribution of calculated values for the temeprature section of the experiment

Figure 2. Distribution of calculated values for the diet section of the experiment

The distribution of the new values calculated appear good. The values for crickets at 41˚C are not visible after day 21 because of a mortality rate above 50%. The values were deemed unfit for analysis below that threshold.The data is clean and tidy. We can continue to the data visualization.

Data visualization

We will now visualize the three different values calculated above to observe how they are affected by temperature and diet over the development of each crickets. In all three figures below, the values are averaged by week and treatment. The error bars represent the standard error for each of those means.

Mass

All crickets started the experiment at the same time, with similar mass. Cricket mass increased with time at all temperature and diets. Final mass was affected by temperature with the smallest recorded at 26˚C and largest at 35˚C. Crickets fed 1:1 and 1:2 proteins to carbohydrates had the largest mass at 6 weeks. Crickets fed higher protein contents ended the experiment smaller. There were then different rates of mass increase at different temperature and diet.

(a) For temperatures between 26 and 41˚C

Figure 3. Average mass for the 6 first weeks of development

Daily Growth

The rate of mass increase is the growth rate, here daily growth. Daily growth increases with time until it peaks and then decreases. This is visible for all temperatures and diets except 26˚C. Daily growth peaked for crickets held at 32, 35 and 38˚C and 29˚C on day 28 and 35 respectively. The highest values of daily growth we measured under diets of 1:1 and 1:2 proteins to carbohydrates. All dietary treatments peaked on day 28. Growth rates are however imperfect. They are only an account of the absolute weight gained per day and does not take into account the initial weight of the animal for the desired period.

Figure 4. Average daily growth for the 6 first weeks of development

Proportional Growth

To take into account the mass of the animal at the beginning of the analyzed period, we need to look at the animal’s proportional growth. Similar to daily growth, proportional growth increases, peaks and than decreases for all temperature except 26 and 41˚C that only decreases with time. Differences in proportional growth are noticeable only on day 7 where it is maximum for diet of 1:4 and minimum for diet of 4:1 proteins to carbohydrates. The peak of proportional growth happens earlier than the one for growth rate at day 14 instead of 28. After day 14, values of proportional growth decrease towards zero.

Figure 5. Average proportional growth for the 6 first weeks of development

Data Wrangling

Temperature

Test for Significance

Since our goal is to identify conditions maximizing growth throughout ontogeny. We first need to identify if significant differences can be identified in growth between the different temperatures to support the visual differences we see on the figures.

Code

# Now for temperature, let's start with a model
sig_temp <- sig |>
  filter(diet == "4") |>
  mutate(week = as.factor(week))
# We will modelize proportional growth, as it is was visualized earlier to be maximized at decreasing temperatures.
sig_temp_model <- lm(propgrowth ~ temperature + week + temperature:week, sig_temp)
anova(sig_temp_model)

Analysis of Variance Table

Response: propgrowth
                   Df   Sum Sq Mean Sq F value    Pr(>F)    
temperature         1  1715719 1715719 119.826 < 2.2e-16 ***
week                6 31398424 5233071 365.477 < 2.2e-16 ***
temperature:week    6  2124625  354104  24.731 < 2.2e-16 ***
Residuals        1198 17153541   14318                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Temperature, week and their interaction are identified as significant in the anova (p < 2.2e-16). We then need to look at the different week separately. To do that, we will perform pairwise comparisons for each week to see which have post-hoc differences between temperatures.

Significant difference between temperatures in the proportional growth are identified at every week. We now need to identify the temperature that maximizes growth for each of them.

Weekly Optimum

To get the optimal temperature maximizing proportional growth at every week, we will use thermal performance curves. The rTPC package contains different models of thermal performance curves to modelize how temperature impacts biological rates (Padfield et al., 2021). The four models chosen (see code) were selected after trials for best fits. Below, we have a series of loops that will first separate the data per week and store them in a list. It then uses a loop in a loop to modelize the weekly proportional growth data and extract the optimal temperature. The average thermal optimum for proportional growth for every week from the 4 models can be found in the table below.

Code

# We first create empty lists
modelw_list <- list() # List to be filled with the data separated by week.
augmented_data_list <- list() # List to store the augmented data from the models.


# Loop to create and process each modeldw data frame from each week of the experiment
for (i in 1:6) {
  modeldw <- sig |>
    mutate(day = factor(day)) |>
    subset(diet == "4") |>
    filter(!(week == "0")) |>
    filter(week == as.character(i)) |> # This will filter weeks from 1 to 6.
    mutate(temperature = as.numeric(as.character(temperature)))
  
  # Store the data frame in the list
  modelw_list[[i]] <- modeldw # This will store data of weeks 1 to 6 in a list in position 1 to 6.
}


# The next rows will provide info on the models to run in the next loop.
models_names <- c("weibull_1995", "thomas_2017", "briere2_1999", "boatman_2017")
params <- c("temp=temperature, a, topt, b, c", "temp=temperature, a, b, c, d, e", "temp=temperature, tmin, tmax, a, b", "temp = temperature, rmax, tmin, tmax, a, b")
iters <- c("4,4,4,4", "3,3,3,3,3", "4,4,4,4", "3,3,3,3,3")

# Create a list to store fit results
fit_list <- list()

# Create a list to store summary data frames
summary_list <- list()

# Loop through each model
for (i in 1:4) { # The first section of the loop prepares the information to be fed to the modelizing function nls_multstart()
  model_name <- models_names[i]
  param <- strsplit(params[i], ",")[[1]]  # Separate params string
  iter <- as.numeric(strsplit(iters[i], ",")[[1]])  # Convert iter to a numeric vector
  
  # Fit the model using nls_multstart
  fit_list[[model_name]] <- list() # Creates a new level in the list fit_list to store this model's outputs.
  for (j in 1:6) { # This second loop modelize the data of the 6 weeks and 4 different models stated above.
    fit_list[[model_name]][[j]] <- nls_multstart(
      formula = as.formula(paste("propgrowth ~", model_name, "(", paste(param, collapse = ","), ")")),
      data = modelw_list[[j]],
      iter = iter,
      start_lower = get_start_vals(modelw_list[[j]]$temperature, modelw_list[[j]]$propgrowth, model_name = model_name) - 10,
      start_upper = get_start_vals(modelw_list[[j]]$temperature, modelw_list[[j]]$propgrowth, model_name = model_name) + 10,
      lower = get_lower_lims(modelw_list[[j]]$temperature, modelw_list[[j]]$propgrowth, model_name = model_name),
      upper = get_upper_lims(modelw_list[[j]]$temperature, modelw_list[[j]]$propgrowth, model_name = model_name),
      supp_errors = 'Y',
      convergence_count = FALSE
    )
    # Generate augmented data using broom::augment. Augmented data create a data frame of extrapolated values from the model created. This allows to create figures down the line.
    new_data <- data.frame(temperature = seq(25, 42, 0.1))
    augmented_data <- augment(fit_list[[model_name]][[j]], newdata = new_data)
    augmented_data <- mutate(augmented_data, week = j)
    augmented_data <- mutate(augmented_data, model = model_name)
    
    # Store augmented data in the list
    augmented_data_list[[model_name]][[j]] <- augmented_data
  }
  
  # After every model, the second loop stop and the first loop has a last action before it goes again with the next model. This step puts in a list of data frames the name of the model, the optimal temperature for every week identified by said model (using get_topt), the AICc score (necessary to rank models) and the week number.
  summary_list[[model_name]] <- tibble(
    model = rep(model_name, 6),
    topt = sapply(fit_list[[model_name]], function(fit) get_topt(fit)),
    AICc = sapply(fit_list[[model_name]], AICc),
    week = 1:6
  )
}

# To continue analysis, we need to combine the summary list of all model operations into one data frame.
combined_summary <- bind_rows(summary_list)

# This allows for the creation of a new column "delta" to rank models for next steps. 
combined_summary <- combined_summary|>
  group_by(week)|>
  arrange(AICc)|>
  mutate(delta= -(AICc - AICc[1]))|>
  ungroup()|>
  group_by(model)
# This table is the one that will be printed in the htlm document.
q <- combined_summary

q <- q |>
  select(topt, week)|>
  group_by(week) |>
  summarise("Optimal temperature" = round(mean(topt), 2))
  
colnames(q) <- c("Week", "Optimal temperature")

kable(q)

Week	Optimal temperature
1	37.34
2	35.86
3	33.56
4	29.21
5	26.41
6	26.78

Plot of Optimum

To Image the table above, we will plot all the optimal temperatures for every week of development identified by the four different models. We will plot those optimums through time to see if a relationship can be identified between temperature and time. Each point is an optimal temperature identified by the model. The size of the point represents its weight according to AICc scores (indicate how well the model fits the data). The larger the point the larger the weight.

Figure 6. Optimal temeprature at every week of development according to four models.

Continuous Optimization

We see a clear decrease in the optimal temperature to maximize proportional growth with time. The goal of this research if to create new temperature profiles to raise crickets to maximize growth. Here I propose to use an inverse logistic regression to modelize the temperature profile.

Code

# We then modelize the data using an inverse logarithm function using nls and store it as "opttemp". This function was chosen after trials with different functions (linear, exponential decay). To lighten the script, those were not included and only the one with the best fit was chosen.
opttemp <- nls(topt ~ a + b / (1 + exp(c * (week - d))), 
               data = combined_summary, 
               start = list(a = 22, b = 10, c = 0.5, d = 4),  # Initial parameter values
               weights = delta,
               algorithm = "port")

# We create a new data frame to extrapolate the data from the model created above.
new_data_opttemp <- data.frame(week = seq(0, 6, 0.5))
prev <- augment(opttemp, newdata = new_data_opttemp)


# With the augmented data, we create a figure with the extrapolated model and the data from all the optimal temperatures for the 4 different models.
ggplot(combined_summary, aes(week, topt)) +
  geom_point(aes(color=model, size=delta)) +
  geom_line(aes(week, .fitted),prev, col = 'blue') +
  theme_classic() +
  theme(
    axis.title = element_text(size = 18),
    legend.box.margin = margin(-15),
    strip.placement = "left",
    axis.text.x = element_text(size = 12, margin = margin(t = 10)),
    axis.text.y = element_text(size = 12, margin = margin(r = 5)),
    legend.text = element_text(size = 10),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10)))+
  ylab("Temperature (˚C)")+
  xlab("Experimental week")+
  labs(color = "Models", size = "Weight")

Figure 7. Optimal temeprature at every week of development according to four models and an optimized inverse logistic model.

This identifies a potential temperature profile to maximize proportional growth throughout the cricket development. This will need to be trialed in the future to see how this temperature profile affect cricket growth and development. Preliminary trials are promising.

Diets

Test for Significance

With our goal to identify conditions maximizing growth throughout ontogeny, we will test if significant differences can be identified in growth between the different diets to support the visual differences we see on the figures.

Code

# Let's continue with identifying diet that maximize growth throughout ontogeny.
# We start by isolating the data for the different diets.

sig_diet <- sig |>
  filter(temperature == "32") |> # All the crickets raised at 32˚C were in the diet trial.
  mutate(week = as.factor(week))

sig_diet_model <- lm(propgrowth ~ diet + week + diet:week, sig_diet)
anova(sig_diet_model)

Analysis of Variance Table

Response: propgrowth
            Df   Sum Sq Mean Sq  F value    Pr(>F)    
diet         1    58746   58746   6.8672  0.008865 ** 
week         6 46186051 7697675 899.8321 < 2.2e-16 ***
diet:week    6   308335   51389   6.0072 3.187e-06 ***
Residuals 1546 13225363    8555                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diet, week and their interaction are identified as significant in the anova (p < 0.05). We then need to look at the different week separately. To do that, we will perform pairwise comparisons for each week to see which week has post-hoc differences.

Code

# Perform pairwise comparisons for each week
pairwise_diet_week <- lapply(levels(sig_diet$week), function(week_level) {
  subset_data <- subset(sig_diet, week == week_level)
  pairwise.test <- pairwise.t.test(subset_data$propgrowth, subset_data$diet, p.adjust.method = "bonferroni")
  return(pairwise.test)
})

# Print results
names(pairwise_diet_week) <- levels(sig_diet$week)
print(pairwise_diet_week$`1`)


    Pairwise comparisons using t tests with pooled SD 

data:  subset_data$propgrowth and subset_data$diet 

  1       2       3       4       5       6      
2 1.00000 -       -       -       -       -      
3 0.14748 1.00000 -       -       -       -      
4 0.00340 1.00000 1.00000 -       -       -      
5 0.00087 0.54104 1.00000 1.00000 -       -      
6 0.00054 0.21633 1.00000 1.00000 1.00000 -      
7 0.03582 1.00000 1.00000 1.00000 1.00000 1.00000

P value adjustment method: bonferroni

Code

We see only differences in proportional growth in the first week of development.

Week 1 Optimum

To modelize optimal diets, the quadratic model is most often used. We will then modelize the relationship between proportional growth and diet at the first week to extract an optimal diet. We will then plot it in Figure 9 below. The black points represent the individual proportional growth for the different diets. The red dots is the mean proportional growth for that diet, the black line is the quadratic model computed above and the vertical red line is the optimal diet for the first week of the trial.

Code

# We will now identify which diet maximized growth on the first week of the trial. 

# Isolation of the first week.
  modeldw1 <- sig %>%
    subset(temperature == "32") %>%
    filter(week == "1")
  
# Fit the model using nls_multstart to a quadratic model. The quadratic model is generally used in this setting.
fit_modeldw1 <- nls_multstart(propgrowth ~ quadratic_2008(temp = dietperc, a, b, c),
                                 data = modeldw1,
                                 iter = c(5, 5, 5),
                                 start_lower = get_start_vals(modeldw1$dietperc, modeldw1$propgrowth, model_name = 'quadratic_2008') - 10,
                                 start_upper = get_start_vals(modeldw1$dietperc, modeldw1$propgrowth, model_name = 'quadratic_2008') + 10,
                                 lower = get_lower_lims(modeldw1$dietperc, modeldw1$propgrowth, model_name = 'quadratic_2008'),
                                 upper = get_upper_lims(modeldw1$dietperc, modeldw1$propgrowth, model_name = 'quadratic_2008'),
                                 supp_errors = 'Y',
                                 convergence_count = FALSE)
  

# The percentage of protein maximizing growth in the first week
  modeldw1topt <- get_topt(fit_modeldw1)

  # 27% of protein maximizes growth in the first week. 
  
# New data frame for the model created above
new_data <- data.frame(dietperc = seq(5, 70, 0.1))
augmented_dataw1 <- augment(fit_modeldw1, newdata = new_data)

meandietperc1 <- modeldw1 |>
  group_by(dietperc) |>
  summarise(propgrowth = mean(propgrowth))

ggplot(data = modeldw1, aes(x = dietperc, y = propgrowth)) +
  geom_point(position = position_jitter(width = 0.7)) +
  geom_point(data = meandietperc1, aes(x = dietperc, y = propgrowth), color = "red", size = 3) +
  geom_line(data = augmented_dataw1, aes(x = dietperc, y = .fitted)) +
  geom_vline(aes(xintercept = modeldw1topt), color = "red", size = 1.2) +
  annotate("text", x = 29, y = 820, 
           label = "Optimal diet: 27%", color = "black", size = 5, angle = 270, hjust = 0) +
  theme_classic() +
  theme(
    axis.title = element_text(size = 18),
    legend.box.margin = margin(-15),
    strip.placement = "left",
    axis.text.x = element_text(size = 12, margin = margin(t = 10)),
    axis.text.y = element_text(size = 12, margin = margin(r = 5)),
    legend.text = element_text(size = 10),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10))
  ) +
  ylab("Proportional growth (%)") +
  xlab("Percentage of protein (%)")

Figure 8. Proportional growth at week 1 at 7 different diets.

Since there are no additional differences in growth between diets for the rest of the weeks, we will identify the diet maximizing mass at the end of the experiment.

General Optimal Diet

To identify a diet maximizing mass at the end of the experiment, we will use a quadratic model to extract the apex. We will then plot it in figure 10. The black points represent the individual mass for the different diets. The red dots is the mean mass for that diet, the black line is the quadratic model computed above and the vertical red line is the optimal diet for the trial overall.

Code

# Isolation of the first week.
modeldw6 <- sig %>%
  subset(temperature == "32") %>%
  filter(week == "6")

# Fit the model using nls_multstart to a quadratic model. The quadratic model is generally used in this setting.
fit_modeldw6 <- nls_multstart(mass ~ quadratic_2008(temp = dietperc, a, b, c),
                              data = modeldw6,
                              iter = c(5, 5, 5),
                              start_lower = get_start_vals(modeldw6$dietperc, modeldw6$mass, model_name = 'quadratic_2008') - 10,
                              start_upper = get_start_vals(modeldw6$dietperc, modeldw6$mass, model_name = 'quadratic_2008') + 10,
                              lower = get_lower_lims(modeldw6$dietperc, modeldw6$mass, model_name = 'quadratic_2008'),
                              upper = get_upper_lims(modeldw6$dietperc, modeldw6$mass, model_name = 'quadratic_2008'),
                              supp_errors = 'Y',
                              convergence_count = FALSE)


# The percentage of protein maximizing mass at the end of the trial
modeldw6topt <- get_topt(fit_modeldw6)

# 38.062% of protein maximizes final mass. 

# New data frame for the model created above
new_data <- data.frame(dietperc = seq(5, 70, 0.1))
augmented_dataw6 <- augment(fit_modeldw6, newdata = new_data)

meandietperc6 <- modeldw6 |>
  group_by(dietperc) |>
  summarise(mass = mean(mass))

ggplot(data = modeldw6, aes(x = dietperc, y = mass))+
  geom_point(position = position_jitter(w=0.7))+
  geom_point(data=meandietperc6, aes(x=dietperc, y=mass), color = "red", size = 3) +
  geom_line(data = augmented_dataw6, aes(dietperc, .fitted))+
  geom_vline(aes(xintercept=modeldw6topt),color = "red", size = 1.2)+
  annotate("text", x = 41, y = 560, 
           label = "Optimal diet: 38%", color = "black", size = 5, angle = 270, hjust = 0) +
  theme_classic() +
  theme(
    axis.title = element_text(size = 18),
    legend.box.margin = margin(-15),
    strip.placement = "left",
    axis.text.x = element_text(size = 12, margin = margin(t = 10)),
    axis.text.y = element_text(size = 12, margin = margin(r = 5)),
    legend.text = element_text(size = 10),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.title.y = element_text(margin = margin(r = 10)))+
  ylab("Mass (mg)")+
  xlab("Percentage of protein (%)")

Figure 9. Final mass of 7 dietary treatments.

The diet maximizing the mass at 6 weeks of development contains 38% proteins according to the apex of the quadratic model shown above. This is close to the 1:1 protein to carbohydrate feed.

Conclusion

Optimizing diets is not providing such a clear way forward. There are indication that feeding crickets a diet of 27% proteins for the first week and 38% proteins for the rest of the growing period may maximize growth and size at adulthood, but this would need to be trialed in the future.

We are clearly able to identify a temperature profile to maximize proportional growth throughout ontogeny. An inverse logistic regression adequately modelizes the optimal temperature. This will be trialed in the coming months.

References

Anderson, R., Bayer, P.E., Edwards, D., 2020. Climate change and the need for agricultural adaptation. Curr. Opin. Plant Biol., Biotic interactions ● AGRI 2019 56, 197–202. https://doi.org/10.1016/j.pbi.2019.12.006

Halloran, A., Flore, R., Vantomme, P., Roos, N. (Eds.), 2018. Edible insects in sustainable food systems. Springer International Publishing, Cham. https://doi.org/10.1007/978-3-319-74011-9

May, J., Lott, B., Simmons, J., 1998. The effect of environmental temperature and body weight on growth rate and feed:gain of male broilers. Poult. Sci. 77, 499–501. https://doi.org/10.1093/ps/77.4.499

NRC, 1994. Nutrient Requirements of Poultry: Ninth Revised Edition, 1994. National Academies Press.

Padfield, D., O’Sullivan, H., Pawar, S., 2021. rTPC and nls.multstart: A new pipeline to fit thermal performance curves in r. Methods in Ecology and Evolution 12, 1138–1143. https://doi.org/10.1111/2041-210X.13585