diff --git a/app/DESCRIPTION b/app/DESCRIPTION
index 4d97921..fb6e4d4 100644
--- a/app/DESCRIPTION
+++ b/app/DESCRIPTION
@@ -18,17 +18,19 @@ Depends:
 Imports:
     bslib,
     config,
-    esquisse (>= 2.0.0),
+    esquisse (>= 2.1.0),
     generics,
     ggplot2,
     golem,
+    mirai (>= 2.5.0),
+    parallel,
     plotly,
     readr,
     shiny,
-    singlehit (>= 0.2.0)
+    singlehit (>= 0.3.0)
 Suggests:
     pkgload,
-    testthat (>= 3.0.0),
+    testthat (>= 3.1.7),
     withr
 Remotes:
     seanthimons/singlehit
diff --git a/app/R/app_server.R b/app/R/app_server.R
index fcb9b9c..9dfdb65 100644
--- a/app/R/app_server.R
+++ b/app/R/app_server.R
@@ -8,6 +8,9 @@ app_server <- function(input, output, session) {
   data <- mod_upload_server("upload")
   analysis <- mod_fit_server("fit", data)
   mod_compare_server("compare", analysis)
-  mod_figures_server("figures", analysis)
-  mod_explore_server("explore", analysis)
+  # Plots and Explore get the analysis with the bootstraps attached once a run
+  # finishes; Compare does not use them.
+  analysis_with_bands <- mod_bootstrap_server("bootstrap", analysis)
+  mod_figures_server("figures", analysis_with_bands)
+  mod_explore_server("explore", analysis_with_bands)
 }
