The structural stability of nested experimental plankton system

packages:

library("readxl")
library("dplyr")
library("ggplot2")
library("patchwork")
library("RColorBrewer")
library("deSolve")

set directions

# raw data direction
if (!dir.exists("../data")) {
  dir.create("../data", recursive = TRUE)
}

# output direction
if (!dir.exists("../output")) {
  dir.create("../output", recursive = TRUE)
}

# plot direction
if (!dir.exists("../plot")) {
  dir.create("../plot", recursive = TRUE)
}

data_dir <- "../data"
output_dir <- "../output"
plot_dir <- "../plot"

Read data and preliminary operation

read data sheet

file_path <- file.path(data_dir, "WP2_2_clean.xlsx")
data<-read_xlsx(file_path, sheet = "Sheet1")

# data <- data %>% select(-date, -grazer_treatment,-diversity_level,-comment)

colnames(data)

 [1] "identifier"               "system"                  
 [3] "experiment"               "date"                    
 [5] "day"                      "grazer_treatment"        
 [7] "nutrient_treatment"       "species_richness"        
 [9] "diversity_level"          "ID_bucket"               
[11] "number_bucket"            "number_bag"              
[13] "name"                     "group"                   
[15] "species"                  "Fix_Area_bag"            
[17] "Chla_bag"                 "Fix_Area_recalc"         
[19] "Fix_Area_accrual"         "Chla_recalc"             
[21] "Fv/Fm"                    "ABS/RC"                  
[23] "Eto/RC"                   "GR"                      
[25] "GR_total"                 "Fix_Area_recalc_expected"
[27] "RYi_Fix_Area"             "RYi_Fix_Area_accrual"    
[29] "PP"                       "filteredVol_PP_ml"       
[31] "PP_µg/L"                  "PP_µmol/L"               
[33] "%N"                       "%C"                      
[35] "filteredVol_CN_ml"        "N_mg/L"                  
[37] "C_mg/L"                   "N_µmol/L"                
[39] "C_µmol/L"                 "C:N"                     
[41] "C:P"                      "N:P"                     
[43] "C_accrual"                "RYi_C"                   
[45] "C:Fix_Area"               "Fix_Area:N"              
[47] "Fix_Area:P"               "cells/L_species"         
[49] "vol/cell_µm3"             "source"                  
[51] "biovolume_µm3/L_species"  "biovolume_mm3/L_species" 
[53] "biovol_perc_species"      "cells/L_total"           
[55] "biovolume_mm3/L_total"    "comment"                 

take necessary columns:

time, treatment class, species, labels and group, population size (Fix_Area_Recalc)

And in each treatment, bucket, we sum the same species up.

df <- data %>% select(day, nutrient_treatment, species_richness, number_bucket, number_bag, group, species, Fix_Area_recalc)
df <- df[df$species!="Medium",]
df$number_bag <- as.factor(df$number_bag)
df$species <- as.factor(df$species)
df$Fix_Area_recalc <- as.numeric(as.character(df$Fix_Area_recalc))

df_total <- df %>%
  group_by(day, nutrient_treatment, species_richness, number_bucket, group, species) %>%
  summarise(Fix_Area_recalc = sum(Fix_Area_recalc, na.rm = TRUE), .groups = "drop")

colnames(df)

[1] "day"                "nutrient_treatment" "species_richness"  
[4] "number_bucket"      "number_bag"         "group"             
[7] "species"            "Fix_Area_recalc"   

Basic plot

We then set a project from species to color:

display.brewer.all(type = "qual")

species <- unique(df$species)

species_colors <- brewer.pal(n = length(species), name = "Set2")
names(species_colors) <- species

par(mar = c(2, 15, 2, 10))

barplot(
  height = rep(1, length(species_colors)),
  names.arg = names(species_colors),
  col = species_colors,
  border = NA,
  horiz = TRUE,    
  las = 1          
)

Now, the species is related to the color shown above.

Then we divide the data by treatment: N pulse, N limited

df_Npulse <- df[df$nutrient_treatment == "N pulse",]
df_Nlimited <- df[df$nutrient_treatment == "N limited",]

Then split the data by species richness as well as bucket.

split_bucket <- function(df) {
  return(lapply(split(df, df$species_richness), function(sub_df){
    split(sub_df, sub_df$number_bucket)
  }))
}

draw_splited_bucket <- function(df){
  
  data_list <- split_bucket(df)
  
  df_total <- df %>%
    group_by(day, nutrient_treatment, species_richness, number_bucket, group, species) %>%
    summarise(Fix_Area_recalc = mean(Fix_Area_recalc, na.rm = TRUE), .groups = "drop")
  
  sum_list <- split_bucket(df_total)
  
  plot_list <- list()
  sum_plot_list <- list()
  
  dummy <- data.frame(
    species = factor(species, levels = species),
    day = NA_real_,
    Fix_Area_recalc = NA_real_,
    number_bag = NA_real_
    )
  
  plot_geoms <- list(
    geom_line(linewidth = 1.2, show.legend = FALSE),
    geom_point(size = 2, show.legend = FALSE)
  )

  dummy_legend <- list(
    geom_line(
        data = dummy,
        aes(x = day, y = Fix_Area_recalc, color = species),
        linewidth = 1.2,
        alpha = 0, 
        show.legend = TRUE,
        na.rm = TRUE
      )
  )

  plot_theme <- list(
    theme_minimal(),
    theme(
      axis.title = element_blank(),
      axis.ticks = element_blank()
    )
  )

  species_scale <- list(
    scale_color_manual(
      name = "Species",
      values = species_colors,
      drop = FALSE,
      guide = guide_legend(
        override.aes = list(alpha = 1, size = 1, 
                            linetype = 1, shape = 16)
      )
    )
  )
  
  for (richness in names(data_list)) {
    p_richness_list <- list()
    s_richness_list <- list()
    for (bucket in names(data_list[[richness]])) {
    
      p1 <- ggplot(data_list[[richness]][[bucket]],
                 aes(x = day, y = Fix_Area_recalc, 
                     group = number_bag, color = species)) +
        plot_geoms+
        plot_theme+
        dummy_legend+
        species_scale+
        labs(title = bucket)
      
      p_richness_list[[bucket]]<-p1
      
      p_sum<- ggplot(sum_list[[richness]][[bucket]], 
               aes(x = day, 
                   y = Fix_Area_recalc, 
                   color = species)) +
        plot_geoms+
        plot_theme+
        dummy_legend+
        species_scale+
        labs(title = bucket)
        
      s_richness_list[[bucket]]<-p_sum
    }
    plot_list[[richness]]<- p_richness_list
    sum_plot_list[[richness]] <- s_richness_list
  }
  return(list(
    data_list,
    sum_list,
    plot_list,
    sum_plot_list
  ))
}

