# Figure 2 — data and code

Supplementary material for *From function to fitness: a tripartite view of sleep*
(G. F. Gilestro). Everything needed to reproduce Figure 2 and every quantity
quoted in its legend.

## Provenance

All primary behavioural data are from Joyce, M. *et al.* Divergent evolution of
sleep in *Drosophila* species. *Nat. Commun.* **15**, 5091 (2024) — reference 135
of the manuscript. The files here are derived products: four-state hidden Markov
models fitted to those recordings, and the bout durations obtained by Viterbi
decoding them. They are not the raw recordings, which remain available through
the data availability statement of that paper.

States are ordered by activity, `activity_rank` 0 being the deepest rest:

| rank | `canonical_label` | HMM state |
|------|-------------------|-----------|
| 0 | Deep sleep   | `s0` |
| 1 | Light sleep  | `s1` |
| 2 | Quiet awake  | `s2` |
| 3 | Full awake   | `s3` |

The ordering by depth was validated in the source paper against a robotically
measured arousal threshold.

## Files

### `metastability_model.py`

The model and the figure, in one documented script. Its header docstring states
the model, its assumptions and its limitations in full. In outline it defines
rest as a state leaking to wakefulness at hazard λ set by the arousal threshold
θ, fixes θ for each state from the measured bout durations, computes the optimal
threshold both numerically and in closed form, and draws all six panels.

Run it from this directory:

```
python3 metastability_model.py
```

It writes `Figure_2.svg` and `Figure_2.png`, and prints the measured hazards,
the optimum, and the panel e model fit. Requires `numpy`, `pandas`, `scipy` and
`matplotlib`; nothing else, and no network access.

### `hmm_bout_durations.csv`

950,651 Viterbi-decoded bouts from 776 individual flies across eight datasets
(seven species plus a wild-caught *melanogaster* strain), about 2,400 fly-days.
One row per bout. This is the uncircular dataset: decoded state durations, not
the geometric dwell times the HMM itself assumes.

| column | meaning |
|--------|---------|
| `species` | dataset label (`D.ere`, `D.mel`, `D.sec`, `D.sims`, `D.vir`, `D.wil`, `D.yak`, `wild_d.mel`) |
| `id` | individual fly (ethoscope recording ID and channel) |
| `state` | HMM state index |
| `canonical_label` | state name, per the table above |
| `activity_rank` | 0–3, 0 = deepest rest |
| `p_move` | emission probability of movement for that state |
| `duration_min` | bout length in minutes |
| `n_epochs` | bout length in 60 s epochs |
| `right_censored` | `True` if the bout was truncated at a recording boundary |

**On the analysed *n*.** The script excludes the two per-species *melanogaster*
model sets, whose deep-sleep dwell collapses to about a minute; the canonical
`get_HMM("M")` model covers that lineage. Panels a, b and f therefore analyse
**516,543** bouts across six datasets, of which 2.8% are right-censored. The
figure legend quotes this analysed count; the file holds all 950,651 so the
exclusion can be inspected or reversed.

### `hmm_lambda_summary.csv`

Per-state exit hazards and mean dwell times for every fitted model, one row per
state per model. `A_ii` is the self-transition probability,
`lambda_per_min = (1 - A_ii) / dt`, and `mean_dwell_min = 1 / lambda_per_min`.
`model_set` distinguishes the two per-species fits (`2`, `new`) from the
canonical *melanogaster* model (`M`). Provided for reference; the script reads
its hazards from the decoded bouts instead.

### `hmm_transition_matrices.json`

The fitted models themselves, keyed `species|model_set` (plus `get_HMM_M`), 17
in total. Each carries `n_states`, `dt_seconds`, `transmat`, `startprob`,
`emissionprob`, `p_move` and `A_ii`. Top-level keys give the epoch length
(60 s), the emission symbols (`0` = not moving, `1` = moving) and the hazard
convention. Provided so the decoding can be reproduced from the source
recordings.

### `sleep_state_24h_summary.csv`

Hourly state occupancies across the 24 h day, used for panel e. One row per
species per Zeitgeber hour. `s0_mean`…`s3_mean` are the mean fractions of flies
in each state with their SEMs, and sum to 1 within rounding; `asleep_mean` is
`s0 + s1`. `n_flies` is the number contributing to that hour. Panel e uses
`s0` (deep sleep) over ZT12–ZT23, the night.

## Reproducing the figure

Extract all files to one directory and run the script there. It resolves the
data by filename in the working directory. Expected output, which matches the
figure legend:

```
MEASURED from hmm_bout_durations.csv (516543 bouts, 2.8% censored):
  Deep sleep   mean bout  27.03 min -> lambda = 0.0370/min
  Light sleep  mean bout   4.06 min -> lambda = 0.2463/min
  Quiet awake  mean bout   2.19 min -> lambda = 0.4572/min
  Full awake   mean bout  36.75 min -> lambda = 0.0272/min

Optimal threshold theta* = 3.83 (numeric) / 3.84 (closed form)
  model fit to the front-loader mean: R^2=0.892
Panel F: monotone ladder (deep > light > quiet) in 6/6 datasets
```

The thresholds in panels a–b and the resettling time in b are set from the
measured bouts. The cost parameters in panels c–d are illustrative rather than
fitted; the qualitative results there require only that genuine threats be more
salient than benign stimuli.
