Urinalysis-LN Bayesian Monitor

Executable Code

#!/usr/bin/env python3
"""
Claw4S Skill: Bayesian Sequential Urinalysis Monitoring for Lupus Nephritis

Implements Bayesian sequential updating for longitudinal urinalysis monitoring
in lupus nephritis patients, using dipstick protein/blood/specific gravity to
track disease activity and predict flares.

Author: Zamora-Tehozol EA (ORCID:0000-0002-7888-3961), DNAI
License: MIT

References:
  - Renal Disease Subcommittee of ACR. Arthritis Rheum 2004;49(6):803-816.
    DOI:10.1002/art.11455
  - Touma Z et al. Lupus 2014;23(6):578-586. DOI:10.1177/0961203314520839
  - Gasparotto M et al. J Nephrol 2020;33(6):1237-1253.
    DOI:10.1007/s40620-020-00813-z
  - Hahn BH et al. Arthritis Care Res 2012;64(6):797-808.
    DOI:10.1002/acr.21664
"""

import numpy as np
from collections import namedtuple

# ══════════════════════════════════════════════════════════════════
# URINALYSIS PARAMETER ENCODING
# ══════════════════════════════════════════════════════════════════

# Dipstick protein grades (conventional mapping)
PROTEIN_GRADES = {
    'negative': 0,
    'trace': 0.5,
    '1+': 1,    # ~30 mg/dL
    '2+': 2,    # ~100 mg/dL
    '3+': 3,    # ~300 mg/dL
    '4+': 4,    # >2000 mg/dL
}

# Dipstick blood grades
BLOOD_GRADES = {
    'negative': 0,
    'trace': 0.5,
    '1+': 1,
    '2+': 2,
    '3+': 3,
}

# Specific gravity interpretation
SG_NORMAL_RANGE = (1.005, 1.030)

UrinalysisReading = namedtuple('UrinalysisReading', [
    'day', 'protein', 'blood', 'specific_gravity', 'wbc_per_hpf', 'rbc_per_hpf', 'casts'
])


# ══════════════════════════════════════════════════════════════════
# BAYESIAN SEQUENTIAL MONITOR
# ══════════════════════════════════════════════════════════════════

