""" roman-lead.py Science Journaling Club -- a simplified reconstruction of the Roman-era lead dose-response chain reported by McConnell et al. (2025), PNAS 122(3), e2419630121, doi:10.1073/pnas.2419630121. WHAT THIS IS ------------ The paper runs Arctic ice-core lead records backwards through a Lagrangian particle dispersion model (FLEXPART) to recover annual European lead emissions, forwards again to map air concentrations over Europe at roughly 20 km resolution, then through two published epidemiological regressions to reach childhood blood lead levels and an IQ decrement. This script is NOT that. It is the Science Journaling Club's own, deliberately crude reconstruction of the same causal chain, built so that a reader can follow every arithmetic step by hand. We have no ice cores, no dispersion model and no cohort data. We take the paper's published annual emission estimate as our INPUT, replace the whole atmospheric model with a single well-mixed box, and then apply the two published regressions exactly as the paper states them. Everything between emissions and air concentration is ours; everything after that is the published epidemiology. THE MODEL, IN FOUR STAGES ------------------------- Stage 1 emissions -> air concentration. A steady-state deposition-limited box over Europe. At steady state everything emitted is removed by deposition inside the box, so E = v_d * C * A => C = E / (v_d * A) with E the annual emission (ug/yr), v_d an effective bulk removal velocity (dry deposition plus wet scavenging, m/yr), A the area of the European domain (m^2), and C the annual-mean near-surface air lead concentration (ug/m^3). No mixing height appears, because at steady state the box's depth cancels out. Stage 2 air concentration -> childhood blood lead. The log-linear air-to-blood regression for 1-5 year olds quoted by McConnell et al. (2025), itself fitted to modern paired air-lead / blood-lead data: dBLL = exp[a + b * ln(C)], a = 1.932 +/- 0.181, b = 0.140 +/- 0.054, C in ug/m^3, dBLL in ug/dL dBLL is the ENHANCEMENT above a non-air baseline. Total blood lead is dBLL plus a Neolithic baseline of about 1.0 ug/dL inferred from tooth enamel. Stage 3 blood lead -> IQ. The log-linear regression the paper applies, derived from the Lanphear et al. (2005) international pooled analysis: IQ deficit = -3.315 (+/- 0.6155) * ln(BLL + 1) The Roman-attributable decrement is the deficit at the total Roman blood lead minus the deficit at the Neolithic baseline, so that the pre-existing background burden is not counted twice. Stage 4 bookkeeping. Sensitivity of the final IQ number to each assumption, one at a time, plus a seeded Monte Carlo over all of them together, plus comparison against modern reference values. ASSUMPTIONS, STATED ------------------- A1 Europe is one well-mixed box of area 1.44e13 m^2 (10W to 40E, 30N to 65N on a sphere of radius 6371 km). Real air lead spanned 0.16 to 157 ng/m^3 across that domain; a box model produces one number and cannot see a smelter. A2 Effective bulk removal velocity 1.5 cm/s, carried over 0.5 to 3.0 cm/s. This lumps dry deposition of fine and coarse smelter aerosol together with wet scavenging, which is the larger term for submicron particles. A3 Steady state. Emissions are treated as constant over the Pax Romana, which is roughly true on the century timescale the ice cores resolve and plainly false year to year. A4 No advective export out of the box. Lead leaving Europe on the westerlies is, in this model, lead that deposited inside Europe. That biases C high. A5 The air-to-blood regression, fitted on modern children breathing air lead of order 0.01 to 1 ug/m^3, is extrapolated down to about 5e-4 ug/m^3. This is the single largest liberty taken in the chain. A6 Population-average exposure. One number for an empire of roughly 60 to 80 million people, farmers and smelter workers alike. LIMITS ------ * A box model cannot reproduce the paper's spatial map, which is most of what the paper is for. Our C is an area average and nothing else. * The two regressions are observational fits with wide confidence bands, applied to a population for which no blood samples exist and never will. * Roman lead exposure also came from pipes, cookware, sapa-sweetened wine, cosmetics and contaminated soil. The chain modelled here is the airborne route alone, which is why the paper calls its numbers minimum estimates. * IQ is a twentieth-century instrument. "2.6 IQ points in the Roman Empire" means "the decrement a modern cohort study would predict from that blood lead", not a score anybody ever sat down and took. Python 3.12. numpy required. Every random number is seeded. """ import math import random import numpy as np SEED = 20250106 # the paper's publication date, 6 January 2025 random.seed(SEED) RNG = np.random.default_rng(SEED) SEC_PER_YEAR = 365.25 * 24 * 3600 # 3.15576e7 s # ---------------------------------------------------------------- constants # Stage 1: the box R_EARTH_KM = 6371.0 LON_W, LON_E = -10.0, 40.0 LAT_S, LAT_N = 30.0, 65.0 def domain_area_m2(lon_w, lon_e, lat_s, lat_n, radius_km=R_EARTH_KM): """Area of a lat-lon rectangle on a sphere, in square metres.""" dlam = math.radians(lon_e - lon_w) band = math.sin(math.radians(lat_n)) - math.sin(math.radians(lat_s)) return (radius_km * 1000.0) ** 2 * dlam * band A_EUROPE = domain_area_m2(LON_W, LON_E, LAT_S, LAT_N) # m^2 VD_CM_S = 1.5 # effective removal velocity E_CENTRAL_KT = 3.5 # kt Pb/yr, midpoint of 3.0-4.2 E_SCEN1_KT = 4.2 # paper scenario 1, Rio Tinto dominant E_SCEN2_KT = 3.0 # paper scenario 2, distributed sources # Stage 2: air -> blood, 1-5 year olds (McConnell et al. 2025) BLL_A, BLL_A_SD = 1.932, 0.181 BLL_B, BLL_B_SD = 0.140, 0.054 BASELINE_BLL = 1.0 # ug/dL, Neolithic tooth enamel BASELINE_BLL_LO, BASELINE_BLL_HI = 0.5, 1.5 # Stage 3: blood -> IQ (Lanphear et al. 2005 pooled analysis, as applied) IQ_SLOPE, IQ_SLOPE_SD = -3.315, 0.6155 # Reference blood lead levels for comparison REF_LEVELS = [ ("Neolithic Europe (tooth enamel)", 1.0, "pre-metallurgy baseline"), ("US children 2011-2016 (NHANES GM)", 0.8, "post-phase-out"), ("CDC blood lead reference value, 2021", 3.5, "triggers case management"), ("US children 1976-1980 (NHANES GM)", 15.2, "leaded-petrol peak"), ("CDC level of concern, 1991-2012", 10.0, "withdrawn as a threshold"), ] # ---------------------------------------------------------------- the chain def air_concentration(emission_kt_yr, vd_cm_s=VD_CM_S, area_m2=A_EUROPE): """Stage 1. Steady-state deposition-limited box. Returns ug/m^3.""" emission_ug_yr = emission_kt_yr * 1e6 * 1e3 * 1e6 # kt -> kg -> g -> ug vd_m_yr = (vd_cm_s / 100.0) * SEC_PER_YEAR return emission_ug_yr / (vd_m_yr * area_m2) def blood_enhancement(c_ug_m3, a=BLL_A, b=BLL_B): """Stage 2. Log-linear air-to-blood. Returns ug/dL above baseline.""" if c_ug_m3 <= 0: return 0.0 return math.exp(a + b * math.log(c_ug_m3)) def iq_deficit(bll_ug_dl, slope=IQ_SLOPE): """Stage 3. Log-linear blood-to-IQ. Returns a negative number of points.""" return slope * math.log(bll_ug_dl + 1.0) def attributable_iq(bll_total, baseline=BASELINE_BLL, slope=IQ_SLOPE): """IQ points lost to the Roman increment alone, as a positive number.""" return -(iq_deficit(bll_total, slope) - iq_deficit(baseline, slope)) def run_chain(emission_kt_yr, vd_cm_s=VD_CM_S, area_m2=A_EUROPE, a=BLL_A, b=BLL_B, baseline=BASELINE_BLL, slope=IQ_SLOPE): """The whole chain, one call. Returns a dict of every intermediate.""" c = air_concentration(emission_kt_yr, vd_cm_s, area_m2) d_bll = blood_enhancement(c, a, b) total = baseline + d_bll return { "E_kt": emission_kt_yr, "C_ug_m3": c, "C_ng_m3": c * 1e3, "dBLL": d_bll, "BLL_total": total, "IQ_total_deficit": iq_deficit(total, slope), "IQ_baseline_deficit": iq_deficit(baseline, slope), "IQ_attributable": attributable_iq(total, baseline, slope), } # ---------------------------------------------------------------- reporting LINES = [] def say(s=""): print(s) LINES.append(s) def rule(title=""): if title: say() say("=" * 74) say(title) say("=" * 74) else: say("-" * 74) rule("SCIENCE JOURNALING CLUB :: ROMAN LEAD DOSE-RESPONSE CHAIN") say("seed = %d numpy %s" % (SEED, np.__version__)) say("Club reconstruction. Emissions taken from McConnell et al. (2025);") say("everything between emissions and IQ is recomputed here from scratch.") # ---- Stage 0: the box itself ---------------------------------------------- rule("STAGE 0 :: THE BOX") say("domain %.0fE to %.0fE, %.0fN to %.0fN" % (LON_W, LON_E, LAT_S, LAT_N)) say("area A %.4g m^2 = %.3f million km^2" % (A_EUROPE, A_EUROPE / 1e12)) say("removal velocity v_d %.2f cm/s = %.4g m/yr" % (VD_CM_S, (VD_CM_S / 100) * SEC_PER_YEAR)) say("scavenging volume v_d * A = %.4g m^3/yr" % ((VD_CM_S / 100) * SEC_PER_YEAR * A_EUROPE)) say("at steady state all emitted lead deposits inside the box, so C = E/(v_d A)") # ---- Stage 1-3 at central values ------------------------------------------- rule("THE CHAIN AT CENTRAL ASSUMPTIONS") central = run_chain(E_CENTRAL_KT) say("1. emissions E = %.2f kt Pb/yr" % central["E_kt"]) say("2. air lead C = %.4e ug/m^3 = %.3f ng/m^3" % (central["C_ug_m3"], central["C_ng_m3"])) say("3. blood lead rise dBLL = %.3f ug/dL" % central["dBLL"]) say(" Neolithic baseline = %.2f ug/dL" % BASELINE_BLL) say(" total child blood lead = %.3f ug/dL" % central["BLL_total"]) say("4. IQ deficit at total BLL = %.3f points" % central["IQ_total_deficit"]) say(" IQ deficit at baseline = %.3f points" % central["IQ_baseline_deficit"]) say(" ROMAN-ATTRIBUTABLE LOSS = %.3f IQ points" % central["IQ_attributable"]) say() say("published for comparison: air 0.49-1.04 ng/m^3 (domain means),") say("dBLL about 2.4 ug/dL, total about 3.4 ug/dL, decrement 2.5-3 IQ points.") # ---- the paper's own two emission scenarios -------------------------------- rule("THE PAPER'S TWO EMISSION SCENARIOS, RUN THROUGH OUR BOX") say("%-34s%9s%10s%9s%9s%9s" % ("scenario", "E kt/yr", "C ng/m3", "dBLL", "BLL tot", "IQ lost")) rule() for name, e in [("1 Rio Tinto dominant", E_SCEN1_KT), ("2 distributed sources", E_SCEN2_KT), ("club midpoint", E_CENTRAL_KT)]: r = run_chain(e) say("%-34s%9.2f%10.3f%9.3f%9.3f%9.3f" % (name, r["E_kt"], r["C_ng_m3"], r["dBLL"], r["BLL_total"], r["IQ_attributable"])) # ---- validation against the paper's own published cross-checks ------------- rule("VALIDATION :: DOES OUR IQ REGRESSION REPRODUCE THEIR PUBLISHED CHECKS?") for label, bll, published in [ ("US children 1976-1980, GM 15.2 ug/dL", 15.2, 9.2), ("US children 2011-2016, GM 0.8 ug/dL", 0.8, 1.9)]: ours = -iq_deficit(bll) say("%-40s ours %6.2f published %4.1f diff %.2f" % (label, ours, published, abs(ours - published))) say("The regression is reproduced to within rounding, so any disagreement") say("downstream is in the atmospheric box, not in the epidemiology.") # ---- how steep is the curve at low dose? ----------------------------------- rule("STEEPNESS :: MARGINAL IQ COST OF ONE MORE MICROGRAM PER DECILITRE") say("d(IQ)/d(BLL) = slope / (BLL + 1), so the curve is steepest at zero dose.") say() say("%11s%13s%15s%13s" % ("BLL ug/dL", "IQ deficit", "marginal pts", "vs first ug")) rule() first_marg = abs(IQ_SLOPE) / 1.0 for bll in [0.0, 0.5, 1.0, 2.0, 3.39, 5.0, 10.0, 15.2, 25.0, 40.0]: marg = abs(IQ_SLOPE) / (bll + 1.0) say("%11.2f%13.3f%15.4f%13.3f" % (bll, iq_deficit(bll), marg, marg / first_marg)) say() say("The first ug/dL costs %.3f IQ points." % first_marg) say("The fortieth costs %.4f, a factor of %.1f less." % (abs(IQ_SLOPE) / 41.0, first_marg / (abs(IQ_SLOPE) / 41.0))) say("Roman children sat on the steep part. That is the whole reason a blood") say("lead level a modern clinic would call unremarkable does this much damage.") # ---- the air-to-blood transfer function is extremely flat ------------------ rule("THE FLATNESS OF THE AIR-TO-BLOOD STEP") say("dBLL = exp(%.3f) * C^%.3f = %.4f * C^%.3f" % (BLL_A, BLL_B, math.exp(BLL_A), BLL_B)) say("A tenfold change in air lead changes blood lead by 10^%.3f = %.4fx, " "i.e. %.1f percent." % (BLL_B, 10 ** BLL_B, 100 * (10 ** BLL_B - 1))) say() say("%12s%12s%13s%10s" % ("C ug/m^3", "C ng/m^3", "dBLL ug/dL", "IQ lost")) rule() for c in [1e-6, 1e-5, 1e-4, 5.134e-4, 1e-3, 1e-2, 1e-1, 1.0]: d = blood_enhancement(c) say("%12.1e%12.4g%13.4f%10.3f" % (c, c * 1e3, d, attributable_iq(BASELINE_BLL + d))) say() pristine = blood_enhancement(1e-6) say("At effectively pristine air (1 pg/m^3) the regression still returns") say("%.3f ug/dL, which lands on the Neolithic tooth-enamel figure of about" % pristine) say("1.0 ug/dL. That is either a good sign or a coincidence, and we cannot") say("tell which from outside the data.") # ---- emission sweep -------------------------------------------------------- rule("SWEEP :: THE CHAIN ACROSS THE PLAUSIBLE EMISSION RANGE") say("%9s%11s%9s%11s%10s%12s" % ("E kt/yr", "C ng/m^3", "dBLL", "BLL total", "IQ lost", "vs CDC 3.5")) rule() for e in [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.2, 5.0, 6.0, 8.0, 12.0]: r = run_chain(e) say("%9.1f%11.3f%9.3f%11.3f%10.3f%12.2f" % (e, r["C_ng_m3"], r["dBLL"], r["BLL_total"], r["IQ_attributable"], r["BLL_total"] / 3.5)) say() lo_r, hi_r = run_chain(1.0), run_chain(12.0) say("Emissions span a factor of 12 across that table.") say("The IQ answer spans a factor of %.2f." % (hi_r["IQ_attributable"] / lo_r["IQ_attributable"])) say("Two compressive steps in series flatten a twelvefold input into a") say("%.2f-fold output." % (hi_r["IQ_attributable"] / lo_r["IQ_attributable"])) # ---- one-at-a-time sensitivity -------------------------------------------- rule("SENSITIVITY :: ONE ASSUMPTION AT A TIME, EVERYTHING ELSE HELD CENTRAL") base_iq = central["IQ_attributable"] def chain_iq(**kw): e = kw.pop("emission_kt_yr", E_CENTRAL_KT) return run_chain(e, **kw)["IQ_attributable"] tests = [ ("emission rate E", "3.0 - 4.2 kt/yr", chain_iq(emission_kt_yr=E_SCEN2_KT), chain_iq(emission_kt_yr=E_SCEN1_KT)), ("emission rate E, wide", "2.0 - 6.0 kt/yr", chain_iq(emission_kt_yr=2.0), chain_iq(emission_kt_yr=6.0)), ("removal velocity v_d", "0.5 - 3.0 cm/s", chain_iq(vd_cm_s=3.0), chain_iq(vd_cm_s=0.5)), ("domain area A", "1.0 - 2.0 e13 m^2", chain_iq(area_m2=2.0e13), chain_iq(area_m2=1.0e13)), ("air-to-blood intercept a", "1.932 +/- 0.181", chain_iq(a=BLL_A - BLL_A_SD), chain_iq(a=BLL_A + BLL_A_SD)), ("air-to-blood exponent b", "0.140 +/- 0.054", chain_iq(b=BLL_B + BLL_B_SD), chain_iq(b=BLL_B - BLL_B_SD)), ("Neolithic baseline BLL", "0.5 - 1.5 ug/dL", chain_iq(baseline=BASELINE_BLL_HI), chain_iq(baseline=BASELINE_BLL_LO)), ("IQ regression slope", "-3.315 +/- 0.6155", chain_iq(slope=IQ_SLOPE + IQ_SLOPE_SD), chain_iq(slope=IQ_SLOPE - IQ_SLOPE_SD)), ] tests.sort(key=lambda t: abs(t[3] - t[2]), reverse=True) say("%-28s%-20s%8s%8s%8s%11s" % ("assumption", "range carried", "low", "high", "span", "% of base")) rule() for name, rng_txt, lo_v, hi_v in tests: lo_v, hi_v = sorted((lo_v, hi_v)) span = hi_v - lo_v say("%-28s%-20s%8.3f%8.3f%8.3f%11.1f" % (name, rng_txt, lo_v, hi_v, span, 100 * span / base_iq)) say() say("central answer = %.3f IQ points" % base_iq) dom = tests[0] say("DOMINANT ASSUMPTION: %s (%s), moving the answer by" % (dom[0], dom[1])) say("%.3f points, which is %.0f percent of the central value." % (abs(dom[3] - dom[2]), 100 * abs(dom[3] - dom[2]) / base_iq)) say() say("The atmospheric assumptions (v_d, A) are the ones the club had the most") say("freedom to choose, and they matter least. The epidemiological transfer") say("functions, which we copied verbatim and cannot adjust, matter most.") # ---- Monte Carlo ----------------------------------------------------------- rule("MONTE CARLO :: ALL ASSUMPTIONS VARIED TOGETHER, 200000 SEEDED DRAWS") N = 200000 e_draw = RNG.uniform(3.0, 4.2, N) # paper scenarios vd_draw = np.exp(RNG.uniform(math.log(0.5), math.log(3.0), N)) # log-uniform area_draw = RNG.uniform(1.0e13, 2.0e13, N) a_draw = RNG.normal(BLL_A, BLL_A_SD, N) b_draw = RNG.normal(BLL_B, BLL_B_SD, N) base_draw = RNG.uniform(BASELINE_BLL_LO, BASELINE_BLL_HI, N) slope_draw = RNG.normal(IQ_SLOPE, IQ_SLOPE_SD, N) e_ug = e_draw * 1e15 vd_m_yr = (vd_draw / 100.0) * SEC_PER_YEAR c_draw = e_ug / (vd_m_yr * area_draw) dbll_draw = np.exp(a_draw + b_draw * np.log(c_draw)) bll_draw = base_draw + dbll_draw iq_draw = -(slope_draw * np.log(bll_draw + 1.0) - slope_draw * np.log(base_draw + 1.0)) qs = [2.5, 16, 50, 84, 97.5] say("%-26s" % "quantity" + "".join("%11s" % ("p" + str(q)) for q in qs)) rule() for label, arr in [("air lead ng/m^3", c_draw * 1e3), ("blood lead rise ug/dL", dbll_draw), ("total blood lead ug/dL", bll_draw), ("IQ points lost", iq_draw)]: vals = np.percentile(arr, qs) say("%-26s" % label + "".join("%11.3f" % v for v in vals)) say() say("median IQ loss %.3f points" % np.median(iq_draw)) say("68%% interval %.3f to %.3f" % (np.percentile(iq_draw, 16), np.percentile(iq_draw, 84))) say("95%% interval %.3f to %.3f" % (np.percentile(iq_draw, 2.5), np.percentile(iq_draw, 97.5))) say("P(loss > 1 point) %.1f percent" % (100 * np.mean(iq_draw > 1.0))) say("P(loss > 2 points) %.1f percent" % (100 * np.mean(iq_draw > 2.0))) say("P(loss in 2.5 - 3.0) %.1f percent" % (100 * np.mean((iq_draw > 2.5) & (iq_draw < 3.0)))) say("P(total BLL > CDC 3.5) %.1f percent" % (100 * np.mean(bll_draw > 3.5))) say() say("The 95% band is wide because we let the two regression coefficients") say("wander over their published standard errors. That is honest. It also") say("means our interval is not the paper's interval, which comes from a") say("different and much more careful error propagation.") # ---- comparison table ------------------------------------------------------ rule("CONTEXT :: WHERE 3.4 ug/dL SITS AMONG NUMBERS WE CAN CHECK") roman = central["BLL_total"] say("%-38s%8s%12s%11s" % ("population / reference", "BLL", "IQ deficit", "vs Roman")) rule() rows = REF_LEVELS + [("Roman Empire, Pax Romana (this model)", roman, "club output")] rows.sort(key=lambda r: r[1]) for name, bll, note in rows: say("%-38s%8.2f%12.2f%11.2f" % (name, bll, -iq_deficit(bll), bll / roman)) say() say("Roman child, total blood lead %.2f ug/dL" % roman) say("CDC reference value since 2021 3.50 ug/dL") say("ratio %.3f" % (roman / 3.5)) say("A Roman child at the empire-wide average would today sit just under the") say("threshold at which an American clinic opens a case file.") say() gas_era = -iq_deficit(15.2) modern = -iq_deficit(0.8) say("Leaded-petrol-era US deficit %.2f points" % gas_era) say("Post-phase-out US deficit %.2f points" % modern) say("Recovered by removing tetraethyl lead from petrol: %.2f points" % (gas_era - modern)) say("Roman airborne burden %.2f points" % base_iq) say("Roman burden as a share of the American petrol burden: %.0f percent" % (100 * base_iq / (gas_era - modern))) # ---- what emission rate would it take? ------------------------------------- rule("INVERSION :: WHAT WOULD IT TAKE TO REACH THE LEADED-PETROL ERA?") target_bll = 15.2 need_dbll = target_bll - BASELINE_BLL need_c = math.exp((math.log(need_dbll) - BLL_A) / BLL_B) vd_m_yr_c = (VD_CM_S / 100.0) * SEC_PER_YEAR need_e_kt = need_c * vd_m_yr_c * A_EUROPE / 1e15 say("target total blood lead %.1f ug/dL" % target_bll) say("required air lead %.4g ug/m^3 = %.4g ng/m^3" % (need_c, need_c * 1e3)) say("required emission rate %.4g kt Pb/yr" % need_e_kt) say("that is %.4g times the Roman rate." % (need_e_kt / E_CENTRAL_KT)) say() say("The exponent of 0.14 is doing all of that work. Because blood lead") say("responds to the 0.14 power of air lead, matching twentieth-century") say("American children on the airborne route alone would need an absurd") say("emission rate. Which is the point: their lead was not mostly airborne.") # ---- the club's schematic emission history --------------------------------- rule("SCHEMATIC EMISSION HISTORY USED BY THE INTERACTIVE COMPANION") say("Shape taken from the narrative published in the paper and its 2018") say("predecessor; absolute scaling is ours, pinned so the Pax Romana plateau") say("sits at the club's 3.5 kt/yr midpoint. This is NOT digitised data.") say() HISTORY = [ (-500, 0.05, "Iron Age background; Greek silver at Laurion beginning"), (-400, 0.10, "Athenian silver coinage at full output"), (-300, 0.20, "Carthaginian working of Iberian ores"), (-206, 0.70, "Rome takes Iberia from Carthage"), (-150, 1.40, "Republican mining expansion in Hispania"), (-120, 1.80, "Republican maximum, late 2nd century BCE"), (-91, 1.00, "Social War"), (-80, 0.60, "Sullan civil wars; output near background"), (-49, 0.40, "Caesar crosses the Rubicon"), (-15, 1.20, "Augustan settlement; sharp rise begins"), (0, 2.40, "Empire consolidating"), (50, 3.20, "Julio-Claudian and Flavian output"), (100, 3.60, "Trajanic peak"), (150, 3.50, "Antonine plateau"), (165, 3.40, "Antonine Plague arrives"), (180, 1.90, "output falls with the plague"), (250, 0.80, "third-century crisis"), (400, 0.50, "late Empire"), (600, 0.20, "post-Roman minimum"), ] say("%7s%10s%12s%10s%9s %s" % ("year", "E kt/yr", "nbPb pg/g", "BLL tot", "IQ lost", "event")) rule() PG_PER_KT = 1.6 / 3.5 # club scaling: 3.5 kt/yr -> 1.6 pg/g non-background Pb ICE_BACKGROUND = 0.55 # pg/g, natural Holocene Greenland lead for yr, e, ev in HISTORY: r = run_chain(e) lbl = "%d %s" % (abs(yr), "BCE" if yr < 0 else "CE") say("%7s%10.2f%12.3f%10.3f%9.3f %s" % (lbl, e, ICE_BACKGROUND + e * PG_PER_KT, r["BLL_total"], r["IQ_attributable"], ev)) say() say("ice calibration: nbPb(pg/g) = %.4f * E(kt/yr), background %.2f pg/g" % (PG_PER_KT, ICE_BACKGROUND)) # ---- figure data dump ------------------------------------------------------ rule("FIGURE DATA :: EXACT NUMBERS PLOTTED IN THE ARTICLE") say("Figure 2a, air-to-blood transfer function, C in ng/m^3 -> dBLL ug/dL") for c_ng in [0.01, 0.03, 0.1, 0.3, 0.513, 1.0, 3.0, 10.0, 30.0, 100.0, 300.0]: say(" %8.3f %7.3f" % (c_ng, blood_enhancement(c_ng * 1e-3))) say() say("Figure 2b, blood-to-IQ transfer function, BLL ug/dL -> deficit points") for bll in [0, 1, 2, 3.39, 5, 8, 12, 15.2, 20, 30, 40]: say(" %8.2f %7.3f" % (bll, iq_deficit(bll))) say() say("Figure 3, sensitivity spans in IQ points (sorted, widest first)") for name, rng_txt, lo_v, hi_v in tests: lo_v, hi_v = sorted((lo_v, hi_v)) say(" %-28s%8.3f%8.3f span %6.3f" % (name, lo_v, hi_v, hi_v - lo_v)) say() say("Figure 4, blood lead ladder") for name, bll, note in rows: say(" %-38s%7.2f deficit %5.2f (%s)" % (name, bll, -iq_deficit(bll), note)) say() say("Figure 1, chain compression factors") say(" emissions 3.0 -> 4.2 kt/yr is a factor of %.3f" % (4.2 / 3.0)) say(" air %.3f -> %.3f ng/m^3, factor %.3f" % (run_chain(3.0)["C_ng_m3"], run_chain(4.2)["C_ng_m3"], 4.2 / 3.0)) say(" blood %.3f -> %.3f ug/dL, factor %.3f" % (run_chain(3.0)["dBLL"], run_chain(4.2)["dBLL"], run_chain(4.2)["dBLL"] / run_chain(3.0)["dBLL"])) say(" IQ %.3f -> %.3f points, factor %.3f" % (run_chain(3.0)["IQ_attributable"], run_chain(4.2)["IQ_attributable"], run_chain(4.2)["IQ_attributable"] / run_chain(3.0)["IQ_attributable"])) rule("END") with open(__file__.replace(".py", "-output.txt"), "w", encoding="utf-8") as fh: fh.write("\n".join(LINES) + "\n")