#!/usr/bin/env python3
"""
HCC-METASCORE: Biomarker-Driven Composite Scoring Framework for
Systemic Therapy Signal Prioritisation in HCC with Extrahepatic Metastatic Spread

Purpose: Research tool to narrow the field of analytical focus and structure
multidisciplinary discussion. NOT for clinical prescribing or treatment decisions.

Weight rationale:
- AFP/AFP-L3 (0.12): Validated vascular invasion proxy (Lok et al. 2010)
- DCP/PIVKA-II (0.08): Portal vein invasion predictor (Imamura et al. 2008)
- CTC/ctDNA (0.10): Active haematogenous dissemination signal (Ye et al. 2021)
- MVI (0.14): Strongest histological metastatic predictor (Roayaie et al. 2004)
- EMT markers (0.10): Direct invasive capacity drivers (Schulze et al. 2015)
- TME/immune (0.14): Checkpoint inhibitor responsiveness candidate (Llovet et al. 2021)
- Angiogenic burden (0.12): Core sorafenib/bevacizumab target (Llovet et al. 2008)
- FGFR signal (0.06): Lenvatinib-specific target (Kudo et al. 2018)
- Genetic drivers (0.10): TP53/PTEN/RB1 metastatic burden
- Wnt/epigenetic (0.04): EMT and stem-cell metastatic programme
"""

from __future__ import annotations
import random
from dataclasses import dataclass, field
from typing import Optional, List

# ─────────────────────────────────────────────
# Patient input dataclass
# ─────────────────────────────────────────────

@dataclass
class HCCPatient:
    # Circulating biomarkers
    afp_ng_ml: float = 20.0                  # AFP in ng/mL; normal <20
    afp_l3_pct: float = 5.0                  # AFP-L3 fraction %; high-risk threshold ~15%
    dcp_mau_ml: float = 40.0                 # DCP/PIVKA-II; high-risk threshold ~40 mAU/mL
    ctc_per_75ml: int = 0                    # CTCs per 7.5 mL blood
    ctdna_vaf_pct: float = 0.0               # ctDNA variant allele fraction %

    # Vascular invasion
    microvascular_invasion: bool = False      # MVI on histology
    macrovascular_invasion: bool = False      # Portal/hepatic vein tumour thrombus

    # EMT markers
    e_cadherin_loss: bool = False            # Reduced E-cadherin expression
    vimentin_positive: bool = False          # Vimentin upregulation
    ctnnb1_mutation: bool = False            # Wnt/beta-catenin activating mutation

    # TME / immune
    nlr: float = 2.5                         # Neutrophil-to-lymphocyte ratio
    pd_l1_tps_pct: float = 0.0              # PD-L1 tumour proportion score %
    til_density: str = "low"                 # "low", "moderate", "high"

    # Angiogenesis
    vegf_pg_ml: float = 100.0               # Serum VEGF in pg/mL
    fgfr_amplification: bool = False         # FGFR1/2/3/4 amplification or mutation

    # Genetic drivers
    tp53_mutation: bool = False
    pten_loss: bool = False
    rb1_loss: bool = False

    # Epigenetic
    epigenetic_dysregulation: bool = False   # Aberrant methylation / chromatin remodelling reported


# ─────────────────────────────────────────────
# Domain scoring functions
# ─────────────────────────────────────────────

def score_afp(afp: float, afp_l3: float) -> tuple[float, str]:
    """
    AFP >400 ng/mL and AFP-L3 >15% are established thresholds for MVI risk.
    Ref: Lok et al. Hepatology 2010.
    """
    s = 0.0
    if afp > 400:
        s += 55
    elif afp > 200:
        s += 35
    elif afp > 100:
        s += 20
    elif afp > 20:
        s += 8
    if afp_l3 > 15:
        s += 30
    elif afp_l3 > 10:
        s += 15
    return min(s, 100), f"AFP={afp:.0f} ng/mL, AFP-L3={afp_l3:.1f}%"


def score_dcp(dcp: float) -> tuple[float, str]:
    """
    DCP/PIVKA-II >40 mAU/mL associated with portal vein invasion.
    Ref: Imamura et al. Hepatology 2008.
    """
    if dcp > 400:
        return 90, f"DCP={dcp:.0f} mAU/mL [markedly elevated]"
    if dcp > 100:
        return 60, f"DCP={dcp:.0f} mAU/mL [elevated]"
    if dcp > 40:
        return 30, f"DCP={dcp:.0f} mAU/mL [above threshold]"
    return 5, f"DCP={dcp:.0f} mAU/mL [normal range]"


