prediction-functions.R

## This script is used for getting deSolve (numeric) predictions for the minimum release width or
## release diameter of underdominance, homing gene drive, TADE modification, or TADE suppression systems.
## It also finds the phat of TADE modification and TADE suppression systems numerically.
if (!requireNamespace("here", quietly = TRUE)) install.packages("here")
library(here)

source(here("deSolve-simulations/underdominance-functions.R"))
source(here("deSolve-simulations/homing-gd-functions.R"))
source(here("deSolve-simulations/tade-modification-functions.R"))
source(here("deSolve-simulations/tade-suppression-functions.R"))

###### Functions for quickly running simulations and returning whether the drive spread ##########

panmictic_tade_spreads_at_p0 = function(s, h, c, pd_overall, max_time, tol, suppression){
  # Arguments:
  #   s: selection coefficient (dd fitness = 1 - s, dw fitness = 1 - hs)
  #   h: dominance coefficient
  #   c: rate at which drive cuts "a" nondisrupted target alleles in the germline
  #   pd_overall: the overall drive allele introduction frequency. Assuming a heterozygous release, pdw = 2*pd_overall.
  #   max_time: number of timesteps to simulate
  #   tol: tolerance to use for comparisons; drive spreads if (pfinal > (pd_overall + tol))
  #   suppression: if T, simulate TADE suppression. If F, simulate TADE modification
  #
  # Returns:
  #   Logical of whether the drive spread
  
  if (!suppression){
    spreads = simulate_panmictic_tade_modification(s = s, h = h, c = c, pd_overall = pd_overall, release_heterozygotes = T,
                                                   max_time = max_time, show_frequencies = F, tol = tol)
  } else {
    spreads = simulate_panmictic_tade_suppression(s = s, h = h, c = c, pd_overall = pd_overall, max_time = max_time,
                                                  show_frequencies = F, tol = tol)
  }
  return(spreads)
}


