H5N1 Avian Influenza in Australia: State Burden, Stochastic Scenarios and Surveillance Indicators

1 Scientific purpose and interpretation

This report creates five publication-quality surveillance figures for H5 high pathogenicity avian influenza (HPAI), including H5N1 clade 2.3.4.4b, in Australia. It combines an official state-level data snapshot with deliberately conservative stochastic scenarios.

Important interpretation: the stochastic results are scenario simulations, not forecasts. They show what could occur under explicitly stated assumptions. They should not be interpreted as predictions of the true number of infected birds, because confirmed detections depend strongly on carcass discovery, reporting, sampling, laboratory confirmation, species distribution, and surveillance effort.

Official snapshot embedded below:

• 123 confirmed H5 detections in wild birds as at 11:00 AEST, 6 August 2026.
• State distribution: WA 10, SA 87, NSW 2, QLD 1, VIC 23; ACT, NT and TAS 0.
• No detections in poultry or the wider agricultural production system at that snapshot.
• The Australian public-health risk was assessed as low.

Primary official sources:

1. Australian Government Department of Agriculture, Fisheries and Forestry, Bird Flu dashboard: https://www.agriculture.gov.au/campaigns/birdflu
2. Australian Centre for Disease Control, Bird flu: https://www.cdc.gov.au/diseases/bird-flu-avian-influenza
3. World Organisation for Animal Health, first Australian HPAI H5N1 wild-bird notification: https://www.woah.org/en/australia-notifies-first-case-of-high-pathogenicity-avian-influenza-h5n1-in-a-wild-bird/

2 Embedded observed data

# -----------------------------------------------------------------------------
# OFFICIAL H5 DATA SNAPSHOT, 6 AUGUST 2026
# The value 123 is deliberately written explicitly in the code, as requested.
# These are confirmed detections in wild birds, not estimated infections.
# -----------------------------------------------------------------------------

national_h5_confirmed_detections <- 123L
national_poultry_detections <- 0L
snapshot_date <- as.Date("2026-08-06")

state_data <- tribble(
  ~state_code, ~state_name,                   ~confirmed_detections, ~confirmed_events,
  "ACT",       "Australian Capital Territory", 0L,                    0L,
  "NSW",       "New South Wales",              2L,                    2L,
  "NT",        "Northern Territory",            0L,                    0L,
  "QLD",       "Queensland",                    1L,                    1L,
  "SA",        "South Australia",              87L,                   22L,
  "TAS",       "Tasmania",                     0L,                    0L,
  "VIC",       "Victoria",                     23L,                   7L,
  "WA",        "Western Australia",            10L,                   10L
) %>%
  mutate(
    share_percent = 100 * confirmed_detections / national_h5_confirmed_detections,
    detections_per_event = if_else(confirmed_events > 0,
                                   confirmed_detections / confirmed_events,
                                   NA_real_)
  )

stopifnot(sum(state_data$confirmed_detections) == national_h5_confirmed_detections)
stopifnot(sum(state_data$confirmed_events) == 42L)

kable(
  state_data,
  digits = 2,
  caption = "Observed Australian H5 wild-bird detections and confirmed events, 6 August 2026"
)

Observed Australian H5 wild-bird detections and confirmed events, 6 August 2026

state_code	state_name	confirmed_detections	confirmed_events	share_percent	detections_per_event
ACT	Australian Capital Territory	0	0	0.00	NA
NSW	New South Wales	2	2	1.63	1.00
NT	Northern Territory	0	0	0.00	NA
QLD	Queensland	1	1	0.81	1.00
SA	South Australia	87	22	70.73	3.95
TAS	Tasmania	0	0	0.00	NA
VIC	Victoria	23	7	18.70	3.29
WA	Western Australia	10	10	8.13	1.00

The difference between detections and events is scientifically important. Multiple positive birds may belong to one epidemiological event. Therefore, neither measure should automatically be treated as incidence in the underlying wild-bird population.

3 Mathematical model

3.1 Observation-aware stochastic detection model