def score_ctc_ctdna(ctc: int, vaf: float) -> tuple[float, str]:
    """
    CTC detection indicates haematogenous dissemination.
    ctDNA VAF tracks tumour burden. Ref: Ye et al. Hepatology 2021.
    """
    s = 0.0
    if ctc >= 5:
        s += 50
    elif ctc >= 2:
        s += 30
    elif ctc == 1:
        s += 15
    if vaf >= 5.0:
        s += 40
    elif vaf >= 1.0:
        s += 25
    elif vaf >= 0.5:
        s += 10
    return min(s, 100), f"CTCs={ctc}/7.5mL, ctDNA VAF={vaf:.1f}%"


def score_mvi(micro: bool, macro: bool) -> tuple[float, str]:
    """
    MVI is the strongest histological predictor of extrahepatic spread.
    Macrovascular = BCLC-C. Ref: Roayaie et al. Ann Surg 2004.
    """
    if macro:
        return 95, "Macrovascular invasion present (portal/hepatic vein)"
    if micro:
        return 50, "Microvascular invasion present on histology"
    return 0, "No vascular invasion identified"


def score_emt(e_cad_loss: bool, vimentin: bool, ctnnb1: bool) -> tuple[float, str]:
    """
    EMT marker composite. E-cadherin loss + vimentin = full mesenchymal shift.
    CTNNB1 mutation activates Wnt/beta-catenin, drives stemness.
    Ref: Schulze et al. Nature Genetics 2015.
    """
    s = 0.0
    details = []
    if e_cad_loss:
        s += 35
        details.append("E-cad loss")
    if vimentin:
        s += 35
        details.append("vimentin+")
    if ctnnb1:
        s += 25
        details.append("CTNNB1 mut")
    return min(s, 100), ", ".join(details) if details else "No EMT markers detected"


def score_tme(nlr: float, pd_l1: float, til: str) -> tuple[float, str]:
    """
    TME immune profile. High PD-L1 + high TILs = immune-inflamed phenotype,
    associated with checkpoint inhibitor responsiveness.
    NLR >5 indicates systemic inflammation suppressing immune response.
    Ref: Llovet et al. Nature Reviews Clinical Oncology 2021.
    """
    s = 0.0
    details = []
    # PD-L1 expression
    if pd_l1 >= 10:
        s += 40
        details.append(f"PD-L1 TPS {pd_l1:.0f}%")
    elif pd_l1 >= 1:
        s += 20
        details.append(f"PD-L1 TPS {pd_l1:.0f}%")
    else:
        details.append("PD-L1 <1%")
    # TIL density
    if til == "high":
        s += 35
        details.append("TIL-high")
    elif til == "moderate":
        s += 15
        details.append("TIL-moderate")
    else:
        details.append("TIL-low")
    # NLR: high NLR is immunosuppressive — penalises the immune signal
    if nlr > 5:
        s = max(s - 20, 0)
        details.append(f"NLR {nlr:.1f} [immunosuppressed systemic milieu]")
    else:
        details.append(f"NLR {nlr:.1f}")
    return min(s, 100), ", ".join(details)


def score_angiogenesis(vegf: float) -> tuple[float, str]:
    """
    VEGF drives HCC vascularisation and is the primary mechanistic target
    of sorafenib and bevacizumab. Ref: Llovet et al. NEJM 2008.
    """
    if vegf > 400:
        return 90, f"VEGF={vegf:.0f} pg/mL [markedly elevated]"
    if vegf > 250:
        return 65, f"VEGF={vegf:.0f} pg/mL [elevated]"
    if vegf > 150:
        return 40, f"VEGF={vegf:.0f} pg/mL [moderately elevated]"
    if vegf > 80:
        return 20, f"VEGF={vegf:.0f} pg/mL [mildly elevated]"
    return 5, f"VEGF={vegf:.0f} pg/mL [near normal]"


