# ild-track

Executable Python skill.

Run: `python3 ild_track.py`

## Code

```python
#!/usr/bin/env python3
"""
ILD-TRACK: Interstitial Lung Disease Progression Tracker
Longitudinal FVC/DLCO Decline Modeling with Monte Carlo Uncertainty
for Autoimmune-Associated ILD (SSc-ILD, RA-ILD, Myositis-ILD)

Authors: Erick Adrián Zamora Tehozol, DNAI, Claw 🦞
License: MIT | RheumaAI × Frutero Club
"""

import math
import random
import json
from dataclasses import dataclass, field, asdict
from typing import List, Optional, Tuple
from datetime import datetime, timedelta

# ── Reference Data ──────────────────────────────────────────────────────────

# Annual FVC decline rates (%predicted/year) from literature
# Distler 2019 (SENSCIS), Tashkin 2006 (SLS-I), Solomon 2016, Flaherty 2017
DECLINE_RATES = {
    "SSc-ILD":      {"mean": -5.0, "sd": 2.5},   # Distler NEJM 2019
    "RA-ILD":       {"mean": -3.5, "sd": 2.0},    # Solomon Chest 2016
    "Myositis-ILD": {"mean": -6.5, "sd": 3.0},    # Rapid progressive subset
    "IPAF":         {"mean": -3.0, "sd": 1.8},    # Interstitial pneumonia w/ autoimmune features
    "ANCA-ILD":     {"mean": -2.5, "sd": 1.5},    # MPA-associated
    "Sjogren-ILD":  {"mean": -2.8, "sd": 1.5},
}

# Treatment modification factors (relative reduction of decline)
TREATMENT_EFFECTS = {
    "nintedanib":       0.44,   # 44% reduction — SENSCIS, Distler 2019
    "mycophenolate":    0.50,   # ~50% — SLS-II, Tashkin 2016
    "tocilizumab":      0.60,   # ~60% FVC stabilization — focuSSced, Khanna 2020
    "rituximab":        0.55,   # ~55% — Md Yusof 2017, Lancet Resp Med
    "cyclophosphamide": 0.40,   # ~40% — SLS-I short-term
    "azathioprine":     0.30,   # modest effect
    "pirfenidone":      0.35,   # off-label SSc-ILD data
    "none":             0.0,
}

# Risk factors that accelerate decline
RISK_MODIFIERS = {
    "UIP_pattern":          1.35,  # UIP worse than NSIP — Bouros 2002
    "extensive_disease":    1.40,  # >20% HRCT extent — Goh 2008
    "anti_MDA5":            1.60,  # rapidly progressive — Moghadam-Kia 2017
    "anti_Scl70":           1.25,  # diffuse SSc marker
    "smoking_current":      1.30,
    "smoking_former":       1.10,
    "age_over_65":          1.15,
    "male_sex":             1.10,  # slightly worse prognosis
    "low_baseline_DLCO":    1.30,  # DLCO <40% predicted
    "pulmonary_hypertension": 1.45,
}

# DLCO/FVC ratio thresholds for pulmonary hypertension screening
# Coghlan 2014, Wells 2018
PH_SCREENING = {
    "high_risk": 0.50,    # DLCO/FVC < 0.50 → echo + RHC
    "moderate_risk": 0.70, # 0.50-0.70 → echo screening
}

# Progression criteria (ATS/ERS 2018, Flaherty 2022 INBUILD)
PROGRESSION_CRITERIA = {
    "definite": {"fvc_decline": -10.0, "period_months": 24},  # ≥10% absolute in 24mo
    "probable": {"fvc_decline": -5.0, "dlco_decline": -15.0, "period_months": 12},
    "marginal": {"fvc_decline": -5.0, "period_months": 12},  # 5-10% needs context
}


@dataclass
class PFTMeasurement:
    """Single pulmonary function test measurement."""
    date: str                    # ISO format YYYY-MM-DD
    fvc_percent: float           # FVC % predicted
    dlco_percent: Optional[float] = None  # DLCO % predicted (may be unavailable)
    fev1_fvc_ratio: Optional[float] = None  # to exclude obstructive pattern
    tlc_percent: Optional[float] = None     # total lung capacity


@dataclass
class PatientProfile:
    """Patient clinical profile for ILD tracking."""
    diagnosis: str               # Key from DECLINE_RATES
    age: int
    sex: str                     # "M" or "F"
    measurements: List[PFTMeasurement] = field(default_factory=list)
    treatments: List[str] = field(default_factory=list)  # active treatments
    risk_factors: List[str] = field(default_factory=list)
    hrct_extent_percent: Optional[float] = None  # % lung involvement on HRCT
    autoantibodies: List[str] = field(default_factory=list)  # anti-Scl70, anti-MDA5, etc.


@dataclass
class ProgressionAssessment:
    """Results of ILD progression analysis."""
    observed_fvc_slope: float         # %predicted per year
    observed_dlco_slope: Optional[float]
    predicted_fvc_slope: float        # model-predicted
    fvc_at_12mo: Tuple[float, float, float]  # (mean, CI_low, CI_high)
    fvc_at_24mo: Tuple[float, float, float]
    progression_status: str           # "Stable", "Marginal", "Progressive", "Rapidly Progressive"
    progression_confidence: float     # 0-1
    ph_risk: str                      # PH screening recommendation
    treatment_response: str           # assessment of current therapy
    recommendations: List[str]
    mc_simulations: int
    risk_multiplier: float


def _parse_date(s: str) -> datetime:
    return datetime.strptime(s, "%Y-%m-%d")


def _slope_years(measurements: List[PFTMeasurement], attr: str) -> Optional[float]:
    """Compute annualized slope via OLS for a PFT attribute."""
    points = []
    for m in measurements:
        val = getattr(m, attr, None)
        if val is not None:
            t = _parse_date(m.date)
            points.append((t, val))
    if len(points) < 2:
        return None
    points.sort(key=lambda x: x[0])
    t0 = points[0][0]
    xs = [(p[0] - t0).days / 365.25 for p in points]
    ys = [p[1] for p in points]
    n = len(xs)
    mx = sum(xs) / n
    my = sum(ys) / n
    num = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    den = sum((x - mx) ** 2 for x in xs)
    if den < 1e-9:
        return None
    return num / den  # %predicted per year


def compute_risk_multiplier(profile: PatientProfile) -> float:
    """Aggregate risk factor multiplier."""
    m = 1.0
    for rf in profile.risk_factors:
        if rf in RISK_MODIFIERS:
            m *= RISK_MODIFIERS[rf]
    if profile.age > 65:
        if "age_over_65" not in profile.risk_factors:
            m *= RISK_MODIFIERS["age_over_65"]
    if profile.sex == "M":
        if "male_sex" not in profile.risk_factors:
            m *= RISK_MODIFIERS["male_sex"]
    # HRCT extent
    if profile.hrct_extent_percent is not None and profile.hrct_extent_percent > 20:
        if "extensive_disease" not in profile.risk_factors:
            m *= RISK_MODIFIERS["extensive_disease"]
    # Autoantibodies
    for ab in profile.autoantibodies:
        key = f"anti_{ab}"
        if key in RISK_MODIFIERS and key not in profile.risk_factors:
            m *= RISK_MODIFIERS[key]
    # DLCO check
    if profile.measurements:
        latest = sorted(profile.measurements, key=lambda x: x.date)[-1]
        if latest.dlco_percent is not None and latest.dlco_percent < 40:
            if "low_baseline_DLCO" not in profile.risk_factors:
                m *= RISK_MODIFIERS["low_baseline_DLCO"]
    return m


def compute_treatment_factor(treatments: List[str]) -> float:
    """Combined treatment effect (best single agent — no additive assumption)."""
    if not treatments:
        return 0.0
    effects = [TREATMENT_EFFECTS.get(t.lower(), 0.0) for t in treatments]
    return max(effects) if effects else 0.0


def monte_carlo_projection(
    baseline_fvc: float,
    annual_decline_mean: float,
    annual_decline_sd: float,
    risk_mult: float,
    treatment_factor: float,
    horizon_years: float,
    n_sim: int = 5000,
    seed: Optional[int] = 42
) -> Tuple[float, float, float]:
    """
    Project FVC at horizon_years with Monte Carlo.
    Returns (mean, 2.5th percentile, 97.5th percentile).
    """
    if seed is not None:
        random.seed(seed)

    adjusted_mean = annual_decline_mean * risk_mult * (1.0 - treatment_factor)
    adjusted_sd = annual_decline_sd * math.sqrt(risk_mult)

    results = []
    for _ in range(n_sim):
        rate = random.gauss(adjusted_mean, adjusted_sd)
        projected = baseline_fvc + rate * horizon_years
        projected = max(projected, 0.0)  # can't go below 0
        results.append(projected)

    results.sort()
    mean_val = sum(results) / len(results)
    ci_low = results[int(0.025 * n_sim)]
    ci_high = results[int(0.975 * n_sim)]
    return (round(mean_val, 1), round(ci_low, 1), round(ci_high, 1))


def assess_ph_risk(measurements: List[PFTMeasurement]) -> str:
    """Screen for pulmonary hypertension via DLCO/FVC ratio."""
    latest = sorted(measurements, key=lambda x: x.date)[-1]
    if latest.dlco_percent is None or latest.fvc_percent < 1:
        return "Unable to assess (DLCO unavailable)"
    ratio = latest.dlco_percent / latest.fvc_percent
    if ratio < PH_SCREENING["high_risk"]:
        return "HIGH RISK — recommend echocardiogram + right heart catheterization (DLCO/FVC ratio < 0.50)"
    elif ratio < PH_SCREENING["moderate_risk"]:
        return "MODERATE RISK — recommend screening echocardiogram (DLCO/FVC ratio 0.50-0.70)"
    return "LOW RISK — routine monitoring (DLCO/FVC ratio > 0.70)"


def classify_progression(
    obs_fvc_slope: Optional[float],
    obs_dlco_slope: Optional[float],
    span_months: float
) -> Tuple[str, float]:
    """
    Classify ILD progression per ATS/ERS criteria.
    Returns (status, confidence).
    """
    if obs_fvc_slope is None:
        return ("Insufficient data", 0.0)

    annualized = obs_fvc_slope  # already annualized

    if annualized <= -10.0:
        return ("Rapidly Progressive", 0.95)
    elif annualized <= -5.0:
        # Check if DLCO also declining
        if obs_dlco_slope is not None and obs_dlco_slope <= -15.0:
            return ("Progressive (FVC + DLCO decline)", 0.90)
        conf = min(0.85, 0.60 + span_months / 48.0)
        return ("Progressive", conf)
    elif annualized <= -2.0:
        return ("Marginal decline — close monitoring", 0.65)
    elif annualized <= 2.0:
        return ("Stable", 0.80)
    else:
        return ("Improving", 0.70)


def assess_treatment_response(
    obs_slope: Optional[float],
    expected_untreated: float,
    treatment_factor: float
) -> str:
    """Evaluate if treatment is working as expected."""
    if obs_slope is None or treatment_factor == 0:
        return "No treatment or insufficient data"
    expected_treated = expected_untreated * (1.0 - treatment_factor)
    if obs_slope > expected_treated + 1.0:
        return "BETTER than expected — treatment appears effective"
    elif obs_slope > expected_treated - 2.0:
        return "AS EXPECTED — treatment response adequate"
    else:
        return "WORSE than expected — consider treatment escalation or switch"


def generate_recommendations(
    profile: PatientProfile,
    progression: str,
    ph_risk: str,
    tx_response: str,
    fvc_12mo: Tuple[float, float, float]
) -> List[str]:
    """Generate evidence-based clinical recommendations."""
    recs = []
    if "Rapidly Progressive" in progression or "Progressive" in progression:
        if not any(t in profile.treatments for t in ["nintedanib", "pirfenidone"]):
            recs.append("Consider antifibrotic therapy (nintedanib — SENSCIS/INBUILD evidence)")
        if "mycophenolate" not in profile.treatments and profile.diagnosis == "SSc-ILD":
            recs.append("Consider mycophenolate mofetil (SLS-II, Tashkin 2016)")
        if profile.diagnosis == "Myositis-ILD" and "rituximab" not in profile.treatments:
            recs.append("Consider rituximab for myositis-ILD (Md Yusof 2017)")
        recs.append("Repeat HRCT in 3-6 months to assess radiographic progression")
        recs.append("PFTs every 3 months during active decline")

    if "Marginal" in progression:
        recs.append("PFTs every 3-4 months — marginal decline requires close monitoring")
        recs.append("Repeat HRCT in 6 months if decline continues")

    if "Stable" in progression:
        recs.append("Continue current management — PFTs every 6 months")

    if "HIGH RISK" in ph_risk:
        recs.append("URGENT: Echocardiogram + RHC for pulmonary hypertension evaluation")
        recs.append("Consider referral to PH specialist if confirmed")
    elif "MODERATE" in ph_risk:
        recs.append("Screening echocardiogram recommended for PH")

    if "WORSE" in tx_response:
        recs.append("Treatment response suboptimal — multidisciplinary discussion recommended")
        if "tocilizumab" not in profile.treatments and profile.diagnosis == "SSc-ILD":
            recs.append("Consider tocilizumab (focuSSced trial, Khanna 2020)")

    if fvc_12mo[0] < 50:
        recs.append("FVC projected <50% — evaluate for lung transplant referral")
        recs.append("6-minute walk test and cardiopulmonary exercise testing recommended")

    if not profile.treatments or profile.treatments == ["none"]:
        recs.append("Patient currently untreated — immunosuppression evaluation needed")

    # Oxygen assessment
    if fvc_12mo[0] < 55:
        recs.append("Assess supplemental oxygen need (resting + exertional SpO2)")

    return recs if recs else ["Continue current management with routine monitoring"]


def analyze(profile: PatientProfile, n_sim: int = 5000) -> ProgressionAssessment:
    """
    Main analysis function.
    Takes patient profile with longitudinal PFTs, returns full assessment.
    """
    # Validate
    if profile.diagnosis not in DECLINE_RATES:
        raise ValueError(f"Unknown diagnosis: {profile.diagnosis}. Options: {list(DECLINE_RATES.keys())}")
    if len(profile.measurements) < 1:
        raise ValueError("At least 1 PFT measurement required")

    # Sort measurements
    profile.measurements.sort(key=lambda m: m.date)

    # Observed slopes
    obs_fvc = _slope_years(profile.measurements, "fvc_percent")
    obs_dlco = _slope_years(profile.measurements, "dlco_percent")

    # Time span
    if len(profile.measurements) >= 2:
        d0 = _parse_date(profile.measurements[0].date)
        d1 = _parse_date(profile.measurements[-1].date)
        span_months = (d1 - d0).days / 30.44
    else:
        span_months = 0

    # Reference decline
    ref = DECLINE_RATES[profile.diagnosis]

    # Risk & treatment
    risk_mult = compute_risk_multiplier(profile)
    tx_factor = compute_treatment_factor(profile.treatments)

    # Latest FVC as baseline for projection
    latest_fvc = profile.measurements[-1].fvc_percent

    # Monte Carlo projections
    fvc_12 = monte_carlo_projection(latest_fvc, ref["mean"], ref["sd"], risk_mult, tx_factor, 1.0, n_sim)
    fvc_24 = monte_carlo_projection(latest_fvc, ref["mean"], ref["sd"], risk_mult, tx_factor, 2.0, n_sim)

    # Predicted slope
    predicted_slope = ref["mean"] * risk_mult * (1.0 - tx_factor)

    # Progression classification
    prog_status, prog_conf = classify_progression(obs_fvc, obs_dlco, span_months)

    # PH risk
    ph = assess_ph_risk(profile.measurements)

    # Treatment response
    tx_resp = assess_treatment_response(obs_fvc, ref["mean"] * risk_mult, tx_factor)

    # Recommendations
    recs = generate_recommendations(profile, prog_status, ph, tx_resp, fvc_12)

    return ProgressionAssessment(
        observed_fvc_slope=round(obs_fvc, 2) if obs_fvc else 0.0,
        observed_dlco_slope=round(obs_dlco, 2) if obs_dlco else None,
        predicted_fvc_slope=round(predicted_slope, 2),
        fvc_at_12mo=fvc_12,
        fvc_at_24mo=fvc_24,
        progression_status=prog_status,
        progression_confidence=round(prog_conf, 2),
        ph_risk=ph,
        treatment_response=tx_resp,
        recommendations=recs,
        mc_simulations=n_sim,
        risk_multiplier=round(risk_mult, 3),
    )


def format_report(profile: PatientProfile, result: ProgressionAssessment) -> str:
    """Format a clinical report."""
    lines = [
        "=" * 70,
        "  ILD-TRACK: Interstitial Lung Disease Progression Report",
        "  RheumaAI × Frutero Club",
        "=" * 70,
        f"\n  Diagnosis: {profile.diagnosis}",
        f"  Age/Sex: {profile.age}{profile.sex}",
        f"  Active treatments: {', '.join(profile.treatments) if profile.treatments else 'None'}",
        f"  Risk factors: {', '.join(profile.risk_factors) if profile.risk_factors else 'None'}",
        f"  HRCT extent: {profile.hrct_extent_percent}%" if profile.hrct_extent_percent else "",
        f"\n  Measurements ({len(profile.measurements)} timepoints):",
    ]
    for m in profile.measurements:
        dlco_str = f", DLCO {m.dlco_percent}%" if m.dlco_percent else ""
        lines.append(f"    {m.date}: FVC {m.fvc_percent}%{dlco_str}")

    lines += [
        f"\n{'─' * 70}",
        f"  OBSERVED TRENDS",
        f"    FVC slope: {result.observed_fvc_slope:+.2f} %predicted/year",
        f"    DLCO slope: {result.observed_dlco_slope:+.2f} %predicted/year" if result.observed_dlco_slope else "    DLCO slope: N/A",
        f"\n  MODEL PREDICTION (risk multiplier: {result.risk_multiplier}x)",
        f"    Expected FVC slope: {result.predicted_fvc_slope:+.2f} %predicted/year",
        f"    FVC at 12 months: {result.fvc_at_12mo[0]}% (95% CI: {result.fvc_at_12mo[1]}–{result.fvc_at_12mo[2]}%)",
        f"    FVC at 24 months: {result.fvc_at_24mo[0]}% (95% CI: {result.fvc_at_24mo[1]}–{result.fvc_at_24mo[2]}%)",
        f"    Monte Carlo simulations: {result.mc_simulations}",
        f"\n{'─' * 70}",
        f"  PROGRESSION STATUS: {result.progression_status}",
        f"  Confidence: {result.progression_confidence:.0%}",
        f"\n  PH Screening: {result.ph_risk}",
        f"  Treatment Response: {result.treatment_response}",
        f"\n{'─' * 70}",
        f"  RECOMMENDATIONS:",
    ]
    for i, r in enumerate(result.recommendations, 1):
        lines.append(f"    {i}. {r}")

    lines += [
        f"\n{'─' * 70}",
        "  References:",
        "    Distler O et al. NEJM 2019;380:2518-28 (SENSCIS)",
        "    Tashkin DP et al. NEJM 2006;354:2655-66 (SLS-I)",
        "    Tashkin DP et al. Lancet Resp Med 2016;4:708-19 (SLS-II)",
        "    Khanna D et al. Lancet 2020;395:1407-18 (focuSSced)",
        "    Flaherty KR et al. NEJM 2019;381:1718-27 (INBUILD)",
        "    Goh NS et al. AJRCCM 2008;177:1248-54 (staging system)",
        "    Coghlan JG et al. Ann Rheum Dis 2014;73:1340-49 (DETECT)",
        "=" * 70,
    ]
    return "\n".join(lines)


# ── Demo Scenarios ──────────────────────────────────────────────────────────

def demo():
    """Run 3 clinical scenarios."""

    # Scenario 1: SSc-ILD on mycophenolate, progressive
    print("\n" + "▶ SCENARIO 1: SSc-ILD with progressive decline on mycophenolate")
    p1 = PatientProfile(
        diagnosis="SSc-ILD",
        age=52, sex="F",
        measurements=[
            PFTMeasurement("2025-01-15", 78.0, 65.0),
            PFTMeasurement("2025-06-20", 74.0, 60.0),
            PFTMeasurement("2025-12-10", 69.0, 55.0),
            PFTMeasurement("2026-03-15", 66.0, 50.0),
        ],
        treatments=["mycophenolate"],
        risk_factors=["UIP_pattern", "extensive_disease"],
        hrct_extent_percent=35,
        autoantibodies=["Scl70"],
    )
    r1 = analyze(p1)
    print(format_report(p1, r1))

    # Scenario 2: RA-ILD stable on nintedanib
    print("\n" + "▶ SCENARIO 2: RA-ILD stable on nintedanib")
    p2 = PatientProfile(
        diagnosis="RA-ILD",
        age=68, sex="M",
        measurements=[
            PFTMeasurement("2025-03-01", 62.0, 55.0),
            PFTMeasurement("2025-09-15", 61.0, 54.0),
            PFTMeasurement("2026-03-10", 60.5, 53.0),
        ],
        treatments=["nintedanib"],
        risk_factors=["smoking_former"],
        hrct_extent_percent=15,
    )
    r2 = analyze(p2)
    print(format_report(p2, r2))

    # Scenario 3: Myositis-ILD anti-MDA5 rapidly progressive, untreated
    print("\n" + "▶ SCENARIO 3: Anti-MDA5 myositis-ILD, rapidly progressive, untreated")
    p3 = PatientProfile(
        diagnosis="Myositis-ILD",
        age=38, sex="F",
        measurements=[
            PFTMeasurement("2026-01-05", 85.0, 70.0),
            PFTMeasurement("2026-02-10", 72.0, 58.0),
            PFTMeasurement("2026-03-18", 60.0, 45.0),
        ],
        treatments=[],
        risk_factors=[],
        autoantibodies=["MDA5"],
    )
    r3 = analyze(p3)
    print(format_report(p3, r3))

    print("\n✅ All 3 scenarios completed successfully")


if __name__ == "__main__":
    demo()

```