Let \(Y_{s,t}\) be the number of newly confirmed detections in scenario \(s\) and week \(t\). The simulation uses a negative-binomial observation process with a time-varying effective reproduction multiplier and density-dependent attenuation:

\[ Y_{s,t+1} \sim \operatorname{NegBin}\left(\mu_{s,t+1}, k_s\right), \]

\[ \mu_{s,t+1} = \left[Y_{s,t}R_{s,t}\left(1-\frac{C_{s,t}}{K_s}\right)_+ + \lambda_s\right]q_{s,t}, \]

\[ R_{s,t}=R_{0,s}\exp(-\delta_s t), \]

\[ C_{s,t+1}=C_{s,t}+Y_{s,t+1}. \]

Here:

• \(R_{0,s}\) is an initial effective propagation multiplier for detected clusters, not a directly estimated biological \(R_0\);
• \(\delta_s\) represents reduction through surveillance, behavioural change, carcass removal and biosecurity;
• \(K_s\) is a soft finite ceiling for additional confirmed detections over the short simulation horizon;
• \(\lambda_s\) is a low external-importation or discovery term;
• \(q_{s,t}\) is a surveillance ascertainment multiplier;
• \(k_s\) is the negative-binomial dispersion parameter, allowing clustering and superspreading-like variability;
• \((x)_+=\max(x,0)\).

This structure avoids mechanically generating unrealistically explosive counts. It is still a simplified surveillance model rather than a mechanistic wild-bird transmission model.

3.2 Scenario parameterisation

scenario_parameters <- tribble(
  ~scenario, ~description, ~R0, ~decline_delta, ~K_additional, ~importation_lambda, ~ascertainment_q, ~dispersion_k,
  "A: Rapid containment",
  "Strong surveillance and effective local interruption",
  0.72, 0.055, 140, 0.15, 0.90, 0.45,

  "B: Sustained low-level detection",
  "Persistent coastal detections with gradual control",
  0.92, 0.035, 200, 0.30, 1.00, 0.40,

  "C: Moderate wildlife expansion",
  "Limited local amplification across affected coastal systems",
  1.08, 0.045, 280, 0.45, 1.05, 0.35,

  "D: Adverse but bounded",
  "Higher short-term amplification followed by intensified control",
  1.22, 0.070, 380, 0.65, 1.10, 0.30
) %>%
  mutate(
    initial_weekly_detections = 18L,
    initial_cumulative = national_h5_confirmed_detections,
    horizon_weeks = params$horizon_weeks,
    simulations = params$n_sim
  )

kable(
  scenario_parameters,
  digits = 3,
  caption = "Parameters for conservative stochastic H5 detection scenarios"
)

Parameters for conservative stochastic H5 detection scenarios

scenario	description	R0	decline_delta	K_additional	importation_lambda	ascertainment_q	dispersion_k	initial_weekly_detections	initial_cumulative	horizon_weeks	simulations
A: Rapid containment	Strong surveillance and effective local interruption	0.72	0.055	140	0.15	0.90	0.45	18	123	12	3000
B: Sustained low-level detection	Persistent coastal detections with gradual control	0.92	0.035	200	0.30	1.00	0.40	18	123	12	3000
C: Moderate wildlife expansion	Limited local amplification across affected coastal systems	1.08	0.045	280	0.45	1.05	0.35	18	123	12	3000
D: Adverse but bounded	Higher short-term amplification followed by intensified control	1.22	0.070	380	0.65	1.10	0.30	18	123	12	3000

Parameter values are transparent scenario assumptions, not fitted estimates. They are chosen to create plausible bounded stress tests rather than extreme projections.

4 Figure 1: State-level burden choropleth

# ozmaps supplies Australian state polygons as an sf object.
aus_states <- ozmaps::ozmap_states

# Harmonise state labels robustly across ozmaps versions.
name_field <- intersect(c("NAME", "name", "STATE_NAME", "state"), names(aus_states))[1]
if (is.na(name_field)) stop("Could not identify the state-name field in ozmap_states.")