def score_fgfr(fgfr_amp: bool) -> tuple[float, str]:
    """
    FGFR amplification/dysregulation is targeted by lenvatinib but not sorafenib.
    Presence is a lenvatinib-differentiating signal. Ref: Kudo et al. Lancet 2018.
    """
    if fgfr_amp:
        return 85, "FGFR amplification/dysregulation detected"
    return 0, "No FGFR amplification detected"


def score_genetic_drivers(tp53: bool, pten: bool, rb1: bool) -> tuple[float, str]:
    """
    Convergent loss of TP53, PTEN, and RB1 is associated with aggressive
    metastatic phenotype in extrahepatic HCC deposits.
    """
    s = 0.0
    details = []
    if tp53:
        s += 40
        details.append("TP53 mut")
    if pten:
        s += 35
        details.append("PTEN loss")
    if rb1:
        s += 25
        details.append("RB1 loss")
    return min(s, 100), ", ".join(details) if details else "No high-risk driver mutations detected"


def score_wnt_epigenetic(ctnnb1: bool, epi: bool) -> tuple[float, str]:
    """
    Wnt/beta-catenin and epigenetic dysregulation contribute to pro-metastatic
    gene programme activation. CTNNB1 mutation also noted in EMT domain.
    Ref: Hoshida et al. New England Journal of Medicine 2008.
    """
    s = 0.0
    details = []
    if ctnnb1:
        s += 50
        details.append("CTNNB1 mutation")
    if epi:
        s += 45
        details.append("Epigenetic dysregulation reported")
    return min(s, 100), ", ".join(details) if details else "No Wnt/epigenetic dysregulation noted"


# ─────────────────────────────────────────────
# Weights
# ─────────────────────────────────────────────

WEIGHTS = {
    "afp":          0.12,
    "dcp":          0.08,
    "ctc_ctdna":    0.10,
    "mvi":          0.14,
    "emt":          0.10,
    "tme":          0.14,
    "angiogenesis": 0.12,
    "fgfr":         0.06,
    "genetic":      0.10,
    "wnt_epi":      0.04,
}
assert abs(sum(WEIGHTS.values()) - 1.0) < 1e-9, "Weights must sum to 1.0"


# ─────────────────────────────────────────────
# Result dataclass
# ─────────────────────────────────────────────

@dataclass
class HCCResult:
    composite_score: float
    ci_lower: float
    ci_upper: float
    score_category: str
    domains: list
    pathway_signal: dict
    interpretive_notes: List[str] = field(default_factory=list)


# ─────────────────────────────────────────────
# Core computation
# ─────────────────────────────────────────────

def compute_domain_scores(p: HCCPatient) -> list:
    return [
        ("afp",          *score_afp(p.afp_ng_ml, p.afp_l3_pct),          WEIGHTS["afp"]),
        ("dcp",          *score_dcp(p.dcp_mau_ml),                         WEIGHTS["dcp"]),
        ("ctc_ctdna",    *score_ctc_ctdna(p.ctc_per_75ml, p.ctdna_vaf_pct), WEIGHTS["ctc_ctdna"]),
        ("mvi",          *score_mvi(p.microvascular_invasion, p.macrovascular_invasion), WEIGHTS["mvi"]),
        ("emt",          *score_emt(p.e_cadherin_loss, p.vimentin_positive, p.ctnnb1_mutation), WEIGHTS["emt"]),
        ("tme",          *score_tme(p.nlr, p.pd_l1_tps_pct, p.til_density), WEIGHTS["tme"]),
        ("angiogenesis", *score_angiogenesis(p.vegf_pg_ml),                WEIGHTS["angiogenesis"]),
        ("fgfr",         *score_fgfr(p.fgfr_amplification),               WEIGHTS["fgfr"]),
        ("genetic",      *score_genetic_drivers(p.tp53_mutation, p.pten_loss, p.rb1_loss), WEIGHTS["genetic"]),
        ("wnt_epi",      *score_wnt_epigenetic(p.ctnnb1_mutation, p.epigenetic_dysregulation), WEIGHTS["wnt_epi"]),
    ]