rde_drive_spreads_at_release_stat = function(system, dimension, drive_params, D,release_stat, p0, N, max_time, tol, verbose){
  # For a 1D RDE, determines whether the system spreads at this release width
  # Arguments:
  #   -system: this is one of "underdominance", "homing_gene_drive", "tade_modification", or "tade_suppression"
  #   -dimension: 1 or 2
  #   -drive_params:
  #       (1) if "underdominance", this is a list of "s", "phat", "cubic_approximation"
  #       (2) if "homing_gene_drive", this is a list of "s", "h", "c", "cubic_approximation"
  #       (3) if "tade_modification", this is a list of "s", "h", "c", "litter_size"
  #       (4) if "tade_suppression", this is a list of "s", "h", "c", "litter_size", "capacity"
  #   -D: diffusion constant
  #   -release_stat: if dimension=1, this is the release_width. 
  #                  if dimension=2, this is the release_diameter.
  #   -p0: the frequency of drive INDIVIDUALS within the release width
  #       note: for "underdominance" and "homing_gene_drive", homozygotes are released (pd = p0 within the release width)
  #             for "tade_modification" and "tade_suppression", heterozygotes are released (pd = 0.5*p0 within the release width)
  #   -N: if dimension=1, the number of slices (dx = 1/N). 
  #       if dimension=2, the number of slices in the x-direction and y-direction (dx = dy = 1/N; N*N cells total)
  #   -max_time: number of time steps to simulate
  #   -tol: for modification systems, drive spreads if (pd.last > (p.d.full.0 + tol)). 
  #         for suppression, drive spreads if (pd.last > (p.d.full.0 + tol)) or if (N.last < (capacity*(1 - release_area) - tol)) 
  #   -verbose: whether to print results
  #
  # Returns:
  #   whether the drive spreads in the RDE (T or F)
  
  if (system == "underdominance"){
    s = drive_params$s
    phat = drive_params$phat
    cubic_approximation = drive_params$cubic_approximation
    
    if (dimension == 1){
      spreads = simulate_underdominance_rde_1D(phat = phat, D = D, s = s, 
                                       release_width = release_stat, 
                                       bubble = F,
                                       p0 = p0,
                                       N = N, max_time = max_time,
                                       cubic_approximation = cubic_approximation, 
                                       plot = F, 
                                       return_out = F,
                                       tol = tol)  
    } else {
      spreads = simulate_underdominance_rde_2D(phat = phat, D = D, s = s, 
                                       release_diameter = release_stat,
                                       p0 = p0, 
                                       N = N,
                                       max_time = max_time,
                                       cubic_approximation = cubic_approximation,
                                       plot = F, return_out = F,
                                       tol = tol)
    }
  } else if (system == "homing_gene_drive"){
    s = drive_params$s
    h = drive_params$h
    c = drive_params$c
    cubic_approximation = drive_params$cubic_approximation
    
    if (dimension == 1){
      spreads = simulate_homing_gene_drive_rde_1D(s = s, h = h, c = c, D = D,
                                                  release_width = release_stat,
                                                  bubble = F, p0 = p0, 
                                                  N = N, max_time = max_time,
                                                  cubic_approximation = cubic_approximation, 
                                                  plot = F, 
                                                  return_out = F,
                                                  tol = tol)  
    } else {
      spreads = simulate_homing_gene_drive_rde_2D(s = s, h = h, c = c, D = D,
                                                  release_diameter = release_stat,
                                                  p0 = p0, 
                                                  N = N, max_time = max_time,
                                                  cubic_approximation = cubic_approximation,
                                                  plot = F, return_out = F, tol = tol)
    }
  } else if (system == "tade_modification"){
    s = drive_params$s
    h = drive_params$h
    c = drive_params$c
    litter_size = drive_params$litter_size
    
    if (dimension == 1){
      spreads = simulate_tade_modification_rde_1D(s = s, h = h, c = c, D = D,
                                                  release_width = release_stat, 
                                                  N = N, 
                                                  max_time = max_time, 
                                                  litter_size = litter_size,
                                                  p0 = p0, 
                                                  release_heterozygotes = T,
                                                  plot = F, return_out = F,
                                                  tol = tol)    
    } else {
      spreads = simulate_tade_modification_rde_2D(s = s, h = h, c = c, D = D,
                                                  release_diameter = release_stat,
                                                  N = N,
                                                  max_time = max_time, litter_size = litter_size,
                                                  p0 = p0, release_heterozygotes = T,
                                                  plot = F, return_out = F, tol = tol)
    }
  } else if (system == "tade_suppression"){
    s = drive_params$s
    h = drive_params$h
    c = drive_params$c
    litter_size = drive_params$litter_size
    capacity = drive_params$capacity
    
    if (dimension == 1){
      spreads = simulate_tade_suppression_rde_1D(s = s, h = h, c = c, D = D,
                                                 release_width = release_stat, 
                                                 N = N, max_time = max_time,
                                                 capacity = capacity, litter_size = litter_size,
                                                 p0 = p0, plot = F, return_out = F, tol = tol)
    } else {
      spreads = simulate_tade_suppression_rde_2D(s = s, h = h, c = c, D = D,
                                                 release_diameter = release_stat, 
                                                 N = N, max_time = max_time,
                                                 capacity = capacity, litter_size = litter_size,
                                                 p0 = p0, plot = F, return_out = F, tol = tol)
    }
    
  } else {
    stop(paste0(system, " is not one of: underdominance, homing_gene_drive, tade_modification, tade_suppression"))
  }
  
  if (verbose){
    print(paste0(system, " system ", ifelse(spreads, "DOES spread", "DOESN'T spread"), 
                 " at release ",ifelse(dimension == 1, "width", "diameter")," = ",release_stat," and p0 = ", p0))
  }
  
  return(spreads)
}



######################################################################################################