diff --git a/app/R/app_ui.R b/app/R/app_ui.R
index 42eb8ec..fe6840f 100644
--- a/app/R/app_ui.R
+++ b/app/R/app_ui.R
@@ -7,6 +7,10 @@
 app_ui <- function(request) {
   tagList(
     golem_add_external_resources(),
+    # A pulse bar across the top while the server is busy, and spinners on
+    # outputs that are recalculating. R is single threaded, so a long fit or a
+    # run starting can make the app unresponsive; this says why.
+    shiny::useBusyIndicators(),
     bslib::page_navbar(
       title = "singlehit",
       bslib::nav_panel(
@@ -18,6 +22,7 @@ app_ui <- function(request) {
         bslib::layout_sidebar(sidebar = bslib::sidebar(mod_fit_ui("fit")), mod_fit_results_ui("fit"))
       ),
       bslib::nav_panel("Compare", mod_compare_ui("compare")),
+      bslib::nav_panel("Bootstrap", mod_bootstrap_ui("bootstrap")),
       bslib::nav_panel("Plots", mod_figures_ui("figures")),
       bslib::nav_panel("Explore", mod_explore_ui("explore"))
     )
diff --git a/app/R/fct_bootstrap.R b/app/R/fct_bootstrap.R
new file mode 100644
index 0000000..d9c9790
--- /dev/null
+++ b/app/R/fct_bootstrap.R
@@ -0,0 +1,172 @@
+#' Replicate limits for the bootstrap controls
+#'
+#' The minimum is 100 because a 99% percentile interval from fewer draws is
+#' meaningless. The exact model is far slower per replicate, so it has a smaller
+#' ceiling: the shared pool is first-in first-out, and a very long exact run
+#' would hold up every other session.
+#'
+#' The ceiling matches singlehit's own default of 10,000. Starting a run prepares
+#' every replicate on the app's single R process, which blocked it for 17 s at
+#' 10,000 replicates (two models) before singlehit drew the replicates in one
+#' vectorized call and stopped shipping the whole task list with every batch; it
+#' is now about a second. See dev/benchmark_bootstrap.R.
+#' @noRd
+bootstrap_limits <- list(min = 100L, max = 10000L, exact_max = 1000L)
+
+#' Gentle defaults for the bootstrap controls
+#' @noRd
+bootstrap_defaults <- list(times = 1000L, exact_times = 200L)
+
+check_count <- function(x, label, lower, upper) {
+  if (length(x) != 1L || !is.numeric(x) || is.na(x) || x != floor(x) || x < lower || x > upper) {
+    stop(
+      sprintf(
+        "%s must be a whole number between %s and %s.",
+        label,
+        format(lower, big.mark = ","),
+        format(upper, big.mark = ",")
+      ),
+      call. = FALSE
+    )
+  }
+  as.integer(x)
+}
+
+#' Work out what to run for each fitted model
+#'
+#' Seeds follow `analyze_dose_response()`: `seed + index - 1` over the fitted
+#' models, so a run can be reproduced from the exported script.
+#'
+#' @param analysis A `qdr_analysis`.
+#' @param times Replicates for the exponential and beta-Poisson models.
+#' @param exact_times Replicates for the exact model; used only if it was fitted.
+#' @param seed Optional whole number; `NULL` or `NA` for none.
+#' @return A named list (by model) of `list(model, times, seed)`.
+#' @noRd
+bootstrap_plan <- function(analysis, times, exact_times = NULL, seed = NULL) {
+  models <- names(analysis$fits)
+  times <- check_count(times, "Replicates", bootstrap_limits$min, bootstrap_limits$max)
+  if ("exact_beta_poisson" %in% models) {
+    if (is.null(exact_times)) {
+      exact_times <- bootstrap_defaults$exact_times
+    }
+    exact_times <- check_count(exact_times, "Exact-model replicates", bootstrap_limits$min, bootstrap_limits$exact_max)
+  }
+  if (is.null(seed) || length(seed) != 1L || is.na(seed)) {
+    seed <- NULL
+  } else {
+    seed <- check_count(seed, "The seed", 0L, .Machine$integer.max - length(models))
+  }
+  plan <- lapply(seq_along(models), function(i) {
+    list(
+      model = models[[i]],
+      times = if (models[[i]] == "exact_beta_poisson") exact_times else times,
+      seed = if (is.null(seed)) NULL else seed + i - 1L
+    )
+  })
+  names(plan) <- models
+  plan
+}
+
+#' Summarize the progress of several jobs
+#'
+#' @param progress Named list of `singlehit::bootstrap_progress()` results.
+#' @return A list with a per-model data frame and the overall fraction.
+#' @noRd
+progress_summary <- function(progress) {
+  done <- vapply(progress, function(p) p$done, integer(1))
+  total <- vapply(progress, function(p) p$total, integer(1))
+  list(
+    models = data.frame(
+      model = names(progress),
+      done = unname(done),
+      total = unname(total),
+      fraction = ifelse(total > 0L, unname(done) / unname(total), 1),
+      stringsAsFactors = FALSE
+    ),
+    fraction = if (sum(total) > 0L) sum(done) / sum(total) else 1,
+    finished = all(done == total)
+  )
+}
+
+#' Split bootstraps into those that can draw bands and those that cannot
+#'
+#' `prediction_curve()` stops when fewer than two replicates converged, which
+#' would break the Plots tab, so such a bootstrap is reported and not attached.
+#'
+#' @param boots Named list of `qdr_bootstrap` objects.
+#' @return A list with `attached` (named list) and `dropped` (model names).
+#' @noRd
+attachable_bootstraps <- function(boots) {
+  ok <- vapply(boots, function(boot) sum(boot$converged, na.rm = TRUE) >= 2L, logical(1))
+  list(attached = boots[ok], dropped = names(boots)[!ok])
+}
+
+#' How many replicates converged, per model
+#' @noRd
+convergence_table <- function(boots) {
+  data.frame(
+    Model = model_label_ui(names(boots)),
+    `Converged replicates` = vapply(
+      boots,
+      function(boot) sprintf("%d of %d", sum(boot$converged, na.rm = TRUE), nrow(boot)),
+      character(1)
+    ),
+    check.names = FALSE,
+    row.names = NULL
+  )
+}
+
+#' Percentile intervals for one bootstrap, one row per quantity
+#'
+#' Reshapes `singlehit::bootstrap_confint()` (long, one row per quantity and
+#' level) into one row per quantity with a column per confidence level.
+#' @noRd
+bootstrap_interval_table <- function(boot) {
+  ci <- as.data.frame(singlehit::bootstrap_confint(boot))
+  quantities <- unique(ci$term)
+  labels <- c(n50 = "N50", ed10 = "ED10", ed50 = "ED50")
+  shown <- ifelse(quantities %in% names(labels), labels[quantities], quantities)
+  out <- data.frame(
+    Model = model_label_ui(ci$model[[1L]]),
+    Quantity = unname(shown),
+    stringsAsFactors = FALSE
+  )
+  for (level in sort(unique(ci$level))) {
+    part <- ci[ci$level == level, , drop = FALSE]
+    idx <- match(quantities, part$term)
+    out[[paste0(format(100 * level, trim = TRUE), "% interval")]] <- paste0(
+      "(",
+      format_number(part$lower[idx]),
+      ", ",
+      format_number(part$upper[idx]),
+      ")"
+    )
+  }
+  out
+}
+
+#' Interval tables for all models, stacked
+#' @noRd
+bootstrap_tables <- function(boots) {
+  do.call(rbind, unname(lapply(boots, bootstrap_interval_table)))
+}
+
+#' Attach finished bootstraps to the analysis they were computed from
+#'
+#' Bands must never describe a different fit than the one on screen, so the
+#' bootstraps are attached only when `base` is identical to the analysis they
+#' were run against.
+#'
+#' @param base The current `qdr_analysis`, or `NULL`.
+#' @param boots Named list of finished `qdr_bootstrap` objects (possibly empty).
+#' @param bootstrapped_base The analysis the bootstraps were run against.
+#' @return `base` with the bootstraps attached, or `base` unchanged.
+#' @noRd
+attach_bands <- function(base, boots, bootstrapped_base) {
+  if (!is.null(base) && length(boots) > 0L && identical(bootstrapped_base, base)) {
+    singlehit::add_bootstrap(base, boots)
+  } else {
+    base
+  }
+}
diff --git a/app/R/fct_explore.R b/app/R/fct_explore.R
index 5be4156..f1d0feb 100644
--- a/app/R/fct_explore.R
+++ b/app/R/fct_explore.R
@@ -57,3 +57,31 @@ explore_datasets <- function(analysis) {
     model_residuals = model_residuals
   )
 }
+
+#' Input ids inside esquisse's module, relative to the Explore module
+#' @noRd
+esquisse_ids <- list(
+  aesthetics = "esquisse-geomaes-aes_1-dragvars",
+  geom = "esquisse-geomaes-geom_1",
+  x_transform = "esquisse-controls-axes-transX"
+)
+
+#' What the "Draw confidence bands" button fills in on the Fitted curves data
+#' @noRd
+bands_preset <- list(
+  aesthetics = list(xvar = "dose", yvar = "estimate", ymin = "lower", ymax = "upper", fill = "level", facet = "model"),
+  geom = "ribbon",
+  x_transform = "log10"
+)
+
+#' Does esquisse's reported drag-and-drop mapping contain the given variables?
+#'
+#' @param target The `target` element of the dragula input value: a list keyed
+#'   by slot name.
+#' @param wanted Named list: slot name to variable name.
+#' @noRd
+mapping_matches <- function(target, wanted) {
+  all(vapply(names(wanted), function(slot) {
+    identical(as.character(unlist(target[[slot]])), wanted[[slot]])
+  }, logical(1)))
+}
diff --git a/app/R/fct_figures.R b/app/R/fct_figures.R
index 3198a13..61e51d5 100644
--- a/app/R/fct_figures.R
+++ b/app/R/fct_figures.R
@@ -51,3 +51,16 @@ bounded_number <- function(x, lower, upper, default) {
   }
   min(max(x, lower), upper)
 }
+
+#' The view to switch to when confidence bands become available
+#'
+#' The overlay never draws bands, so someone looking at it right after a
+#' bootstrap would see no change and conclude the bands did not apply.
+#'
+#' @param type The currently selected view.
+#' @param has_bands Whether the analysis now carries bootstraps.
+#' @return The view to switch to, or `NULL` to leave the choice alone.
+#' @noRd
+view_for_bands <- function(type, has_bands) {
+  if (isTRUE(has_bands) && identical(type, "overlay")) "panels" else NULL
+}
diff --git a/app/R/fct_workers.R b/app/R/fct_workers.R
new file mode 100644
index 0000000..d499340
--- /dev/null
+++ b/app/R/fct_workers.R
@@ -0,0 +1,54 @@
+#' mirai compute profile shared by every bootstrap run in the app
+#'
+#' One named pool serves all sessions and all models. Without it each async
+#' bootstrap job would start its own pool sized to 75% of the cores, and only
+#' `collect_bootstrap()` shuts a job's own pool down.
+#' @noRd
+bootstrap_compute_profile <- "singlehitapp"
+
+#' Default number of bootstrap workers
+#'
+#' 75% of the physical cores, capped at 8 so a hosted deployment on a many-core
+#' machine does not claim all of its memory. Override with
+#' `options(singlehitapp.workers = n)`.
+#'
+#' @param cores Detected core count; injectable for tests.
+#' @return A whole number of at least 1.
+#' @noRd
+default_bootstrap_workers <- function(cores = parallel::detectCores(logical = FALSE)) {
+  chosen <- getOption("singlehitapp.workers")
+  if (!is.null(chosen)) {
+    return(max(1L, as.integer(chosen)))
+  }
+  if (length(cores) != 1L || is.na(cores) || cores < 1L) {
+    cores <- parallel::detectCores(logical = TRUE)
+  }
+  if (length(cores) != 1L || is.na(cores) || cores < 1L) {
+    return(1L)
+  }
+  max(1L, min(8L, as.integer(floor(0.75 * cores))))
+}
+
+#' Start the shared bootstrap workers if they are not running
+#'
+#' Safe to call before every run. Starting the pool ourselves means no job
+#' "owns" it, so one session finishing a job cannot shut the pool down under
+#' another session's job.
+#'
+#' @param workers Number of daemons to start.
+#' @noRd
+start_bootstrap_workers <- function(workers = default_bootstrap_workers()) {
+  if (!mirai::daemons_set(.compute = bootstrap_compute_profile)) {
+    mirai::daemons(workers, .compute = bootstrap_compute_profile)
+  }
+  invisible(TRUE)
+}
+
+#' Stop the shared bootstrap workers if they are running
+#' @noRd
+stop_bootstrap_workers <- function() {
+  if (mirai::daemons_set(.compute = bootstrap_compute_profile)) {
+    mirai::daemons(0L, .compute = bootstrap_compute_profile)
+  }
+  invisible(TRUE)
+}
diff --git a/app/R/mod_bootstrap.R b/app/R/mod_bootstrap.R
new file mode 100644
index 0000000..36669c3
--- /dev/null
+++ b/app/R/mod_bootstrap.R
@@ -0,0 +1,336 @@
+#' Bootstrap module UI
+#'
+#' @param id Module id.
+#' @noRd
+mod_bootstrap_ui <- function(id) {
+  ns <- NS(id)
+  # fillable = FALSE: the results are tall, so let the page scroll and give each
+  # card its natural height instead of squeezing the cards and plots to fit.
+  bslib::layout_sidebar(
+    fillable = FALSE,
+    sidebar = bslib::sidebar(
+      numericInput(
+        ns("times"),
+        "Replicates",
+        value = bootstrap_defaults$times,
+        min = bootstrap_limits$min,
+        max = bootstrap_limits$max,
+        step = 100
+      ),
+      uiOutput(ns("exact_input")),
+      radioButtons(
+        ns("resample"),
+        "Resampling",
+        choices = c("Observed response probabilities" = "observed", "Fitted model (parametric)" = "fitted")
+      ),
+      numericInput(ns("seed"), "Random seed (optional)", value = NA, min = 0, step = 1),
+      bslib::input_task_button(ns("run"), "Run bootstrap", label_busy = "Preparing replicates..."),
+      actionButton(ns("cancel"), "Cancel", class = "btn-outline-secondary"),
+      p(
+        class = "text-muted small mt-3",
+        "Runs in the background on a shared pool of workers, so you can keep using the other tabs.",
+        "When it finishes, confidence bands appear on the Plots and Explore tabs.",
+        "Starting a run briefly keeps the app busy, and the button shows when it is."
+      )
+    ),
+    uiOutput(ns("status")),
+    uiOutput(ns("results"))
+  )
+}
+
+#' One progress bar per model
+#'
+#' Per-model bars matter because the exact model's batches are far slower, so a
+#' single overall bar would stall at the end.
+#' @noRd
+progress_bars <- function(progress, elapsed) {
+  summary <- progress_summary(progress)
+  tagList(
+    p(class = "mb-2", sprintf("Running: %.0f s elapsed", elapsed)),
+    lapply(seq_len(nrow(summary$models)), function(i) {
+      row <- summary$models[i, ]
+      percent <- round(100 * row$fraction)
+      div(
+        class = "mb-2",
+        div(class = "small", model_label_ui(row$model)),
+        div(
+          class = "progress",
+          role = "progressbar",
+          `aria-valuenow` = percent,
+          `aria-valuemin` = 0,
+          `aria-valuemax` = 100,
+          div(class = "progress-bar", style = sprintf("width: %d%%;", percent), sprintf("%d%%", percent))
+        )
+      )
+    })
+  )
+}
+
+#' Bootstrap module server
+#'
+#' One job per fitted model runs on the shared mirai pool. Progress is polled
+#' without blocking; results are collected only once every job has finished. A
+#' run is cancelled when the analysis changes, when Cancel is pressed, and when
+#' the session ends, so no workers keep computing for a result nobody can see.
+#'
+#' @param id Module id.
+#' @param analysis Reactive returning a `qdr_analysis`, or `NULL`.
+#' @return A reactive returning the analysis with the completed bootstraps
+#'   attached (so Plots and Explore draw bands), or the base analysis otherwise.
+#' @noRd
+mod_bootstrap_server <- function(id, analysis) {
+  moduleServer(id, function(input, output, session) {
+    state <- reactiveValues(
+      status = "idle",
+      base = NULL,
+      progress = NULL,
+      boots = list(),
+      dropped = character(),
+      warnings = character(),
+      message = NULL,
+      started = NULL,
+      elapsed = NULL,
+      tick = 0L
+    )
+    jobs <- reactiveVal(list())
+
+    cancel_jobs <- function() {
+      for (job in isolate(jobs())) {
+        try(mirai::stop_mirai(job$map), silent = TRUE)
+      }
+      jobs(list())
+    }
+
+    reset <- function() {
+      state$status <- "idle"
+      state$base <- NULL
+      state$progress <- NULL
+      state$boots <- list()
+      state$dropped <- character()
+      state$warnings <- character()
+      state$message <- NULL
+    }
+
+    fail <- function(message) {
+      cancel_jobs()
+      state$status <- "error"
+      state$message <- message
+    }
+
+    # A different analysis means different fits: stop any run and drop results.
+    observeEvent(analysis(), {
+      cancel_jobs()
+      reset()
+    }, ignoreNULL = FALSE)
+
+    session$onSessionEnded(function() isolate(cancel_jobs()))
+
+    # Starts one job per model. The shared pool is (re)started first, so no job
+    # ever owns it and finishing one cannot shut it down under another session.
+    launch <- function(base, plan, resample) {
+      started <- list()
+      tryCatch(
+        {
+          start_bootstrap_workers()
+          for (item in plan) {
+            started[[item$model]] <- singlehit::bootstrap_dose_response_async(
+              base$fits[[item$model]],
+              times = item$times,
+              resample = resample,
+              seed = item$seed,
+              compute = bootstrap_compute_profile
+            )
+          }
+          started
+        },
+        error = function(e) {
+          for (job in started) {
+            try(mirai::stop_mirai(job$map), silent = TRUE)
+          }
+          e
+        }
+      )
+    }
+
+    observeEvent(input$run, {
+      base <- analysis()
+      if (is.null(base)) {
+        showNotification("Run the analysis on the Fit tab first.", type = "warning")
+      } else if (identical(state$status, "running")) {
+        showNotification("A bootstrap is already running.", type = "warning")
+      } else {
+        plan <- tryCatch(
+          bootstrap_plan(base, input$times, input$exact_times, input$seed),
+          error = function(e) e
+        )
+        if (inherits(plan, "error")) {
+          reset()
+          state$status <- "error"
+          state$message <- conditionMessage(plan)
+        } else {
+          # Preparing the replicates blocks the app for a few seconds on larger
+          # runs, so say so instead of looking frozen.
+          started <- withProgress(
+            message = "Preparing bootstrap replicates",
+            detail = "Large runs take a few seconds to start.",
+            value = 0.5,
+            launch(base, plan, input$resample)
+          )
+          if (inherits(started, "error")) {
+            reset()
+            state$status <- "error"
+            state$message <- paste("The bootstrap could not be started:", conditionMessage(started))
+          } else {
+            reset()
+            jobs(started)
+            state$base <- base
+            state$started <- Sys.time()
+            state$progress <- lapply(started, singlehit::bootstrap_progress)
+            state$status <- "running"
+          }
+        }
+      }
+    })
+
+    observeEvent(input$cancel, {
+      if (identical(state$status, "running")) {
+        cancel_jobs()
+        state$status <- "cancelled"
+      }
+    })
+
+    finish <- function() {
+      collected <- with_captured_warnings(lapply(jobs(), singlehit::collect_bootstrap))
+      jobs(list())
+      state$warnings <- collected$warnings
+      if (!is.null(collected$error)) {
+        state$status <- "error"
+        state$message <- paste("The bootstrap run failed:", collected$error)
+      } else {
+        split <- attachable_bootstraps(collected$value)
+        state$boots <- split$attached
+        state$dropped <- split$dropped
+        state$elapsed <- as.numeric(difftime(Sys.time(), state$started, units = "secs"))
+        state$status <- "done"
+      }
+    }
+
+    # Polling reads progress only; it never waits on a job. A cancelled run has
+    # status "cancelled", so its stopped jobs (which report as finished) are
+    # never collected.
+    observe({
+      req(identical(state$status, "running"))
+      invalidateLater(500, session)
+      isolate({
+        progress <- tryCatch(lapply(jobs(), singlehit::bootstrap_progress), error = function(e) e)
+        if (inherits(progress, "error")) {
+          fail(conditionMessage(progress))
+        } else {
+          state$progress <- progress
+          # Re-render the status every poll so the elapsed time keeps moving even
+          # when no batch has finished since the last one.
+          state$tick <- state$tick + 1L
+          if (all(vapply(progress, function(p) p$finished, logical(1)))) {
+            finish()
+          }
+        }
+      })
+    })
+
+    output$exact_input <- renderUI({
+      req(analysis())
+      if ("exact_beta_poisson" %in% names(analysis()$fits)) {
+        tagList(
+          numericInput(
+            session$ns("exact_times"),
+            "Exact-model replicates",
+            value = bootstrap_defaults$exact_times,
+            min = bootstrap_limits$min,
+            max = bootstrap_limits$exact_max,
+            step = 50
+          ),
+          p(class = "text-muted small", "The exact model is much slower per replicate; start low.")
+        )
+      }
+    })
+
+    output$status <- renderUI({
+      if (is.null(analysis())) {
+        return(p(class = "text-muted", "Run the analysis on the Fit tab first."))
+      }
+      switch(
+        state$status,
+        idle = p(class = "text-muted", "Choose the settings and click Run bootstrap."),
+        running = {
+          state$tick
+          progress_bars(state$progress, as.numeric(difftime(Sys.time(), state$started, units = "secs")))
+        },
+        done = tagList(
+          div(
+            class = "alert alert-success",
+            sprintf(
+              "Bootstrap complete in %.1f seconds. Confidence bands now appear on the Plots and Explore tabs.",
+              state$elapsed
+            )
+          ),
+          if (length(state$dropped) > 0L) {
+            div(
+              class = "alert alert-warning",
+              sprintf(
+                "Too few replicates converged to draw bands for: %s.",
+                paste(model_label_ui(state$dropped), collapse = ", ")
+              )
+            )
+          },
+          if (length(state$warnings) > 0L) {
+            div(class = "alert alert-warning", tags$ul(lapply(state$warnings, tags$li)))
+          }
+        ),
+        cancelled = div(class = "alert alert-warning", "Bootstrap cancelled. No confidence bands were added."),
+        error = div(class = "alert alert-danger", state$message)
+      )
+    })
+
+    output$results <- renderUI({
+      req(identical(state$status, "done"), length(state$boots) > 0L)
+      ns <- session$ns
+      tagList(
+        bslib::card(
+          fill = FALSE,
+          bslib::card_header("Percentile intervals"),
+          tableOutput(ns("intervals")),
+          tableOutput(ns("convergence"))
+        ),
+        bslib::card(
+          fill = FALSE,
+          bslib::card_header("Parameter uncertainty"),
+          lapply(names(state$boots), function(model) {
+            plotOutput(ns(paste0("plot_", model)), height = "320px", fill = FALSE)
+          })
+        )
+      )
+    })
+    output$intervals <- renderTable(bootstrap_tables(state$boots))
+    output$convergence <- renderTable(convergence_table(state$boots))
+
+    observeEvent(state$boots, {
+      for (model in names(state$boots)) {
+        local({
+          chosen <- model
+          output[[paste0("plot_", chosen)]] <- renderPlot(ggplot2::autoplot(state$boots[[chosen]]), res = 96)
+        })
+      }
+    })
+
+    # Held in a reactiveVal, which ignores a set to an identical value, so Plots
+    # and Explore only re-render when the analysis they show really changes, not
+    # on every status change (Run, Cancel, errors). The exact model's curves are
+    # slow to redraw.
+    shown <- reactiveVal(NULL)
+    observe({
+      boots <- if (identical(state$status, "done")) state$boots else list()
+      shown(attach_bands(analysis(), boots, state$base))
+    })
+    reactive(shown())
+  })
+}
diff --git a/app/R/mod_explore.R b/app/R/mod_explore.R
index 54f6514..7de397d 100644
--- a/app/R/mod_explore.R
+++ b/app/R/mod_explore.R
@@ -25,14 +25,27 @@ mod_explore_ui <- function(id) {
       ns = ns,
       div(
         class = "mb-3",
-        selectInput(ns("dataset"), "Data to explore", choices = explore_dataset_choices, width = "320px"),
+        div(
+          class = "d-flex flex-wrap align-items-end gap-3 mb-2",
+          selectInput(ns("dataset"), "Data to explore", choices = explore_dataset_choices, width = "320px"),
+          conditionalPanel(
+            "output.has_bands == 'yes'",
+            ns = ns,
+            actionButton(
+              ns("bands_preset"),
+              "Draw confidence bands",
+              icon = icon("chart-area"),
+              class = "btn-outline-primary mb-3"
+            )
+          )
+        ),
         p(
           class = "text-muted small mb-0",
           paste(
             "Drag variables onto the slots to build a plot. For a log dose axis, open Axes and set the X-axis",
-            "transform to log10. Confidence bands (Fitted curves, once a bootstrap has been run): map lower to ymin",
-            "and upper to ymax, choose the Ribbon geometry, and split by model and level, for example facet by model",
-            "and fill by level."
+            "transform to log10. Confidence bands (Fitted curves, once a bootstrap has been run): use the Draw",
+            "confidence bands button, or map lower to ymin and upper to ymax, choose the Ribbon geometry, and split",
+            "by model and level, for example facet by model and fill by level."
           )
         )
       ),
@@ -73,6 +86,13 @@ mod_explore_server <- function(id, analysis) {
     outputOptions(output, "has_analysis", suspendWhenHidden = FALSE)
 
     all_datasets <- reactive(explore_datasets(req(analysis())))
+    # Judged on the data, not on whether a bootstrap is attached: a bootstrap
+    # with too few converged replicates gives bare curves with no bands to draw.
+    output$has_bands <- renderText({
+      has_bands <- !is.null(analysis()) && any(!is.na(all_datasets()$curves$lower))
+      if (has_bands) "yes" else "no"
+    })
+    outputOptions(output, "has_bands", suspendWhenHidden = FALSE)
     dataset <- reactive({
       req(input$dataset %in% names(all_datasets()))
       all_datasets()[[input$dataset]]
@@ -87,6 +107,56 @@ mod_explore_server <- function(id, analysis) {
       import_from = NULL
     )
 
+    # Pre-fill esquisse's slots. Its drag-and-drop, geometry and axis pickers are
+    # ordinary shiny inputs, so they can be driven with their public update
+    # functions, but each step has to wait for esquisse to react to the one before
+    # it: choosing the dataset rebuilds the variable list, and every change to the
+    # mapping makes esquisse re-pick the geometry. `preset` records which step is
+    # being waited for. Messages go out from onFlushed() so they follow whatever
+    # esquisse sent in the same flush.
+    preset <- reactiveVal(NULL)
+    reported_mapping <- reactive(input[[esquisse_ids$aesthetics]]$target)
+
+    set_geometry_and_axis <- function() {
+      preset(NULL)
+      session$onFlushed(function() {
+        esquisse::updateDropInput(session, esquisse_ids$geom, selected = bands_preset$geom)
+        updateSelectizeInput(session, esquisse_ids$x_transform, selected = bands_preset$x_transform)
+      }, once = TRUE)
+    }
+    set_aesthetics <- function() {
+      if (mapping_matches(reported_mapping(), bands_preset$aesthetics)) {
+        return(set_geometry_and_axis())
+      }
+      preset("aesthetics")
+      session$onFlushed(function() {
+        esquisse::updateDragulaInput(session, esquisse_ids$aesthetics, selected = bands_preset$aesthetics)
+      }, once = TRUE)
+    }
+
+    observeEvent(input$bands_preset, {
+      if (identical(input$dataset, "curves")) {
+        set_aesthetics()
+      } else {
+        preset("data")
+        updateSelectInput(session, "dataset", selected = "curves")
+      }
+    })
+    # Drop a pending preset whose conditions no longer hold, so it cannot fire
+    # later when the user happens to build the same mapping by hand.
+    observeEvent(analysis(), preset(NULL), ignoreInit = TRUE)
+    observeEvent(input$dataset, if (!identical(input$dataset, "curves")) preset(NULL))
+    observeEvent(dataset(), {
+      if (identical(preset(), "data") && identical(input$dataset, "curves")) set_aesthetics()
+    })
+    observeEvent(reported_mapping(), {
+      # The dataset change also reports its own (stale) mapping first; only act
+      # once esquisse reports the mapping that was sent.
+      if (identical(preset(), "aesthetics") && mapping_matches(reported_mapping(), bands_preset$aesthetics)) {
+        set_geometry_and_axis()
+      }
+    })
+
     dataset
   })
 }
diff --git a/app/R/mod_figures.R b/app/R/mod_figures.R
index 16f6bff..0754c01 100644
--- a/app/R/mod_figures.R
+++ b/app/R/mod_figures.R
@@ -7,6 +7,7 @@ mod_figures_ui <- function(id) {
   bslib::layout_sidebar(
     sidebar = bslib::sidebar(
       radioButtons(ns("type"), "View", choices = figure_types),
+      uiOutput(ns("band_hint")),
       conditionalPanel(
         "input.type == 'single'",
         ns = ns,
@@ -64,6 +65,43 @@ mod_figures_server <- function(id, analysis) {
     })
     output$plot <- renderPlot(figure(), res = 96)
 
+    # When bands first become available, leave the overlay (which never draws
+    # them) for the panel view, so the new bands are visible straight away.
+    # Held in a reactiveVal, which ignores a set to an identical value, so the
+    # handler below runs only when bands really appear or disappear and does not
+    # depend on how often the upstream analysis is recomputed.
+    has_bands <- reactiveVal(FALSE)
+    observe(has_bands(length(analysis()$bootstraps) > 0L))
+    observeEvent(
+      has_bands(),
+      {
+        target <- view_for_bands(input$type, has_bands())
+        if (!is.null(target)) {
+          updateRadioButtons(session, "type", selected = target)
+          showNotification(
+            "Confidence bands are drawn in the panel and single-model views, so this view was switched to one panel per model.",
+            type = "message",
+            duration = 10
+          )
+        }
+      },
+      ignoreInit = TRUE
+    )
+
+    output$band_hint <- renderUI({
+      req(analysis(), input$type)
+      hint <- function(text) div(class = "alert alert-info small py-2 mb-2", text)
+      if (identical(input$type, "overlay")) {
+        hint("The overlay never shows confidence bands. Choose One panel per model or Single model to see them.")
+      } else if (identical(input$type, "single") && !is.null(input$model)) {
+        if (!input$model %in% names(analysis()$bootstraps)) {
+          hint("No confidence bands for this model yet: run the bootstrap.")
+        }
+      } else if (length(analysis()$bootstraps) == 0L) {
+        hint("Run the bootstrap to add confidence bands to this view.")
+      }
+    })
+
     output$download_button <- renderUI({
       req(analysis())
       downloadButton(session$ns("download"), "Download PNG")
diff --git a/app/R/mod_fit.R b/app/R/mod_fit.R
index bd3b0b7..2105109 100644
--- a/app/R/mod_fit.R
+++ b/app/R/mod_fit.R
@@ -18,7 +18,9 @@ mod_fit_ui <- function(id) {
       selected = c("exponential", "beta_poisson")
     ),
     numericInput(ns("alpha"), "Significance level (alpha)", value = 0.05, min = 0.001, max = 0.5, step = 0.01),
-    actionButton(ns("run"), "Run analysis", class = "btn-primary")
+    # A task button turns busy the moment it is clicked, in the browser, so it
+    # still shows something while the server is too busy to answer.
+    bslib::input_task_button(ns("run"), "Run analysis", label_busy = "Fitting models...")
   )
 }
 
diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R
index 42f03e0..98fcb6c 100644
--- a/app/R/mod_upload.R
+++ b/app/R/mod_upload.R
@@ -10,7 +10,8 @@ mod_upload_ui <- function(id) {
   auto <- c("Auto-detect" = "")
   tagList(
     fileInput(ns("file"), "Dose-response file", accept = c(".csv", ".tsv", ".txt")),
-    actionButton(ns("example"), "Use example data"),
+    actionButton(ns("example"), "Use Ward rotavirus example data"),
+    p(class = "text-muted small mt-1", "Ward et al. (1986) human rotavirus challenge study."),
     bslib::accordion(
       id = ns("mapping_panel"),
       open = FALSE,
@@ -71,7 +72,7 @@ mod_upload_server <- function(id) {
     observeEvent(input$file, load_source(input$file$datapath))
     observeEvent(
       input$example,
-      load_source(system.file("extdata", "dose-response-example.csv", package = "singlehit"))
+      load_source(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit"))
     )
 
     # Same reader settings as singlehit::read_dose_response(), so the offered
diff --git a/app/R/run_app.R b/app/R/run_app.R
index 687fb88..3ae8dc9 100644
--- a/app/R/run_app.R
+++ b/app/R/run_app.R
@@ -8,21 +8,30 @@
 #' @importFrom shiny shinyApp
 #' @importFrom golem with_golem_options
 run_app <- function(
-	onStart = NULL,
-	options = list(),
-	enableBookmarking = NULL,
-	uiPattern = "/",
-	...
+  onStart = NULL,
+  options = list(),
+  enableBookmarking = NULL,
+  uiPattern = "/",
+  ...
 ) {
-	with_golem_options(
-		app = shinyApp(
-			ui = app_ui,
-			server = app_server,
-			onStart = onStart,
-			options = options,
-			enableBookmarking = enableBookmarking,
-			uiPattern = uiPattern
-		),
-		golem_opts = list(...)
-	)
+  user_on_start <- onStart
+  with_golem_options(
+    app = shinyApp(
+      ui = app_ui,
+      server = app_server,
+      # The bootstrap workers start lazily on the first run, so a deployment that
+      # never bootstraps never claims them; this only makes sure they stop with
+      # the app.
+      onStart = function() {
+        shiny::onStop(stop_bootstrap_workers)
+        if (is.function(user_on_start)) {
+          user_on_start()
+        }
+      },
+      options = options,
+      enableBookmarking = enableBookmarking,
+      uiPattern = uiPattern
+    ),
+    golem_opts = list(...)
+  )
 }
diff --git a/app/README.md b/app/README.md
index 0492b73..5c94388 100644
--- a/app/README.md
+++ b/app/README.md
@@ -8,7 +8,7 @@ depends on the app.
 
 ## What it does today
 
-- **Data:** upload a delimited file or use the bundled example, with optional
+- **Data:** upload a delimited file or use the bundled Ward rotavirus example, with optional
   column mapping when names are not recognized automatically.
 - **Fit:** choose at least two models (exponential, beta-Poisson, exact
   beta-Poisson), see the dose-trend test, parameter estimates, and goodness of
@@ -17,23 +17,38 @@ depends on the app.
   deviance comparison), the assessment table, and the experimental consensus
   vote, labelled as such.
 
+- **Bootstrap:** percentile confidence intervals, run in the background on a
+  shared pool of workers (a progress bar per model, a Cancel button, and a
+  default of 1,000 replicates; the exact model has its own smaller budget, as
+  it is much slower). When a run finishes, confidence bands appear on the
+  Plots and Explore tabs.
 - **Plots:** the model overlay, one panel per model, or a single model, with PNG
-  download. Confidence bands will appear in the panel and single-model views
-  once the planned bootstrap support lands; the overlay never draws bands.
+  download. Confidence bands appear in the panel and single-model views once a
+  bootstrap has been run; the overlay never draws bands.
 
 - **Explore:** build your own plots interactively with
   [esquisse](https://dreamrs.github.io/esquisse/) by dragging variables onto
   aesthetics. Choose the observed dose groups, the fitted curves, or the
-  residuals. Turn on "Use {plotly} to render plot" under Options for hover and zoom.
+  residuals. Once a bootstrap has been run, "Draw confidence bands" fills in a
+  ribbon plot of the fitted curves on a log dose axis (it drives esquisse's
+  inputs, so it was written against esquisse 2.1.0; if a later esquisse renames
+  them it switches to the fitted curves but leaves the plot unfilled). Turn on
+  "Use {plotly} to render plot" under Options for hover and zoom.
   The interface font is overridden in `inst/app/www/custom.css`.
 
-Bootstrap uncertainty, pooling, and export are planned.
+Pooling and export are planned.
 
 ## License
 
 The app is licensed under GPL-3 (see `LICENSE.md`) because it imports esquisse,
 which is GPL-3. The `singlehit` package itself stays MIT.
 
+## Requires singlehit 0.3.0 or later
+
+The Bootstrap and Plots tabs use `singlehit` functions (`add_bootstrap()`,
+`bootstrap_progress()`) that were added in 0.3.0, so install `singlehit` 0.3.0
+or later (for example from GitHub with `pak::pak("seanthimons/singlehit")`).
+
 ## Running during development
 
 From the repository root:
diff --git a/app/dev/benchmark_bootstrap.R b/app/dev/benchmark_bootstrap.R
new file mode 100644
index 0000000..9df2a2a
--- /dev/null
+++ b/app/dev/benchmark_bootstrap.R
@@ -0,0 +1,75 @@
+# Benchmark how long a bootstrap run blocks the app and how long it takes in total.
+#
+# "Blocked" is the time the calling R process is busy starting the run, which is
+# when a Shiny app cannot respond to anyone. "Total" is until every batch is
+# collected. Run from the app/ directory:
+#
+#   source("dev/benchmark_bootstrap.R")
+#   results <- benchmark_bootstrap()                      # cheap models
+#   results <- benchmark_bootstrap(exact_times = c(200))  # add the exact model
+#
+# It loads the in-repo singlehit, so it measures the code in your checkout.
+pkgload::load_all("..", export_all = FALSE, helpers = FALSE, attach_testthat = FALSE, quiet = TRUE)
+
+benchmark_bootstrap <- function(
+  times = c(1000L, 2500L, 5000L, 10000L),
+  reps = 3L,
+  workers = 8L,
+  exact_times = integer(),
+  models = c("exponential", "beta_poisson")
+) {
+  data <- singlehit::read_dose_response(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit"))
+  all_models <- union(models, if (length(exact_times) > 0L) "exact_beta_poisson")
+  analysis <- singlehit::analyze_dose_response(
+    data,
+    models = all_models,
+    bootstrap_times = 0L,
+    exact_bootstrap_times = 0L,
+    check_trend = FALSE
+  )
+  profile <- "benchmark"
+  mirai::daemons(workers, .compute = profile)
+  on.exit(mirai::daemons(0L, .compute = profile), add = TRUE)
+
+  run_once <- function(fits, n) {
+    started <- Sys.time()
+    jobs <- lapply(fits, function(fit) {
+      singlehit::bootstrap_dose_response_async(fit, times = n, seed = 1, compute = profile)
+    })
+    blocked <- as.numeric(difftime(Sys.time(), started, units = "secs"))
+    while (!all(vapply(jobs, function(job) singlehit::bootstrap_progress(job)$finished, logical(1)))) {
+      Sys.sleep(0.2)
+    }
+    lapply(jobs, singlehit::collect_bootstrap)
+    c(blocked = blocked, total = as.numeric(difftime(Sys.time(), started, units = "secs")))
+  }
+
+  cases <- rbind(
+    data.frame(models = paste(models, collapse = " + "), replicates = times, exact = FALSE),
+    if (length(exact_times) > 0L) data.frame(models = "exact_beta_poisson", replicates = exact_times, exact = TRUE)
+  )
+  rows <- lapply(seq_len(nrow(cases)), function(i) {
+    fits <- if (cases$exact[[i]]) analysis$fits["exact_beta_poisson"] else analysis$fits[models]
+    n_reps <- if (cases$exact[[i]]) 1L else reps
+    runs <- t(vapply(seq_len(n_reps), function(r) run_once(fits, cases$replicates[[i]]), numeric(2)))
+    cat(sprintf(
+      "%-28s n=%5d: blocked %5.2fs (max %5.2fs) | total %6.1fs (%d run%s)\n",
+      cases$models[[i]],
+      cases$replicates[[i]],
+      mean(runs[, "blocked"]),
+      max(runs[, "blocked"]),
+      mean(runs[, "total"]),
+      n_reps,
+      if (n_reps == 1L) "" else "s"
+    ))
+    data.frame(
+      cases[i, c("models", "replicates")],
+      workers = workers,
+      blocked_mean = mean(runs[, "blocked"]),
+      blocked_max = max(runs[, "blocked"]),
+      total_mean = mean(runs[, "total"]),
+      runs = n_reps
+    )
+  })
+  invisible(do.call(rbind, rows))
+}
diff --git a/app/tests/testthat/helper-data.R b/app/tests/testthat/helper-data.R
index a162129..9a3ca66 100644
--- a/app/tests/testthat/helper-data.R
+++ b/app/tests/testthat/helper-data.R
@@ -11,3 +11,14 @@ example_analysis <- function(models = c("exponential", "beta_poisson")) {
     check_trend = FALSE
   )
 }
+
+ward_data <- function() {
+  singlehit::read_dose_response(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit"))
+}
+
+# Small, fast bootstrap (sequential, in-process) attached to the tiny example.
+example_analysis_with_bands <- function(models = c("exponential", "beta_poisson")) {
+  analysis <- example_analysis(models)
+  boots <- singlehit::bootstrap_dose_response_models(analysis$fits, times = 20, seed = 1, backend = "sequential")
+  singlehit::add_bootstrap(analysis, boots)
+}
diff --git a/app/tests/testthat/test-fct_bootstrap.R b/app/tests/testthat/test-fct_bootstrap.R
new file mode 100644
index 0000000..249d99d
--- /dev/null
+++ b/app/tests/testthat/test-fct_bootstrap.R
@@ -0,0 +1,95 @@
+test_that("the plan gives each model the right budget and seed offset", {
+  analysis <- example_analysis(c("exponential", "beta_poisson"))
+  plan <- bootstrap_plan(analysis, times = 500, exact_times = 300, seed = 10)
+  expect_named(plan, c("exponential", "beta_poisson"))
+  expect_equal(vapply(plan, function(p) p$times, integer(1)), c(exponential = 500L, beta_poisson = 500L))
+  # Same offsets as analyze_dose_response(): seed + index - 1.
+  expect_equal(vapply(plan, function(p) p$seed, integer(1)), c(exponential = 10L, beta_poisson = 11L))
+})
+
+test_that("the exact model uses its own budget, and is ignored when not fitted", {
+  with_exact <- example_analysis(c("exponential", "exact_beta_poisson"))
+  plan <- bootstrap_plan(with_exact, times = 1000, exact_times = 200, seed = NULL)
+  expect_equal(plan$exponential$times, 1000L)
+  expect_equal(plan$exact_beta_poisson$times, 200L)
+  expect_null(plan$exact_beta_poisson$seed)
+
+  # A missing exact input falls back to the default budget.
+  expect_equal(bootstrap_plan(with_exact, 1000, NULL)$exact_beta_poisson$times, bootstrap_defaults$exact_times)
+
+  # Not fitted: a bad exact value is irrelevant and not validated.
+  expect_no_error(bootstrap_plan(example_analysis(), times = 1000, exact_times = 99999))
+})
+
+test_that("a blank seed means no seed", {
+  analysis <- example_analysis()
+  expect_null(bootstrap_plan(analysis, 1000, NULL, NA)$exponential$seed)
+  expect_null(bootstrap_plan(analysis, 1000, NULL, NULL)$exponential$seed)
+})
+
+test_that("invalid settings are refused with a clear message", {
+  analysis <- example_analysis()
+  expect_error(bootstrap_plan(analysis, 50), "between 100 and 10,000")
+  expect_error(bootstrap_plan(analysis, 20000), "between 100 and 10,000")
+  expect_no_error(bootstrap_plan(analysis, 10000))
+  expect_error(bootstrap_plan(analysis, 1000.5), "whole number")
+  expect_error(bootstrap_plan(analysis, NA_real_), "whole number")
+  expect_error(bootstrap_plan(analysis, c(100, 200)), "whole number")
+  expect_error(bootstrap_plan(analysis, 1000, seed = -1), "seed")
+  expect_error(bootstrap_plan(analysis, 1000, seed = 1.5), "seed")
+  expect_error(bootstrap_plan(analysis, 1000, seed = .Machine$integer.max), "seed")
+
+  with_exact <- example_analysis(c("exponential", "exact_beta_poisson"))
+  expect_error(bootstrap_plan(with_exact, 1000, exact_times = 5000), "between 100 and 1,000")
+  expect_error(bootstrap_plan(with_exact, 1000, exact_times = 10), "Exact-model replicates")
+})
+
+test_that("progress is summarized per model and overall", {
+  progress <- list(
+    exponential = list(done = 50L, total = 50L, finished = TRUE),
+    exact_beta_poisson = list(done = 1L, total = 10L, finished = FALSE)
+  )
+  summary <- progress_summary(progress)
+  expect_equal(summary$models$model, c("exponential", "exact_beta_poisson"))
+  expect_equal(summary$models$fraction, c(1, 0.1))
+  expect_equal(summary$fraction, 51 / 60)
+  expect_false(summary$finished)
+  expect_true(progress_summary(list(a = list(done = 3L, total = 3L, finished = TRUE)))$finished)
+})
+
+test_that("bootstraps with fewer than two converged replicates are not attached", {
+  make <- function(converged) structure(list(converged = converged), class = "data.frame", row.names = seq_along(converged))
+  boots <- list(exponential = make(c(TRUE, TRUE, FALSE)), exact_beta_poisson = make(c(TRUE, FALSE, FALSE)))
+  split <- attachable_bootstraps(boots)
+  expect_equal(names(split$attached), "exponential")
+  expect_equal(split$dropped, "exact_beta_poisson")
+})
+
+test_that("the interval and convergence tables describe a real bootstrap", {
+  analysis <- example_analysis()
+  boots <- singlehit::bootstrap_dose_response_models(analysis$fits, times = 100, seed = 1, backend = "sequential")
+  tables <- bootstrap_tables(boots)
+  expect_true(all(c("Model", "Quantity", "95% interval", "99% interval") %in% names(tables)))
+  expect_setequal(tables$Model, c("Exponential", "Beta-Poisson"))
+  expect_true("N50" %in% tables$Quantity)
+  expect_true(all(grepl("^\\(.+, .+\\)$", tables[["95% interval"]])))
+
+  convergence <- convergence_table(boots)
+  expect_equal(convergence$Model, c("Exponential", "Beta-Poisson"))
+  expect_match(convergence[["Converged replicates"]], "^\\d+ of 100$")
+})
+
+test_that("bands are attached only to the analysis they were computed from", {
+  analysis <- example_analysis()
+  boots <- singlehit::bootstrap_dose_response_models(analysis$fits, times = 100, seed = 1, backend = "sequential")
+  other <- example_analysis(c("exponential", "exact_beta_poisson"))
+
+  attached <- attach_bands(analysis, boots, analysis)
+  expect_setequal(names(attached$bootstraps), c("exponential", "beta_poisson"))
+
+  # A different analysis (or none) must never receive another fit's bands.
+  expect_length(attach_bands(other, boots, analysis)$bootstraps, 0L)
+  expect_null(attach_bands(NULL, boots, analysis))
+  expect_identical(attach_bands(analysis, list(), analysis), analysis)
+  expect_identical(attach_bands(analysis, boots, NULL), analysis)
+})
diff --git a/app/tests/testthat/test-fct_explore.R b/app/tests/testthat/test-fct_explore.R
index e94a3d7..9c6d21e 100644
--- a/app/tests/testthat/test-fct_explore.R
+++ b/app/tests/testthat/test-fct_explore.R
@@ -50,3 +50,17 @@ test_that("residuals pair each model with fitted values and deviance residuals",
   expect_true(all(c("dose", "response", ".fitted", ".resid") %in% names(residuals)))
   expect_equal(nrow(residuals), 6L)
 })
+
+test_that("the bands preset names columns the Fitted curves data actually has", {
+  curves <- explore_datasets(analysis_with_bands())$curves
+  expect_true(all(unlist(bands_preset$aesthetics) %in% names(curves)))
+  expect_true(bands_preset$geom %in% c("ribbon", "area"))
+})
+
+test_that("mapping_matches() compares esquisse's slot mapping with the wanted variables", {
+  wanted <- list(xvar = "dose", yvar = "estimate")
+  expect_true(mapping_matches(list(xvar = "dose", yvar = "estimate", fill = character()), wanted))
+  expect_false(mapping_matches(list(xvar = "dose", yvar = "other"), wanted))
+  expect_false(mapping_matches(list(xvar = "dose"), wanted))
+  expect_false(mapping_matches(NULL, wanted))
+})
diff --git a/app/tests/testthat/test-fct_figures.R b/app/tests/testthat/test-fct_figures.R
index 1acfdba..7b5f826 100644
--- a/app/tests/testthat/test-fct_figures.R
+++ b/app/tests/testthat/test-fct_figures.R
@@ -48,3 +48,11 @@ test_that("bounded_number clamps and falls back", {
   expect_equal(bounded_number(NA_real_, 2, 20, 7), 7)
   expect_equal(bounded_number(NULL, 2, 20, 7), 7)
 })
+
+test_that("new bands move someone off the overlay, which never draws them, and nothing else", {
+  expect_equal(view_for_bands("overlay", TRUE), "panels")
+  expect_null(view_for_bands("overlay", FALSE))
+  expect_null(view_for_bands("panels", TRUE))
+  expect_null(view_for_bands("single", TRUE))
+  expect_null(view_for_bands(NULL, TRUE))
+})
diff --git a/app/tests/testthat/test-fct_workers.R b/app/tests/testthat/test-fct_workers.R
new file mode 100644
index 0000000..f6d572b
--- /dev/null
+++ b/app/tests/testthat/test-fct_workers.R
@@ -0,0 +1,35 @@
+test_that("the default worker count is 75% of the cores, capped at 8, and at least 1", {
+  withr::local_options(singlehitapp.workers = NULL)
+  expect_equal(default_bootstrap_workers(cores = 16L), 8L)
+  expect_equal(default_bootstrap_workers(cores = 8L), 6L)
+  expect_equal(default_bootstrap_workers(cores = 4L), 3L)
+  expect_equal(default_bootstrap_workers(cores = 1L), 1L)
+})
+
+test_that("an unknown core count falls back to a usable number", {
+  withr::local_options(singlehitapp.workers = NULL)
+  workers <- default_bootstrap_workers(cores = NA_integer_)
+  expect_true(is.integer(workers) && workers >= 1L && workers <= 8L)
+})
+
+test_that("the singlehitapp.workers option overrides the default", {
+  withr::local_options(singlehitapp.workers = 2)
+  expect_equal(default_bootstrap_workers(cores = 16L), 2L)
+  withr::local_options(singlehitapp.workers = 0)
+  expect_equal(default_bootstrap_workers(cores = 16L), 1L)
+})
+
+test_that("starting and stopping the shared workers is idempotent", {
+  skip_if_not_installed("mirai")
+  withr::defer(stop_bootstrap_workers())
+  expect_false(mirai::daemons_set(.compute = bootstrap_compute_profile))
+
+  start_bootstrap_workers(2)
+  expect_true(mirai::daemons_set(.compute = bootstrap_compute_profile))
+  expect_no_error(start_bootstrap_workers(2))
+  expect_true(mirai::daemons_set(.compute = bootstrap_compute_profile))
+
+  stop_bootstrap_workers()
+  expect_false(mirai::daemons_set(.compute = bootstrap_compute_profile))
+  expect_no_error(stop_bootstrap_workers())
+})
diff --git a/app/tests/testthat/test-mod_bootstrap.R b/app/tests/testthat/test-mod_bootstrap.R
new file mode 100644
index 0000000..ef12df2
--- /dev/null
+++ b/app/tests/testthat/test-mod_bootstrap.R
@@ -0,0 +1,142 @@
+# These tests run real jobs on a tiny shared worker pool. 100 replicates is the
+# smallest allowed run: 5 batches per model.
+wait_for <- function(session, condition, tries = 150) {
+  for (i in seq_len(tries)) {
+    session$elapse(600)
+    if (condition()) {
+      return(TRUE)
+    }
+    Sys.sleep(0.2)
+  }
+  FALSE
+}
+
+local_workers <- function(workers = 2, env = parent.frame()) {
+  skip_if_not_installed("mirai")
+  withr::local_options(singlehitapp.workers = workers, .local_envir = env)
+  withr::defer(stop_bootstrap_workers(), envir = env)
+}
+
+test_that("a run completes and the returned analysis carries the bootstraps", {
+  local_workers()
+  analysis <- shiny::reactiveVal(example_analysis())
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = analysis), {
+    expect_length(session$returned()$bootstraps, 0L)
+
+    session$setInputs(times = 100, resample = "observed", seed = 1, run = 1)
+    expect_equal(state$status, "running")
+    expect_true(wait_for(session, function() state$status != "running"))
+    expect_equal(state$status, "done")
+
+    returned <- session$returned()
+    expect_setequal(names(returned$bootstraps), c("exponential", "beta_poisson"))
+    expect_equal(nrow(returned$bootstraps$exponential), 100L)
+    expect_equal(state$progress$exponential$done, state$progress$exponential$total)
+    expect_s3_class(singlehit::bootstrap_confint(returned$bootstraps$beta_poisson), "data.frame")
+  })
+})
+
+test_that("the same seed gives the same replicates", {
+  local_workers()
+  run_once <- function() {
+    deviance <- NULL
+    shiny::testServer(
+      mod_bootstrap_server,
+      args = list(analysis = shiny::reactiveVal(example_analysis())),
+      {
+        session$setInputs(times = 100, resample = "observed", seed = 42, run = 1)
+        wait_for(session, function() state$status != "running")
+        deviance <<- session$returned()$bootstraps$beta_poisson$deviance
+      }
+    )
+    deviance
+  }
+  first <- run_once()
+  second <- run_once()
+  expect_length(first, 100L)
+  expect_equal(first, second)
+})
+
+test_that("invalid settings show an error and start nothing", {
+  local_workers()
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = shiny::reactiveVal(example_analysis())), {
+    session$setInputs(times = 10, resample = "observed", seed = NA, run = 1)
+    expect_equal(state$status, "error")
+    expect_match(state$message, "between 100 and 10,000")
+    expect_length(jobs(), 0L)
+    expect_false(mirai::daemons_set(.compute = bootstrap_compute_profile))
+  })
+})
+
+test_that("running with no analysis is refused", {
+  local_workers()
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = shiny::reactiveVal(NULL)), {
+    session$setInputs(times = 100, resample = "observed", seed = NA, run = 1)
+    expect_equal(state$status, "idle")
+    expect_length(jobs(), 0L)
+    expect_null(session$returned())
+  })
+})
+
+test_that("a second Run while running is ignored", {
+  # One worker and 1000 replicates per model keeps the run going long enough.
+  local_workers(1)
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = shiny::reactiveVal(example_analysis())), {
+    session$setInputs(times = 1000, resample = "observed", seed = 1, run = 1)
+    expect_equal(state$status, "running")
+    first <- jobs()
+    session$setInputs(run = 2)
+    expect_identical(jobs(), first)
+    session$setInputs(cancel = 1)
+    expect_equal(state$status, "cancelled")
+  })
+})
+
+test_that("cancelling stops the run and adds no bands", {
+  local_workers(1)
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = shiny::reactiveVal(example_analysis())), {
+    session$setInputs(times = 1000, resample = "observed", seed = 1, run = 1)
+    expect_equal(state$status, "running")
+
+    session$setInputs(cancel = 1)
+    expect_equal(state$status, "cancelled")
+    expect_length(jobs(), 0L)
+    # Stopped jobs report as finished; the poll must not collect them.
+    session$elapse(1200)
+    expect_equal(state$status, "cancelled")
+    expect_length(session$returned()$bootstraps, 0L)
+  })
+})
+
+test_that("a new analysis cancels a run and discards finished results", {
+  local_workers()
+  analysis <- shiny::reactiveVal(example_analysis())
+  shiny::testServer(mod_bootstrap_server, args = list(analysis = analysis), {
+    session$setInputs(times = 100, resample = "observed", seed = 1, run = 1)
+    expect_true(wait_for(session, function() state$status != "running"))
+    expect_equal(state$status, "done")
+    expect_length(session$returned()$bootstraps, 2L)
+
+    analysis(example_analysis(c("exponential", "exact_beta_poisson")))
+    session$flushReact()
+    expect_equal(state$status, "idle")
+    expect_length(session$returned()$bootstraps, 0L)
+  })
+})
+
+test_that("the Run buttons are task buttons that turn busy as soon as they are clicked", {
+  bootstrap_html <- as.character(mod_bootstrap_ui("bootstrap"))
+  expect_match(bootstrap_html, "bslib-task-button")
+  expect_match(bootstrap_html, "Preparing replicates...", fixed = TRUE)
+
+  fit_html <- as.character(mod_fit_ui("fit"))
+  expect_match(fit_html, "bslib-task-button")
+  expect_match(fit_html, "Fitting models...", fixed = TRUE)
+})
+
+test_that("the app turns on Shiny's busy indicators", {
+  # useBusyIndicators() works by setting these attributes on the page element.
+  html <- as.character(app_ui(NULL))
+  expect_match(html, "dataset.shinyBusyPulse = 'true'", fixed = TRUE)
+  expect_match(html, "dataset.shinyBusySpinners = 'true'", fixed = TRUE)
+})
diff --git a/app/tests/testthat/test-mod_explore.R b/app/tests/testthat/test-mod_explore.R
index f54aae4..9d25f86 100644
--- a/app/tests/testthat/test-mod_explore.R
+++ b/app/tests/testthat/test-mod_explore.R
@@ -52,3 +52,143 @@ test_that("the Explore tab explains the log axis and how to draw confidence band
   expect_match(html, "map lower to ymin and upper to ymax", fixed = TRUE)
   expect_match(html, "Ribbon", fixed = TRUE)
 })
