summary-statistics.Rmd

---
title: "Summarizing SLiM output"
output: 
  github_document:
    toc: true
    toc_depth: 5
date: "2025-10-21"
---

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

## Note

* Make sure to open the zip files in the `data/` subdirectories to access the SLiM output data

## Load in functions and packages
```{r eval=F, message = F}
packages = c("openxlsx", "here") 
installed = packages %in% installed.packages()[, "Package"]
if (any(!installed)) {
  install.packages(packages[!installed])
}
lapply(packages, library, character.only = TRUE)

source(here("plotting/plotting-functions.R")) # this calls the other scripts
```

## 1. Critical bubble in SLiM

* Data directory: `Containment/data/underdominance-alleles/bubble`
* Summary Excel sheet: `Containment/data/summaries/critical-bubble-underdominance-alleles.xlsx`

* Create initial Excel sheet with paramset, param_no, D, phat, s. Also record:
  + Critical bubble AUC 
  + Corresponding release size under full critical bubble
  + AUC of bubble within SLiM boundaries
```{r eval=F, message=F}
param_no = c(3,6,8,9,11,13,14,18,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48)
paramset = paste0("paramset", param_no)
D = c(rep(1e-5,3),1e-4, rep(1e-5,2), 1e-4,1e-5,1e-4,1e-3,1e-5, 1e-5,1e-4, 1e-4,1e-3,1e-4,1e-2,1e-5,1e-2,1e-3,1e-3,1e-2,1e-4, 1e-4, 1e-3, 1e-3, 1e-2, 1e-4, 1e-3, 1e-2, 1e-3, 1e-2, 1e-2, 1e-3,1e-2, 1e-2)
phat = c(0.05, 0.05, 0.20, 0.05, 0.05, 0.20, 0.05, 0.40, 0.20, 0.05, 0.40, 0.20, 0.05, 0.20, 0.05, 0.40, 0.05, 0.40, 0.05, 0.20, 0.05, 0.05, 0.40, 0.20, 0.20, 0.40, 0.20, 0.40, 0.20, 0.20, 0.40, 0.20, 0.40, 0.40, 0.40, 0.40)
s = c(0.80, 0.30, 0.80, 0.80, 0.05, 0.30, 0.30, 0.80, 0.80, 0.80, 0.30, 0.05, 0.05, 0.30, 0.30, 0.80, 0.80, 0.05, 0.30, 0.80, 0.05, 0.05, 0.30, 0.05, 0.30, 0.80, 0.80, 0.05, 0.05, 0.30, 0.30, 0.05, 0.80, 0.05, 0.30, 0.05)
data = tibble(paramset, param_no, D, phat, s)

n = 36
full_cb_auc = rep(-1, n)
full_cb_size = rep(-1, n)
slim_cb_auc = rep(-1,n)
K = 100000
for (i in 1:n){
  phat = data$phat[i]; D = data$D[i]; s = data$s[i]; param.no = data$param_no[i]
  cb = create_bubble_csv(param_no = param.no, s=s, D=D, phat=phat, write_out = F)
  slim_cb_auc[i] = AUC(x = cb$x, y = cb$p)
  full_cb_auc[i] = underdominance_width(phat = phat, D=D, s=s)  
  full_cb_size[i] = full_cb_auc[i]*K
}

data2 = data %>% add_column(full_cb_auc, full_cb_size, slim_cb_auc) %>%
  mutate(slim_bubble_div_full_bubble = slim_cb_auc/full_cb_auc)

# add columns with placeholders
slim_cb_q0.05 = rep(-1,n)
slim_cb_q0.95 = rep(-1,n)
slim_cb_critical_release_size = rep(-1,n)

data3 = data2 %>% add_column(slim_cb_critical_release_size,
                             slim_cb_q0.05,
                             slim_cb_q0.95)


write.xlsx(data3,  here("data/summaries/critical-bubble-underdominance-alleles.xlsx"))
```

* Run SLiM on the cluster (make sure to output critical bubble csvs for SLiM using the `critical-bubble-frequencies.R` script)

* Loop through SLiM output data and get:
 + SLiM transition range
 + SLiM critical release size