aus_map <- aus_states %>%
  mutate(state_name_map = .data[[name_field]]) %>%
  left_join(state_data, by = c("state_name_map" = "state_name")) %>%
  mutate(
    confirmed_detections = replace_na(confirmed_detections, 0L),
    burden_class = cut(
      confirmed_detections,
      breaks = c(-Inf, 0, 2, 10, 25, Inf),
      labels = c("0", "1–2", "3–10", "11–25", ">25"),
      right = TRUE
    )
  )

centroids <- suppressWarnings(st_point_on_surface(aus_map))
centroid_coords <- st_coordinates(centroids)
centroids <- centroids %>%
  mutate(x = centroid_coords[, 1], y = centroid_coords[, 2])

p1 <- ggplot(aus_map) +
  geom_sf(aes(fill = confirmed_detections), colour = "white", linewidth = 0.55) +
  geom_text(
    data = st_drop_geometry(centroids),
    aes(x = x, y = y, label = paste0(state_code, "\n", confirmed_detections)),
    fontface = "bold", size = 3.5
  ) +
  scale_fill_viridis_c(
    option = "C", trans = "sqrt", begin = 0.08, end = 0.95,
    breaks = c(0, 1, 2, 10, 25, 50, 87),
    name = "Confirmed\ndetections"
  ) +
  coord_sf(crs = 3577, expand = FALSE) +
  labs(
    title = "State burden of confirmed H5 detections in Australian wild birds",
    subtitle = "Laboratory-confirmed detections; square-root colour scale preserves low-count visibility",
    caption = paste0(
      "Data: Australian Government DAFF, 11:00 AEST, 6 August 2026. ",
      "National total = ", national_h5_confirmed_detections,
      ". Counts represent detected positive birds, not estimated infections."
    )
  ) +
  theme_surveillance() +
  theme(axis.text = element_blank(), axis.title = element_blank(), panel.grid = element_blank())

p1

ggsave("figures/Figure_1_state_H5_burden_map.png", p1, width = 10, height = 7, dpi = 400)

5 Figure 2: Detection burden versus epidemiological events

plot_state <- state_data %>%
  filter(confirmed_detections > 0) %>%
  arrange(confirmed_detections) %>%
  mutate(state_code = factor(state_code, levels = state_code))

p2a <- ggplot(plot_state, aes(x = confirmed_detections, y = state_code)) +
  geom_segment(aes(x = 0, xend = confirmed_detections, yend = state_code),
               linewidth = 1.2, colour = "grey75") +
  geom_point(aes(size = confirmed_events), shape = 21, stroke = 0.8) +
  geom_text(aes(label = confirmed_detections), hjust = -0.65, size = 3.8, fontface = "bold") +
  scale_size_continuous(range = c(4, 12), name = "Confirmed events") +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title = "Confirmed detections and event structure",
    subtitle = "Point size represents independent confirmed events",
    x = "Confirmed positive wild birds", y = NULL
  ) +
  theme_surveillance()

p2b <- plot_state %>%
  filter(confirmed_events > 0) %>%
  ggplot(aes(x = reorder(state_code, detections_per_event), y = detections_per_event)) +
  geom_col(width = 0.68) +
  geom_hline(yintercept = 1, linetype = "dashed", linewidth = 0.6) +
  geom_text(aes(label = sprintf("%.1f", detections_per_event)), vjust = -0.45, fontface = "bold") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.14))) +
  labs(
    title = "Within-event detection intensity",
    subtitle = "Higher values may indicate clustered mortality or intensified sampling",
    x = "State", y = "Detections per confirmed event"
  ) +
  theme_surveillance()

p2 <- p2a + p2b + plot_annotation(
  title = "H5 surveillance burden: positive birds versus confirmed events",
  caption = "An event may contain multiple positive birds; differences can reflect both epidemiology and surveillance practice."
)

p2

ggsave("figures/Figure_2_detection_event_structure.png", p2, width = 13, height = 6.5, dpi = 400)

6 Stochastic simulation engine

