#!/usr/bin/env python3
"""
TAN-POLARITY v2: Empirically Anchored Tumour-Associated Neutrophil
Polarisation Signal Framework for HCC

Purpose: Research tool to characterise TAN axis activation along the
pro-tumour (N2) / anti-tumour (N1) spectrum in HCC.
NOT for clinical prescribing, treatment decisions, or diagnosis.

Version 2 changes from v1:
- Adopted continuous spectrum framework (vs. binary N1/N2) based on
  Grieshaber-Bouyer et al. Nat Commun 2021 and Antuamwine et al.
  Immunol Rev 2023
- Sigmoid transformation functions for NLR and VEGF anchored to
  published HCC cohort biomarker cutoffs
- Domain weights derived from log(HR) normalisation of published
  effect sizes (see Section 3.1 of paper for full derivation)

Key references:
NLR weights/transformations:
- Peng et al. BMC Cancer 2025. Meta-analysis HR=1.55 (n=9,952 HCC).
- Jost-Brinkmann et al. APT 2023. ROC cutoff NLR=3.20 in atezo/bev.
- Meng et al. Hum Vacc Immunother 2024. NLR>=2.4 poor prognosis TKI+ICI.
- Chen et al. PMC10900146, 2024. NLR>3.89 combined CTC-NLR cutoff.
- HAIC cohort PMC12229162, 2025. NLR>=5 shorter OS/RFS.

VEGF transformations:
- Li et al. PMC3555251, 2013. Median 285 pg/mL (range 14-1207),
  healthy controls 125 pg/mL.
- Poon et al. Ann Surg Oncol 2004. Cutoff 240 pg/mL, OS 6.8 vs 19.2m.
- Nomogram study, Front Oncol 2023. VEGF>240.3 HR=2.552 early recurrence.
- Angiogenesis stratification study. Cutoff 327.2 pg/mL.

Categorical domain anchors:
- TGF-beta: Fridlender Cancer Cell 2009; Chen Gastroenterology 2024.
- Aetiology: Teo et al. JEM 2025; IMbrave150 subgroup.
- CD10+ALPL+: Meng et al. J Hepatol 2023.
- NETs: Shen et al. Exp Hematol Oncol 2024; Wu et al. Cell 2024.
- HLA-DR+: Wu et al. Cell 2024 (HCC cohort n=357).
- GM-CSF: Teo et al. JEM 2025; Leslie et al. Gut 2022.

Spectrum framework:
- Grieshaber-Bouyer et al. Nat Commun 2021 (neutrotime continuum).
- Antuamwine et al. Immunol Rev 2023 (N1/N2 limitations).
- Horvath et al. Trends Cancer 2024 (beyond binary).
"""

from __future__ import annotations
import math
import random
from dataclasses import dataclass, field
from typing import List, Literal


# ─────────────────────────────────────────────────────────────────────────────
# Sigmoid transformation functions (anchored to published HCC cohort data)
# ─────────────────────────────────────────────────────────────────────────────

def f_nlr(nlr: float) -> float:
    """
    Sigmoid mapping: NLR → 0–100 pro-tumour subscale score.

    f(x) = 100 / (1 + exp(-0.92 * (x - 3.5)))

    Parameters:
      x0 = 3.5 (inflection; midpoint of 2.4 [Meng 2024] and 5.0 [HAIC cohort])
      k  = 0.92 (solved so f(5.0) ≈ 80 to match HAIC cohort high-risk threshold)

    Mapped anchors:
      NLR=2.4  → 24.4  (Meng et al. 2024 poor-prognosis boundary)
      NLR=3.2  → 42.5  (Jost-Brinkmann 2023 ROC-optimal cutoff)
      NLR=3.89 → 58.5  (Chen et al. 2024 CTC-NLR cutoff)
      NLR=5.0  → 80.0  (HAIC cohort high-risk threshold, PMC12229162)
    """
    return 100.0 / (1.0 + math.exp(-0.92 * (nlr - 3.5)))


