#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ iguana-raft.py Science Journaling Club - club-built companion calculation for the field note "Eight Thousand Kilometres on a Raft of Weeds: How Iguanas Reached Fiji". ANCHOR PAPER (not ours): Scarpetta, S. G., et al. (2025). Iguanas rafted more than 8,000 km from North America to Fiji. PNAS 122(12), e2318622122. doi:10.1073/pnas.2318622122 The paper is a phylogenomics paper. It sequenced >4,000 loci and >2 million base pairs, recovered the Fijian genus Brachylophus as sister to the North American desert iguana Dipsosaurus, dated that split to the late Paleogene, and ran biogeographic model selection that puts the ancestral range in North America. None of that is reproduced here and none of it is ours. WHAT THIS SCRIPT IS ------------------- This is the club's own simplified feasibility model. It asks a question the paper does not try to answer numerically: given ordinary published numbers for Pacific surface currents and for lizard physiology, is a drifting iguana a physically coherent object, and how often would you expect one to make landfall? Four blocks: BLOCK 1 DRIFT TIME. Turn published North and South Equatorial Current speed ranges into crossing times for 8,000 km, with a path-sinuosity factor because a raft does not travel on a great circle. BLOCK 2 ENERGY. Allometric resting metabolic rate for a squamate of mass M, temperature-corrected with a Q10, run against a fat store. Gives survivable days with no food. Then the same budget with the raft itself treated as food, which is the case that matters for a herbivore. BLOCK 3 WATER. A physical evaporative-water-loss budget: cutaneous loss through a skin resistance driven by the vapour-density deficit between skin and marine air, respiratory loss tied to oxygen uptake, excretory loss tied to catabolism, metabolic water as a credit. Gives survivable days with no fresh water. Then a seeded Monte Carlo over rainy days on the raft, crossed with current speed, to find the survivable window. BLOCK 4 DEEP TIME. Rafts launched per century times survival probability per voyage gives a Poisson rate of successful colonisations. Expected waiting time, and probability of at least one success inside the window the phylogeny allows. ASSUMPTIONS (all of them, so you can break them) ------------------------------------------------ A1 Great-circle distance North America to Fiji is taken as 8,000 km, the paper's own figure. Real drift paths are longer; handled by SINUOSITY. A2 A raft moves at the speed of the surface current it sits in, with a small windage credit folded into the "fast" case. Rafts with sail-like canopies move faster than the water; rafts that are mostly waterlogged move slower. A3 Body mass defaults to 70 g, a mid-sized adult Dipsosaurus. The actual Paleogene disperser is unknown and could be several times heavier. A4 Resting metabolic rate is 40% of the standard reptile field metabolic rate allometry, because a fasting animal on a raft is not foraging or territorial. Q10 = 2.5 away from the 30 C reference. A5 Fat store is 10% of body mass, 80% of it mobilisable, 39.3 kJ per gram. A6 Total body water is 70% of body mass and death comes at a 25% loss of that water. Published lethal dehydration limits for xeric lizards straddle this. A7 Skin resistance to water vapour is 300 s/cm, at the permeable end of the range measured for desert lizards, so the water budget is pessimistic. A8 Body temperature on the raft sits 1.5 C above air temperature. A raft is wet, often overcast, and a lizard on it cannot thermoregulate the way it can in a desert. A9 Rain is an independent Bernoulli draw each day. Real tropical rainfall is clustered, which makes long dry runs MORE likely than this model says. A10 A drinking bout refills at most 12% of body mass in water. A11 In Block 4, per-voyage survival is a product of four independent probabilities. They are not independent in reality. A12 The raft's vegetation supplies a fraction phi(t) = exp(-t/TAU_LEAF) of the animal's daily energy requirement, because foliage dies, is eaten, and becomes poorer food as the voyage goes on. TAU_LEAF is a free parameter and it is the single softest number in the whole script. Block 2b solves for the value at which the animal exactly starves, rather than assuming one and hoping. LIMITS (things this model simply does not contain) -------------------------------------------------- L1 Raft integrity. A vegetation raft waterlogs, is eaten by shipworm, and comes apart. This is very probably the dominant term and it enters only as a tunable probability in Block 4, never as physics. L2 Storms, waves washing the animal off, and predation from above and below. L3 Salt. Sea spray on skin and in drinking water; iguanid nasal salt glands help, and are not modelled. L4 Thermal extremes. Day-night swings on an open raft are averaged away. L5 Landfall geometry. Whether a drifting object actually intersects Fiji is an ocean-circulation question this script does not attempt. L6 Everything about the founding population: one gravid female, two animals, inbreeding, and the demography of a colonisation that nearly fails. Every number printed by this script is reproducible: python iguana-raft.py Seeds are fixed at the top of each stochastic block. """ from __future__ import annotations import math import random try: import numpy as np HAVE_NUMPY = True except Exception: # pragma: no cover np = None HAVE_NUMPY = False MASTER_SEED = 20250317 # PNAS publication date of the anchor paper random.seed(MASTER_SEED) if HAVE_NUMPY: np.random.seed(MASTER_SEED) W = 78 def rule(ch="-"): print(ch * W) def head(title): print() rule("=") print(title) rule("=") def sub(title): print() print(title) rule("-") # ============================================================================ # BLOCK 1. DRIFT TIME ACROSS 8,000 KM # ============================================================================ DISTANCE_KM = 8000.0 SINUOSITY = 1.35 # A1: path length / great-circle distance # Published surface-current speed ranges for the tropical Pacific zonal system. # The North Equatorial Current and the South Equatorial Current are the two # westward conveyors a raft leaving western North America can end up in. CURRENTS = [ ("slow", 0.15, "weak NEC limb, near the northern flank"), ("median", 0.35, "typical NEC / SEC surface speed"), ("fast", 0.70, "SEC core, plus windage on a canopied raft"), ] def drift_days(speed_ms, distance_km=DISTANCE_KM, sinuosity=SINUOSITY): """Days to cover distance_km at speed_ms, inflated by path sinuosity.""" metres = distance_km * 1000.0 * sinuosity return metres / (speed_ms * 86400.0) def block1(): head("BLOCK 1 DRIFT TIME: HOW LONG IS THE VOYAGE?") print(f"Great-circle distance assumed : {DISTANCE_KM:,.0f} km") print(f"Path sinuosity factor (A1) : {SINUOSITY:.2f}") print(f"Actual water distance travelled : {DISTANCE_KM*SINUOSITY:,.0f} km") sub("Crossing time by current speed") print(f"{'case':<8}{'m/s':>7}{'km/day':>10}{'straight d':>13}" f"{'with detour d':>15}{' months':>9} note") out = {} for name, v, note in CURRENTS: straight = DISTANCE_KM * 1000.0 / (v * 86400.0) d = drift_days(v) out[name] = d print(f"{name:<8}{v:>7.2f}{v*86.4:>10.1f}{straight:>13.1f}" f"{d:>15.1f}{d/30.44:>9.1f} {note}") sub("Comparison against a published drift estimate (not a validation)") ref_lo, ref_hi = 80.0, 120.0 ratio = out['fast'] / ref_hi print("Montenegro et al. (2008) model sweet-potato drift trajectories across") print(f"the tropical Pacific and report transit times of {ref_lo:.0f}-{ref_hi:.0f} days for") print("comparable distances in their fastest current. That is a different") print("object (a tuber, not a raft carrying a lizard) on a different route,") print("so this is a plausibility comparison, not an independent check of") print("this model: nothing here is computed from a second, unrelated method.") print(f"Our fast case lands at {out['fast']:.0f} days, which is " f"{ratio:.2f}x the top of that range, not inside it.") print("We report the gap honestly rather than rounding it away. Plausible") print("reasons a raft could run slower than a semi-submerged tuber include") print("more wind and wave drag, a longer effective path than the great") print("circle sinuosity factor assumes, and current speed varying along") print("the route rather than sitting at the single fastest value used here.") return out # ============================================================================ # BLOCK 2. ENERGY BUDGET # ============================================================================ MASS_G = 70.0 # A3 FAT_FRACTION = 0.10 # A5 FAT_USABLE = 0.80 # A5 FAT_KJ_PER_G = 39.3 # A5 FMR_A, FMR_B = 0.196, 0.889 # reptile field metabolic rate, kJ/day, mass in g RESTING_FRACTION = 0.40 # A4 Q10 = 2.5 # A4 T_REF_C = 30.0 # A4 T_AIR_C = 27.0 # tropical Pacific marine boundary layer along the route RH = 0.80 # A7 companion: marine relative humidity TB_OFFSET = 1.5 # A8 def body_temp(t_air=T_AIR_C): return t_air + TB_OFFSET def rmr_kj_per_day(mass_g=MASS_G, t_body=None): """Resting metabolic rate, kJ/day, from the reptile FMR allometry (A4).""" if t_body is None: t_body = body_temp() fmr = FMR_A * mass_g ** FMR_B rest = RESTING_FRACTION * fmr return rest * Q10 ** ((t_body - T_REF_C) / 10.0) def fat_store_kj(mass_g=MASS_G): return mass_g * FAT_FRACTION * FAT_USABLE * FAT_KJ_PER_G def fasting_days(mass_g=MASS_G, t_body=None): return fat_store_kj(mass_g) / rmr_kj_per_day(mass_g, t_body) def block2(drift): head("BLOCK 2 ENERGY: HOW LONG CAN IT GO WITHOUT EATING?") tb = body_temp() fmr = FMR_A * MASS_G ** FMR_B rest = rmr_kj_per_day() store = fat_store_kj() print(f"Body mass (A3) : {MASS_G:.0f} g") print(f"Air temperature along route : {T_AIR_C:.1f} C") print(f"Body temperature on raft (A8) : {tb:.1f} C") print(f"Reptile FMR allometry {FMR_A}*M^{FMR_B} : {fmr:.2f} kJ/day at 30 C") print(f"Resting fraction (A4) : {RESTING_FRACTION:.2f}") print(f"Q10 correction 30 C -> {tb:.1f} C : " f"{Q10 ** ((tb - T_REF_C) / 10.0):.4f}") print(f"Resting metabolic rate on the raft : {rest:.3f} kJ/day") print() print(f"Fat store {FAT_FRACTION:.0%} of mass, {FAT_USABLE:.0%} usable, " f"{FAT_KJ_PER_G} kJ/g") print(f"Usable energy reserve : {store:.1f} kJ") print(f"SURVIVABLE DAYS WITH NO FOOD : {store/rest:.1f} days") print(f"Fat burned per day : {rest/FAT_KJ_PER_G:.4f} g/day") sub("Fasting endurance against the three crossing times") fd = store / rest print(f"{'case':<8}{'voyage d':>11}{'fast-only d':>14}{'shortfall d':>14}" f"{' verdict':>12}") for name, _, _ in CURRENTS: d = drift[name] short = d - fd verdict = "STARVES" if short > 0 else "survives" print(f"{name:<8}{d:>11.1f}{fd:>14.1f}{short:>14.1f}{verdict:>12}") sub("The raft as food: a herbivore is eating its own vehicle") LEAF_KJ_PER_G_DRY = 7.0 # gross energy of foliage, dry mass ASSIM = 0.50 # digestible fraction for a hindgut-fermenting lizard for name, _, _ in CURRENTS: d = drift[name] need_kj = rest * d need_leaf = need_kj / (LEAF_KJ_PER_G_DRY * ASSIM) print(f"{name:<8} {d:>7.1f} d needs {need_kj:>8.1f} kJ = " f"{need_leaf:>7.1f} g dry foliage " f"({need_leaf/MASS_G:>5.2f} x its own body mass)") print() print("A single uprooted tree carries kilograms of leaf. Energy is not the") print("wall, provided the vegetation is alive and the animal is a herbivore.") print("A carnivore on the same raft has nothing to eat but the raft's other") print("passengers, and runs out.") sub("Block 2b: how long does the raft have to keep feeding it?") print("Mass of foliage is not the constraint. Freshness is. Model the raft as") print("supplying phi(t) = exp(-t/TAU) of the daily requirement (A12), so the") print("animal draws burn*(1 - phi) from fat. Solve for the TAU at which fat") print("runs out exactly at landfall.") print() print(f"{'case':>8}{'voyage d':>11}{'critical TAU d':>17}" f"{' meaning':>10}") crit = {} for name, _, _ in CURRENTS: d = drift[name] lo, hi = 1.0, 1e6 for _ in range(200): mid = 0.5 * (lo + hi) used = rest * (d - mid * (1.0 - math.exp(-d / mid))) if used > store: lo = mid # raft spoils too fast, animal starves else: hi = mid crit[name] = 0.5 * (lo + hi) print(f"{name:>8}{d:>11.1f}{crit[name]:>17.1f}" f" raft must stay edible on this timescale") print() print("Read that carefully, because it is the most awkward number we found.") print("Even in the fast current the raft has to go on being worth eating on a") print(f"{crit['fast']:.0f}-day timescale, most of the length of the crossing. On the slow") print(f"current it has to be {crit['slow']:.0f} days, which is not a raft, it is a") print("small floating island with a working food web. The vegetation cannot be") print("ballast. It has to be alive, or dead in a way a hindgut fermenter can") print("still use.") return fd, crit # ============================================================================ # BLOCK 3. WATER BUDGET # ============================================================================ TBW_FRACTION = 0.70 # A6 LETHAL_LOSS = 0.25 # A6 SKIN_RESISTANCE = 300.0 # s/cm, A7 O2_EXTRACTION = 0.15 # fraction of inspired O2 taken up KJ_PER_L_O2 = 19.8 # oxycalorific equivalent for lipid catabolism EXCRETORY_G_PER_KJ = 0.045 METABOLIC_WATER_G_PER_G_FAT = 1.07 DRINK_CAP_FRACTION = 0.12 # A10 def sat_vapour_density(t_c): """Saturation water-vapour density in g per cubic metre (Magnus + ideal gas).""" es = 610.94 * math.exp(17.625 * t_c / (t_c + 243.04)) # Pa return es * 0.018016 / (8.3145 * (t_c + 273.15)) * 1000.0 def surface_area_cm2(mass_g=MASS_G): """Standard lizard surface-area allometry, cm^2.""" return 10.0 * mass_g ** (2.0 / 3.0) def water_budget(mass_g=MASS_G, t_air=T_AIR_C, rh=RH): """Return a dict of water fluxes in grams per day, plus the survivable days.""" tb = t_air + TB_OFFSET rho_skin = sat_vapour_density(tb) rho_air = rh * sat_vapour_density(t_air) deficit = max(rho_skin - rho_air, 0.0) # g/m^3 area = surface_area_cm2(mass_g) # cm^2 # deficit g/m^3 -> g/cm^3 is 1e-6; area cm^2 / resistance s/cm -> cm^3/s cut = area * (deficit * 1e-6) / SKIN_RESISTANCE * 86400.0 rest = rmr_kj_per_day(mass_g, tb) vo2_l = rest / KJ_PER_L_O2 air_l = vo2_l / (0.209 * O2_EXTRACTION) resp = air_l * (deficit * 1e-3) # g/L * L excr = EXCRETORY_G_PER_KJ * rest fat_g = rest / FAT_KJ_PER_G metab = METABOLIC_WATER_G_PER_G_FAT * fat_g net = cut + resp + excr - metab tbw = TBW_FRACTION * mass_g tolerable = LETHAL_LOSS * tbw days = tolerable / net if net > 0 else float("inf") return { "t_air": t_air, "t_body": tb, "rho_skin": rho_skin, "rho_air": rho_air, "deficit": deficit, "area": area, "cut": cut, "resp": resp, "excr": excr, "metab": metab, "net": net, "tbw": tbw, "tolerable": tolerable, "days": days, "rmr": rest, } def block3a(): head("BLOCK 3 WATER: HOW LONG WITHOUT A DRINK?") b = water_budget() print(f"Air temperature : {b['t_air']:.1f} C") print(f"Relative humidity : {RH:.0%}") print(f"Skin temperature : {b['t_body']:.1f} C") print(f"Saturation density at skin : {b['rho_skin']:.2f} g/m3") print(f"Vapour density of marine air : {b['rho_air']:.2f} g/m3") print(f"Vapour density deficit : {b['deficit']:.2f} g/m3") print(f"Skin area 10*M^(2/3) : {b['area']:.1f} cm2") print(f"Skin resistance (A7) : {SKIN_RESISTANCE:.0f} s/cm") sub("Daily water ledger, grams") print(f" cutaneous evaporation -{b['cut']:.4f}") print(f" respiratory evaporation -{b['resp']:.4f}") print(f" urine and faeces -{b['excr']:.4f}") print(f" metabolic water (credit) +{b['metab']:.4f}") rule(".") print(f" NET LOSS -{b['net']:.4f} g/day") print() print(f"Total body water {TBW_FRACTION:.0%} of mass : {b['tbw']:.2f} g") print(f"Tolerable deficit at {LETHAL_LOSS:.0%} loss : {b['tolerable']:.2f} g") print(f"SURVIVABLE DAYS WITH NO WATER : {b['days']:.1f} days") print() print("Cutaneous loss is " f"{100*b['cut']/(b['cut']+b['resp']+b['excr']):.0f}% of gross loss. The animal") print("is leaking through its skin, slowly, and nothing else is close.") return b def block3b(): sub("Sensitivity: the water clock against air temperature") print(f"{'T air C':>9}{'T skin C':>10}{'deficit g/m3':>14}" f"{'net g/day':>12}{'days dry':>11}") rows = [] for t in (24, 26, 27, 28, 30, 32, 34, 36): b = water_budget(t_air=t) rows.append((t, b)) print(f"{t:>9.0f}{b['t_body']:>10.1f}{b['deficit']:>14.2f}" f"{b['net']:>12.4f}{b['days']:>11.1f}") print() print("Warm air is what kills, and it kills through the vapour deficit rather") print("than through metabolism. Between 27 C and 36 C the dry-survival window") print(f"falls from {rows[2][1]['days']:.1f} days to {rows[-1][1]['days']:.1f} days.") return rows def block3c(): sub("Sensitivity: body mass") print(f"{'mass g':>8}{'RMR kJ/d':>11}{'fat kJ':>9}{'fast d':>9}" f"{'net H2O g/d':>13}{'dry d':>8}") rows = [] for m in (20, 40, 70, 120, 200, 350): b = water_budget(mass_g=m) fd = fasting_days(m) rows.append((m, fd, b["days"])) print(f"{m:>8.0f}{rmr_kj_per_day(m):>11.3f}{fat_store_kj(m):>9.1f}" f"{fd:>9.1f}{b['net']:>13.4f}{b['days']:>8.1f}") print() print("Bigger is better on both clocks, because reserves scale with mass and") print("losses scale with surface area. A 350 g animal is a better boat than a") print("20 g one by roughly a factor of two on water and three on fuel.") return rows # ---- Monte Carlo: rain on the raft -------------------------------------- MC_SEED = 611 MC_TRIALS = 20000 def simulate_voyage(rng, days_total, p_rain, mass_g=MASS_G, t_air=T_AIR_C, raft_food=True): """One voyage. Returns (survived, day_of_death_or_arrival, cause).""" b = water_budget(mass_g=mass_g, t_air=t_air) net = b["net"] tolerable = b["tolerable"] drink_cap = DRINK_CAP_FRACTION * mass_g energy = fat_store_kj(mass_g) burn = rmr_kj_per_day(mass_g, t_air + TB_OFFSET) deficit = 0.0 n = int(math.ceil(days_total)) for day in range(1, n + 1): deficit += net if rng.random() < p_rain: deficit = max(0.0, deficit - drink_cap) if deficit >= tolerable: return False, day, "dehydration" if not raft_food: energy -= burn if energy <= 0.0: return False, day, "starvation" return True, n, "landfall" def block3d(drift): sub("Monte Carlo: does rain on the raft rescue the voyage?") print(f"seed = {MC_SEED}, trials per cell = {MC_TRIALS}, " f"raft treated as food = yes") print("Rain is an independent daily Bernoulli draw (A9). A drinking bout") print(f"restores at most {DRINK_CAP_FRACTION:.0%} of body mass in water (A10).") print() p_grid = [0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.40, 0.50] header = f"{'p(rain)':>8}" + "".join(f"{name:>13}" for name, _, _ in CURRENTS) print(header) print(f"{'':>8}" + "".join(f"{drift[n]:>10.0f} d " for n, _, _ in CURRENTS)) rule(".") table = {} for p in p_grid: row = [] for name, _, _ in CURRENTS: rng = random.Random(MC_SEED + int(p * 1000) * 7 + len(name)) wins = 0 for _ in range(MC_TRIALS): ok, _d, _c = simulate_voyage(rng, drift[name], p) wins += ok row.append(wins / MC_TRIALS) table[p] = row print(f"{p:>8.2f}" + "".join(f"{v:>13.4f}" for v in row)) print() print("Read the zero row first. With no rain at all, every crossing kills the") print("animal, at every current speed, without exception. The voyage is only") print("possible because the route runs through the wettest air on the planet.") return p_grid, table def block3e(drift): sub("Where the survivable window actually is") print("Minimum rainfall frequency for a 50% and a 90% chance of arrival,") print("by current speed, at 27 C.") print(f"{'case':>8}{'voyage d':>11}{'p for 50%':>12}{'p for 90%':>12}" f"{'rain days/yr at p50':>22}") out = {} for name, _, _ in CURRENTS: d = drift[name] p50 = p90 = None p = 0.0 while p <= 0.60001: rng = random.Random(MC_SEED + 991 + int(p * 1000)) wins = 0 trials = 4000 for _ in range(trials): ok, _d, _c = simulate_voyage(rng, d, p) wins += ok f = wins / trials if p50 is None and f >= 0.50: p50 = p if p90 is None and f >= 0.90: p90 = p break p += 0.01 out[name] = (p50, p90) s50 = f"{p50:.2f}" if p50 is not None else " > .60" s90 = f"{p90:.2f}" if p90 is not None else " > .60" rd = f"{p50*365:.0f}" if p50 is not None else "-" print(f"{name:>8}{d:>11.0f}{s50:>12}{s90:>12}{rd:>22}") print() print("For scale: the Intertropical Convergence Zone and the South Pacific") print("Convergence Zone, which sit across this route, deliver rain on a large") print("fraction of days year-round. A raft crossing the tropical Pacific is") print("not crossing a desert. It is crossing a rainforest with no floor.") return out def block3f(drift): sub("Deterministic voyage log at the median current (no rain, worst case)") d = drift["median"] b = water_budget() burn = rmr_kj_per_day() store = fat_store_kj() print(f"{'day':>6}{'km':>9}{'route %':>10}{'water left %':>14}" f"{'fat left %':>12} state") marks = [0, 30, 60, 90, 120, 180, 240, 300, int(round(d))] log = [] for day in marks: km = DISTANCE_KM * min(day / d, 1.0) water_left = max(0.0, 1.0 - (b["net"] * day) / b["tolerable"]) fat_left = max(0.0, 1.0 - (burn * day) / store) if water_left <= 0: state = "dead of thirst" elif fat_left <= 0: state = "dead of hunger" elif day >= d: state = "landfall" else: state = "adrift" log.append((day, km, 100 * min(day / d, 1.0), 100 * water_left, 100 * fat_left, state)) print(f"{day:>6}{km:>9.0f}{100*min(day/d,1.0):>10.1f}" f"{100*water_left:>14.1f}{100*fat_left:>12.1f} {state}") print() print("This is the pessimistic log: no rain, nothing edible. The water clock") print("stops the voyage before the first tenth of the route is behind it.") return log TAU_LEAF = 900.0 # A12, generous; Block 2b shows the critical value def block3g(drift): sub("Voyage log: one seeded crossing, median current, 15% of days rainy") print(f"TAU_LEAF = {TAU_LEAF:.0f} d (A12). Seed for this single voyage = 7711.") print("Water reserve is what remains of the tolerable dehydration deficit.") print("Energy reserve is what remains of the usable fat store.") print() d = drift["median"] b = water_budget() burn = rmr_kj_per_day() store = fat_store_kj() drink_cap = DRINK_CAP_FRACTION * MASS_G rng = random.Random(7711) n = int(math.ceil(d)) deficit = 0.0 energy = store series = [] dry_run = 0 worst_dry = 0 rain_days = 0 for day in range(1, n + 1): deficit += b["net"] if rng.random() < 0.15: deficit = max(0.0, deficit - drink_cap) rain_days += 1 dry_run = 0 else: dry_run += 1 worst_dry = max(worst_dry, dry_run) phi = math.exp(-day / TAU_LEAF) energy -= burn * (1.0 - phi) alive = deficit < b["tolerable"] and energy > 0.0 series.append((day, 100 * max(0.0, 1 - deficit / b["tolerable"]), 100 * max(0.0, energy / store), alive)) if not alive: break print(f"{'day':>6}{'km':>9}{'route %':>10}{'water %':>10}{'energy %':>11}" f" log entry") marks = [1, 36, 72, 108, 143, 179, 215, 250, 286, 322, n] for day in marks: if day > len(series): break _, w, e, alive = series[day - 1] frac = min(day / d, 1.0) km = DISTANCE_KM * frac note = "landfall" if day >= n else ("adrift" if alive else "dead") print(f"{day:>6}{km:>9.0f}{100*frac:>10.1f}{w:>10.1f}{e:>11.1f}" f" {note}") print() print(f"rainy days on this crossing : {rain_days} of {len(series)} " f"({100*rain_days/len(series):.1f}%)") print(f"longest dry run : {worst_dry} days " f"(lethal at {b['days']:.1f})") print(f"water reserve at landfall : {series[-1][1]:.1f}%") print(f"energy reserve at landfall : {series[-1][2]:.1f}%") print(f"outcome : " f"{'ARRIVED' if series[-1][3] else 'DIED'}") return series # ============================================================================ # BLOCK 4. DEEP TIME LOTTERY # ============================================================================ WINDOW_MYR = 34.0 # paper's median Brachylophus-Dipsosaurus split WINDOW_MAX_MYR = 42.0 # upper reach of the credible intervals reported def waiting_time(launches_per_century, p_success): """Expected centuries and Myr between successful colonisations.""" rate = launches_per_century * p_success # successes per century if rate <= 0: return float("inf"), float("inf") return 1.0 / rate, 1.0 / rate / 10000.0 def block4(mc_survival): head("BLOCK 4 DEEP TIME: NEAR-IMPOSSIBLE, DIVIDED BY THIRTY-FOUR MILLION YEARS") print("Per-voyage success is a product of four terms (A11).") print() P_RAFT = 0.02 P_ANIMAL = mc_survival P_LANDFALL = 0.01 P_PROPAGULE = 0.05 q = P_RAFT * P_ANIMAL * P_LANDFALL * P_PROPAGULE print(f" raft still afloat and vegetated at the far end {P_RAFT:>10.4f}") print(f" animal survives water and energy budget {P_ANIMAL:>10.4f}") print(f" raft actually intersects Fiji rather than ocean {P_LANDFALL:>10.4f}") print(f" founding propagule viable (gravid female, or 2+) {P_PROPAGULE:>10.4f}") rule(".") print(f" PER-VOYAGE SUCCESS PROBABILITY q {q:>10.3e}") print(f" that is one voyage in {1/q:,.0f}") sub("Expected waiting time for one successful colonisation") print(f"{'launches/century':>18}{'successes/century':>20}" f"{'wait, centuries':>18}{'wait, Myr':>13}") waits = {} for lam in (0.1, 1.0, 5.0, 10.0, 50.0, 100.0): c, myr = waiting_time(lam, q) waits[lam] = myr print(f"{lam:>18.1f}{lam*q:>20.3e}{c:>18,.0f}{myr:>13.2f}") sub(f"Probability of at least one arrival inside the phylogenetic window") print(f"Window taken as {WINDOW_MYR:.0f} Myr (median split) and " f"{WINDOW_MAX_MYR:.0f} Myr (upper interval).") print(f"{'launches/century':>18}{'expected arrivals':>20}" f"{'P(>=1) 34 Myr':>16}{'P(>=1) 42 Myr':>16}") for lam in (0.1, 1.0, 5.0, 10.0, 50.0, 100.0): n34 = lam * q * WINDOW_MYR * 1e6 / 100.0 n42 = lam * q * WINDOW_MAX_MYR * 1e6 / 100.0 print(f"{lam:>18.1f}{n34:>20.3f}{1-math.exp(-n34):>16.4f}" f"{1-math.exp(-n42):>16.4f}") sub("The same fact said three ways") lam = 10.0 total = lam * WINDOW_MYR * 1e6 / 100.0 print(f"At {lam:.0f} launches per century over {WINDOW_MYR:.0f} Myr, " f"{total:,.0f} rafts set out.") print(f"At q = {q:.3e}, {total*q:.2f} of them are expected to succeed.") print(f"Probability at least one does: {1-math.exp(-total*q):.4f}") print() print(f"A voyage that fails {100*(1-q):.5f}% of the time is a near-certainty if") print("you run it three million times. Deep time is not a long wait. It is a") print("large sample.") sub("How hard is the conclusion? Sweep q over three orders of magnitude") print(f"{'q':>12}{'one in':>14}{'arrivals in 34 Myr':>22}{'P(>=1)':>10}" f" (10 launches/century)") for mult in (0.01, 0.1, 0.5, 1.0, 2.0, 10.0, 100.0): qq = q * mult n = 10.0 * qq * WINDOW_MYR * 1e6 / 100.0 print(f"{qq:>12.3e}{1/qq:>14,.0f}{n:>22.3f}{1-math.exp(-n):>10.4f}") print() print("A ten-fold pessimism about q still leaves an 85% chance of at least one") print("arrival. A hundred-fold pessimism drops it to 17%, and that is where the") print("conclusion finally becomes fragile. Everything hangs on q to within") print("about two orders of magnitude, and q is a guess.") return q, waits # ============================================================================ def main(): print("=" * W) print("IGUANAS ON A RAFT: THE CLUB'S OWN FEASIBILITY CALCULATION") print("Companion to Scarpetta et al. (2025) PNAS 122(12) e2318622122") print(f"master seed = {MASTER_SEED} numpy = {HAVE_NUMPY}") print("Everything below is the club's simplified model, not the paper's") print("analysis. The paper did the phylogenomics. We did the arithmetic.") print("=" * W) drift = block1() fast_days, crit_tau = block2(drift) block3a() block3b() block3c() p_grid, table = block3d(drift) block3e(drift) block3f(drift) vlog = block3g(drift) # survival probability used in Block 4: median current, 15% of days rainy, # a deliberately conservative reading of tropical Pacific rainfall. rng = random.Random(MC_SEED + 4242) wins = 0 N = 40000 for _ in range(N): ok, _d, _c = simulate_voyage(rng, drift["median"], 0.15) wins += ok p_animal = wins / N print() print(f"[Block 4 input] survival at median current ({drift['median']:.0f} d), " f"p(rain)=0.15, {N} trials: {p_animal:.4f}") block4(p_animal) head("SUMMARY OF THE NUMBERS THE ARTICLE QUOTES") b = water_budget() print(f"voyage, slow current {drift['slow']:.0f} days") print(f"voyage, median current {drift['median']:.0f} days") print(f"voyage, fast current {drift['fast']:.0f} days") print(f"fasting endurance, 70 g, 28.5 C {fast_days:.1f} days") print(f"dry endurance, 70 g, 27 C {b['days']:.1f} days") print(f"net water loss {b['net']:.4f} g/day") print(f"resting metabolic rate {rmr_kj_per_day():.3f} kJ/day") print(f"usable fat store {fat_store_kj():.1f} kJ") print(f"survival, median + 15% rain {p_animal:.4f}") print(f"critical raft-food TAU, median {crit_tau['median']:.0f} days") print(f"critical raft-food TAU, fast {crit_tau['fast']:.0f} days") print(f"voyage-log landfall water/energy {vlog[-1][1]:.1f}% / {vlog[-1][2]:.1f}%") print() print("Reproduce with: python analysis/iguana-raft.py") if __name__ == "__main__": main()