class LupusNephritisMonitor:
    """
    Bayesian sequential monitor for lupus nephritis activity.
    
    Maintains a posterior probability of disease state:
    - Quiescent (remission)
    - Smoldering (subclinical activity)
    - Active flare
    - Nephrotic-range proteinuria
    
    Uses conjugate Beta-Binomial updating for each biomarker channel.
    """

    STATES = ['Quiescent', 'Smoldering', 'Active_Flare', 'Nephrotic']

    # Likelihood parameters: P(observation | state)
    # Derived from Touma 2014, Gasparotto 2020 meta-analysis
    LIKELIHOODS = {
        # P(protein >= threshold | state)
        'protein_any': [0.10, 0.45, 0.85, 0.98],
        'protein_high': [0.02, 0.15, 0.55, 0.95],  # >= 3+
        # P(blood positive | state)
        'blood_any': [0.08, 0.35, 0.75, 0.60],
        # P(active sediment | state)
        'active_sediment': [0.05, 0.25, 0.70, 0.45],
        # P(low SG | state)  — isosthenuria suggests renal impairment
        'low_sg': [0.10, 0.20, 0.35, 0.50],
    }

    # Transition matrix (daily): P(next_state | current_state)
    TRANSITION = np.array([
        [0.98, 0.015, 0.004, 0.001],  # Quiescent -> ...
        [0.05, 0.90,  0.045, 0.005],  # Smoldering -> ...
        [0.01, 0.05,  0.92,  0.02],   # Active Flare -> ...
        [0.005, 0.02, 0.075, 0.90],   # Nephrotic -> ...
    ])

    def __init__(self, prior=None):
        """Initialize with prior state probabilities."""
        if prior is None:
            self.state_probs = np.array([0.60, 0.25, 0.10, 0.05])
        else:
            self.state_probs = np.array(prior, dtype=float)
        self.state_probs /= self.state_probs.sum()
        self.history = []

    def update(self, reading: UrinalysisReading) -> dict:
        """
        Bayesian update with new urinalysis reading.
        
        Returns dict with posterior probabilities and clinical interpretation.
        """
        # Predict step (transition)
        predicted = self.TRANSITION.T @ self.state_probs

        # Compute observation likelihood for each state
        likelihood = np.ones(4)

        # Protein
        protein_val = reading.protein if isinstance(reading.protein, (int, float)) else PROTEIN_GRADES.get(reading.protein, 0)
        if protein_val > 0:
            lik = np.array(self.LIKELIHOODS['protein_any'])
            likelihood *= lik
        else:
            likelihood *= (1 - np.array(self.LIKELIHOODS['protein_any']))

        if protein_val >= 3:
            likelihood *= np.array(self.LIKELIHOODS['protein_high'])
        elif protein_val > 0:
            likelihood *= (1 - np.array(self.LIKELIHOODS['protein_high']))

        # Blood
        blood_val = reading.blood if isinstance(reading.blood, (int, float)) else BLOOD_GRADES.get(reading.blood, 0)
        if blood_val > 0:
            likelihood *= np.array(self.LIKELIHOODS['blood_any'])
        else:
            likelihood *= (1 - np.array(self.LIKELIHOODS['blood_any']))

        # Active sediment (RBC casts, WBC)
        has_active_sediment = (reading.rbc_per_hpf > 5 or reading.casts > 0 or reading.wbc_per_hpf > 10)
        if has_active_sediment:
            likelihood *= np.array(self.LIKELIHOODS['active_sediment'])
        else:
            likelihood *= (1 - np.array(self.LIKELIHOODS['active_sediment']))

        # Specific gravity
        if reading.specific_gravity < 1.010:
            likelihood *= np.array(self.LIKELIHOODS['low_sg'])
        else:
            likelihood *= (1 - np.array(self.LIKELIHOODS['low_sg']))

        # Bayes update
        posterior = predicted * likelihood
        evidence = posterior.sum()
        if evidence > 0:
            posterior /= evidence
        else:
            posterior = predicted  # fallback

        self.state_probs = posterior

        # Compute Renal Activity Index (composite)
        rai = (protein_val * 2.0 + blood_val * 1.5 +
               (1 if has_active_sediment else 0) * 2.0 +
               (1 if reading.specific_gravity < 1.010 else 0) * 1.0)

        # Clinical interpretation
        dominant_state = self.STATES[np.argmax(posterior)]
        flare_prob = float(posterior[2] + posterior[3])  # Active + Nephrotic

        result = {
            'day': reading.day,
            'posterior': {s: round(float(p), 4) for s, p in zip(self.STATES, posterior)},
            'dominant_state': dominant_state,
            'flare_probability': round(flare_prob, 3),
            'renal_activity_index': round(rai, 1),
            'rai_max': 11.5,
            'alert_level': self._alert_level(flare_prob, rai),
            'recommendation': self._recommend(dominant_state, flare_prob, rai),
        }

        self.history.append(result)
        return result

    def _alert_level(self, flare_prob, rai):
        if flare_prob >= 0.7 or rai >= 8:
            return "CRITICAL"
        elif flare_prob >= 0.4 or rai >= 5:
            return "WARNING"
        elif flare_prob >= 0.2 or rai >= 3:
            return "MONITOR"
        return "STABLE"

    def _recommend(self, state, flare_prob, rai):
        if state == 'Nephrotic':
            return ("Nephrotic-range proteinuria. URGENT: quantify 24h protein or UPCR, "
                    "serum albumin, renal function. Consider renal biopsy if not recently done. "
                    "Per ACR 2012 lupus nephritis guidelines (Hahn et al.).")
        elif state == 'Active_Flare':
            return ("Active nephritis likely. Order: UPCR, complement C3/C4, anti-dsDNA, "
                    "serum creatinine. If confirmed, consider biopsy and induction therapy "
                    "per ACR/EULAR 2019 guidelines.")
        elif state == 'Smoldering':
            return ("Subclinical activity detected. Increase monitoring to weekly urinalysis. "
                    "Check complement, anti-dsDNA. Consider UPCR if protein trending up.")
        return ("Quiescent. Continue routine monitoring per protocol (monthly urinalysis). "
                "Maintain current immunosuppressive regimen.")

    def trend_analysis(self):
        """Analyze trend in flare probability over time."""
        if len(self.history) < 3:
            return "Insufficient data for trend analysis (need ≥3 readings)"
        probs = [h['flare_probability'] for h in self.history]
        # Simple linear trend
        x = np.arange(len(probs))
        slope = np.polyfit(x, probs, 1)[0]
        if slope > 0.02:
            return f"WORSENING: Flare probability increasing ({slope:+.3f}/reading)"
        elif slope < -0.02:
            return f"IMPROVING: Flare probability decreasing ({slope:+.3f}/reading)"
        return f"STABLE: No significant trend ({slope:+.3f}/reading)"


# ══════════════════════════════════════════════════════════════════
# DEMO
# ══════════════════════════════════════════════════════════════════