draw_Npulse <- draw_splited_bucket(df_Npulse)
Npulse_list <- draw_Npulse[[1]]
Npulse_sum <- draw_Npulse[[2]]
Npulse_plots <- draw_Npulse[[3]]
Npulse_sumplots <- draw_Npulse[[4]]

draw_Nlimited <- draw_splited_bucket(df_Nlimited)
Nlimited_list <- draw_Nlimited[[1]]
Nlimited_sum <- draw_Nlimited[[2]]
Nlimited_plots <- draw_Nlimited[[3]]
Nlimited_sumplots <- draw_Nlimited[[4]]

Now, each experiment is organized in a table. _list is the list of all experiment, _sum is the sum of each species, plots is the population size change plot, as well as sumplot is the plot of sum species population (or average). For example:

(Npulse_plots[["1"]][["D1_1"]]|Npulse_plots[["8"]][["D8_1"]])+
  plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

We first show the single species case.

draw_all_species <- function(draw_list){
  (draw_list[["1"]][["D1_1"]] | draw_list[["1"]][["D1_2"]] | draw_list[["1"]][["D1_3"]] | draw_list[["1"]][["D1_4"]])/
  (draw_list[["1"]][["D1_5"]] | draw_list[["1"]][["D1_6"]] | draw_list[["1"]][["D1_7"]] | draw_list[["1"]][["D1_8"]]) +
  plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")
}
draw_all_species(Npulse_plots)& 
  plot_annotation(
    title = "Npulse single species"
  )

The 8 curves in each graph is the data of 8 bags. We could note that each species shows the similar curve. However, species 1, 4, 6 first decrease then increase, species 3 shows two convex curve, as well as species 7 first increase to a peak then decrease. Only 2, 5, 8 might follows a logistic increase.

Compare to the Nlimited case:

draw_all_species(Nlimited_plots)& 
  plot_annotation(
    title = "Nlimited single species"
  )

Shows the similar trend, but with more noise in D1_4, D1_7, D1_8.

Then we plot the population change in co-cultivate cases:

draw_multi_species <- function(draw_list){
  (draw_list[["2"]][["D2_1"]] | draw_list[["2"]][["D2_2"]] | draw_list[["2"]][["D2_3"]] | draw_list[["2"]][["D2_4"]]) /
    (draw_list[["4"]][["D4_1"]] | draw_list[["4"]][["D4_2"]] | draw_list[["4"]][["D4_3"]] | draw_list[["4"]][["D4_4"]]) /
    (draw_list[["8"]][["D8_1"]] | draw_list[["8"]][["D8_2"]] | draw_list[["8"]][["D8_3"]] | draw_list[["8"]][["D8_4"]])+
  plot_layout(guides = 'collect') &
  theme(legend.position = "right")
}
draw_multi_species(Npulse_plots)& 
  plot_annotation(
    title = "Npulse multi-species"
  )

Species Chlamydomonas (7) and Monoraphidium (8) shows a much higher population size compared to others, and their changing trend seems keep the same. We may need to rescale this two species later.

draw_multi_species(Nlimited_plots)& 
  plot_annotation(
    title = "Nlimited multi-species"
  )

Nlimited case shows more complex results.

We may also worried that species do not interact to each other. If this is true, then the population curve would keep similar in all cases. We plot each curve for single species:

ggplot(df[df$species==species[3] & df$nutrient_treatment=="N pulse",],aes(x = day, y = Fix_Area_recalc, group = paste(number_bucket, number_bag), colour = number_bucket))+
  geom_line()+
  geom_point()+
  theme_minimal()

Species shows similar change in each combinations, indicate the interactions might very weak.

To make the graph more clear, we show the sum plot below, for Npulse treatment:

draw_all_species(Npulse_sumplots)& 
  plot_annotation(
    title = "Npulse single species sum"
  )

draw_multi_species(Npulse_sumplots)& 
  plot_annotation(
    title = "Npulse multi-species sum"
  )

draw_all_species(Nlimited_sumplots)& 
  plot_annotation(
    title = "Nlimited single species sum"
  )

draw_multi_species(Nlimited_sumplots)& 
  plot_annotation(
    title = "Npulse multi-species sum"
  )

ggplot(df_total[df_total$species==species[6] & df_total$nutrient_treatment=="N pulse",],aes(x = day, y = Fix_Area_recalc*species_richness, color = number_bucket))+
  geom_line()+
  geom_point()

Model selection and parameter inference

We write the Logistic model:

\[ \frac{\mathrm{d}N(t)}{\mathrm{d}t} = r N ( 1 - \frac{N}{K} ) \]

Here are only two parameters to be determined, intrinsic growth rate \(r\) and environment capacity \(K\). We could infer these parameters from the 5 points data. However, logistic model could only show very limited dynamics. As what we show below, population size either grow to the environment capacity, or decrease to it.

Ks <- 10
t <- seq(0, 10, by = 0.05)

N0_1 <- 0.1
rs_1 <- c(1, 2, 5, 10)

N0_2 <- 20
rs_2 <- c(0.1, 0.2, 0.5, 1)

logistic <- function(K, N0, r, t) {
  K / (1 + ((K - N0)/N0) * exp(-r * t))
}

df1 <- do.call(rbind, lapply(rs_1, function(r) {
  data.frame(
    t = t,
    N = logistic(Ks, N0_1, r, t),
    r = factor(paste0("r=", r)),
    group = "N0=0.1"
  )
}))

df2 <- do.call(rbind, lapply(rs_2, function(r) {
  data.frame(
    t = t,
    N = logistic(Ks, N0_2, r, t),
    r = factor(paste0("r=", r)),
    group = "N0=15"
  )
}))

df_logi <- rbind(df1, df2)

p <- ggplot(df_logi, aes(x = t, y = N, color = r, linetype = group)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = Ks, linetype = "dashed") +
  geom_hline(yintercept = Ks/2, linetype = "dotted", color = "grey60") +
  labs(
    title = "Logistic Growth Curves",
    x = "t", 
    y = "N(t)",
    color = "Growth rate r",
    linetype = "Initial N0"
  ) +
  theme_minimal() +
  theme(
    legend.position = "right"
  )

print(p)

The increase curve is the general logistic curve where increase rate first increase then decrease. the decrease curve could only decrease with a decreasing rate. As a result, only species 2, 5, 8 could somehow fit to the logistic model. For the others, we may need

1. model the initial decline. which could be modeled by the change of intrinsic rate, which represent the adaption to a new environment

\[ \frac{\mathrm{d}N(t)}{\mathrm{d}t} = r(t) N ( 1 - \frac{N}{K} )\\ r(t) = r_{\infty} - (r_{\infty} - r_{0})e^{-\alpha t} \]

Here \(r(0) = r_{0}\) and \(r(t)_{t \rightarrow \infty} \rightarrow r_{ \infty }\) . But this model introduced 3 more parameters, makes total parameters 5, which is the same number of the observed data. This case we may not be able to give a efficient prediction.

2. model the increase-decrease trend. This could be modeled by lead in an inhibitor

\[ \frac{\mathrm{d}N(t)}{\mathrm{d}t} = r N ( 1 - \frac{N}{K} ) - \gamma N W\\ \frac{\mathrm{d}W(t)}{\mathrm{d}t} = \mu N - \kappa W \]

Here \(W\) could be some toxic substances. But again, here need 5 parameters.

As a result, we first try to infer the parameters for logistic model of species 2,5,8.

We could solve the logistic model as

\[ N(t) = \frac{K}{1 + A e^{-rt}} \]

where \(A\) depends on the initial state \(A = \frac{K}{N_0} - 1\). If we view \(N_0\) as a known parameter, which is the population size at day \(=0\). But if we think the initial state also an observed data, we actually need to infer 3 parameters here.

we plot these three species again

(Npulse_sumplots[["1"]][["D1_2"]]|Npulse_sumplots[["1"]][["D1_5"]]|Npulse_sumplots[["1"]][["D1_8"]])+
plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

nonlinear least squares

We first use the nonlinear least squares method to estimate the parameters.

3 parameters:

df_N_t_2 <- Npulse_sum[["1"]][["D1_2"]][c("day","Fix_Area_recalc")]
fit2 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = df_N_t_2,
  start = list(
    K = 42000,
    r = 0.5,
    A = -0.45
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

df_N_t_5 <- Npulse_sum[["1"]][["D1_5"]][c("day","Fix_Area_recalc")]
fit5 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = df_N_t_5,
  start = list(
    K = 310000,
    r = 0.36,
    A = 5.2
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

df_N_t_8 <- Npulse_sum[["1"]][["D1_8"]][c("day","Fix_Area_recalc")]
fit8 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = df_N_t_8,
  start = list(
    K = 1400000,
    r = 0.9,
    A = 65
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
coeff_pulse <- rbind(coef(fit2),coef(fit5),coef(fit8))
coeff_pulse

              K         r          A
[1,]   5215.622 0.5079697 -0.4557603
[2,]  39175.577 0.3614746  5.2230213
[3,] 177963.965 0.9050188 67.8913708

summary(fit2)

Formula: Fix_Area_recalc ~ K/(1 + A * exp(-r * day))

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
K 5215.62184  170.11494  30.659  0.00106 **
r    0.50797    0.11304   4.494  0.04612 * 
A   -0.45576    0.02048 -22.250  0.00201 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 191.7 on 2 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 1.674e-06

summary(fit5)

Formula: Fix_Area_recalc ~ K/(1 + A * exp(-r * day))

Parameters:
   Estimate Std. Error t value Pr(>|t|)  
K 3.918e+04  6.798e+03   5.763   0.0288 *
r 3.615e-01  7.626e-02   4.740   0.0417 *
A 5.223e+00  9.090e-01   5.746   0.0290 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1245 on 2 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 5.103e-06

summary(fit8)

Formula: Fix_Area_recalc ~ K/(1 + A * exp(-r * day))

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
K 1.780e+05  7.550e+03  23.570   0.0018 **
r 9.050e-01  1.032e-01   8.769   0.0128 * 
A 6.789e+01  2.762e+01   2.458   0.1332   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4783 on 2 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 2.352e-06

fit2_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/df_N_t_2$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = df_N_t_2,
  start = list(
    K = 41000,
    r = 0.5
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

fit5_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/df_N_t_5$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = df_N_t_5,
  start = list(
    K = 350000,
    r = 0.3
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

fit8_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/df_N_t_8$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = df_N_t_8,
  start = list(
    K = 1400000,
    r = 0.9
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
coeff_pulse_l <- rbind(coef(fit2_l),coef(fit5_l),coef(fit8_l))
coeff_pulse_l

             K           r
[1,] 223951.10 -0.07770653
[2,]  44746.51  0.29877416
[3,] 192622.86  0.69485278

coeff_pulse

              K         r          A
[1,]   5215.622 0.5079697 -0.4557603
[2,]  39175.577 0.3614746  5.2230213
[3,] 177963.965 0.9050188 67.8913708

coeff_pulse_l

             K           r
[1,] 223951.10 -0.07770653
[2,]  44746.51  0.29877416
[3,] 192622.86  0.69485278

summary(fit2_l)

Formula: Fix_Area_recalc ~ K/(1 + (K/df_N_t_2$Fix_Area_recalc[1] - 1) * 
    exp(-r * day))

Parameters:
    Estimate Std. Error t value Pr(>|t|)
K  2.240e+05  2.933e+07   0.008    0.994
r -7.771e-02  3.627e-01  -0.214    0.844

Residual standard error: 1560 on 3 degrees of freedom

Number of iterations till stop: 32 
Achieved convergence tolerance: 7.503
Reason stopped: 算法的步因素0.000488281的大小被减少到小于0.000976563的'minFactor'值

summary(fit5_l)

Formula: Fix_Area_recalc ~ K/(1 + (K/df_N_t_5$Fix_Area_recalc[1] - 1) * 
    exp(-r * day))

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
K 4.475e+04  9.524e+03   4.698  0.01824 * 
r 2.988e-01  4.239e-02   7.048  0.00587 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1249 on 3 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 5.27e-06

summary(fit8_l)

Formula: Fix_Area_recalc ~ K/(1 + (K/df_N_t_8$Fix_Area_recalc[1] - 1) * 
    exp(-r * day))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K 1.926e+05  1.364e+04   14.12 0.000770 ***
r 6.948e-01  3.721e-02   18.68 0.000335 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7732 on 3 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 5.035e-06

plot(df_N_t_5$day, df_N_t_5$Fix_Area_recalc)
curve( coef(fit5)[["K"]] / (1 + coef(fit5)[["A"]] * exp(- coef(fit5)[["r"]] * x )),
      add = TRUE, lwd = 2, col = "blue")

curve( coef(fit5_l)[["K"]] / (1 + (coef(fit5_l)[["K"]]/df_N_t_5$Fix_Area_recalc[1] - 1) * exp(- coef(fit5_l)[["r"]] * x )),
      add = TRUE, lwd = 2, col = "green")

(Nlimited_sumplots[["1"]][["D1_2"]]|Nlimited_sumplots[["1"]][["D1_5"]]|Nlimited_sumplots[["1"]][["D1_8"]])+
plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

d_N_t_2 <- Nlimited_sum[["1"]][["D1_2"]][c("day","Fix_Area_recalc")]
fit_limit2 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = d_N_t_2,
  start = list(
    K = 38000,
    r = 1.6,
    A = -0.22
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
d_N_t_5 <- Nlimited_sum[["1"]][["D1_5"]][c("day","Fix_Area_recalc")]
fit_limit5 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = d_N_t_5,
  start = list(
    K = 58000,
    r = 0.3,
    A = 0.6
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
d_N_t_8 <- Nlimited_sum[["1"]][["D1_8"]][c("day","Fix_Area_recalc")]
fit_limit8 <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = d_N_t_8,
  start = list(
    K = 220000,
    r = 0.45,
    A = 3.5
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
coeff_limit <- rbind(coef(fit_limit2),coef(fit_limit5),coef(fit_limit8))
coeff_limit

             K         r          A
[1,]  4786.891 1.6346296 -0.2162845
[2,]  7302.432 0.3120379  0.6700429
[3,] 28109.213 0.4463273  3.5334045

fit_limit2_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/d_N_t_2$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = d_N_t_2,
  start = list(
    K = 38000,
    r = 1.6
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

fit_limit5_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/d_N_t_5$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = d_N_t_5,
  start = list(
    K = 58000,
    r = 0.3
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

fit_limit8_l <- nls(
  Fix_Area_recalc ~ K / (1 + (K/d_N_t_8$Fix_Area_recalc[1] - 1) * exp(-r * day )),
  data = d_N_t_8,
  start = list(
    K = 220000,
    r = 0.45
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)
coeff_limit_l <- rbind(coef(fit_limit2_l),coef(fit_limit5_l),coef(fit_limit8_l))
coeff_limit_l

             K         r
[1,]  4786.891 1.6346229
[2,]  7355.162 0.2939016
[3,] 28454.925 0.4247108

plot(d_N_t_5$day, d_N_t_5$Fix_Area_recalc)
curve( coef(fit_limit5)[["K"]] / (1 + coef(fit_limit5)[["A"]] * exp(- coef(fit_limit5)[["r"]] * x )),
      add = TRUE, lwd = 2, col = "blue")

curve( coef(fit_limit5_l)[["K"]] / (1 + (coef(fit_limit5_l)[["K"]]/d_N_t_5$Fix_Area_recalc[1] - 1) * exp(- coef(fit_limit5_l)[["r"]] * x )),
      add = TRUE, lwd = 2, col = "green")

((Npulse_sumplots[["1"]][["D1_2"]] + ylim(0,40000))|(Npulse_sumplots[["2"]][["D2_3"]] + ylim(0,40000))|(Npulse_sumplots[["1"]][["D1_5"]]) + ylim(0,40000))+
plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

((Npulse_plots[["1"]][["D1_2"]] + ylim(0,40000))|(Npulse_plots[["2"]][["D2_3"]] + ylim(0,40000))|(Npulse_plots[["1"]][["D1_5"]]) + ylim(0,40000))+
plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

((Nlimited_sumplots[["1"]][["D1_2"]] + ylim(0,10000))|(Nlimited_sumplots[["2"]][["D2_3"]] + ylim(0,10000))|(Nlimited_sumplots[["1"]][["D1_5"]]) + ylim(0,10000))+
plot_layout(guides = 'collect') &
  theme(legend.position = "bottom")

df_2_total <- Npulse_list[["1"]][["D1_2"]][c("day","Fix_Area_recalc")]
fit2_all <- nls(
  Fix_Area_recalc ~ K / (1 + A * exp(-r * day )),
  data = df_2_total,
  start = list(
    K = 42000,
    r = 0.5,
    A = -0.45
    ), 
  control = nls.control(maxiter = 100, warnOnly = TRUE)
)

summary(fit2)

Formula: Fix_Area_recalc ~ K/(1 + A * exp(-r * day))

Parameters:
    Estimate Std. Error t value Pr(>|t|)   
K 5215.62184  170.11494  30.659  0.00106 **
r    0.50797    0.11304   4.494  0.04612 * 
A   -0.45576    0.02048 -22.250  0.00201 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 191.7 on 2 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 1.674e-06

summary(fit2_all)

Formula: Fix_Area_recalc ~ K/(1 + A * exp(-r * day))

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K 5215.62650  338.26445   15.42  < 2e-16 ***
r    0.50797    0.22477    2.26   0.0298 *  
A   -0.45576    0.04073  -11.19 1.96e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1078 on 37 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 2.85e-06

Interaction parameters

We have the data of co-cultivate of species 2 and 5. we consider the generalized Lotka-Volterra model:

\[ \frac{\mathrm{d}N_{1}(t)}{\mathrm{d}t} = N_{1} ( r_{1} - \alpha_{1,1}N_{1} - \alpha_{1,2}N_{2} )\\ \frac{\mathrm{d}N_{2}(t)}{\mathrm{d}t} = N_{2} ( r_{2} - \alpha_{2,2}N_{2} - \alpha_{2,1}N_{1} ) \]

We already get the \(r_1\) and \(r_2\) from the former inference of single species. For paramter \(\alpha_{i,i}\) it is actually \(r_i/K_i\). As a result, we only need to infer the parameter \(\alpha_{1,2}\) and \(\alpha_{2,1}\), which are the inter-species interaction strength.

We could also use the nonlinear least squares method. But this time, we should first use linear regression on the difference approximation of the data, to estimate a rough region of parameters. This step is done by hand in the former part, since parameters in logistic curve is quite obvious.

To simplify the model, as well as make the computation more efficient, we rescale the logistic model:

\[ \frac{\mathrm{d}N}{\mathrm{d}t} = r N ( 1 - \frac{N}{K} ) \]

Because of the significant different of \(K\) for each species, directly using these values in gLV will makes the interaction a very small value, which is hard to evaluate. As a result, we rewrite it as:

\[ \frac{\mathrm{d}N_{K}}{\mathrm{d}t} = r N_{K} ( 1 - N_{K} )\\ \frac{\mathrm{d}N_{rK}}{\mathrm{d}t} = N_{rK} ( r- N_{rK} ) \]

where \(N_{K} = N / K\),or , \(N_{rK} = r N /K\), while \(r\) kept the same. now all species are comparable. We use the rescaled population size to write the gLV:

\[ \frac{\mathrm{d}N_{1}}{\mathrm{d}t} = N_{1} ( r_{1} - N_{1} - \alpha_{1,2}N_{2} )\\ \frac{\mathrm{d}N_{2}}{\mathrm{d}t} = N_{2} ( r_{2} - N_{2} - \alpha_{2,1}N_{1} ) \]

###
#observed data:
### 
df_inter <- Npulse_sum[["2"]][["D2_3"]]
#t  <- c(0, 2, 4, 6, 8) 
t <- df_inter[df_inter$species == species[2],]$day
N1 <- df_inter[df_inter$species == species[2],]$Fix_Area_recalc
N2 <- df_inter[df_inter$species == species[5],]$Fix_Area_recalc

r1 <- coef(fit2)[["r"]]
r2 <- coef(fit5)[["r"]]
K1 <- coef(fit2)[["K"]]
K2 <- coef(fit5)[["K"]]
N1 <- r1 * N1 / K1
N2 <- r2 * N2 / K2


###############################################
### difference linear regression
###############################################
dN1 <- diff(N1) / diff(t)
#dN1 <- c(dN1, dN1[length(dN1)])
M1 <- N1[-length(N1)] + diff(N1)/2

dN2 <- diff(N2) / diff(t)
#dN2 <- c(dN2, dN2[length(dN2)])
M2 <- N2[-length(N2)] + diff(N2)/2

x  <- M1 * M2
y1 <- r1*M1 - M1^2 - dN1
y2 <- r2*M2 - M2^2 - dN2 

# linear regression
fit12 <- lm(y1 ~ 0 + x)
fit21 <- lm(y2 ~ 0 + x)

a12_init <- as.numeric(coef(fit12)[1])
a21_init <- as.numeric(coef(fit21)[1])

summary(fit12)

Call:
lm(formula = y1 ~ 0 + x)

Residuals:
       1        2        3        4 
 0.09948  0.01199 -0.02928 -0.02034 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)
x -0.08902    0.48658  -0.183    0.867

Residual standard error: 0.0614 on 3 degrees of freedom
Multiple R-squared:  0.01103,   Adjusted R-squared:  -0.3186 
F-statistic: 0.03347 on 1 and 3 DF,  p-value: 0.8665

summary(fit21)

Call:
lm(formula = y2 ~ 0 + x)

Residuals:
        1         2         3         4 
 0.006637 -0.004077 -0.006187  0.003453 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)  
x  0.12502    0.04817   2.595   0.0807 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.006079 on 3 degrees of freedom
Multiple R-squared:  0.6918,    Adjusted R-squared:  0.5891 
F-statistic: 6.734 on 1 and 3 DF,  p-value: 0.08072

estimate_initial_interaction <- function(t,N1,N2,r1,r2){

  ###############################################
  ### difference linear regression
  ###############################################
  dN1 <- diff(N1) / diff(t)
  M1 <- N1[-length(N1)] + diff(N1)/2
  dN2 <- diff(N2) / diff(t)
  M2 <- N2[-length(N2)] + diff(N2)/2

  x  <- M1 * M2
  y1 <- r1*M1 - M1^2 - dN1
  y2 <- r2*M2 - M2^2 - dN2 

  # linear regression
  fit12 <- lm(y1 ~ 0 + x)
  fit21 <- lm(y2 ~ 0 + x)
  
  return(list(fit12,fit21))
}

gLV <- function(t, y, params) {
    n1 <- y[1]; n2 <- y[2]
    a12 <- params["a12"]
    a21 <- params["a21"]

    dn1 <- r1*n1 - n1*n1 - a12*n1*n2
    dn2 <- r2*n2 - n2*n2 - a21*n1*n2

    list(c(dn1, dn2))
}

estimate_nls_interaction <- function(t,N1,N2,r1,r2, a12_init, a21_init){

  ###############################################
  ### nls + ODE（deSolve）
  ###############################################

  simulate_gLV <- function(a12, a21) {
    params <- c(a12=a12, a21=a21)
    y0 <- c(N1[1], N2[1])
    out <- ode(y=y0, times=t, func=gLV, parms=params)
    return(out[,2:3])
  }

  residual_fn <- function(p) {
    a12 <- p[1]; a21 <- p[2]
    pred <- simulate_gLV(a12, a21)
    res <- c(pred[,1] - N1, pred[,2] - N2)
    return(res)
  }

  obj_fn <- function(p) sum(residual_fn(p)^2)

  p0 <- c(a12_init, a21_init)

  fit <- optim(p0, obj_fn)
  
  return(fit)
}

df_inter<-Npulse_sum[["2"]][["D2_3"]]
df_inter1 <- df_inter[df_inter$species == species[2],]
df_inter2 <- df_inter[df_inter$species == species[5],]

t <- df_inter1$day
N1 <- df_inter1$Fix_Area_recalc
N2 <- df_inter2$Fix_Area_recalc
r1 <- coef(fit2)[["r"]]
r2 <- coef(fit5)[["r"]]
K1 <- coef(fit2)[["K"]]
K2 <- coef(fit5)[["K"]]
N1 <- r1 * N1 / K1
N2 <- r2 * N2 / K2

pulse_fit_inter <- estimate_initial_interaction(t,N1,N2,r1,r2)

pulse_a12_init <- as.numeric(coef(pulse_fit_inter[[1]])[1])
pulse_a21_init <- as.numeric(coef(pulse_fit_inter[[2]])[1])

pulse_nls_inter <- estimate_nls_interaction(t,N1,N2,r1,r2, pulse_a12_init,pulse_a21_init)

summary(pulse_fit_inter[[1]])

Call:
lm(formula = y1 ~ 0 + x)

Residuals:
       1        2        3        4 
 0.09948  0.01199 -0.02928 -0.02034 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)
x -0.08902    0.48658  -0.183    0.867

Residual standard error: 0.0614 on 3 degrees of freedom
Multiple R-squared:  0.01103,   Adjusted R-squared:  -0.3186 
F-statistic: 0.03347 on 1 and 3 DF,  p-value: 0.8665

summary(pulse_fit_inter[[2]])

Call:
lm(formula = y2 ~ 0 + x)

Residuals:
        1         2         3         4 
 0.006637 -0.004077 -0.006187  0.003453 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)  
x  0.12502    0.04817   2.595   0.0807 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.006079 on 3 degrees of freedom
Multiple R-squared:  0.6918,    Adjusted R-squared:  0.5891 
F-statistic: 6.734 on 1 and 3 DF,  p-value: 0.08072

pulse_nls_inter$par

[1] -0.001995256  0.130470828

pulse_a12_hat <- pulse_nls_inter$par[1]
pulse_a21_hat <- pulse_nls_inter$par[2]

simulate_gLV <- function(a12, a21) {
    params <- c(a12=a12, a21=a21)
    y0 <- c(N1[1], N2[1])
    out <- ode(y=y0, times=t, func=gLV, parms=params)
    return(out[,2:3])
}

pred <- simulate_gLV(pulse_a12_hat, pulse_a21_hat)

df <- data.frame(
  time = rep(t, 4),
  abundance = c(N1, pred[,1], N2, pred[,2]),
  type = factor(rep(c("N1 obs", "N1 fit", "N2 obs", "N2 fit"), 
                    each = length(t))),
  species = rep(c("N1", "N1", "N2", "N2"), each = length(t))
)

# 为 obs/fit 区分点和线
df$mode <- ifelse(grepl("obs", df$type), "obs", "fit")

# 绘图
p_inter25_pulse <- ggplot(df, aes(x = time, y = abundance, color = species)) +
  # 实测点
  geom_point(data = subset(df, mode == "obs"), size = 2) +
  # 拟合线
  geom_line(data = subset(df, mode == "fit"), linewidth = 1) +
  labs(
    title = "gLV NLS Npulse",
    x = "time",
    y = "abundance",
    color = "Species"
  ) +
  scale_color_manual(values = c("N1" = "blue", "N2" = "red")) +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )

print(p_inter25_pulse)

df_inter<-Nlimited_sum[["2"]][["D2_3"]]
df_inter1 <- df_inter[df_inter$species == species[2],]
df_inter2 <- df_inter[df_inter$species == species[5],]
t <- df_inter1$day
N1 <- df_inter1$Fix_Area_recalc
N2 <- df_inter2$Fix_Area_recalc
r1 <- coef(fit_limit2)[["r"]]
r2 <- coef(fit_limit5)[["r"]]
K1 <- coef(fit_limit2)[["K"]]
K2 <- coef(fit_limit5)[["K"]]
N1 <- r1 * N1 / K1
N2 <- r2 * N2 / K2

limited_fit_inter <- estimate_initial_interaction(t,N1,N2,r1,r2)

limit_a12_init <- as.numeric(coef(limited_fit_inter[[1]])[1])
limit_a21_init <- as.numeric(coef(limited_fit_inter[[2]])[1])

limit_nls_inter <- estimate_nls_interaction(t,N1,N2,r1,r2, limit_a12_init,limit_a21_init)

summary(limited_fit_inter[[1]])

Call:
lm(formula = y1 ~ 0 + x)

Residuals:
      1       2       3       4 
-0.7864  0.3549  0.3253  0.2978 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)
x  -1.0974     0.4793  -2.289    0.106

Residual standard error: 0.5594 on 3 degrees of freedom
Multiple R-squared:  0.636, Adjusted R-squared:  0.5147 
F-statistic: 5.242 on 1 and 3 DF,  p-value: 0.106

summary(limited_fit_inter[[2]])

Call:
lm(formula = y2 ~ 0 + x)

Residuals:
        1         2         3         4 
-0.017142 -0.012457  0.001443  0.034461 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)
x 0.0006912  0.0200245   0.035    0.975

Residual standard error: 0.02337 on 3 degrees of freedom
Multiple R-squared:  0.000397,  Adjusted R-squared:  -0.3328 
F-statistic: 0.001192 on 1 and 3 DF,  p-value: 0.9746

limit_nls_inter$par

[1] -0.51443631 -0.00219877

limit_a12_hat <- limit_nls_inter$par[1]
limit_a21_hat <- limit_nls_inter$par[2]

simulate_gLV <- function(a12, a21) {
    params <- c(a12=a12, a21=a21)
    y0 <- c(N1[1], N2[1])
    out <- ode(y=y0, times=t, func=gLV, parms=params)
    return(out[,2:3])
}

pred <- simulate_gLV(limit_a12_hat, limit_a21_hat)

df <- data.frame(
  time = rep(t, 4),
  abundance = c(N1, pred[,1], N2, pred[,2]),
  type = factor(rep(c("N1 obs", "N1 fit", "N2 obs", "N2 fit"), 
                    each = length(t))),
  species = rep(c("N1", "N1", "N2", "N2"), each = length(t))
)

# 为 obs/fit 区分点和线
df$mode <- ifelse(grepl("obs", df$type), "obs", "fit")

# 绘图
p_inter25_limited <- ggplot(df, aes(x = time, y = abundance, color = species)) +
  # 实测点
  geom_point(data = subset(df, mode == "obs"), size = 2) +
  # 拟合线
  geom_line(data = subset(df, mode == "fit"), linewidth = 1) +
  labs(
    title = "gLV NLS Nlimited",
    x = "time",
    y = "abundance",
    color = "Species"
  ) +
  scale_color_manual(values = c("N1" = "blue", "N2" = "red")) +
  theme_minimal() +
  theme(
    legend.position = "bottom"
  )

print(p_inter25_limited)

we then change back:

now the functions are

\[ \frac{\mathrm{d}N_{1}}{\mathrm{d}t} = N_{1} ( r_{1} - N_{1} - \alpha_{1,2}N_{2} )\\ \frac{\mathrm{d}N_{2}}{\mathrm{d}t} = N_{2} ( r_{2} - N_{2} - \alpha_{2,1}N_{1} ) \]

where \(N1,N2\) are \(N_{rK} = r N_{K}\).

\[ \frac{\mathrm{d}N_{1}/r_1}{\mathrm{d}t} = N_{1}/r_1 ( r_{1} - r_1 N_{1}/r_1 - r_1 r_2 \alpha_{1,2} N_{2}/(r_1 r_2) )\\ \Rightarrow \frac{\mathrm{d}N_{1}}{\mathrm{d}t} = r_{1}N_{1} ( 1 - N_{1} - r_{2}/r_{1} \alpha_{1,2} N_{2} ) \\ \frac{\mathrm{d}N_{2}}{\mathrm{d}t} = r_{2} N_{2} ( 1 - N_{2} - r_{1}/r_{2} \alpha_{2,1}N_{1} ) \]

As a result, we write \(\alpha_{1,2} = \frac{r_2}{r_1}\alpha_{1,2}\), \(\alpha_{2,1} = \frac{r_1}{r_2}\alpha_{2,1}\),

At the positive stable equilibrium, we have:

\[ r_{1} - N^*_{1} - \alpha_{1,2} N^*_{2} = 0\\ r_{2} - N^*_{2} - \alpha_{2,1} N^*_{1} = 0 \]

As a result, we could solve:

\[ N^*_{1} = \frac{r_{1} - \alpha_{1,2}r_{2}}{1 - \alpha_{1,2} \alpha_{2,1}}\\ N^*_{2} = \frac{r_{2} - \alpha_{2,1}r_{1}}{1 - \alpha_{1,2} \alpha_{2,1}} \]

To make sure \(N^*>0\), with the condition that \(1-\alpha_{1,2}\alpha{2,1} >0\), we need to make sure

\[ r_{1} - \alpha_{1,2}r_{2}>0\\ r_{2} - \alpha_{2,1}r_{1}>0 \]

This region is the feasible domain of \(r_1,r_2\) here.

a12 <- pulse_a12_hat
a21 <- pulse_a21_hat

# 网格分辨率（越小越精细，但越慢）
step <- 0.01

# r1,r2 范围
r1_seq <- seq(-2, 2, by = step)
r2_seq <- seq(-2, 2, by = step)

# 生成网格并计算不等式
grid <- expand.grid(r1 = r1_seq, r2 = r2_seq)
grid$cond1 <- with(grid, r1 > a12 * r2)      # r1 - a12*r2 > 0
grid$cond2 <- with(grid, r2 > a21 * r1)      # r2 - a21*r1 > 0

# 目标点
pt <- data.frame(r1 = as.numeric(coef(fit2)["r"]), r2 = as.numeric(coef(fit5)["r"]))

# 绘图
p_pulse <- ggplot() +
  # 条件1 的区域（半透明蓝）
  geom_tile(
    data = subset(grid, cond1),
    aes(x = r1, y = r2),
    fill = "steelblue",
    alpha = 0.28
  ) +
  # 条件2 的区域（半透明红）
  geom_tile(
    data = subset(grid, cond2),
    aes(x = r1, y = r2),
    fill = "firebrick",
    alpha = 0.28
  ) +
  # 边界直线： r1 = a12 * r2 -> r2 = r1 / a12
  geom_abline(slope = 1 / a12, intercept = 0, color = "steelblue4", linetype = "dashed") +
  # 边界直线： r2 = a21 * r1 -> slope = a21
  geom_abline(slope = a21, intercept = 0, color = "firebrick4", linetype = "dashed") +
  # 标注目标点
  geom_point(data = pt, aes(x = r1, y = r2), color = "black", size = 3) +
  geom_text(data = pt, aes(x = r1, y = r2, label = "(r1,r2)"),
            nudge_x = 0.12, nudge_y = 0.05, size = 3.5, hjust = 0) +
  coord_equal(xlim = c(-2, 2), ylim = c(-2, 2), expand = FALSE) +
  labs(
    #title = expression("Regions: " ~ r[1] - alpha[1,2] * r[2] > 0 ~ "  and  " ~ r[2] - alpha[2,1] * r[1] > 0),
    title = "Npulse",
    x = expression(r[1]),
    y = expression(r[2])
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_line(color = "grey90"),
    legend.position = "none"
  )

print(p_pulse)

a12 <- limit_a12_hat
a21 <- limit_a21_hat

# 网格分辨率（越小越精细，但越慢）
step <- 0.01

# r1,r2 范围
r1_seq <- seq(-2, 2, by = step)
r2_seq <- seq(-2, 2, by = step)

# 生成网格并计算不等式
grid <- expand.grid(r1 = r1_seq, r2 = r2_seq)
grid$cond1 <- with(grid, r1 > a12 * r2)      # r1 - a12*r2 > 0
grid$cond2 <- with(grid, r2 > a21 * r1)      # r2 - a21*r1 > 0

# 目标点
pt <- data.frame(r1 = as.numeric(coef(fit_limit2)["r"]), r2 = as.numeric(coef(fit_limit5)["r"]))

# 绘图
p_limited <- ggplot() +
  # 条件1 的区域（半透明蓝）
  geom_tile(
    data = subset(grid, cond1),
    aes(x = r1, y = r2),
    fill = "steelblue",
    alpha = 0.28
  ) +
  # 条件2 的区域（半透明红）
  geom_tile(
    data = subset(grid, cond2),
    aes(x = r1, y = r2),
    fill = "firebrick",
    alpha = 0.28
  ) +
  # 边界直线： r1 = a12 * r2 -> r2 = r1 / a12
  geom_abline(slope = 1 / a12, intercept = 0, color = "steelblue4", linetype = "dashed") +
  # 边界直线： r2 = a21 * r1 -> slope = a21
  geom_abline(slope = a21, intercept = 0, color = "firebrick4", linetype = "dashed") +
  # 标注目标点
  geom_point(data = pt, aes(x = r1, y = r2), color = "black", size = 3) +
  geom_text(data = pt, aes(x = r1, y = r2, label = "(r1,r2)"),
            nudge_x = -0.2, nudge_y = 0.15, size = 3.5, hjust = 0) +
  coord_equal(xlim = c(-2, 2), ylim = c(-2, 2), expand = FALSE) +
  labs(
    #title = expression("Regions: " ~ r[1] - alpha[1,2] * r[2] > 0 ~ "  and  " ~ r[2] - alpha[2,1] * r[1] > 0),
    title = "Nlimited",
    x = expression(r[1]),
    y = expression(r[2])
  ) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_line(color = "grey90"),
    legend.position = "none"
  )

print(p_limited)

p_pulse | p_limited