```{r eval=F, message = F}
data3 = read.xlsx(here("data/summaries/critical-bubble-underdominance-alleles.xlsx")) %>% as_tibble()
n = 36 # paramsets

slim.results.dir = here("data/underdominance-alleles/bubble/")
for (i in 1:n){
  param.no = data2$param_no[i]; predicted = data2$full_cb_size[i]
  csv = read_csv(paste0(slim.results.dir, "paramset", param.no,"-cb-probabilistic-compiled.csv"))
  results = fit.L2.logistic.wrapper(csv = csv, predicted_size = predicted)
  transition = results$transition
  data3$slim_cb_q0.05[i] = transition$x.low
  data3$slim_cb_q0.95[i] = transition$x.high
  data3$slim_cb_critical_release_size[i] = results$critical_release_size
}

write.xlsx(data3,  here("data/summaries/critical-bubble-underdominance-alleles.xlsx"))
```

## 2. Rectangular release of underdominance alleles

* Data directory: `Containment/data/underdominance-alleles/rectangle1D`
* Summary Excel sheet: `Containment/data/summaries/rectangle1D-underdominance-alleles.xlsx`

Get critical bubble AUC & corresponding release size from `Containment/data/summaries/critical-bubble-underdominance-alleles.xlsx`
```{r eval=F, message=F}
data = read.xlsx(here("data/summaries/critical-bubble-underdominance-alleles.xlsx")) %>% as_tibble() %>% 
  dplyr::select(paramset, param_no, D, phat, s, full_cb_size)

# Add placeholder columns for deSolve rectangular prediction, SLiM transition range bounds, and SLiM critical release size
n = nrow(data)
deSolve_rectangle_critical_size = rep(-1,n)
slim_rectangle_q0.05 = rep(-1,n)
slim_rectangle_q0.95 = rep(-1,n)
slim_rectangle_critical_release_size = rep(-1,n)

data2 = data %>% add_column(deSolve_rectangle_critical_size,
                            slim_rectangle_q0.05,
                            slim_rectangle_q0.95,
                            slim_rectangle_critical_release_size)

write.xlsx(data2,  here("data/summaries/rectangle1D-underdominance-alleles.xlsx"))
```

* Get deSolve predictions using Rscript (note: this will take a long time to run)
* Call this from the command line in `Containment/deSolve-simulations`. deSolve predictions will directly go into the Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-predictions.R underdominance > underdominance-out.txt &
```

* Run SLiM on the cluster

* Loop through SLiM output data and get:
  + SLiM transition range
  + SLiM critical release size

```{r eval=F, message=F}
data = read.xlsx(here("data/summaries/rectangle1D-underdominance-alleles.xlsx")) %>% as_tibble() 
slim.results.dir = here("data/underdominance-alleles/rectangle1D/")

n = nrow(data)
for (i in 1:n){
  param.no = data$param_no[i]; predicted = data$deSolve_rectangle_critical_size[i]
  name = paste0("paramset", param.no,"-release-width-1D-compiled.csv")
  csv = read_csv(paste0(slim.results.dir, name))
  results = fit.L2.logistic.wrapper(csv = csv, predicted_size = predicted)
  transition = results$transition
  data$slim_rectangle_q0.05[i] = transition$x.low
  data$slim_rectangle_q0.95[i] = transition$x.high
  data$slim_rectangle_critical_release_size[i] = results$critical_release_size
}

write.xlsx(data,  here("data/summaries/rectangle1D-underdominance-alleles.xlsx"))
```

## 3. Rectangular release of homing underdominance gene drives

* Data directory: `Containment/data/homing-gene-drive/rectangle1D`
* Summary Excel sheet: `Containment/data/summaries/rectangle1D-homing-gene-drive.xlsx`

* Create initial Excel sheet with paramset, param_no, D, s, h, c, phat
```{r eval=F, message = F}
param_no = c(5,7,11, 13, 17, 19, 21, 23, 25, 27, 29, 31, 35, 36, 37, 38, 39, 40, 43, 44, 45, 46, 47, 48, 51, 52, 53, 54, 55, 56, 59, 60, 61, 62, 63, 64)
paramset = paste0("paramset", param_no)
D = c(1e-05, 1e-05, 1e-04, 1e-04, 1e-03, 1e-03, 1e-02, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-05, 1e-05, 1e-05, 1e-05,
      1e-05, 1e-04, 1e-04, 1e-04, 1e-04, 1e-04, 1e-04,1e-03, 1e-03, 1e-03, 1e-03, 1e-03, 1e-03, 1e-02, 1e-02, 1e-02, 1e-02, 1e-02, 1e-02)