find_phat_tade = function(s, h, c, suppression = F,
                          low = 0.01, high = 0.5, 
                          spreads_low = NULL, spreads_high=NULL,
                          phat.tol = 1e-5, condition.tol = 0,
                          max_time = 100, step = 0.01){
  
  # Finds phat numerically through a binary search algorithm
  # Note: "p0" in this context refers to the overall introduction frequency of the drive. 2*p0 is the starting frequency of heterozygotes.
  #
  # Arguments:
  #   s: selection coefficient (dd fitness = 1 - s, dw fitness = 1 - hs)
  #   h: dominance coefficient
  #   c: rate at which drive cuts "a" nondisrupted target alleles in the germline
  #   suppression: if T, this is for TADE suppression. If F, this is for TADE modification.
  #   low: the lower guess for phat
  #   high: the higher guess for phat
  #   spreads_low: whether the drive spreads at p0 = low, or NULL if not tested yet
  #   spreads_high: whether the drive spreads at p0 = high, or NULL if not tested yet
  #   phat.tol: if the difference in p0 values at which the drive doesn't spread and does spread is less than this, return the p0 at which it does spread
  #   condition.tol: for drive to have "spread", pfinal > (p0 + condition.tol)
  #   max_time: number of timesteps for the panmictic simulations to determine if drive spreads at a p0
  #   step: if drive doesn't spread at both low and high, next low = high and high = low - step
  #         and if drive does spread at both low and high, next high = low and low = high - step.
  #
  # Returns:
  #   numeric guess of phat, the minimum p0 at which the drive spreads 
  #
  # Assumptions:
  #   -Random mating / no spatial structure
  #   -Drive is released in heterozygotes, so p0 cannot exceed 0.5
  
  if (is.null(spreads_low)){
    spreads_low = panmictic_tade_spreads_at_p0(s = s, h = h, c = c, pd_overall = low, max_time = max_time, 
                                               tol = condition.tol, suppression = suppression)
  }
  
  if (is.null(spreads_high)){
    spreads_high = panmictic_tade_spreads_at_p0(s = s, h = h, c = c, pd_overall = high, max_time = max_time, 
                                               tol = condition.tol, suppression = suppression)
  }
  
  
  if (spreads_low & spreads_high){
    high = low
    low = high - step
    if (low < 0){
      low = 0 # never spreads at 0, so if spreads_low = F and spreads_high = T, should use midpoint
    }
    
    phat = find_phat_tade(low = low, high = high,
                          spreads_low = NULL, spreads_high=1,
                          step = step,
                          s = s, h = h, c = c, suppression = suppression,
                          phat.tol = phat.tol, condition.tol = condition.tol, max_time = max_time)
    
    
  } else if (!spreads_low & !spreads_high){
    
    if (near(high, 0.5)){
      return(NA) # drive can't spread
    }
    
    low = high
    high = low + step
    
    if (high > 0.5){
      high = 0.5 # max introduction; if drive can't spread here, it can never spread
    }
    
    phat = find_phat_tade(low = low, high = high, 
                          spreads_low = 0, spreads_high = NULL,
                          step = step,
                          s = s, h = h, c = c, suppression = suppression,
                          phat.tol = phat.tol, condition.tol = condition.tol, max_time = max_time)
    
  } else {
    # doesn't spread low but does spread high
    midpt = (low+high)/2
    
    spreads_midpt = panmictic_tade_spreads_at_p0(s = s, h = h, c = c, pd_overall = midpt, max_time = max_time, 
                                                tol = condition.tol, suppression = suppression)
    
    if (spreads_midpt){
      if (high - midpt < phat.tol){
        return(midpt)
      }
      
      phat = find_phat_tade(low = low, high = midpt, 
                            spreads_low = 0, spreads_high = 1,
                            step = step,
                            s = s, h = h, c = c, suppression = suppression,
                            phat.tol = phat.tol, condition.tol = condition.tol, max_time = max_time)
      
    } else {
      if (high-midpt < phat.tol){
        return(high)
      }
      
      phat = find_phat_tade(low = midpt, high = high,
                            spreads_low = 0, spreads_high = 1,
                            step = step,
                            s = s, h = h, c = c, suppression = suppression,
                            phat.tol = phat.tol, condition.tol = condition.tol, max_time = max_time)
      
    }
  }
}

reminders = function(system){
  # Arguments:
  #   -system: this is one of "underdominance", "homing_gene_drive", "tade_modification", or "tade_suppression"
  # Returns:
  #   No returns, just prints to the screen (1) the parameters that need to be in drive_params when calling find_critical_rde_release_stat
  #   and (2) the default N used in 1D D=1e-5 simulations
  
  if (system == "underdominance"){
    print("drive_params should be a list of s, phat, cubic_approximation")
    print("dd are always released")
    print("baseline N = 25000")
  } else if (system == "homing_gene_drive"){
    print("drive_params should be a list of s, h, c, cubic_approximation")
    print("dd are always released")
    print("baseline N = 25000")
  } else if (system == "tade_modification"){
    print("drive_params should be a list of s, h, c, litter_size")
    print("dwAa are always released")
    print("baseline N = 10000")
  } else if (system == "tade_suppression") {
    print("drive_params should be a list of s, h, c, litter_size, capacity")
    print("dwAa are always released")
    print("baseline N = 25000")
  } else {
    stop(paste0(system, " is not one of: underdominance, homing_gene_drive, tade_modification, tade_suppression"))
  }
  
}