simulate_one_path <- function(par, horizon = params$horizon_weeks) {
  y <- integer(horizon + 1)
  cumulative <- numeric(horizon + 1)
  Rt <- numeric(horizon + 1)

  y[1] <- par$initial_weekly_detections
  cumulative[1] <- par$initial_cumulative
  Rt[1] <- par$R0

  for (tt in seq_len(horizon)) {
    Rt[tt + 1] <- par$R0 * exp(-par$decline_delta * tt)
    additional_so_far <- cumulative[tt] - par$initial_cumulative
    saturation <- max(0, 1 - additional_so_far / par$K_additional)

    mu <- (y[tt] * Rt[tt + 1] * saturation + par$importation_lambda) * par$ascertainment_q
    mu <- max(mu, 0.02)

    # R's rnbinom uses Var(Y)=mu + mu^2/size.
    y[tt + 1] <- rnbinom(1, size = par$dispersion_k, mu = mu)
    cumulative[tt + 1] <- cumulative[tt] + y[tt + 1]
  }

  tibble(
    week = 0:horizon,
    weekly_detections = y,
    cumulative_detections = cumulative,
    Rt = Rt
  )
}

run_scenario <- function(par_row, n_sim = params$n_sim) {
  par <- as.list(par_row)
  map_dfr(seq_len(n_sim), function(i) {
    simulate_one_path(par) %>%
      mutate(simulation = i, scenario = par$scenario)
  })
}

simulation_raw <- scenario_parameters %>%
  split(.$scenario) %>%
  map_dfr(run_scenario)

simulation_summary <- simulation_raw %>%
  group_by(scenario, week) %>%
  summarise(
    weekly_median = median(weekly_detections),
    weekly_q10 = quantile(weekly_detections, 0.10),
    weekly_q90 = quantile(weekly_detections, 0.90),
    cumulative_median = median(cumulative_detections),
    cumulative_q025 = quantile(cumulative_detections, 0.025),
    cumulative_q10 = quantile(cumulative_detections, 0.10),
    cumulative_q90 = quantile(cumulative_detections, 0.90),
    cumulative_q975 = quantile(cumulative_detections, 0.975),
    Rt_median = median(Rt),
    .groups = "drop"
  )

write.csv(simulation_summary, "outputs/stochastic_scenario_summary.csv", row.names = FALSE)

7 Figure 3: Weekly stochastic scenario envelopes

p3 <- ggplot(simulation_summary, aes(x = week, y = weekly_median)) +
  geom_ribbon(aes(ymin = weekly_q10, ymax = weekly_q90, fill = scenario),
              alpha = 0.22, colour = NA) +
  geom_line(aes(colour = scenario), linewidth = 1.0) +
  geom_hline(yintercept = 0, linewidth = 0.4) +
  facet_wrap(~scenario, ncol = 2, scales = "free_y") +
  scale_x_continuous(breaks = 0:params$horizon_weeks) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.08))) +
  guides(fill = "none", colour = "none") +
  labs(
    title = "Stochastic weekly H5 detection scenarios",
    subtitle = paste0(
      "Median and 10th–90th percentile envelopes from ",
      format(params$n_sim, big.mark = ","), " simulations per scenario"
    ),
    x = "Weeks after 6 August 2026",
    y = "New confirmed detections per week",
    caption = "Scenario simulations—not forecasts. Negative-binomial variation represents clustered detection and surveillance uncertainty."
  ) +
  theme_surveillance()

p3

ggsave("figures/Figure_3_weekly_stochastic_scenarios.png", p3, width = 12, height = 8, dpi = 400)

8 Figure 4: Cumulative detection trajectories with uncertainty

