---
title: "Checking for alternative stable states"
description: "Testing alternative stable states in R: why a bimodal distribution is not proof, why a driver mimics one, and how perturbation and hysteresis give the evidence."
date: "2026-08-03 09:00"
categories: [ecology tutorial, R, community ecology, model diagnostics, theoretical ecology]
image: thumbnail.png
image-alt: "Two bimodal state distributions that look identical, one from a driver threshold and one from true bistability, distinguished only by a perturbation test."
---
```{r setup}
#| include: false
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
dev = "png", dpi = 120, fig.path = "figures/",
fig.width = 7, fig.height = 4.4)
library(ggplot2)
theme_te <- theme_minimal(base_size = 12) +
theme(panel.background = element_rect(fill = "white", colour = NA),
plot.background = element_rect(fill = "white", colour = NA),
panel.grid.minor = element_blank(),
plot.title = element_text(face = "bold", size = 12))
theme_set(theme_te)
col_a <- "#3a6b52"; col_b <- "#b5651d"; col_ref <- "#555555"
set.seed(169)
env <- runif(500, 0, 4)
state_mono <- 9/(1 + exp(4*(env - 2))) + rnorm(500, 0, 0.5) # single state tracking a driver
z <- rbinom(500, 1, 0.5)
state_bi <- ifelse(z == 1, 8.5, 0.6) + rnorm(500, 0, 0.5) # two states under the same driver
cor_mono <- cor(env, state_mono); cor_bi <- cor(env, state_bi)
# bistable competition model for the perturbation test
r1 <- r2 <- 0.5; K1 <- K2 <- 100; a12 <- a21 <- 1.5
run_lv <- function(N) { for (t in 1:2000) {
N <- pmax(N + c(r1*N[1]*(1-(N[1]+a12*N[2])/K1), r2*N[2]*(1-(N[2]+a21*N[1])/K2))*0.1, 0) }; round(N) }
push_hi <- run_lv(c(30, 70)); push_lo <- run_lv(c(70, 30))
```
Alternative stable states are one of the most consequential ideas in ecology and one of the hardest to prove. The trouble is that the pattern they leave, a system that sits in one of two very different conditions, is also produced by things that are not alternative stable states at all: a sharp response to an environmental gradient, a threshold in a driver, plain measurement of two different places. A field data set almost never comes with the manipulation that would settle the question, so most claims rest on patterns that have innocent explanations. These three checks are what a claim of alternative stable states has to survive.
## Check 1: a bimodal distribution is not bistability
The most common evidence offered is a histogram with two peaks: many sites are either high or low, few in between, so surely there are two states. But a single-valued response to a driver produces exactly the same histogram. If a system tracks one stable state that happens to change steeply with some environmental variable, then wherever the driver is common at its low and high ends, the states pile up at two values with a gap in the middle, and it is not bistable at all.
```{r bimodal}
set.seed(169)
env <- runif(500, 0, 4)
state_mono <- 9/(1 + exp(4*(env - 2))) + rnorm(500, 0, 0.5) # ONE state, steep function of a driver
state_bi <- ifelse(rbinom(500, 1, 0.5) == 1, 8.5, 0.6) + rnorm(500, 0, 0.5) # TWO states
c(mono_gap = mean(state_mono > 2 & state_mono < 6), # both are bimodal: few in the middle
bi_gap = mean(state_bi > 2 & state_bi < 6))
```
Both distributions are bimodal, with only a small fraction of sites in the middle, and by eye the two histograms are the same. One comes from a single state bent steeply by a driver, the other from genuine bistability. The histogram cannot tell them apart, which means a bimodal distribution is a starting point for a question, never an answer to it.
```{r fig-bimodal}
#| fig-cap: "Two state distributions, both bimodal. The left comes from a single stable state responding steeply to a driver; the right from true bistability. The shapes are indistinguishable, so bimodality alone decides nothing."
#| fig-alt: "Two histograms side by side, both with a low peak and a high peak and few values in between, looking essentially the same."
db <- rbind(data.frame(x = state_mono, src = "one state, steep driver"),
data.frame(x = state_bi, src = "two stable states"))
ggplot(db, aes(x, fill = src)) +
geom_histogram(bins = 30, colour = "white", linewidth = 0.1) +
facet_wrap(~src) +
scale_fill_manual(values = c("one state, steep driver" = col_b, "two stable states" = col_a), guide = "none") +
labs(x = "State", y = "Count", title = "Same bimodal pattern, different cause") +
theme(strip.text = element_text(size = 9))
```
## Check 2: does a driver explain the two groups?
The first thing to try is cheap: see whether a measured driver accounts for which state each site is in. If it does, the parsimonious reading is a threshold response, not two states; if the states are unrelated to every driver you can measure, alternative stable states become more plausible. It is a screen, not a proof, but it removes the easy explanations.
```{r driver}
c(driver_explains_mono = round(cor(env, state_mono), 2),
driver_explains_bi = round(cor(env, state_bi), 2))
```
In the single-state data the driver and the state are tightly related, a correlation of `r sprintf("%.2f", cor_mono)`: the environment tells you almost exactly which mode a site is in, so there is no need to invoke two states. In the bistable data the same driver is essentially uncorrelated with the state, `r sprintf("%.2f", cor_bi)`: sites in opposite states sit at the same environmental value. That is the signature worth chasing, though it is still not conclusive, because an unmeasured driver could always be doing the work of the one you measured.
## Check 3: the perturbation and the hysteresis
The only decisive evidence is dynamical, and it takes a manipulation. Alternative stable states make a specific promise: push the system hard enough to cross the boundary and it stays in the new state after the push ends, under the same conditions as before. A monostable system, perturbed, slides back. So the definitive test is to displace the system and watch whether it returns or holds. In the bistable model, starting the community on either side of the separatrix sends it to a different, permanent state.
```{r perturb}
# same parameters, two starting points either side of the boundary: the state that persists differs
rbind(pushed_toward_2 = run_lv(c(30, 70)),
pushed_toward_1 = run_lv(c(70, 30)))
```
A community displaced toward species 2 settles at `r sprintf("(%.0f, %.0f)", push_hi[1], push_hi[2])` and stays there; displaced toward species 1 it settles at `r sprintf("(%.0f, %.0f)", push_lo[1], push_lo[2])` and stays there. The persistence of the new state after the perturbation is removed is the thing a single stable state cannot do, and it is what the classic field experiments on rocky shores and lakes were designed to show. The comparable signature along a driver is [hysteresis](../hysteresis-and-tipping-points/): if pushing the driver up collapses the system at one value and pulling it back only restores it at a much lower value, the loop is hard to explain without two states. Both tests need the system to be driven, not just observed.
```{r fig-perturb}
#| fig-cap: "A perturbation test. Displaced across the boundary toward one species or the other, the community settles into a different state and stays: the state depends on the perturbation, not just the conditions, which a single stable state cannot do."
#| fig-alt: "A bar chart of the final abundances of two species under two perturbations, showing species 1 dominant in one case and species 2 dominant in the other under identical conditions."
dp <- rbind(data.frame(perturbation = "pushed toward species 2", sp = c("species 1", "species 2"), N = push_hi),
data.frame(perturbation = "pushed toward species 1", sp = c("species 1", "species 2"), N = push_lo))
ggplot(dp, aes(sp, N, fill = sp)) +
geom_col(width = 0.6) + facet_wrap(~perturbation) +
scale_fill_manual(values = c("species 1" = col_a, "species 2" = col_b), guide = "none") +
labs(x = NULL, y = "Final abundance", title = "The state persists after the push") +
theme(strip.text = element_text(size = 9))
```
## The weight of the evidence
Three checks, and they climb a ladder of difficulty. Bimodality is free and proves nothing; a driver screen is cheap and can rule out the easy alternatives; only a perturbation or a demonstrated hysteresis loop, both of which need the system actually manipulated, can carry the claim. This is why rigorously confirmed alternative stable states in nature are far rarer than the idea's popularity suggests, and why careful reviews keep returning most field examples to the "suggestive, not established" pile. None of this is an argument that alternative stable states are unreal; the theory is sound and the best experimental cases are convincing. It is an argument that the bar is a manipulation, and that a bimodal histogram from a snapshot survey, however striking, has not cleared it.
## Related tutorials
- [Priority effects and alternative stable states](../priority-effects-and-alternative-states/)
- [Basins of attraction](../basins-of-attraction/)
- [Hysteresis and tipping points](../hysteresis-and-tipping-points/)
- [Checking an Allee analysis](../checking-an-allee-analysis/)
## References
- Schroder A, Persson L, De Roos AM 2005. Oikos 110(1):3-19 (10.1111/j.0030-1299.2005.13962.x).
- Petraitis PS, Latham RE 1999. Ecology 80(2):429-442 (10.1890/0012-9658(1999)080[0429:TIOASS]2.0.CO;2).