def pathway_signal_profile(domains: list) -> dict:
    """
    Maps domain scores to three therapy signal axes.
    Returns descriptive signal levels, not recommendations.
    """
    domain_map = {d[0]: d[1] for d in domains}

    sorafenib_signal = (
        domain_map["angiogenesis"] * 0.50 +
        domain_map["afp"] * 0.30 +
        domain_map["mvi"] * 0.20
    )
    lenvatinib_signal = (
        domain_map["angiogenesis"] * 0.40 +
        domain_map["fgfr"] * 0.30 +
        domain_map["afp"] * 0.20 +
        domain_map["mvi"] * 0.10
    )
    atezo_bev_signal = (
        domain_map["tme"] * 0.55 +
        domain_map["angiogenesis"] * 0.30 +
        domain_map["ctc_ctdna"] * 0.15
    )

    def level(score):
        if score >= 60: return "VERY HIGH"
        if score >= 40: return "HIGH"
        if score >= 20: return "MODERATE"
        return "LOW"

    return {
        "Sorafenib pathway signal": (round(sorafenib_signal, 1), level(sorafenib_signal)),
        "Lenvatinib pathway signal": (round(lenvatinib_signal, 1), level(lenvatinib_signal)),
        "Atezo/Bev pathway signal": (round(atezo_bev_signal, 1), level(atezo_bev_signal)),
    }


def compute_hcc_score(patient: HCCPatient, n_simulations: int = 5000, seed: int = 42) -> HCCResult:
    domains = compute_domain_scores(patient)
    composite = min(sum(score * weight for _, score, _, weight in domains), 100.0)

    # Monte Carlo: perturb continuous inputs only
    rng = random.Random(seed)
    sims = []
    for _ in range(n_simulations):
        def perturb(val, cv=0.12):
            return max(0.0, val * (1 + rng.gauss(0, cv)))

        noisy = HCCPatient(
            afp_ng_ml=perturb(patient.afp_ng_ml),
            afp_l3_pct=min(100, perturb(patient.afp_l3_pct)),
            dcp_mau_ml=perturb(patient.dcp_mau_ml),
            ctc_per_75ml=patient.ctc_per_75ml,  # integer, not perturbed
            ctdna_vaf_pct=perturb(patient.ctdna_vaf_pct, cv=0.15),
            microvascular_invasion=patient.microvascular_invasion,
            macrovascular_invasion=patient.macrovascular_invasion,
            e_cadherin_loss=patient.e_cadherin_loss,
            vimentin_positive=patient.vimentin_positive,
            ctnnb1_mutation=patient.ctnnb1_mutation,
            nlr=perturb(patient.nlr, cv=0.10),
            pd_l1_tps_pct=min(100, perturb(patient.pd_l1_tps_pct)),
            til_density=patient.til_density,
            vegf_pg_ml=perturb(patient.vegf_pg_ml),
            fgfr_amplification=patient.fgfr_amplification,
            tp53_mutation=patient.tp53_mutation,
            pten_loss=patient.pten_loss,
            rb1_loss=patient.rb1_loss,
            epigenetic_dysregulation=patient.epigenetic_dysregulation,
        )
        nd = compute_domain_scores(noisy)
        sims.append(min(sum(s * w for _, s, _, w in nd), 100.0))

    sims.sort()
    ci_lower = round(sims[int(0.025 * n_simulations)], 1)
    ci_upper = round(sims[int(0.975 * n_simulations)], 1)

    if composite < 20:
        category = "LOW"
    elif composite < 40:
        category = "MODERATE"
    elif composite < 60:
        category = "HIGH"
    else:
        category = "VERY HIGH"

    pathway = pathway_signal_profile(domains)

    # Interpretive notes
    notes = []
    domain_map = {d[0]: d[1] for d in domains}
    if domain_map["tme"] >= 40 and domain_map["angiogenesis"] >= 40:
        notes.append("Elevated immune and angiogenic signals co-occur — consistent with a profile where combined anti-PD-L1 / anti-VEGF mechanisms may be relevant to explore.")
    if domain_map["fgfr"] >= 70:
        notes.append("FGFR amplification/dysregulation detected — the only domain in this framework that mechanistically differentiates lenvatinib from sorafenib.")
    if domain_map["tme"] < 20 and domain_map["angiogenesis"] >= 50:
        notes.append("Non-inflamed TME with high angiogenic burden — profile features are more aligned with anti-angiogenic monotherapy mechanisms in the current framework.")
    if patient.ctnnb1_mutation and domain_map["tme"] < 30:
        notes.append("CTNNB1 mutation with low immune signal — consistent with Wnt-activated HCC phenotype, which has been reported to associate with reduced immune infiltration (Ruiz de Galarreta et al. 2019). This does not exclude immune-based strategies but is noted for research awareness.")
    if domain_map["genetic"] >= 70:
        notes.append("High convergent genetic driver burden (TP53/PTEN/RB1) — associated with aggressive metastatic phenotype; warrants close monitoring regardless of therapeutic direction.")

    return HCCResult(
        composite_score=round(composite, 1),
        ci_lower=ci_lower,
        ci_upper=ci_upper,
        score_category=category,
        domains=[{"domain": d[0], "raw_score": round(d[1], 1), "detail": d[2], "weight": d[3], "weighted": round(d[1]*d[3], 2)} for d in domains],
        pathway_signal=pathway,
        interpretive_notes=notes,
    )


