"""Part 1 inputs: dated classic transitions, dated environmental gates, physiological oxygen thresholds and the parametric oxygen-history ensemble. All ages in Ga. Every number carries its source; 'ver' says how it was checked in this session (A = abstract or full text read via Europe PMC; T = taken from this project's timing study events table, which records its own verification; M = standard value recalled, not re-verified). """ import numpy as np import pandas as pd from config import DATA from curves import pal_to_mr # -------------------------------------------------------------------------------------------------------- # Transition dates. kind 'u': uniform on [lo, hi] (used when the source gives a range across analyses); # kind 'n': two-piece normal with the given centre and a 95 per cent interval [lo, hi]. # -------------------------------------------------------------------------------------------------------- TRANSITIONS = { "life": dict(name="Life (LUCA)", estimates=[ ("Moody2024", "n", 4.20, 4.09, 4.33, "Moody et al. 2024 Nat Ecol Evol 8:1654: LUCA 4.2 Ga (4.09-4.33)", "T"), ]), "oxyphoto": dict(name="Oxygenic photosynthesis", estimates=[ ("Fournier2021", "n", 2.960, 2.780, 3.275, "Fournier et al. 2021 Proc R Soc B 288:20210675: crown Cyanobacteria mean 2960 Ma (95% CI 3275-2780)", "A"), ("Moody2024", "n", 2.99, 2.95, 3.09, "Moody et al. 2024 repository: total-group Oxyphotobacteria 2.99 Ga (2.95-3.09)", "T"), ("Shih2017", "u", None, 2.43, 2.60, "Shih et al. 2017 Geobiology 15:19: Oxyphotobacteria/Melainabacteria split 2.5-2.6 Ga is an upper " "limit; crown group ~2.0 Ga; origin placed between the split and the onset of the GOE", "A"), ("Fournier2021euk", "u", None, 2.392, 2.733, "Fournier et al. 2021: crown Cyanobacteria 2733-2392 Ma when only eukaryotic calibrations are used", "A"), ]), "eukaryotes": dict(name="Eukaryotes (LECA)", estimates=[ ("Parfrey2011", "u", None, 1.679, 1.866, "Parfrey et al. 2011 PNAS 108:13624: LECA 1866-1679 Ma", "A"), ("Betts2018", "u", None, 1.21, 1.84, "Betts et al. 2018 Nat Ecol Evol 2:1556: crown eukaryotes < 1.84 Ga (1.21-1.84 as tabulated in this " "project's convergence study)", "A"), ("Strassert2021", "u", None, 1.958, 2.386, "Strassert et al. 2021 Nat Commun 12:1879: LECA between 2386 and 1958 Ma", "A"), ("Kay2026", "n", 1.735, 1.67, 1.80, "Kay et al. 2026 Nature 650:129: LECA 1.80-1.67 Ga", "A"), ]), "mitochondria": dict(name="Mitochondrial endosymbiosis", estimates=[ ("Kay2026", "n", 2.20, 2.08, 2.32, "Kay et al. 2026: endosymbiosis aligned to ~2.2 Ga (mFECA 2.37-2.13 Ga; width taken from mFECA)", "A"), ("Betts2018", "u", None, 1.21, 2.053, "Betts et al. 2018: mitochondrial symbiosis 2.053-1.21 Ga", "A"), ]), "animals": dict(name="Animals (crown Metazoa)", estimates=[ ("dosReis2015", "u", None, 0.650, 0.833, "dos Reis et al. 2015 Curr Biol 25:2939: 833-650 Ma", "A"), ("LoleDurbin2026", "u", None, 0.700, 0.800, "Lole Durbin et al. 2026 Sci Adv 12:eaeg6289: ~800-700 Ma", "A"), ("Carlisle2024", "u", None, 0.5934, 0.6132, "Carlisle et al. 2024 Sci Adv 10:eadp7161: 613.2-593.4 Ma averaged over analyses", "A"), ]), "landplants": dict(name="Land plants (crown Embryophyta)", estimates=[ ("Morris2018", "u", None, 0.470, 0.515, "Morris et al. 2018 PNAS 115:E2274 (Cambrian origin most probable); cryptospores 0.473-0.471 Ga " "(Rubinstein et al. 2010)", "A"), ]), "humans": dict(name="Homo sapiens / technology", estimates=[ ("now", "u", None, 0.0, 0.0003, "Richter et al. 2017 Nature 546:293 (315 ka); zero on this scale", "T"), ]), } # -------------------------------------------------------------------------------------------------------- # Dated environmental events (uniform ranges, Ga) # -------------------------------------------------------------------------------------------------------- ENV = { "habitable": (4.404, 4.51, "Jack Hills zircon 4.404 Ga to Moon-forming impact 4.51 Ga (timing study)", "T"), "goe_onset": (2.426, 2.50, "S-MIF to S-MDF transition between 2501.5 and 2434 Ma (Warke et al. 2020 PNAS " "117:13314); 2460-2426 Ma (Gumsley et al. 2017 PNAS 114:1811)", "A"), "goe_permanent": (2.20, 2.25, "Permanent oxygenation ~2.22 Ga (Poulton et al. 2021 Nature 592:232); " "oscillations until <= 2250-2240 Ma (Gumsley et al. 2017)", "A"), "lomagundi_end": (2.04, 2.08, "End of the Lomagundi excursion ca. 2060 Ma (Gumsley et al. 2017)", "A"), "mesoprot_pulse": (1.36, 1.44, "Oxygenation episode ~1.4 Ga, Xiamaling Fm (Zhang et al. 2016 PNAS 113:1731); " "contested", "A"), "tonian_onset": (0.75, 0.85, "Onset of Neoproterozoic change ~800 Ma: Cr-isotope shift (Planavsky et al. " "2014), 'Tonian transformation' (Brocks et al. 2023 Nature 618:767)", "A"), "noe_window": (0.539, 0.80, "Neoproterozoic oxygenation treated as one protracted window. This is this project's " "construction (timing study), NOT a gate named by Mills et al. 2025, who call the animal " "threshold 'controversial and essentially unknown' and commit only to animals being excluded " "before ~2.2 Ga", "T"), "t800": (0.80, 0.80, "Onset of the Neoproterozoic oxygenation taken as exactly 0.80 Ga (timing study " "sensitivity case)", "T"), "sturtian_onset": (0.715, 0.719, "Sturtian onset ~717 Ma (Rooney et al. 2015 Geology 43:459)", "M"), "sturtian_end": (0.658, 0.662, "Sturtian termination ~660 Ma (Rooney et al. 2015)", "M"), "marinoan_end": (0.634, 0.636, "Marinoan termination ~635 Ma (Rooney et al. 2015)", "M"), "ediacaran_ox": (0.55, 0.56, "Ocean oxygenation ca. 550-560 Ma from Mo (Dahl et al. 2010 PNAS 107:17911)", "A"), "poe": (0.40, 0.42, "Palaeozoic oxygenation: O2 near modern by 420-400 Ma (Lenton et al. 2016 PNAS 113:9704; " "Mills et al. 2025)", "A"), "permocarb_peak": (0.28, 0.31, "Permo-Carboniferous O2 maximum (GEOCARBSULFOR and COPSE curves of this " "study)", "A"), "end_permian": (0.2515, 0.2525, "End-Permian extinction 251.9 Ma (ICS chart)", "M"), "kpg": (0.0655, 0.0665, "Cretaceous-Palaeogene impact 66.0 Ma (ICS chart)", "M"), "eocene_oligocene": (0.0335, 0.0345, "Antarctic glaciation 33.9 Ma (ICS chart)", "M"), } # pool of candidate gates for the look-elsewhere analysis (fixed list, chosen for being the step changes that # a geology textbook would list, not for fitting the transitions) POOL_OXYGEN = ["goe_onset", "goe_permanent", "lomagundi_end", "mesoprot_pulse", "tonian_onset", "ediacaran_ox", "poe", "permocarb_peak"] POOL_ALL = POOL_OXYGEN + ["sturtian_onset", "sturtian_end", "marinoan_end", "end_permian", "kpg", "eocene_oligocene"] # -------------------------------------------------------------------------------------------------------- # Physiological oxygen thresholds fixed in advance (fractions of PAL); log-uniform between lo and hi # -------------------------------------------------------------------------------------------------------- THRESH = { "euk": (1e-5, 4e-3, "aerobic eukaryotes 0.001-0.4% PAL (Mills et al. 2025; nanomolar O2 suffices for " "sterol synthesis, Waldbauer et al. 2011, and aerobic growth, Stolper et al. 2010)"), "animal_sponge": (0.005, 0.04, "sponge-grade animals 0.5-4% PAL (Mills DB et al. 2014 PNAS 111:4168)"), "animal_mills25": (0.001, 0.01, "animals 0.1-1% PAL (Mills et al. 2025)"), "animal_catling": (0.047, 0.047, "diffusion-limited organisms of ~1 mm: > ~1e3 Pa = 4.7% PAL " "(Catling et al. 2005 Astrobiology 5:415)"), "animal_catling_circ": (0.047, 0.47, "cm-scale animals with circulation: 1e3-1e4 Pa = 4.7-47% PAL " "(Catling et al. 2005)"), "human": (0.53, 0.59, "human-type physiology 53-59% PAL (Mills et al. 2025)"), "fire": (0.878, 0.878, "open-air combustion 87.8% PAL (Mills et al. 2025; ~18.5 per cent O2)"), } # mid-Proterozoic oxygen level scenarios (fractions of PAL, log-uniform) MIDPROT = { "low": (0.001, 0.01, "<= 0.1% PAL (Planavsky et al. 2014 Science 346:635) up to 1% PAL"), "high": (0.01, 0.10, ">= 1% PAL (Canfield et al. 2018 Nat Commun 9:2871), >= 4% PAL (Zhang et al. 2016); " "'usually ~1 to 10% PAL' (Mills et al. 2025)"), "all": (0.001, 0.10, "full published range 0.1-10% PAL"), } NEOPROT = (0.01, 0.50, "Neoproterozoic O2 oscillating between ~1 and ~50% PAL (Krause et al. 2022 Sci Adv " "8:eabm8191)") def two_piece_normal(rng, mode, lo, hi, size): s_lo = (mode - lo) / 1.96 s_hi = (hi - mode) / 1.96 w = s_hi / (s_lo + s_hi) up = rng.uniform(size=size) < w z = np.abs(rng.standard_normal(size)) return np.where(up, mode + s_hi * z, mode - s_lo * z) def draw_transition(rng, key, size, which=None): """Return ages (Ga) and the index of the estimate each draw came from (equal weights). `which` may be an estimate label, or a tuple (lo, hi) for a uniform range not in the table (index -1).""" if isinstance(which, (tuple, list)): return rng.uniform(which[0], which[1], size=size), np.full(size, -1) est = TRANSITIONS[key]["estimates"] if which is None: idx = rng.integers(0, len(est), size=size) else: idx = np.full(size, [e[0] for e in est].index(which)) out = np.empty(size) for k, (lab, kind, mode, lo, hi, _, _) in enumerate(est): m = idx == k n = int(m.sum()) if n == 0: continue if kind == "u": out[m] = rng.uniform(lo, hi, size=n) else: out[m] = two_piece_normal(rng, mode, lo, hi, n) return out, idx def draw_env(rng, key, size): """key: a name in ENV, a tuple (lo, hi) in Ga, or {'mix': [keys...]} for an equal-weight mixture of definitions (one definition per draw).""" if isinstance(key, dict): parts = key["mix"] which = rng.integers(0, len(parts), size=size) out = np.empty(size) for k, part in enumerate(parts): m = which == k out[m] = draw_env(rng, part, int(m.sum())) return out lo, hi = (key[0], key[1]) if isinstance(key, (tuple, list)) else (ENV[key][0], ENV[key][1]) return rng.uniform(lo, hi, size=size) def draw_luca_before(rng, start_age, margin=0.01): """LUCA ages (Ga) from the Moody et al. posterior, truncated to be at least `margin` Gyr younger than the start of the habitable window (rejection sampling, i.e. the joint constraint that life cannot predate habitability). Returns the ages and the fraction of first-try draws that had to be redrawn. (The first pass clipped such draws to exactly `margin` after the start, which piles probability on a 10-Myr lag; the referee showed that this inflated the 'late window start' Bayes factor.)""" size = len(start_age) luca, _ = draw_transition(rng, "life", size) bad = luca > start_age - margin frac = float(bad.mean()) for _ in range(200): if not bad.any(): break new, _ = draw_transition(rng, "life", int(bad.sum())) luca[bad] = new bad = luca > start_age - margin luca = np.minimum(luca, start_age - margin) # safety net; unreachable in practice return luca, frac def loguniform(rng, lo, hi, size): if lo == hi: return np.full(size, lo) return np.exp(rng.uniform(np.log(lo), np.log(hi), size=size)) def permanent_crossing_age(curve_percent, theta_pal): """curve_percent: (n_age, n_draw) O2 in per cent with row index = age in Ma (0 = present). Returns, per draw, the oldest age (Ma) of the final uninterrupted interval with O2 >= theta (np.inf if the curve never drops below theta inside the record).""" thr = 100.0 * pal_to_mr(theta_pal) # per cent, per draw below = curve_percent < thr[None, :] first_below = np.argmax(below, axis=0).astype(float) # first age index at which it is below never = ~below.any(axis=0) age = first_below - 1.0 # last age still above age[age < 0] = 0.0 age[never] = np.inf return age def physio_gates(rng, size, gcso, midprot="all", animal="animal_sponge", human="human", euk_first=False): """Gate ages (Ga) for eukaryotes, animals and humans from thresholds fixed in advance and the oxygen history ensemble. Returns dict of arrays plus the sampled levels.""" g_on = draw_env(rng, "goe_onset", size) g_perm = draw_env(rng, "goe_permanent", size) t_neo = draw_env(rng, "tonian_onset", size) L_mp = loguniform(rng, MIDPROT[midprot][0], MIDPROT[midprot][1], size) L_neo = np.maximum(loguniform(rng, NEOPROT[0], NEOPROT[1], size), L_mp) cols = rng.integers(0, gcso.shape[1], size=size) curve = gcso[:, cols] def gate_for(theta, allow_first=False): ph = permanent_crossing_age(curve, theta) / 1000.0 # Ga, inf if always above within 550 Ma ph = np.where(np.isinf(ph), 0.55, ph) base = g_on if allow_first else g_perm return np.where(theta <= L_mp, base, np.where(theta <= L_neo, t_neo, ph)) th_e = loguniform(rng, *THRESH["euk"][:2], size) th_a = loguniform(rng, *THRESH[animal][:2], size) th_h = loguniform(rng, *THRESH[human][:2], size) return dict(eukaryotes=gate_for(th_e, euk_first), mitochondria=gate_for(th_e, euk_first), animals=gate_for(th_a), humans=gate_for(th_h), L_mp=L_mp, L_neo=L_neo, th_a=th_a, th_e=th_e, th_h=th_h) def write_tables(): rows = [] for k, v in TRANSITIONS.items(): for (lab, kind, mode, lo, hi, src, ver) in v["estimates"]: rows.append(dict(table="transition", key=k, name=v["name"], estimate=lab, distribution="uniform" if kind == "u" else "two-piece normal, 95% interval", centre_Ga=mode, lo_Ga=lo, hi_Ga=hi, source=src, verified=ver)) for k, (lo, hi, src, ver) in ENV.items(): rows.append(dict(table="environment", key=k, name=k, estimate="", distribution="uniform", centre_Ga=None, lo_Ga=lo, hi_Ga=hi, source=src, verified=ver)) for k, (lo, hi, src) in THRESH.items(): rows.append(dict(table="threshold_PAL", key=k, name=k, estimate="", distribution="log-uniform", centre_Ga=None, lo_Ga=lo, hi_Ga=hi, source=src, verified="A")) for k, (lo, hi, src) in MIDPROT.items(): rows.append(dict(table="midproterozoic_PAL", key=k, name=k, estimate="", distribution="log-uniform", centre_Ga=None, lo_Ga=lo, hi_Ga=hi, source=src, verified="A")) rows.append(dict(table="neoproterozoic_PAL", key="neo", name="neo", estimate="", distribution="log-uniform", centre_Ga=None, lo_Ga=NEOPROT[0], hi_Ga=NEOPROT[1], source=NEOPROT[2], verified="A")) df = pd.DataFrame(rows) df.to_csv(DATA + "/part1_inputs.csv", index=False) return df if __name__ == "__main__": print(write_tables().to_string())