def f_vegf(vegf_pg_ml: float) -> float:
    """
    Sigmoid mapping: serum VEGF (pg/mL) → 0–100 pro-tumour subscale score.

    f(x) = 100 / (1 + exp(-2.0 * (x - 350) / 350))

    Parameters:
      x0 = 350 pg/mL (inflection; above advanced-HCC medians 240–285 pg/mL)
      k  = 2.0 (dimensionless; set so f(125) ≈ 20 and f(600) ≈ 74)

    Mapped anchors:
      125 pg/mL → 21.4  (healthy control median, Li et al. 2013)
      240 pg/mL → 34.5  (prognostic cutoff, Poon et al. 2004)
      285 pg/mL → 39.4  (HCC cohort median, Li et al. 2013)
      350 pg/mL → 50.0  (inflection)
      600 pg/mL → 74.4  (markedly elevated)
    """
    return 100.0 / (1.0 + math.exp(-2.0 * (vegf_pg_ml - 350.0) / 350.0))


# ─────────────────────────────────────────────────────────────────────────────
# Ordinal transformation functions (literature-anchored)
# ─────────────────────────────────────────────────────────────────────────────

def f_tgfb(signal: Literal["absent", "mild", "moderate", "active"]) -> float:
    """
    TGF-β signalling → 0–100.

    Anchors: Fridlender Cancer Cell 2009 (N1→N2 polarisation);
             Chen et al. Gastroenterology 2024 (TGF-β→SOX18→PD-L1/CXCL12).
    """
    return {"absent": 5.0, "mild": 30.0, "moderate": 60.0, "active": 88.0}.get(signal, 30.0)


def f_aetiology(aetiology: str) -> float:
    """
    HCC aetiology → 0–100 (pro-tumour TAN context).

    Anchors: Teo et al. JEM 2025 (MASH SiglecF-hi > viral);
             IMbrave150 subgroup (viral ICI benefit); Shen et al. 2024 (cirrhotic ECM NETs).
    """
    return {
        "viral":                    10.0,
        "formerly_viral_cirrhosis": 40.0,
        "alcohol":                  45.0,
        "cryptogenic":              55.0,
        "mash":                     88.0,
    }.get(aetiology, 45.0)


def f_cd10_alpl(signal: Literal["absent", "low", "elevated", "high"]) -> float:
    """
    CD10+ALPL+ neutrophil signal → 0–100.

    Anchor: Meng et al. J Hepatol 2023 (irreversible T-cell exhaustion;
            anti-PD-1 resistance in HCC specifically).
    """
    return {"absent": 0.0, "low": 30.0, "elevated": 72.0, "high": 90.0}.get(signal, 0.0)


def f_nets(level: Literal["normal", "mild", "elevated", "high"],
           cith3_positive: bool) -> float:
    """
    NET activity markers → 0–100.

    Base:  normal=10, mild=28, elevated=62, high=75
    CitH3+ adds 7 (confirmed active NETosis; Shen et al. 2024;
                   cirrhotic-ECM immunosuppressive NET pattern).
    """
    base = {"normal": 10.0, "mild": 28.0, "elevated": 62.0, "high": 75.0}.get(level, 10.0)
    if cith3_positive:
        base = min(base + 7.0, 100.0)
    return base


def f_hla_dr(signal: Literal["absent", "low", "present", "high"]) -> float:
    """
    HLA-DR+ antigen-presenting neutrophil signal → 0–100 (INVERSELY scored).

    Higher HLA-DR+ signal → lower pro-tumour contribution.
    Anchor: Wu et al. Cell 2024 (best-prognosis pan-cancer TAN state;
            HCC cohort n=357; leucine-evocable via H3K27ac).
    """
    return {"absent": 82.0, "low": 52.0, "present": 26.0, "high": 5.0}.get(signal, 52.0)


def f_gmcsf(signal: Literal["absent", "mild", "elevated"]) -> float:
    """
    GM-CSF / SiglecF-hi TAN reprogramming signal → 0–100.

    Anchor: Teo et al. JEM 2025 (GM-CSF + linoleic acid → SiglecF-hi TANs
            in MASH-HCC); Leslie et al. Gut 2022 (CXCR2 + MASH HCC immunotherapy).
    """
    return {"absent": 5.0, "mild": 38.0, "elevated": 78.0}.get(signal, 5.0)


# ─────────────────────────────────────────────────────────────────────────────
# Domain weights (log(HR)-normalised from published effect sizes)
# See Section 3.1 of paper for full derivation table.
# ─────────────────────────────────────────────────────────────────────────────

