deSolve-simulations

deSolve simulations README

For details on the deSolve R package, visit: https://cran.r-project.org/web/packages/deSolve/index.html

Setting the number of grid cells in ode.1D and ode.2D

• To simulate reaction-diffusion waves in deSolve, you must specify N, the number of grid cells along the x or y axis (each axis has a length of 1).

  • Note: in 2D, the total number of grid cells is N*N

• The higher N is, the more "local" dynamics become. However, setting N very high has two major drawbacks:

  1. The simulation runs more slowly.
  2. The diffusion term scales with D/(dx^2) where dx = 1/N, so a high N and high D creates a very large multiplier, potentially leading to numerical instability or non-convergence issues.

• To deal with these problems, we first find a good baseline N for our lowest D (1e-5) in 1D. This is high enough to capture local dynamics but low enough for the simulation to run in a reasonable amount of time. Then we calculate the baseline "flux" coefficient: (1e-5)*(baseline N)^2. For higher values of D, we lower N such that the "flux" coefficient remains equal to the baseline.

• In 2D, we set N along each axis such that the total number of grid cells (N^2) matches the 1D case for that D.

• The baseline N values used for D = 1e-5 1D simulations were:

  • 25000 for underdominance alleles, homing underdominance gene drives, and TADE suppression drives
  • 10000 for TADE modification drives (the slowest system to simulate)

• To see this N assignment in practice, check out get-deSolve-predictions.R

Finding the critical width or diameter with deSolve simulations

• Use the find_critical_rde_release_stat function in prediction-functions.R to find the critical width or diameter for one of the following systems:

  • "underdominance" (then the drive_params parameter should be a list of s, phat, and cubic_approximation = T or F)
  • "homing_gene_drive" (then the drive_params parameter should be a list of s, h, c, and cubic_approximation = T or F)
  • "tade_modification" (then the drive_params parameter should be a list of s, h, c, and litter_size)
  • "tade_suppression" (then the drive_params parameter should be a list of s, h, c, litter_size, and capacity)
  • Note: for reminders on what must be included inside the drive_params list, check out the reminders function in prediction-functions.R
  • Check out the specs of find_critical_rde_release_stat for more details on its inputs and outputs

• find_critical_rde_release_stat operates using a binary search algorithm. We are looking for the minimum width or diameter at which the system spreads. See the figure below:

• We start by guessing a "low" and "high" width or diameter. We check whether the drive spreads at "low" and "high" using the system-specific deSolve functions. There's 3 possible outcomes:

  1. Both produce drive spread (ie both points are on the upper plateau in the plot)
     • Action: go down -- set the new "high" equal to this "low". Set the new "low" equal to the "high" minus a small step. Rerun the function.
  2. Neither produce drive spread (ie both points are on the lower plateau in the plot)
     • Action: go up -- set the new "low" equal to this "high". Set the new "high" equal to "low" plus a small step. Rerun the function.
  3. "low" does not produce drive spread but "high" does produce drive spread (ie "low" is on the lower plateau and "right" is on the upper plateau in the plot)
     • Action: calculate the midpoint between "high" and "low". Check if the drive spreads here.
       • If it doesn't spread at the midpoint, keep "high" the same as before, but set new "low" equal to the midpoint. Rerun the function.
       • If it does spread at the midpoint, keep "low" the same as before, but set the new "high" equal to the midpoint. Rerun the function.

• The end conditions are:

  • If the drive can't spread from a width or diameter of 1, it can never spread, so return NA
  • If the distance between the midpoint and high or low is less than the tolerance (stat.tol):
    • Return the midpoint if it allowed for spread.
    • Return the "high" point if the midpoint didn't allow for spread.

Example

• To find the critical width of an underdominance allele with phat = 0.2, s = 0.3, and D = 1e-5, you could do the following:
  • width = find_critical_rde_release_stat(system = "underdominance", dimension = 1, N = 5000, drive_params = list(s = 0.3, phat = 0.2, cubic_approximation = F), D = 1e-5, verbose = T, low = 0.0001, high = 0.1)
  • (N is set low here for the function to run fast)
• Then you can plot the spread of the underdominance allele at the critical width with:
  • simulate_underdominance_rde_1D(phat = 0.2, D = 1e-5, s = 0.3, release_width = width, N = 5000, cubic_approximation = F, plot =T)