find_critical_rde_release_stat = function(system, 
                                          dimension, N,
                                          drive_params, D, 
                                          low = 0.01, high = 0.5, 
                                          spreads_low = NULL, spreads_high=NULL,
                                          step = 0.01,
                                          p0 = 1, max_time = 100, condition.tol = 0,
                                          stat.tol = 1e-5, verbose = T){
  # Finds the critical RDE release width (if dimension = 1) or release diameter (if dimension = 2) numerically through a binary search algorithm
  #
  # Arguments:
  #   -system: this is one of "underdominance", "homing_gene_drive", "tade_modification", or "tade_suppression"
  #   -dimension: 1 or 2
  #   -N: if dimension=1, the number of slices (dx = 1/N). 
  #       if dimension=2, the number of slices in the x-direction and y-direction (dx = dy = 1/N; N*N cells total)
  #   -drive_params:
  #       (1) if "underdominance", this is a list of "s", "phat", "cubic_approximation"
  #       (2) if "homing_gene_drive", this is a list of "s", "h", "c", "cubic_approximation"
  #       (3) if "tade_modification", this is a list of "s", "h", "c", "litter_size"
  #       (4) if "tade_suppression", this is a list of "s", "h", "c", "litter_size", "capacity"
  #   -D: diffusion constant
  #   -low: the lower guess for the critical release width (if dimension=1) or diameter (if dimension=2)
  #   -high: the higher guess for the critical release width (if dimension=1) or diameter (if dimension=2)
  #   -spreads_low: whether the drive spreads at low, or NULL if not tested yet
  #   -spreads_high: whether the drive spreads at high, or NULL if not tested yet
  #   -step: if drive doesn't spread at both low and high, next low = high and high = low - step
  #         and if drive does spread at both low and high, next high = low and low = high - step.
  #   -p0: the frequency of drive INDIVIDUALS within the release width
  #       note: for "underdominance" and "homing_gene_drive", homozygotes are released (pd = p0 within the release width)
  #             for "tade_modification" and "tade_suppression", heterozygotes are released (pd = 0.5*p0 within the release width)
  #   -max_time: number of time steps to simulate
  #   -tol: for modification systems, drive spreads if (pd.last > (p.d.full.0 + tol)). 
  #         for suppression, drive spreads if (pd.last > (p.d.full.0 + tol)) or if (N.last < (capacity*(1 - release_area) - tol)) 
  #   -condition.tol: for drive to have "spread", pfinal > (p0 + condition.tol) or, for TADE suppression, Nfinal < (capacity*(1 - release_area) - tol)
  #   -stat.tol: if the difference in release_stat values at which the drive doesn't spread and does spread is less than this value, 
  #             then return the minimum release_stat at which it does spread
  #   -verbose: whether to print results for each release_stat tested  
  #
  # Returns:
  #   numeric guess of the critical release width (if dimension=1) or diameter (if dimension=2) for the system, based on its reaction-diffusion model
  #
  # Assumptions:
  #   -See the system Rscript for system-specific assumptions
  
  if (is.null(spreads_low)){
    spreads_low = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                    D = D, release_stat = low, p0 = p0, N = N, max_time= max_time,
                                                    tol = condition.tol, verbose = verbose)
  }
  
  if (is.null(spreads_high)){
    spreads_high = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                    D = D, release_stat = high, p0 = p0, N = N, max_time= max_time,
                                                    tol = condition.tol, verbose = verbose)
  }
  
  
  if (spreads_low & spreads_high){
    high = low
    low = high - step
    if (low < 0){
      low = 0 # never spreads at 0, so if spreads_low = F and spreads_high = T, should use midpoint
    }
    
    critical_stat = find_critical_rde_release_stat(low = low, high = high, 
                                                   spreads_low = NULL, spreads_high=1,
                                                   system = system,
                                                   dimension = dimension, 
                                                   N = N,
                                                   drive_params = drive_params, 
                                                   D = D, 
                                                   step = step,
                                                   p0 = p0, 
                                                   max_time = max_time, 
                                                   condition.tol = condition.tol,
                                                   stat.tol = stat.tol, 
                                                   verbose = verbose)
    
  } else if (!spreads_low & !spreads_high){
    
    if (near(high, 1)){
      return(NA) # drive can't spread since width or diameter can't exceed arena length
    }
    
    low = high
    high = low + step
    
    if (high > 1){
      high = 1 # max introduction; if drive can't spread here, it can never spread
    }
    
    critical_stat = find_critical_rde_release_stat(low = low, high = high, 
                                                   spreads_low = 0, spreads_high=NULL,
                                                   system = system,
                                                   dimension = dimension, 
                                                   N = N,
                                                   drive_params = drive_params, 
                                                   D = D, 
                                                   step = step,
                                                   p0 = p0, 
                                                   max_time = max_time, 
                                                   condition.tol = condition.tol,
                                                   stat.tol = stat.tol, 
                                                   verbose = verbose)
    
    
  } else {
    # doesn't spread low but does spread high
    midpt = (low+high)/2
    
    spreads_midpt = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                     D = D, release_stat = midpt, p0 = p0, N = N, max_time= max_time,
                                                     tol = condition.tol, verbose = verbose)
    
    if (spreads_midpt){
      if (high - midpt < stat.tol){
        return(midpt)
      }
      
      critical_stat = find_critical_rde_release_stat(low = low, high = midpt, 
                                                     spreads_low = 0, spreads_high=1,
                                                     system = system,
                                                     dimension = dimension, 
                                                     N = N,
                                                     drive_params = drive_params, 
                                                     D = D, 
                                                     step = step,
                                                     p0 = p0, 
                                                     max_time = max_time, 
                                                     condition.tol = condition.tol,
                                                     stat.tol = stat.tol, 
                                                     verbose = verbose)      
      
      
    } else {
      if (high-midpt < stat.tol){
        return(high)
      }
      
      critical_stat = find_critical_rde_release_stat(low = midpt, high = high, 
                                                     spreads_low = 0, spreads_high=1,
                                                     system = system,
                                                     dimension = dimension, 
                                                     N = N,
                                                     drive_params = drive_params, 
                                                     D = D, 
                                                     step = step,
                                                     p0 = p0, 
                                                     max_time = max_time, 
                                                     condition.tol = condition.tol,
                                                     stat.tol = stat.tol, 
                                                     verbose = verbose)       
      
    }
  }
}