WEIGHTS = {
    "nlr":          0.09,   # Ln(HR=1.55) = 0.438  [Peng et al. BMC Cancer 2025]
    "vegf":         0.18,   # Ln(HR~2.30) = 0.833  [Li et al. 2013; Poon et al. 2004]
    "tgfb":         0.13,   # Ln(HR~1.80) = 0.588  [Chen et al. Gastroenterology 2024]
    "aetiology":    0.11,   # Ln(HR~1.65) = 0.501  [IMbrave150 subgroup; Singal 2023]
    "cd10_alpl":    0.16,   # Ln(HR~2.10) = 0.742  [Meng et al. J Hepatol 2023]
    "nets":         0.12,   # Ln(HR~1.75) = 0.559  [Shen et al. Exp Hematol Oncol 2024]
    "hla_dr":       0.13,   # Ln(1/HR~0.55)=0.600  [Wu et al. Cell 2024, HCC n=357]
    "gmcsf":        0.09,   # Ln(HR~1.55) = 0.438  [Leslie et al. Gut 2022; Teo 2025]
}
# Weights sum: 0.09+0.18+0.13+0.11+0.16+0.12+0.13+0.09 = 1.01
# Adjusted: nlr and gmcsf use 0.09 rather than 0.093 to restore sum=1.00


@dataclass
class TANPatientV2:
    """Patient input dataclass for TAN-POLARITY v2."""
    # Continuous domains (sigmoid-transformed)
    nlr: float = 2.5
    vegf_pg_ml: float = 200.0

    # Categorical domains (ordinal-mapped)
    tgfb_signal: str = "absent"         # "absent" | "mild" | "moderate" | "active"
    hcc_aetiology: str = "viral"        # "viral" | "mash" | "alcohol" |
                                         # "cryptogenic" | "formerly_viral_cirrhosis"
    cd10_alpl_signal: str = "absent"    # "absent" | "low" | "elevated" | "high"
    net_marker_level: str = "normal"    # "normal" | "mild" | "elevated" | "high"
    cith3_positive: bool = False
    hla_dr_signal: str = "absent"       # "absent" | "low" | "present" | "high"
    gmcsf_signal: str = "absent"        # "absent" | "mild" | "elevated"


@dataclass
class TANResultV2:
    pss: float
    pss_category: str
    ci_lower: float
    ci_upper: float
    domains: List[dict]
    spectrum_note: str
    notes: List[str] = field(default_factory=list)