# ─────────────────────────────────────────────
# Output printer
# ─────────────────────────────────────────────

def print_result(result: HCCResult, label: str):
    SEP = "=" * 72
    print(f"\n{SEP}\n{label}\n{SEP}")
    print(f"Composite score: {result.composite_score}/100 [{result.score_category}]")
    print(f"95% CI: [{result.ci_lower}, {result.ci_upper}]  (reflects input measurement variability)")
    print("\nDomain breakdown:")
    for d in result.domains:
        print(f"  {d['domain']:15s} raw={d['raw_score']:5.1f}  weight={d['weight']:.2f}  weighted={d['weighted']:.2f}  | {d['detail']}")
    print("\nPathway signal profile:")
    for agent, (score, level) in result.pathway_signal.items():
        print(f"  {agent}: {score:.1f} [{level}]")
    if result.interpretive_notes:
        print("\nInterpretive notes (hypothesis-generating only):")
        for note in result.interpretive_notes:
            print(f"  * {note}")
    print(f"\n{'─'*72}")
    print("REMINDER: This output is for research focus and discussion only.")
    print("It does not recommend, rank, or eliminate any treatment.")
    print(f"{'─'*72}")


# ─────────────────────────────────────────────
# Demo
# ─────────────────────────────────────────────

def demo():
    scenarios = [
        (
            "Scenario 1 — Early extrahepatic HCC, predominantly angiogenic profile",
            HCCPatient(
                afp_ng_ml=850, afp_l3_pct=22, dcp_mau_ml=180,
                ctc_per_75ml=0, ctdna_vaf_pct=0.8,
                microvascular_invasion=True, macrovascular_invasion=False,
                e_cadherin_loss=True, vimentin_positive=True, ctnnb1_mutation=True,
                nlr=3.2, pd_l1_tps_pct=2.0, til_density="low",
                vegf_pg_ml=310, fgfr_amplification=False,
                tp53_mutation=True, pten_loss=False, rb1_loss=False,
                epigenetic_dysregulation=False,
            ),
        ),
        (
            "Scenario 2 — Advanced HCC with immune-inflamed microenvironment",
            HCCPatient(
                afp_ng_ml=320, afp_l3_pct=8, dcp_mau_ml=95,
                ctc_per_75ml=3, ctdna_vaf_pct=3.2,
                microvascular_invasion=True, macrovascular_invasion=True,
                e_cadherin_loss=True, vimentin_positive=True, ctnnb1_mutation=False,
                nlr=2.1, pd_l1_tps_pct=18.0, til_density="high",
                vegf_pg_ml=190, fgfr_amplification=False,
                tp53_mutation=False, pten_loss=True, rb1_loss=False,
                epigenetic_dysregulation=False,
            ),
        ),
        (
            "Scenario 3 — Aggressive multi-route dissemination, mixed signals",
            HCCPatient(
                afp_ng_ml=12400, afp_l3_pct=41, dcp_mau_ml=880,
                ctc_per_75ml=9, ctdna_vaf_pct=8.7,
                microvascular_invasion=True, macrovascular_invasion=True,
                e_cadherin_loss=True, vimentin_positive=True, ctnnb1_mutation=True,
                nlr=6.8, pd_l1_tps_pct=5.0, til_density="low",
                vegf_pg_ml=520, fgfr_amplification=True,
                tp53_mutation=True, pten_loss=True, rb1_loss=True,
                epigenetic_dysregulation=True,
            ),
        ),
    ]

    for label, patient in scenarios:
        result = compute_hcc_score(patient)
        print_result(result, label)


if __name__ == "__main__":
    demo()