# for finding the critical p0 given a width or diameter
find_critical_rde_p0 = function(system, dimension,
                                release_stat,
                                N,
                                drive_params, D, 
                                low = 0.01, high = 1, 
                                spreads_low = NULL, spreads_high=NULL,
                                step = 0.01, max_time = 100, condition.tol = 0,
                                stat.tol = 1e-5, verbose = T){
  # Finds the critical RDE release p0 for a given release width (if dimension = 1) or 
  # release diameter (if dimension = 2) numerically through a binary search algorithm
  #
  # Note: For the underdominance system and homing_gene_drive system, p0 corresponds to pd.0 since d/d are released.
  #       But for tade_modification and tade_suppression, heterozygotes are released. p0 = frequency of heterozygotes within
  #       the release width. This means that p0 = 1 corresponds to pd.0 = 0.5 within the release shape. 
  #
  # Arguments:
  #   -system: this is one of "underdominance", "homing_gene_drive", "tade_modification", or "tade_suppression"
  #   -dimension: 1 or 2
  #   -release_stat: the fixed release width if 1D or release diameter if 2D
  #   -N: if dimension=1, the number of slices (dx = 1/N). 
  #       if dimension=2, the number of slices in the x-direction and y-direction (dx = dy = 1/N; N*N cells total)
  #   -drive_params:
  #       (1) if "underdominance", this is a list of "s", "phat", "cubic_approximation"
  #       (2) if "homing_gene_drive", this is a list of "s", "h", "c", "cubic_approximation"
  #       (3) if "tade_modification", this is a list of "s", "h", "c", "litter_size"
  #       (4) if "tade_suppression", this is a list of "s", "h", "c", "litter_size", "capacity"
  #   -D: diffusion constant
  #   -low: the lower guess for the critical p0 
  #   -high: the higher guess for the critical p0
  #   -spreads_low: whether the drive spreads at low, or NULL if not tested yet
  #   -spreads_high: whether the drive spreads at high, or NULL if not tested yet
  #   -step: if drive doesn't spread at both low and high, next low = high and high = low - step
  #         and if drive does spread at both low and high, next high = low and low = high - step.
  #   -max_time: number of time steps to simulate
  #   -tol: for modification systems, drive spreads if (pd.last > (p.d.full.0 + tol)). 
  #         for suppression, drive spreads if (pd.last > (p.d.full.0 + tol)) or if (N.last < (capacity*(1 - release_area) - tol)) 
  #   -condition.tol: for drive to have "spread", pfinal > (p0 + condition.tol) or, for TADE suppression, Nfinal < (capacity*(1 - release_area) - tol)
  #   -stat.tol: if the difference in release_stat values at which the drive doesn't spread and does spread is less than this value, 
  #             then return the minimum release_stat at which it does spread
  #   -verbose: whether to print results for each release_stat tested  
  #
  # Returns:
  #   numeric guess of the critical release p0 for the system, based on its reaction-diffusion model. If no p0 can allow spread, returns NA. 
  #
  # Assumptions:
  #   -See the system Rscript for system-specific assumptions
  
  if (is.null(spreads_low)){
    spreads_low = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                    D = D, release_stat = release_stat, p0 = low, N = N, max_time= max_time,
                                                    tol = condition.tol, verbose = verbose)
  }
  
  if (is.null(spreads_high)){
    spreads_high = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                     D = D, release_stat = release_stat, p0 = high, N = N, max_time= max_time,
                                                     tol = condition.tol, verbose = verbose)
  }
  
  
  if (spreads_low & spreads_high){
    high = low
    low = high - step
    if (low < 0){
      low = 0 # never spreads at 0, so if spreads_low = F and spreads_high = T, should use midpoint
    }
    
    critical_stat = find_critical_rde_p0(low = low, high = high,
                                         spreads_low = NULL, spreads_high=1,
                                         system = system, dimension = dimension,
                                         release_stat = release_stat,
                                         N = N,
                                         drive_params = drive_params, D = D, 
                                         step = step, max_time = max_time, condition.tol = condition.tol,
                                         stat.tol = stat.tol, verbose = verbose)
    
  } else if (!spreads_low & !spreads_high){
    
    if (near(high, 1)){
      return(NA) # drive can't spread at all
    }
    
    low = high
    high = low + step
    
    if (high > 1){
      high = 1 # max introduction; if drive can't spread here, it can never spread
    }
    
    critical_stat = find_critical_rde_p0(low = low, high = high,
                                         spreads_low = 0, spreads_high=NULL,
                                         system = system, dimension = dimension,
                                         release_stat = release_stat,
                                         N = N,
                                         drive_params = drive_params, D = D, 
                                         step = step, max_time = max_time, condition.tol = condition.tol,
                                         stat.tol = stat.tol, verbose = verbose)
    
    
  } else {
    # doesn't spread low but does spread high
    midpt = (low+high)/2
    
    spreads_midpt = rde_drive_spreads_at_release_stat(system = system, dimension = dimension, drive_params = drive_params,
                                                      D = D, release_stat = release_stat, p0 = midpt, N = N, max_time= max_time,
                                                      tol = condition.tol, verbose = verbose)
    
    
    if (spreads_midpt){
      if (high - midpt < stat.tol){
        return(midpt)
      }
      
      critical_stat = find_critical_rde_p0(low = low, high = midpt,
                                           spreads_low = 0, spreads_high=1,
                                           system = system, dimension = dimension,
                                           release_stat = release_stat,
                                           N = N,
                                           drive_params = drive_params, D = D, 
                                           step = step, max_time = max_time, condition.tol = condition.tol,
                                           stat.tol = stat.tol, verbose = verbose)
      
      
    } else {
      if (high-midpt < stat.tol){
        return(high)
      }
      
      critical_stat = find_critical_rde_p0(low = midpt, high = high,
                                           spreads_low = 0, spreads_high=NULL,
                                           system = system, dimension = dimension,
                                           release_stat = release_stat,
                                           N = N,
                                           drive_params = drive_params, D = D, 
                                           step = step, max_time = max_time, condition.tol = condition.tol,
                                           stat.tol = stat.tol, verbose = verbose)
      
    }
  }
}