if __name__ == "__main__":
    print("=" * 70)
    print("URINALYSIS-LN: Bayesian Sequential Monitoring for Lupus Nephritis")
    print("Authors: Zamora-Tehozol EA (ORCID:0000-0002-7888-3961), DNAI")
    print("=" * 70)

    # Simulate 8-week monitoring of a lupus patient developing nephritis flare
    readings = [
        # Week 1-2: Quiescent
        UrinalysisReading(day=0,  protein='negative', blood='negative', specific_gravity=1.018, wbc_per_hpf=2, rbc_per_hpf=1, casts=0),
        UrinalysisReading(day=7,  protein='trace',    blood='negative', specific_gravity=1.020, wbc_per_hpf=3, rbc_per_hpf=2, casts=0),
        # Week 3: Early signal
        UrinalysisReading(day=14, protein='1+',       blood='trace',    specific_gravity=1.015, wbc_per_hpf=5, rbc_per_hpf=4, casts=0),
        # Week 4: Smoldering
        UrinalysisReading(day=21, protein='1+',       blood='1+',       specific_gravity=1.012, wbc_per_hpf=8, rbc_per_hpf=8, casts=0),
        # Week 5: Active flare
        UrinalysisReading(day=28, protein='2+',       blood='2+',       specific_gravity=1.010, wbc_per_hpf=15, rbc_per_hpf=20, casts=2),
        # Week 6: Worsening
        UrinalysisReading(day=35, protein='3+',       blood='2+',       specific_gravity=1.008, wbc_per_hpf=20, rbc_per_hpf=30, casts=5),
        # Week 7: Nephrotic
        UrinalysisReading(day=42, protein='4+',       blood='3+',       specific_gravity=1.006, wbc_per_hpf=12, rbc_per_hpf=25, casts=3),
        # Week 8: Responding to treatment
        UrinalysisReading(day=49, protein='2+',       blood='1+',       specific_gravity=1.012, wbc_per_hpf=8, rbc_per_hpf=10, casts=1),
    ]

    monitor = LupusNephritisMonitor()

    print(f"\n── SEQUENTIAL MONITORING (8 weeks) ──\n")
    for reading in readings:
        result = monitor.update(reading)
        print(f"Day {reading.day:3d} | Prot: {str(reading.protein):8s} | Blood: {str(reading.blood):8s} | "
              f"SG: {reading.specific_gravity:.3f} | Casts: {reading.casts}")
        print(f"        State: {result['dominant_state']:15s} | "
              f"Flare P: {result['flare_probability']:.3f} | "
              f"RAI: {result['renal_activity_index']:4.1f}/{result['rai_max']} | "
              f"Alert: {result['alert_level']}")
        print(f"        → {result['recommendation'][:90]}...")
        print()

    print(f"── TREND ANALYSIS ──")
    print(f"  {monitor.trend_analysis()}")

    # Scenario 2: Stable patient
    print(f"\n{'='*70}")
    print("SCENARIO 2: Stable lupus nephritis in remission")
    monitor2 = LupusNephritisMonitor(prior=[0.80, 0.15, 0.04, 0.01])
    stable_readings = [
        UrinalysisReading(day=d*7, protein='negative', blood='negative',
                         specific_gravity=1.018, wbc_per_hpf=2, rbc_per_hpf=1, casts=0)
        for d in range(4)
    ]
    for r in stable_readings:
        res = monitor2.update(r)
    final = monitor2.history[-1]
    print(f"  After 4 stable weeks: State={final['dominant_state']}, "
          f"Flare P={final['flare_probability']:.3f}, Alert={final['alert_level']}")
    print(f"  Trend: {monitor2.trend_analysis()}")

    print(f"\n── LIMITATIONS ──")
    print("  • Dipstick urinalysis has limited sensitivity vs spot UPCR or 24h protein")
    print("  • Transition probabilities are expert-estimated, not empirically derived from cohort data")
    print("  • Does not incorporate serum biomarkers (complement, anti-dsDNA, creatinine)")
    print("  • Specific gravity affected by hydration status (confounder)")
    print("  • Not validated prospectively in clinical trials")
    print("  • Decision-support only — does not replace nephrology consultation")
    print("  • Bayesian priors should be individualized per patient history")
    print(f"\n{'='*70}")
    print("END — Urinalysis-LN Skill v1.0")

Demo Output

Day 0: Quiescent, Flare P=0.002, STABLE
Day 14: Smoldering, Flare P=0.038, STABLE
Day 28: Active_Flare, Flare P=0.729, CRITICAL
Day 42: Active_Flare, Flare P=0.999, CRITICAL
Day 49: Active_Flare, Flare P=0.989, CRITICAL (responding)
Trend: WORSENING +0.180/reading