s = c(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.10, 0.30, 0.10, 0.30, 0.10, 0.30, 0.10, 0.30, 
      0.10, 0.30, 0.10, 0.30, 0.10, 0.30, 0.10, 0.30, 0.10,0.30, 0.10, 0.30, 0.10, 0.30, 0.10, 0.30)
phat = c(0.10, 0.20, 0.10, 0.20, 0.10, 0.20, 0.10, 0.20, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.10, 0.10, 
         0.20, 0.20, 0.05, 0.05, 0.10, 0.10, 0.20, 0.20, 0.05, 0.05, 0.10, 0.10, 0.20, 0.20, 0.05, 0.05, 0.10, 0.10, 0.20, 0.20)
c = rep(1, 36)
h = solve_for_h_gene_drive(s = s, c = c, pg_hat = phat)

data = tibble(paramset, param_no, D, s, h, c, phat)

# create placeholder columns for deSolve prediction, SLiM transition, and SLiM critical release size
n = nrow(data)
deSolve_rectangle_critical_size = rep(-1,n)
slim_rectangle_q0.05 = rep(-1,n)
slim_rectangle_q0.95 = rep(-1,n)
slim_rectangle_critical_release_size = rep(-1,n)

data2 = data %>% add_column(deSolve_rectangle_critical_size,
                            slim_rectangle_q0.05,
                            slim_rectangle_q0.95,
                            slim_rectangle_critical_release_size)

write.xlsx(data2,  here("data/summaries/rectangle1D-homing-gene-drive.xlsx"))
```

* Get deSolve predictions using Rscript (note: this will take a long time to run)
* Call this from the command line in `Containment/deSolve-simulations`. deSolve predictions will directly go into the Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-predictions.R homing_gene_drive > homing_gene_drive_out.txt &
```

* Run SLiM on the cluster

* Loop through SLiM output data and get:
  + SLiM transition range
  + SLiM critical release size
```{r eval=F, message=F}
data = read.xlsx(here("data/summaries/rectangle1D-homing-gene-drive.xlsx")) %>% as_tibble() 
slim.results.dir = here("data/homing-gene-drive/rectangle1D/")

n = nrow(data)
for (i in 1:n){
  param.no = data$param_no[i]; predicted = data$deSolve_rectangle_critical_size[i]
  name = paste0("paramset", param.no,"-gd-release-width-compiled.csv")
  csv = read_csv(paste0(slim.results.dir, name))
  results = fit.L2.logistic.wrapper(csv = csv, predicted_size = predicted)
  transition = results$transition
  data$slim_rectangle_q0.05[i] = transition$x.low
  data$slim_rectangle_q0.95[i] = transition$x.high
  data$slim_rectangle_critical_release_size[i] = results$critical_release_size
}

write.xlsx(data,  here("data/summaries/rectangle1D-homing-gene-drive.xlsx"))
```

## 4. Rectangular release of TADE modification

* Data directory: `Containment/data/tade-modification/rectangle1D`
* Summary Excel sheet: `Containment/data/summaries/rectangle1D-tade-modification.xlsx`

* Create initial Excel sheet with paramset, param_no, D, s, h, c, phat
```{r eval=F, message = F}
param_no = c(5,  6,  7,  8, 13, 14, 15, 16, 21, 22, 23, 24, 26, 28, 30, 32)
paramset = paste0("paramset", param_no)


D = c(1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02)
s = c(0.05, 0.05, 0.05, 0.05, 0.30, 0.30, 0.30, 0.30, 0.50, 0.50, 0.50, 0.50, 0.15, 0.15, 0.15, 0.15)
n = length(D)
h = rep(0.5, n)
c = rep(1, n)

# get phat numerically
phat = rep(-1,n)
data = tibble(paramset, param_no, D, s, h, c, phat)

for (i in 1:n){
  s = data$s[i]; h = data$h[i]; c = data$c[i]
  this.phat = find_phat_tade(s = s, h = h, c = c, suppression = F,
                          max_time = 100, phat.tol = 1e-8)
  data$phat[i] = this.phat
}

# create placeholder columns for deSolve prediction, SLiM transition, and SLiM critical release size
deSolve_rectangle_critical_size = rep(-1,n)
slim_rectangle_q0.05 = rep(-1,n)
slim_rectangle_q0.95 = rep(-1,n)
slim_rectangle_critical_release_size = rep(-1,n)

data2 = data %>% add_column(deSolve_rectangle_critical_size,
                            slim_rectangle_q0.05,
                            slim_rectangle_q0.95,
                            slim_rectangle_critical_release_size)

write.xlsx(data2,  here("data/summaries/rectangle1D-tade-modification.xlsx"))
```