+
+# The preset drives esquisse's own inputs; these tests check the sequence of
+# messages. That esquisse then draws the bands was checked in a browser.
+local_preset_recorder <- function(env = parent.frame()) {
+  calls <- new.env()
+  calls$log <- list()
+  # esquisse makes the same update calls for its own inputs; keep only ours.
+  ours <- c("dataset", unlist(esquisse_ids))
+  record <- function(what) function(session, inputId, ...) {
+    if (inputId %in% ours) {
+      calls$log[[length(calls$log) + 1L]] <- list(what = what, id = inputId, args = list(...))
+    }
+  }
+  testthat::local_mocked_bindings(updateSelectInput = record("dataset"), updateSelectizeInput = record("axis"),
+    .package = "singlehitapp", .env = env)
+  testthat::local_mocked_bindings(updateDragulaInput = record("aesthetics"), updateDropInput = record("geom"),
+    .package = "esquisse", .env = env)
+  calls
+}
+sent <- function(calls) vapply(calls$log, `[[`, "", "what")
+mapping_input <- function() {
+  c(list(xvar = "dose", yvar = "estimate"), list(ymin = "lower", ymax = "upper", fill = "level", facet = "model"))
+}
+mapping_id <- "esquisse-geomaes-aes_1-dragvars"
+
+test_that("the bands button only appears once a bootstrap is attached", {
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis())), {
+    expect_equal(output$has_bands, "no")
+  })
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    expect_equal(output$has_bands, "yes")
+  })
+  html <- as.character(mod_explore_ui("explore"))
+  expect_match(html, "explore-bands_preset", fixed = TRUE)
+  expect_match(html, "output.has_bands == &#39;yes&#39;", fixed = TRUE)
+})
+
+test_that("the bands preset switches the dataset, then the mapping, then the geometry and axis", {
+  calls <- local_preset_recorder()
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(dataset = "observed")
+    session$setInputs(bands_preset = 1)
+    expect_equal(sent(calls), "dataset")
+    expect_equal(calls$log[[1]]$args$selected, "curves")
+
+    session$setInputs(dataset = "curves")
+    session$flushReact()
+    expect_equal(sent(calls), c("dataset", "aesthetics"))
+    expect_equal(calls$log[[2]]$id, mapping_id)
+    expect_equal(calls$log[[2]]$args$selected, bands_preset$aesthetics)
+
+    # esquisse first reports a stale mapping; nothing more is sent for it.
+    stale <- list(target = list(xvar = "dose"))
+    do.call(session$setInputs, setNames(list(stale), mapping_id))
+    session$flushReact()
+    expect_equal(sent(calls), c("dataset", "aesthetics"))
+
+    do.call(session$setInputs, setNames(list(list(target = mapping_input())), mapping_id))
+    session$flushReact()
+    expect_equal(sent(calls), c("dataset", "aesthetics", "geom", "axis"))
+    expect_equal(calls$log[[3]]$args$selected, "ribbon")
+    expect_equal(calls$log[[4]]$args$selected, "log10")
+  })
+})
+
+test_that("the bands preset skips straight to the mapping when the curves are already shown", {
+  calls <- local_preset_recorder()
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(dataset = "curves")
+    session$setInputs(bands_preset = 1)
+    session$flushReact()
+    expect_equal(sent(calls), "aesthetics")
+  })
+})
+
+test_that("pressing the bands preset again only resets the geometry and axis", {
+  calls <- local_preset_recorder()
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(dataset = "curves")
+    do.call(session$setInputs, setNames(list(list(target = mapping_input())), mapping_id))
+    session$setInputs(bands_preset = 1)
+    session$flushReact()
+    expect_equal(sent(calls), c("geom", "axis"))
+  })
+})
+
+test_that("changing the mapping by hand does not trigger the preset", {
+  calls <- local_preset_recorder()
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(dataset = "curves")
+    do.call(session$setInputs, setNames(list(list(target = mapping_input())), mapping_id))
+    session$flushReact()
+    expect_length(calls$log, 0L)
+  })
+})
+
+test_that("a pending preset is dropped when the analysis is re-run", {
+  calls <- local_preset_recorder()
+  analysis <- shiny::reactiveVal(example_analysis_with_bands())
+  shiny::testServer(mod_explore_server, args = list(analysis = analysis), {
+    session$setInputs(dataset = "curves")
+    session$setInputs(bands_preset = 1)
+    session$flushReact()
+    expect_equal(sent(calls), "aesthetics")
+
+    # A different analysis: a reactiveVal ignores an identical value.
+    rerun <- example_analysis_with_bands()
+    rerun$note <- "rerun"
+    analysis(rerun)
+    session$flushReact()
+    do.call(session$setInputs, setNames(list(list(target = mapping_input())), mapping_id))
+    session$flushReact()
+    expect_equal(sent(calls), "aesthetics")
+  })
+})
+
+test_that("a pending preset is dropped when the user leaves the fitted curves", {
+  calls <- local_preset_recorder()
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(dataset = "curves")
+    session$setInputs(bands_preset = 1)
+    session$flushReact()
+    session$setInputs(dataset = "observed")
+    # esquisse warns about a mapping whose columns the observed data lacks.
+    suppressWarnings({
+      do.call(session$setInputs, setNames(list(list(target = mapping_input())), mapping_id))
+      session$flushReact()
+    })
+    expect_equal(sent(calls), "aesthetics")
+  })
+})
+
+test_that("the bands button is hidden when the bootstrap gives no bands", {
+  analysis <- example_analysis_with_bands()
+  analysis$bootstraps$exponential$converged <- FALSE
+  analysis$bootstraps$beta_poisson$converged <- FALSE
+  shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactiveVal(analysis)), {
+    expect_equal(output$has_bands, "no")
+  })
+})
diff --git a/app/tests/testthat/test-mod_figures.R b/app/tests/testthat/test-mod_figures.R
index c7165c3..c8e011c 100644
--- a/app/tests/testthat/test-mod_figures.R
+++ b/app/tests/testthat/test-mod_figures.R
@@ -39,3 +39,76 @@ test_that("the download writes a non-empty PNG and tolerates bad dimensions", {
     expect_equal(signature, as.raw(c(0x89, 0x50, 0x4e, 0x47, 0x0d, 0x0a, 0x1a, 0x0a)))
   })
 })