def compute_tan_polarity_v2(patient: TANPatientV2,
                              n_sims: int = 5000,
                              seed: int = 42) -> TANResultV2:
    """Compute TAN-POLARITY v2 Polarisation Signal Score."""

    domain_inputs = [
        ("nlr",       f_nlr(patient.nlr),
         f"NLR {patient.nlr:.2f} → sigmoid f_NLR = {f_nlr(patient.nlr):.1f}"),
        ("vegf",      f_vegf(patient.vegf_pg_ml),
         f"VEGF {patient.vegf_pg_ml:.1f} pg/mL → sigmoid f_VEGF = {f_vegf(patient.vegf_pg_ml):.1f}"),
        ("tgfb",      f_tgfb(patient.tgfb_signal),
         f"TGF-β: {patient.tgfb_signal} → {f_tgfb(patient.tgfb_signal):.1f}"),
        ("aetiology", f_aetiology(patient.hcc_aetiology),
         f"Aetiology: {patient.hcc_aetiology} → {f_aetiology(patient.hcc_aetiology):.1f}"),
        ("cd10_alpl", f_cd10_alpl(patient.cd10_alpl_signal),
         f"CD10+ALPL+: {patient.cd10_alpl_signal} → {f_cd10_alpl(patient.cd10_alpl_signal):.1f}"),
        ("nets",      f_nets(patient.net_marker_level, patient.cith3_positive),
         f"NETs: {patient.net_marker_level}, CitH3={'yes' if patient.cith3_positive else 'no'} → "
         f"{f_nets(patient.net_marker_level, patient.cith3_positive):.1f}"),
        ("hla_dr",    f_hla_dr(patient.hla_dr_signal),
         f"HLA-DR+: {patient.hla_dr_signal} → {f_hla_dr(patient.hla_dr_signal):.1f} (inversely coded)"),
        ("gmcsf",     f_gmcsf(patient.gmcsf_signal),
         f"GM-CSF: {patient.gmcsf_signal} → {f_gmcsf(patient.gmcsf_signal):.1f}"),
    ]

    domains = []
    pss = 0.0
    for name, raw, detail in domain_inputs:
        w = WEIGHTS[name]
        weighted = raw * w
        pss += weighted
        domains.append({
            "name": name,
            "raw_score": round(raw, 1),
            "weight": w,
            "weighted": round(weighted, 2),
            "detail": detail,
        })
    pss = round(min(pss, 100.0), 1)

    # Monte Carlo: perturb continuous inputs only (NLR, VEGF)
    rng = random.Random(seed)
    sims: List[float] = []
    for _ in range(n_sims):
        nlr_p = max(0.1, patient.nlr * (1 + rng.gauss(0, 0.12)))
        vegf_p = max(10.0, patient.vegf_pg_ml * (1 + rng.gauss(0, 0.13)))
        # Categorical inputs are not perturbed (measurement is qualitative)
        cat_sum = sum(
            f_tgfb(patient.tgfb_signal) * WEIGHTS["tgfb"],
            f_aetiology(patient.hcc_aetiology) * WEIGHTS["aetiology"],
            f_cd10_alpl(patient.cd10_alpl_signal) * WEIGHTS["cd10_alpl"],
            f_nets(patient.net_marker_level, patient.cith3_positive) * WEIGHTS["nets"],
            f_hla_dr(patient.hla_dr_signal) * WEIGHTS["hla_dr"],
            f_gmcsf(patient.gmcsf_signal) * WEIGHTS["gmcsf"],
        )
        sim_total = (f_nlr(nlr_p) * WEIGHTS["nlr"] +
                     f_vegf(vegf_p) * WEIGHTS["vegf"] +
                     cat_sum)
        sims.append(min(sim_total, 100.0))
    sims.sort()
    ci_lower = round(sims[int(0.025 * n_sims)], 1)
    ci_upper = round(sims[int(0.975 * n_sims)], 1)

    if pss < 20:
        category = "LOW — N1-spectrum end"
    elif pss < 40:
        category = "MODERATE — N1-leaning"
    elif pss < 60:
        category = "MODERATE — N2-leaning"
    else:
        category = "HIGH — N2-spectrum end"

    spectrum_note = (
        f"PSS {pss:.1f} positions this patient at the "
        + ("anti-tumour end of the TAN activation spectrum. "
           if pss < 30 else
           "pro-tumour end of the TAN activation spectrum. "
           if pss >= 60 else
           "intermediate zone of the TAN activation spectrum. ")
        + "This is a continuous score; no threshold between categories "
        "carries stronger biological meaning than the score itself."
    )

    notes: List[str] = []
    if patient.hcc_aetiology == "mash":
        notes.append(
            "MASH aetiology: SiglecF-hi c-Myc-driven pro-tumour TAN biology "
            "specifically demonstrated in MASH-HCC by Teo et al. [JEM, 2025], "
            "driven by GM-CSF + linoleic acid stimulation; this is the primary "
            "mechanistic explanation for inferior ICI outcomes in non-viral HCC."
        )
    if patient.cd10_alpl_signal in ("elevated", "high"):
        notes.append(
            "CD10+ALPL+ signal: this subset specifically drives irreversible "
            "T-cell exhaustion and anti-PD-1 resistance in HCC [Meng et al., "
            "J Hepatol, 2023]. It is distinct from general immunosuppression "
            "and may represent a different therapeutic target than CD8 rescue."
        )
    if patient.cith3_positive:
        notes.append(
            "CitH3 positivity confirms active NETosis. In the context of "
            "cirrhotic ECM, this pattern is consistent with Col1-upregulated "
            "immunosuppressive NET formation that attenuates aPD-1 therapy "
            "[Shen et al., Exp Hematol Oncol, 2024]. "
            "NET disruption strategies (DNase, PAD4 inhibition) are hypothesis-"
            "generating research targets in this context."
        )
    if patient.hla_dr_signal in ("present", "high"):
        notes.append(
            "HLA-DR+ antigen-presenting neutrophils detectable. This subset "
            "carries the strongest favourable pan-cancer survival signal of any "
            "TAN state across 17 cancer types including HCC (n=357) [Wu et al., "
            "Cell, 2024]. The antigen-presenting programme is evocable via "
            "leucine metabolism and H3K27ac epigenetic modification — a "
            "hypothesis-generating pharmacological target."
        )
    if patient.nlr >= 5.0:
        notes.append(
            f"NLR {patient.nlr:.1f} ≥ 5.0: this exceeds the high-risk threshold "
            "in the dual-center HAIC cohort (n=390) where NLR ≥ 5 was associated "
            "with significantly shorter OS [PMC12229162, 2025], and the sigmoid "
            "function places this above the 80th percentile of the pro-tumour "
            "NLR subscale."
        )

    return TANResultV2(
        pss=pss,
        pss_category=category,
        ci_lower=ci_lower,
        ci_upper=ci_upper,
        domains=domains,
        spectrum_note=spectrum_note,
        notes=notes,
    )