p4 <- ggplot(simulation_summary, aes(x = week, y = cumulative_median, colour = scenario)) +
  geom_ribbon(
    aes(ymin = cumulative_q025, ymax = cumulative_q975, fill = scenario),
    alpha = 0.12, colour = NA
  ) +
  geom_ribbon(
    aes(ymin = cumulative_q10, ymax = cumulative_q90, fill = scenario),
    alpha = 0.20, colour = NA
  ) +
  geom_line(linewidth = 1.15) +
  geom_point(data = filter(simulation_summary, week == 0), size = 2.7) +
  geom_hline(yintercept = national_h5_confirmed_detections,
             linetype = "dashed", linewidth = 0.65, colour = "grey30") +
  annotate("text", x = 0.3, y = national_h5_confirmed_detections + 8,
           label = paste0("Observed baseline = ", national_h5_confirmed_detections),
           hjust = 0, size = 3.6) +
  scale_x_continuous(breaks = 0:params$horizon_weeks) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.08))) +
  labs(
    title = "Bounded cumulative H5 detection scenarios",
    subtitle = "Dark ribbons: 80% interval; light ribbons: 95% interval",
    x = "Weeks after 6 August 2026",
    y = "Cumulative confirmed detections",
    colour = "Scenario", fill = "Scenario",
    caption = "The finite-capacity term limits short-horizon growth; it is a modelling safeguard, not a biological population ceiling."
  ) +
  theme_surveillance() +
  theme(legend.position = "bottom")

p4

ggsave("figures/Figure_4_cumulative_scenario_trajectories.png", p4, width = 11.5, height = 7, dpi = 400)

9 Figure 5: Threshold-based surveillance decision heatmap

The following plot transforms the simulations into operational quantities. For each scenario and week, it estimates the probability that cumulative detections exceed selected surveillance thresholds. These thresholds are illustrative and should be replaced by thresholds agreed with wildlife-health and biosecurity authorities.

thresholds <- c(150, 175, 200, 250, 300, 400)

threshold_probability <- simulation_raw %>%
  crossing(threshold = thresholds) %>%
  group_by(scenario, week, threshold) %>%
  summarise(
    exceedance_probability = mean(cumulative_detections >= threshold),
    .groups = "drop"
  ) %>%
  mutate(
    scenario_short = str_remove(scenario, ":.*$"),
    threshold_label = paste0("≥", threshold)
  )

p5 <- ggplot(
  threshold_probability,
  aes(x = week, y = threshold_label, fill = exceedance_probability)
) +
  geom_tile(colour = "white", linewidth = 0.35) +
  geom_text(
    aes(label = ifelse(exceedance_probability >= 0.05,
                       sprintf("%.0f%%", 100 * exceedance_probability), "")),
    size = 2.7
  ) +
  facet_wrap(~scenario, ncol = 2) +
  scale_fill_viridis_c(
    option = "B", limits = c(0, 1),
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    labels = scales::percent_format(accuracy = 1),
    name = "Probability"
  ) +
  scale_x_continuous(breaks = 0:params$horizon_weeks) +
  labs(
    title = "Probability of crossing cumulative surveillance thresholds",
    subtitle = "Decision-oriented summary of stochastic uncertainty",
    x = "Weeks after 6 August 2026",
    y = "Cumulative-detection threshold",
    caption = "Illustrative thresholds only. Operational triggers should incorporate species, location, mortality clustering and poultry proximity."
  ) +
  theme_surveillance() +
  theme(panel.grid = element_blank())

p5

ggsave("figures/Figure_5_threshold_exceedance_heatmap.png", p5, width = 12, height = 8, dpi = 400)

10 Numerical scenario summary

final_week_summary <- simulation_summary %>%
  filter(week == params$horizon_weeks) %>%
  transmute(
    Scenario = scenario,
    `Median cumulative detections` = round(cumulative_median),
    `80% interval` = paste0(round(cumulative_q10), "–", round(cumulative_q90)),
    `95% interval` = paste0(round(cumulative_q025), "–", round(cumulative_q975)),
    `Median final-week detections` = round(weekly_median, 1),
    `Median final-week R multiplier` = round(Rt_median, 3)
  )

kable(
  final_week_summary,
  caption = paste0("Scenario outcomes at week ", params$horizon_weeks,
                   "; simulations are not forecasts")
)

Scenario outcomes at week 12; simulations are not forecasts