+
+test_that("the Plots tab says where confidence bands come from", {
+  shiny::testServer(mod_figures_server, args = list(analysis = shiny::reactiveVal(example_analysis())), {
+    session$setInputs(type = "overlay")
+    expect_match(as.character(output$band_hint$html), "never shows confidence bands")
+    session$setInputs(type = "panels")
+    expect_match(as.character(output$band_hint$html), "Run the bootstrap")
+  })
+  shiny::testServer(mod_figures_server, args = list(analysis = shiny::reactiveVal(example_analysis_with_bands())), {
+    session$setInputs(type = "panels")
+    expect_null(output$band_hint$html)
+  })
+})
+
+test_that("the single-model hint covers a model that has no bands", {
+  analysis <- example_analysis_with_bands()
+  analysis$bootstraps$exponential <- NULL
+  shiny::testServer(mod_figures_server, args = list(analysis = shiny::reactiveVal(analysis)), {
+    session$setInputs(type = "single", model = "exponential")
+    expect_match(as.character(output$band_hint$html), "No confidence bands for this model")
+    session$setInputs(type = "single", model = "beta_poisson")
+    expect_null(output$band_hint$html)
+  })
+})
+
+test_that("bands arriving while on the overlay switch the view and say so", {
+  updates <- list()
+  notices <- character()
+  testthat::local_mocked_bindings(
+    updateRadioButtons = function(session, inputId, ..., selected = NULL) {
+      updates[[length(updates) + 1L]] <<- list(inputId = inputId, selected = selected)
+    },
+    showNotification = function(ui, ...) {
+      notices <<- c(notices, as.character(ui))
+    },
+    .package = "singlehitapp"
+  )
+  analysis <- shiny::reactiveVal(example_analysis())
+  shiny::testServer(mod_figures_server, args = list(analysis = analysis), {
+    session$setInputs(type = "overlay")
+    session$flushReact()
+    expect_length(updates, 0L)
+
+    analysis(example_analysis_with_bands())
+    session$flushReact()
+    expect_length(updates, 1L)
+    expect_equal(updates[[1L]], list(inputId = "type", selected = "panels"))
+    expect_match(notices, "switched to one panel per model")
+
+    # The same bands again, or a rerun that still has bands, changes nothing.
+    analysis(example_analysis_with_bands())
+    session$flushReact()
+    expect_length(updates, 1L)
+  })
+})
+
+test_that("bands arriving while on the panel or single view leave the choice alone", {
+  for (view in c("panels", "single")) {
+    updates <- 0L
+    testthat::local_mocked_bindings(
+      updateRadioButtons = function(...) updates <<- updates + 1L,
+      showNotification = function(...) NULL,
+      .package = "singlehitapp"
+    )
+    analysis <- shiny::reactiveVal(example_analysis())
+    shiny::testServer(mod_figures_server, args = list(analysis = analysis), {
+      session$setInputs(type = view, model = "exponential")
+      analysis(example_analysis_with_bands())
+      session$flushReact()
+    })
+    expect_equal(updates, 0L)
+  }
+})
diff --git a/app/tests/testthat/test-mod_fit.R b/app/tests/testthat/test-mod_fit.R
index 96156f7..5cddb84 100644
--- a/app/tests/testthat/test-mod_fit.R
+++ b/app/tests/testthat/test-mod_fit.R
@@ -58,3 +58,15 @@ test_that("changing the data clears the previous analysis", {
     expect_null(session$returned())
   })
 })