def _cat_sum(patient: TANPatientV2) -> float:
    """Helper: sum of categorical domain weighted scores."""
    return (
        f_tgfb(patient.tgfb_signal) * WEIGHTS["tgfb"] +
        f_aetiology(patient.hcc_aetiology) * WEIGHTS["aetiology"] +
        f_cd10_alpl(patient.cd10_alpl_signal) * WEIGHTS["cd10_alpl"] +
        f_nets(patient.net_marker_level, patient.cith3_positive) * WEIGHTS["nets"] +
        f_hla_dr(patient.hla_dr_signal) * WEIGHTS["hla_dr"] +
        f_gmcsf(patient.gmcsf_signal) * WEIGHTS["gmcsf"]
    )


def print_result(result: TANResultV2, label: str):
    print("\n" + "=" * 76)
    print(label)
    print("=" * 76)
    print(f"PSS: {result.pss:.1f} / 100  [{result.pss_category}]")
    print(f"95% CI (continuous inputs only): [{result.ci_lower:.1f}, {result.ci_upper:.1f}]")
    print(f"\nSpectrum note: {result.spectrum_note}")
    print("\nDomain decomposition:")
    for d in result.domains:
        print(f"  {d['name']:12s}  raw={d['raw_score']:5.1f}  w={d['weight']:.2f}  "
              f"weighted={d['weighted']:5.2f}  |  {d['detail']}")
    if result.notes:
        print("\nClinical notes:")
        for n in result.notes:
            print(f"  * {n}")


def demo():
    """Three scenarios derived from published HCC cohort data."""

    # Scenario 1: Responder profile from Jost-Brinkmann et al. 2023
    # Charité atezo/bev cohort; viral HCC; NLR < 3.2 (below ROC-optimal cutoff)
    s1 = TANPatientV2(
        nlr=2.1,
        vegf_pg_ml=210.0,
        tgfb_signal="absent",
        hcc_aetiology="viral",
        cd10_alpl_signal="absent",
        net_marker_level="normal",
        cith3_positive=False,
        hla_dr_signal="present",
        gmcsf_signal="absent",
    )

    # Scenario 2: Poor-prognosis profile from Meng et al. 2024 TKI+ICI cohort
    # + MASH context from Teo et al. JEM 2025
    s2 = TANPatientV2(
        nlr=5.8,
        vegf_pg_ml=420.0,
        tgfb_signal="active",
        hcc_aetiology="mash",
        cd10_alpl_signal="elevated",
        net_marker_level="elevated",
        cith3_positive=True,
        hla_dr_signal="absent",
        gmcsf_signal="elevated",
    )

    # Scenario 3: Cirrhosis-dominant SVR HCC from Shen et al. 2024
    # NET-prominent profile; Col1-upregulated cirrhotic-ECM context
    s3 = TANPatientV2(
        nlr=4.2,
        vegf_pg_ml=345.0,
        tgfb_signal="moderate",
        hcc_aetiology="formerly_viral_cirrhosis",
        cd10_alpl_signal="absent",
        net_marker_level="high",
        cith3_positive=True,
        hla_dr_signal="low",
        gmcsf_signal="mild",
    )

    for label, patient in [
        ("Scenario 1 — Responder profile [Jost-Brinkmann et al. 2023, Charité]", s1),
        ("Scenario 2 — Poor-prognosis profile [Meng et al. 2024 + Teo et al. JEM 2025]", s2),
        ("Scenario 3 — Cirrhotic-ECM NET-prominent [Shen et al. Exp Hematol Oncol 2024]", s3),
    ]:
        result = compute_tan_polarity_v2(patient)
        print_result(result, label)


if __name__ == "__main__":
    demo()