* Get deSolve predictions using Rscript (note: this will take a long time to run)
* Call this from the command line in `Containment/deSolve-simulations`. deSolve predictions will directly go into the Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-predictions.R tade_modification 1 > tade_modification1D_out.txt &
```

* Run SLiM on the cluster

* Loop through SLiM output data and get:
  + SLiM transition range
  + SLiM critical release size
```{r eval=F, message=F}
data = read.xlsx(here("data/summaries/rectangle1D-tade-modification.xlsx")) %>% as_tibble() 
slim.results.dir = here("data/tade-modification/rectangle1D/")

n = nrow(data)
for (i in 1:n){
  param.no = data$param_no[i]; predicted = data$deSolve_rectangle_critical_size[i]
  name = paste0("paramset", param.no,"-tade-mod-release-width-compiled.csv")
  csv = read_csv(paste0(slim.results.dir, name))
  results = fit.L2.logistic.wrapper(csv = csv, predicted_size = predicted)
  transition = results$transition
  data$slim_rectangle_q0.05[i] = transition$x.low
  data$slim_rectangle_q0.95[i] = transition$x.high
  data$slim_rectangle_critical_release_size[i] = results$critical_release_size
}

write.xlsx(data,  here("data/summaries/rectangle1D-tade-modification.xlsx"))
```

## 5. Circular release of TADE modification

* Data directory: `Containment/data/tade-modification/circle2D`
* Summary Excel sheet: `Containment/data/summaries/circle2D-tade-modification.xlsx`

* Get paramset, param_no, D, s, h, c, phat from `rectangle1D-tade-modification.xlsx`
```{r eval=F, message = F}
data = read.xlsx(here("data/summaries/rectangle1D-tade-modification.xlsx")) %>% as_tibble() %>%
  dplyr::select(paramset, param_no, D, s, h, c, phat)

# Create placeholder columns for 2D deSolve prediction, SLiM transition, and SLiM critical release size
deSolve_circle_critical_size = rep(-1,n)
slim_circle_q0.05 = rep(-1,n)
slim_circle_q0.95 = rep(-1,n)
slim_circle_critical_release_size = rep(-1,n)

data2 = data %>% add_column(deSolve_circle_critical_size,
                            slim_circle_q0.05,
                            slim_circle_q0.95,
                            slim_circle_critical_release_size)

write.xlsx(data2, here("data/summaries/circle2D-tade-modification.xlsx"))
```

* Get deSolve predictions using Rscript (note: this will take a long time to run)
* Call this from the command line in `Containment/deSolve-simulations`. deSolve predictions will directly go into the Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-predictions.R tade_modification 2 > tade_modification2D_out.txt &
```

* Run SLiM on the cluster

* Loop through SLiM output data and get:
  + SLiM transition range
  + SLiM critical release size
```{r eval=F, message=F}
data = read.xlsx(here("data/summaries/circle2D-tade-modification.xlsx")) %>% as_tibble() 
slim.results.dir = here("data/tade-modification/circle2D/")

n = nrow(data)
for (i in 1:n){
  param.no = data$param_no[i]; predicted = data$deSolve_circle_critical_size[i]
  name = paste0("paramset", param.no,"-tade-mod-release-diameter-compiled.csv")
  csv = read_csv(paste0(slim.results.dir, name))
  results = fit.L2.logistic.wrapper(csv = csv, predicted_size = predicted)
  transition = results$transition
  data$slim_circle_q0.05[i] = transition$x.low
  data$slim_circle_q0.95[i] = transition$x.high
  data$slim_circle_critical_release_size[i] = results$critical_release_size
}

write.xlsx(data,  here("data/summaries/circle2D-tade-modification.xlsx"))
```


## 6. Circular release of TADE suppression