+
+test_that("the Ward example: exponential fits badly and beta-Poisson is recommended", {
+  shiny::testServer(mod_fit_server, args = list(data = shiny::reactiveVal(ward_data())), {
+    session$setInputs(models = c("exponential", "beta_poisson"), alpha = 0.05, run = 1)
+    analysis <- session$returned()
+    expect_s3_class(analysis, "qdr_analysis")
+    expect_true(analysis$trend$passes)
+    expect_equal(headline_decision(analysis$assessment)$model, "Beta-Poisson")
+    exponential <- analysis$goodness_of_fit[analysis$goodness_of_fit$model == "exponential", ]
+    expect_false(exponential$good_fit)
+  })
+})
diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R
index 312b66c..1118733 100644
--- a/app/tests/testthat/test-mod_upload.R
+++ b/app/tests/testthat/test-mod_upload.R
@@ -43,6 +43,16 @@ test_that("loading a new source discards the previous column mapping", {
     # in a real session; the example must still parse with auto-detection.
     session$setInputs(example = 1)
     expect_true(all(c("dose", "positive", "negative") %in% names(session$returned())))
-    expect_equal(nrow(session$returned()), 3L)
+    expect_equal(nrow(session$returned()), 8L)
+  })
+})
+
+test_that("the example button loads the Ward rotavirus data", {
+  shiny::testServer(mod_upload_server, {
+    session$setInputs(example = 1)
+    data <- session$returned()
+    expect_equal(nrow(data), 8L)
+    expect_equal(range(data$dose), c(0.009, 90000))
+    expect_equal(sum(data$total), 59)
   })
 })