Scenario	Median cumulative detections	80% interval	95% interval	Median final-week detections	Median final-week R multiplier
A: Rapid containment	132	124–179	123–226	0	0.372
B: Sustained low-level detection	143	126–230	124–318	0	0.604
C: Moderate wildlife expansion	149	127–272	124–401	0	0.629
D: Adverse but bounded	154	129–306	125–515	0	0.527

11 Surveillance interpretation

The figures are designed for prevention and monitoring:

1. State map: identifies where confirmed burden is concentrated and where zeros may reflect either absence or insufficient detection.
2. Event-intensity panel: distinguishes isolated detections from clusters containing multiple positive birds.
3. Weekly scenario envelopes: visualise near-term workload and uncertainty for laboratories, wildlife responders and carcass-management teams.
4. Cumulative trajectories: support bounded resource-planning exercises without presenting an unchecked exponential curve.
5. Threshold heatmap: translates uncertainty into probabilities that can be linked to pre-agreed escalation actions.

A responsible operational trigger should not use counts alone. It should combine:

\[ \text{Priority}_{i,t}=w_1Z(\text{mortality cluster})+w_2Z(\text{new species})+ w_3Z(\text{poultry proximity})+w_4Z(\text{geographic spread})+w_5Z(\text{mammal involvement}). \]

Weights \(w_j\) should be set by epidemiologists, wildlife ecologists, veterinary authorities and public-health agencies.

12 Limitations

• Confirmed detections substantially undercount all infected wild birds.
• Surveillance intensity differs by state, coastline accessibility, species and time.
• A detected carcass is not equivalent to a new transmission event.
• The initial weekly value is a scenario starting condition, not a fitted incidence estimate.
• The model does not include migration networks, species-specific susceptibility, environmental persistence, weather, poultry-farm connectivity or genomic linkage.
• The finite-capacity term is included to keep simulations conservative; it should be re-estimated when longer time series become available.
• Because the official situation changes quickly, update state_data, national_h5_confirmed_detections, and snapshot_date before publication.

13 Reproducibility information

cat("Report generated with seed:", params$seed, "\n")

## Report generated with seed: 20260806

cat("Simulations per scenario:", params$n_sim, "\n")

## Simulations per scenario: 3000

cat("Simulation horizon:", params$horizon_weeks, "weeks\n\n")

## Simulation horizon: 12 weeks

sessionInfo()

## R version 4.4.3 (2025-02-28)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.7.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Australia/Sydney
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.3.2 ozmaps_0.4.5    sf_1.0-21       knitr_1.50     
##  [5] stringr_1.5.1   tibble_3.2.1    purrr_1.0.4     tidyr_1.3.1    
##  [9] dplyr_1.1.4     ggplot2_4.0.3  
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.9         generics_0.1.3     class_7.3-23       KernSmooth_2.23-26
##  [5] stringi_1.8.7      digest_0.6.37      magrittr_2.0.3     evaluate_1.0.3    
##  [9] grid_4.4.3         RColorBrewer_1.1-3 oz_1.0-22          fastmap_1.2.0     
## [13] jsonlite_1.9.1     e1071_1.7-16       DBI_1.2.3          viridisLite_0.4.2 
## [17] scales_1.4.0       textshaping_1.0.0  jquerylib_0.1.4    cli_3.6.4         
## [21] rlang_1.3.0        units_0.8-7        withr_3.0.2        cachem_1.1.0      
## [25] yaml_2.3.10        tools_4.4.3        vctrs_0.6.5        R6_2.6.1          
## [29] proxy_0.4-27       lifecycle_1.0.4    classInt_0.4-11    ragg_1.4.0        
## [33] pkgconfig_2.0.3    pillar_1.10.1      bslib_0.9.0        gtable_0.3.6      
## [37] glue_1.8.0         Rcpp_1.0.14        systemfonts_1.3.2  xfun_0.52         
## [41] tidyselect_1.2.1   rstudioapi_0.17.1  farver_2.1.2       htmltools_0.5.8.1 
## [45] labeling_0.4.3     rmarkdown_2.29     compiler_4.4.3     S7_0.2.2