* Data directory: `Containment/data/tade-suppression/circle2D`
* Summary Excel sheet: `Containment/data/summaries/circle2D-tade-suppression.xlsx`

* Create initial Excel sheet with paramset, param_no, D, s, h, c, phat
```{r eval=F, message = F}
param_no = c(5, 6,  7, 8, 13, 14, 15, 16, 26, 28, 30, 32)
paramset = paste0("paramset", param_no)

D = c(1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02, 1e-05, 1e-04, 1e-03, 1e-02)

s = c(0.05, 0.05, 0.05, 0.05, 0.30, 0.30, 0.30, 0.30, 0.15, 0.15, 0.15, 0.15)
n = length(D)
h = rep(0.5, n)
c = rep(1, n)

# get phat numerically
phat = rep(-1,n)
data = tibble(paramset, param_no, D, s, h, c, phat)

for (i in 1:n){
  s = data$s[i]; h = data$h[i]; c = data$c[i]
  this.phat = find_phat_tade(s = s, h = h, c = c, suppression = T,
                          max_time = 100, phat.tol = 1e-8)
  data$phat[i] = this.phat
}

# create placeholder columns for deSolve prediction
deSolve_circle_critical_size = rep(-1,n)

data2 = data %>% add_column(deSolve_circle_critical_size)

write.xlsx(data2,  here("data/summaries/circle2D-tade-suppression.xlsx"))
```

* Get deSolve predictions
  + Note: the "spread" condition is whether the final drive frequency exceeds the initial drive frequency OR if the final population size is less than (initial population size*(1 - release area)). The latter condition indicates that the wave expanded from the release area.

* Get deSolve predictions using Rscript (note: this will take a long time to run)
* Call this from the command line in `Containment/deSolve-simulations`. deSolve predictions will directly go into the Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-predictions.R tade_suppression 2 > tade_suppression2D_out.txt & # 36380
```

* After analyzing the SLiM results and running SLiM on the SLiMgui, we realized that "drive spread" was a lot more complex in the SLiM model. It's hard to say whether chasing effects constitute "drive spread" as this is fundamentally different to the deSolve "drive spread". Thus, we described SLiM results qualitatively, with series of simulation GUI screenshots. These panels are in `Containment/data/slim-screenshots/`
  + `paramset13` (subpanel A of Figure 7; D = 1e-5, s = 0.3)
  + `paramset5` (subpanel B of Figure 7; D = 1e-5, s = 0.05)
  + `paramset6` (subpanel C of Figure 7; D = 1e-4, s = 0.05)
  + `paramset15` (subpanel D of Figure 7; D = 1e-3, s = 0.3)
 
* Panels were arranged using Adobe Illustrator 2025

## 7. (Supplemental) Finding the critical p0 of an underdominance allele with a known release width

* This is for an underdominance allele in 1D with phat = 0.2, s = 0.8, and D = 1e-5, 1e-4, or 1e-3
* Summaries are located in: `Containment/data/summaries/critical-p0-underdominance/`

* Create release_width iterables for each value of D

* D = 1e-5
* Summary: `critical-p0-D1e-5.xlsx`
```{r eval=F, message=F}
this.s = 0.8
this.phat = 0.2
this.D = 1e-5

release_width_iterable = c(0.004589844, 0.008026572, 0.011311337, 0.014596102, 0.017880868, 0.021165633, 0.024450398, 0.027735163, 0.031019929, 0.034304694, 0.037589459, 0.040874225,
                           0.044158990, 0.047443755, 0.050728521, 0.054013286, 0.057298051, 0.060582816, 0.063867582, 0.067152347, 0.070437112, 0.073721878, 0.077006643, 0.080291408,
                           0.083576174, 0.086860939, 0.090145704, 0.093430469, 0.096715235, 0.100000000, 0.131034483, 0.162068966, 0.193103448, 0.224137931, 0.255172414, 0.286206897,
                           0.317241379, 0.348275862, 0.379310345, 0.410344828, 0.441379310, 0.472413793, 0.503448276, 0.534482759, 0.565517241, 0.596551724, 0.627586207, 0.658620690,
                           0.689655172, 0.720689655, 0.751724138, 0.782758621, 0.813793103, 0.844827586, 0.875862069, 0.906896552, 0.937931034, 0.968965517, 1.000000000)
n = length(release_width_iterable)

data = tibble(D = rep(this.D, n),s = rep(this.s, n), phat = rep(this.phat, n),
              release_width = release_width_iterable)

# Add placeholder for critical p0
critical_p0 = rep(-1, n)
data2 = data %>% add_column(critical_p0)

write.xlsx(data2,  here("data/summaries/critical-p0-underdominance/critical-p0-D1e-5.xlsx"))
```

* D = 1e-4
* Summary: `critical-p0-D1e-4.xlsx`
```{r eval=F, message=F}
this.s = 0.8
this.phat = 0.2
this.D = 1e-4

release_width_iterable = c(0.014511528, 0.01756417, 0.02050831, 0.02345245, 0.02639658, 0.02934072, 0.03228486, 0.03522899, 0.03817313, 0.04111727, 0.04406140, 0.04700554, 0.04994968, 0.05289381, 0.05583795, 0.05878209,
                           0.06172622, 0.06467036, 0.06761450, 0.07055863, 0.07350277, 0.07644691, 0.07939104, 0.08233518, 0.08527932, 0.08822345, 0.09116759, 0.09411173, 0.09705586, 0.10000000, 0.13103448,
                           0.16206897, 0.19310345, 0.22413793, 0.25517241, 0.28620690, 0.31724138, 0.34827586, 0.37931034, 0.41034483, 0.44137931, 0.47241379, 0.50344828, 0.53448276,
                           0.56551724, 0.59655172, 0.62758621, 0.65862069, 0.68965517, 0.72068966, 0.75172414, 0.78275862, 0.81379310, 0.84482759, 0.87586207, 0.90689655, 0.93793103, 0.96896552, 1.00000000)

n = length(release_width_iterable)

data = tibble(D = rep(this.D, n),s = rep(this.s, n), phat = rep(this.phat, n),
              release_width = release_width_iterable)

# Add placeholder for critical p0
critical_p0 = rep(-1, n)
data2 = data %>% add_column(critical_p0)

write.xlsx(data2,  here("data/summaries/critical-p0-underdominance/critical-p0-D1e-4.xlsx"))
```

* D = 1e-3
* Summary: `critical-p0-D1e-3.xlsx`
```{r eval=F, message=F}
this.s = 0.8
this.phat = 0.2
this.D = 1e-3

release_width_iterable = c(0.045883749, 0.04837280, 0.05021663, 0.05206046, 0.05390429, 0.05574812, 0.05759195, 0.05943577, 0.06127960, 0.06312343, 0.06496726, 0.06681109,
                           0.06865492, 0.07049874, 0.07234257, 0.07418640, 0.07603023, 0.07787406, 0.07971789, 0.08156172, 0.08340554, 0.08524937,
                           0.08709320, 0.08893703, 0.09078086, 0.09262469, 0.09446851, 0.09631234, 0.09815617, 0.10000000, 0.13103448, 0.16206897,
                           0.19310345, 0.22413793, 0.25517241, 0.28620690, 0.31724138, 0.34827586, 0.37931034, 0.41034483, 0.44137931, 0.47241379,
                           0.50344828, 0.53448276, 0.56551724, 0.59655172, 0.62758621, 0.65862069, 0.68965517, 0.72068966, 0.75172414, 0.78275862,
                           0.81379310, 0.84482759, 0.87586207, 0.90689655, 0.93793103, 0.96896552, 1.00000000)

n = length(release_width_iterable)

data = tibble(D = rep(this.D, n),s = rep(this.s, n), phat = rep(this.phat, n),
              release_width = release_width_iterable)

# Add placeholder for critical p0
critical_p0 = rep(-1, n)
data2 = data %>% add_column(critical_p0)


write.xlsx(data2,  here("data/summaries/critical-p0-underdominance/critical-p0-D1e-3.xlsx"))
```

* Get deSolve critical width prediction (note: this is slow)
* This writes the critical p0 directly to the corresponding Excel sheet
```{bash eval=F}
nohup Rscript get-deSolve-critical-p0-predictions.R 0.00001 > critical-p0-D1e-5-out.txt & # 49072
nohup Rscript get-deSolve-critical-p0-predictions.R 0.0001 > critical-p0-D1e-4-out.txt & # 49095
nohup Rscript get-deSolve-critical-p0-predictions.R 0.001 > critical-p0-D1e-3-out.txt & # 45458
```