bioavailability.py

"""
bioavailability.py — SimBiology BA 계산 로직 이식 모듈 (문헌 5모델 정렬판)
========================================================================
SimBiology + MATLAB 피팅 워크플로를 순수 Python으로 재현.
대리모델(MLP/Neural ODE)이 예측한 림프·혈관 흡수곡선으로부터
생체이용률(Bioavailability, BA)을 계산한다.

이 버전은 5개 ODE 모델의 전신(systemic) 구획 구조를 Kagan (2014,
DMD 42:1890) 종설의 그림과 "동일하게" 맞춘 것이다.
  model 1 = Fig 1A  단일경로, 전신 1구획
  model 2 = Fig 2A  이중경로(림프 구획 없음), 전신 1구획
  model 3 = Fig 2B  이중경로(별도 림프 구획), 전신 2구획  ← 사용자 PPT/SimBiology 모델
  model 4 = Fig 3A  재분포(말초조직→림프), 전신 2구획
  model 5 = Fig 3B  7구획 비선형(Dahlberg 2014, trastuzumab)

흡수 구동 방식(사용자 선택 = 곡선 구동):
  models 1–4 는 대리모델이 예측한 누적 흡수곡선(c_lymph, c_vessel)을
  이중지수로 피팅→미분한 시간의존 질량유입률 kle(t)/kve(t)를 전신 ODE에
  소스항으로 주입한다. 문헌의 ka·SC(t)(예: Fig 2A/2B eq.5,9)에서 SC
  저장고를 명시적으로 풀지 않고 흡수를 곡선으로 미리 계산해 넣은 등가 형태다.
  presystemic 손실(문헌 k_loss / PPT k_decay)은 흡수곡선에 이미 반영된다.

  model 5 (Fig 3B)는 흡수 자체가 명시적 SC 구획에서의 k14 + Michaelis–
  Menten(Vmax,Km) 수송이라 곡선 구동과 호환되지 않는다. 따라서 문헌 구조
  그대로 명시적 SC 구획 + 속도상수 구동으로 적분한다(아래 주석 참조).

워크플로 (원본 SimBiology 검증: BA=96.8% 재현 확인):
  1. 누적 흡수곡선(c_lymph, c_vessel)을 이중지수로 피팅
     y(t) = y0 + (y00-y0)*F*(1-exp(-k_fast*t)) + (y00-y0)*(1-F)*(1-exp(-k_slow*t))
  2. mass rate 유도: [A*exp(-B*t) + C*exp(-D*t)]*dose*0.01
     A=(y00-y0)*F*k_fast, B=k_fast, C=(y00-y0)*(1-F)*k_slow, D=k_slow
  3. SC / IV PK ODE 적분 (모델별 구획 구조)
  4. BA = AUC_SC(central) / AUC_IV(central) * 100

SimBiology 모델 파라미터 (피하주사 SC 모델, PPT '241127 to SG'):
  구획: Central(3.557 L) / Peripheral(1.807 L) / Central lymph(0.312 L) / SC(0.0031 L)
  k_12=0.0992  중앙→말초
  k_21=0.3448  말초→중앙
  k_input=0.1920  림프→중앙
  k_10=0.0043  중앙→제거(청소)
  k_decay      SC 저장고 분해(=문헌 k_loss); 곡선 구동에서는 곡선에 내포
  흐름: 혈관흡수→중앙 직접 / 림프흡수→림프구획→(k_input)→중앙
"""

import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import odeint

# numpy 버전 호환: 2.0+는 trapezoid, 이전은 trapz
_trapz = np.trapezoid if hasattr(np, "trapezoid") else np.trapz


# ── SimBiology 모델 파라미터 (PPT '241127 to SG') ──
K_12     = 0.0992    # 중앙 → 말초 (1/hr)
K_21     = 0.3448    # 말초 → 중앙 (1/hr)
K_INPUT  = 0.1920    # 림프 → 중앙 (1/hr)
K_10     = 0.0043    # 중앙 → 제거 (1/hr)
DEFAULT_DOSE = 40.0  # mg
T_MAX_HR = 1000.0    # 적분 구간 (hr, ≈42일)
N_STEPS  = 20000

# ── model 4 (Fig 3A) 재분포 파라미터 ──
# 전신 말초조직 → 림프계로의 재배액(redistribution) 1차 상수.
# Kagan et al. (2007)이 제안한 경로. 값은 가정값(미검증).
K_RL_DEFAULT = 0.001

# ── model 5 (Fig 3B, Dahlberg 2014) 파라미터 ──
# 논문 본문은 개별 수치를 표로 제공하지 않으므로 아래는 구조 시연용
# 가정 기본값(placeholder, 미검증)이다. 논문의 대칭 제약을 그대로 반영:
#   k45 = k67,  k36 = k34,  k24 = k26,  k52 = k72
M5_DEFAULTS = dict(
    Vmax=5.0,     # SC→중앙 Michaelis–Menten 최대속도 (mg/hr) [가정]
    Km=1.0,       # MM 상수 (mg) [가정]
    k14=0.05,     # SC → 후방 말초림프(comp4) (1/hr) [가정]
    CLd=0.10,     # 중앙 ↔ 말초 분포청소 상당 1차상수 (1/hr) [가정]
    k24=0.02,     # 중앙 → 말초림프(후방 k24 = 전방 k26) (1/hr) [가정]
    k34=0.02,     # 말초 → 말초림프(후방 k34 = 전방 k36) (1/hr) [가정]
    k45=0.30,     # 말초림프 → 중앙림프(후방 k45 = 전방 k67) (1/hr) [가정]
    k52=0.20,     # 중앙림프 → 중앙구획(후방 k52 = 전방 k72) (1/hr) [가정]
)


def biexp_cumulative(t, y0, y00, F, k_fast, k_slow):
    """이중지수 누적 흡수곡선 (MATLAB fittingpara_v3.m와 동일)."""
    return (y0 + (y00 - y0) * F * (1 - np.exp(-k_fast * t))
            + (y00 - y0) * (1 - F) * (1 - np.exp(-k_slow * t)))


def fit_biexp(t_hr, y):
    """
    누적곡선을 이중지수로 피팅. 여러 초기값을 시도해 robust하게.
    MATLAB 제약: Lower=[-1,0,0,0,0] (y0,y00,F,k_fast,k_slow), F∈[0,1].
    반환: (y0, y00, F, k_fast, k_slow) 또는 None(실패).
    """
    y = np.clip(np.asarray(y, float), 0, None)
    lower = [-1, 0, 0, 0, 0]
    upper = [np.inf, np.inf, 1, np.inf, np.inf]
    best, best_sse = None, np.inf
    for k_fast0 in (0.1, 1.0, 8.0):
        for y00_0 in (max(y.max(), 1.0), max(y[-1], 1.0)):
            try:
                p0 = [0, y00_0, 0.01, k_fast0, 0.01]
                popt, _ = curve_fit(biexp_cumulative, t_hr, y, p0=p0,
                                    bounds=(lower, upper), maxfev=5000)
                sse = np.sum((y - biexp_cumulative(t_hr, *popt)) ** 2)
                if sse < best_sse:
                    best_sse, best = sse, popt
            except Exception:
                continue
    return best


def rate_params(popt):
    """피팅 파라미터 → mass rate 계수 (A,B,C,D)."""
    y0, y00, F, k_fast, k_slow = popt
    A = (y00 - y0) * F * k_fast
    B = k_fast
    C = (y00 - y0) * (1 - F) * k_slow
    D = k_slow
    return A, B, C, D


# ══════════════════════════════════════════════════════════════════
#  SC 흡수 모델 ODE (문헌 5모델과 동일 구조)
# ══════════════════════════════════════════════════════════════════

def _sc_ode_m1(y, t, ka, k10):
    """
    모델 1 = 논문 Fig 1A (single-pathway).
    림프·혈관 미구분, 단일 흡수. 전신 1구획(회색 박스).
    문헌 Fig 1A: 흡수는 ka·F로 정성 기술(식 번호 없음).
      dC/dt = (흡수 유입) − k10·C
    곡선 구동: 흡수 유입 = ka(t) = (림프+혈관 합산 곡선)의 유입률.
    """
    (c_cen,) = y
    return [ka(t) - k10 * c_cen]


def _sc_ode_m2(y, t, kle, kve, k10):
    """
    모델 2 = 논문 Fig 2A (dual-pathway, 림프 구획 없음).
    혈관·림프 두 경로가 둘 다 전신으로 직행. 전신 1구획(회색 박스).
    문헌 Fig 2A(McLennan 2003, 렙틴): dC/dt = (k_blood+k_lymph)·SC/Vc ± disp,
    k_loss는 주사부위 presystemic 분해(곡선에 내포).
      dC/dt = kve(t) + kle(t) − k10·C
    """
    (c_cen,) = y
    return [kve(t) + kle(t) - k10 * c_cen]


def _sc_ode_m3(y, t, kle, kve, k10, k12, k21, k_input):
    """
    모델 3 = 논문 Fig 2B (별도 림프 구획) = 사용자 PPT/SimBiology 모델.
    혈관→중앙 직접, 림프→림프구획→(k_input)→중앙. 말초⇄중앙(k12,k21). 전신 2구획.

    문헌 Fig 2B (McLennan 2005a, Kota 2007) eq.9–11 및 PPT '241127 to SG':
      dm_SC/dt = −ṁ_le − ṁ_ve − k_decay·m_SC          (곡선 구동에서는 생략)
      dm_CL/dt = ṁ_le − k_input·m_CL
      dm_CC/dt = ṁ_ve − (k12·m_CC − k21·m_PC) + k_input·m_CL − k10·m_CC
      dm_PC/dt = k12·m_CC − k21·m_PC
    여기서 ṁ_ve = kve(t), ṁ_le = kle(t) (곡선 유입률).
    상태 순서: [중앙 CC, 말초 PC, 중앙림프 CL].
    """
    c_cen, c_per, c_lym = y
    dc_cen = kve(t) + k_input * c_lym - k10 * c_cen - k12 * c_cen + k21 * c_per
    dc_per = k12 * c_cen - k21 * c_per
    dc_lym = kle(t) - k_input * c_lym
    return [dc_cen, dc_per, dc_lym]


# 하위호환 별칭 (구버전 이름)
_sc_ode = _sc_ode_m3
_sc_ode_dyn = _sc_ode_m3


def _sc_ode_m4(y, t, kle, kve, k10, k12, k21, k_input, k_rl=K_RL_DEFAULT):
    """
    모델 4 = 논문 Fig 3A (redistribution, Kagan 2007).
    Fig 2B 구조 + '전신 말초조직 → 림프계' 재배액(재순환).
    문헌 Fig 3A 도식: SC→림프계·혈장, 림프계→혈장(k_input),
    혈장 ↔ 말초조직(k12,k21), **말초조직 → 림프계(k_rl)**, 혈장 소실(k10).
    전신 2구획(혈장 + 말초조직).

      dc_lym = kle(t) − k_input·c_lym + k_rl·c_per      (말초조직→림프 재분포)
      dc_cen = kve(t) + k_input·c_lym − k10·c_cen − k12·c_cen + k21·c_per
      dc_per = k12·c_cen − k21·c_per − k_rl·c_per
    상태 순서: [중앙 c_cen, 말초조직 c_per, 림프계 c_lym].
    """
    c_cen, c_per, c_lym = y
    dc_cen = kve(t) + k_input * c_lym - k10 * c_cen - k12 * c_cen + k21 * c_per
    dc_per = k12 * c_cen - k21 * c_per - k_rl * c_per
    dc_lym = kle(t) - k_input * c_lym + k_rl * c_per
    return [dc_cen, dc_per, dc_lym]


def _sc_ode_m5(y, t, k10, p):
    """
    모델 5 = 논문 Fig 3B (Dahlberg 2014, trastuzumab, 7구획 비선형).
    구획 번호(논문과 동일):
      1: SC 주사부위
      2: Central(혈장)     3: Peripheral(말초조직)
      4: 후방 말초림프      5: 후방 중앙림프
      6: 전방 말초림프      7: 전방 중앙림프
    흡수: SC(1)→후방 말초림프(4) 1차(k14) + SC(1)→중앙(2) Michaelis–Menten(Vmax,Km).
    재분포: 중앙(2)→말초림프(k24=k26), 말초(3)→말초림프(k34=k36).
    림프 배출: 말초림프→중앙림프(k45=k67), 중앙림프→중앙(k52=k72).
    전신 분포: 중앙 ↔ 말초 (CLd), 중앙 제거 k10(=CL/V).
    대칭 제약(논문): k45=k67, k36=k34, k24=k26, k52=k72.

      dm1 = −k14·m1 − Vmax·m1/(Km+m1)
      dm2 = Vmax·m1/(Km+m1) + k52·m5 + k72·m7 − k10·m2 − CLd·m2 + CLd·m3 − (k24+k26)·m2
      dm3 = CLd·m2 − CLd·m3 − (k34+k36)·m3
      dm4 = k14·m1 + k24·m2 + k34·m3 − k45·m4
      dm5 = k45·m4 − k52·m5
      dm6 = k26·m2 + k36·m3 − k67·m6
      dm7 = k67·m6 − k72·m7

    주의: 이 모델의 흡수는 명시적 SC 구획 + MM 수송이라 대리모델 곡선
    구동과 호환되지 않아 속도상수 구동으로 적분한다. 또한 대리모델은 림프
    출력이 단일 채널뿐이라 전방/후방 림프계를 실측 분리할 수 없다(식별 불가).
    아래 M5_DEFAULTS는 구조 시연용 가정값(placeholder)이다.
    """
    m1, m2, m3, m4, m5, m6, m7 = y
    Vmax, Km = p["Vmax"], p["Km"]
    k14, CLd = p["k14"], p["CLd"]
    k24 = k26 = p["k24"]      # 대칭 제약
    k34 = k36 = p["k34"]
    k45 = k67 = p["k45"]
    k52 = k72 = p["k52"]
    mm = Vmax * m1 / (Km + m1)
    dm1 = -k14 * m1 - mm
    dm2 = mm + k52 * m5 + k72 * m7 - k10 * m2 - CLd * m2 + CLd * m3 - (k24 + k26) * m2
    dm3 = CLd * m2 - CLd * m3 - (k34 + k36) * m3
    dm4 = k14 * m1 + k24 * m2 + k34 * m3 - k45 * m4
    dm5 = k45 * m4 - k52 * m5
    dm6 = k26 * m2 + k36 * m3 - k67 * m6
    dm7 = k67 * m6 - k72 * m7
    return [dm1, dm2, dm3, dm4, dm5, dm6, dm7]


# ══════════════════════════════════════════════════════════════════
#  IV 모델 ODE (SC 모델의 전신 구획 구조와 일치)
# ══════════════════════════════════════════════════════════════════

def _iv_ode_1pool(y, t, k10):
    """1구획 IV (모델 1,2용). t=0 전량 투여 후 1차 제거."""
    (c_cen,) = y
    return [-k10 * c_cen]


def _iv_ode_2pool(y, t, k10, k12, k21):
    """2구획 IV (모델 3,4용): 중앙 ↔ 말초, 중앙 제거."""
    c_cen, c_per = y
    return [-k10 * c_cen - k12 * c_cen + k21 * c_per,
            k12 * c_cen - k21 * c_per]


def _iv_ode_m5(y, t, k10, p):
    """
    모델 5 IV: SC 구획(1) 없이 IV dose를 중앙(2)에 직접 투여.
    나머지 6구획(2,3,4,5,6,7)의 전신 분포·재분포·림프 배출은 SC와 동일.
    상태 순서: [m2, m3, m4, m5, m6, m7].
    """
    m2, m3, m4, m5, m6, m7 = y
    CLd = p["CLd"]
    k24 = k26 = p["k24"]
    k34 = k36 = p["k34"]
    k45 = k67 = p["k45"]
    k52 = k72 = p["k52"]
    dm2 = k52 * m5 + k72 * m7 - k10 * m2 - CLd * m2 + CLd * m3 - (k24 + k26) * m2
    dm3 = CLd * m2 - CLd * m3 - (k34 + k36) * m3
    dm4 = k24 * m2 + k34 * m3 - k45 * m4
    dm5 = k45 * m4 - k52 * m5
    dm6 = k26 * m2 + k36 * m3 - k67 * m6
    dm7 = k67 * m6 - k72 * m7
    return [dm2, dm3, dm4, dm5, dm6, dm7]


# 하위호환 별칭
_iv_ode = lambda y, t: _iv_ode_2pool(y, t, K_10, K_12, K_21)
_iv_ode_dyn = _iv_ode_2pool


def compute_bioavailability(c_lymph, c_vessel, t_min, dose=DEFAULT_DOSE,
                            drug="IgG", model=3, m5_params=None,
                            return_detail=False):
    """
    림프·혈관 누적 흡수곡선 → 생체이용률(%).

    인자:
      c_lymph, c_vessel : 누적 % mass 곡선 (대리모델 예측)
      t_min             : 시간축 (분)
      dose              : 투여량 (mg)
      drug              : 약물 종류 ("IgG","INS","ALB") 또는 MW(숫자)
      model             : 1=Fig1A, 2=Fig2A, 3=Fig2B(PPT), 4=Fig3A, 5=Fig3B
      m5_params         : model 5 속도상수 dict (없으면 M5_DEFAULTS)
    반환:
      BA (%)  또는  return_detail=True 시 dict
    """
    # 약물별 전신 PK 파라미터 로드
    try:
        from drug_pk_params import get_drug_params
        pk = get_drug_params(drug)
        k12, k21 = pk["k12"], pk["k21"]
        k_input, k10 = pk["k_input"], pk["k10"]
    except Exception:
        k12, k21, k_input, k10 = K_12, K_21, K_INPUT, K_10
    t_hr = np.asarray(t_min, float) / 60.0

    pl = fit_biexp(t_hr, c_lymph)
    pv = fit_biexp(t_hr, c_vessel)
    if pl is None or pv is None:
        raise RuntimeError("이중지수 피팅 실패 — 곡선을 확인하세요.")

    Al, Bl, Cl, Dl = rate_params(pl)
    Av, Bv, Cv, Dv = rate_params(pv)

    def kle(t):
        return (Al * np.exp(-Bl * t) + Cl * np.exp(-Dl * t)) * dose * 0.01

    def kve(t):
        return (Av * np.exp(-Bv * t) + Cv * np.exp(-Dv * t)) * dose * 0.01

    # 모델 1은 림프+혈관 합산 흡수를 따로 피팅
    pt = fit_biexp(t_hr, np.asarray(c_lymph, float) + np.asarray(c_vessel, float))
    if pt is not None:
        At, Bt, Ct, Dt = rate_params(pt)
        def ka(t):
            return (At * np.exp(-Bt * t) + Ct * np.exp(-Dt * t)) * dose * 0.01
    else:
        ka = lambda t: kle(t) + kve(t)

    tt = np.linspace(0, T_MAX_HR, N_STEPS)

    # ── 모델 선택 (Kagan 2014 기준) ──
    if model == 1:
        # Fig 1A: 단일흡수, 전신 1구획
        sc_full = odeint(_sc_ode_m1, [0], tt, args=(ka, k10))
        sc = np.column_stack([sc_full[:, 0], np.zeros(len(tt)), np.zeros(len(tt))])
        iv = odeint(_iv_ode_1pool, [dose], tt, args=(k10,))
    elif model == 2:
        # Fig 2A: dual, 림프구획 없음, 전신 1구획
        sc_full = odeint(_sc_ode_m2, [0], tt, args=(kle, kve, k10))
        sc = np.column_stack([sc_full[:, 0], np.zeros(len(tt)), np.zeros(len(tt))])
        iv = odeint(_iv_ode_1pool, [dose], tt, args=(k10,))
    elif model == 4:
        # Fig 3A: redistribution(말초조직→림프), 전신 2구획
        sc = odeint(_sc_ode_m4, [0, 0, 0], tt,
                    args=(kle, kve, k10, k12, k21, k_input, K_RL_DEFAULT))
        iv = odeint(_iv_ode_2pool, [dose, 0], tt, args=(k10, k12, k21))
    elif model == 5:
        # Fig 3B: 7구획 비선형. 속도상수 구동(곡선 구동과 비호환).
        p = dict(M5_DEFAULTS)
        if m5_params:
            p.update(m5_params)
        sc7 = odeint(_sc_ode_m5, [dose, 0, 0, 0, 0, 0, 0], tt, args=(k10, p))
        # 중앙(comp2)만 BA에 사용
        sc = np.column_stack([sc7[:, 1], sc7[:, 2], np.zeros(len(tt))])
        iv6 = odeint(_iv_ode_m5, [dose, 0, 0, 0, 0, 0], tt, args=(k10, p))
        iv = np.column_stack([iv6[:, 0], iv6[:, 1]])
    else:  # model == 3 (Fig 2B, PPT/SimBiology, 기본)
        sc = odeint(_sc_ode_m3, [0, 0, 0], tt,
                    args=(kle, kve, k10, k12, k21, k_input))
        iv = odeint(_iv_ode_2pool, [dose, 0], tt, args=(k10, k12, k21))

    auc_sc = _trapz(sc[:, 0], tt)
    auc_iv = _trapz(iv[:, 0], tt)
    BA = auc_sc / auc_iv * 100.0

    if return_detail:
        sc_plasma = sc[:, 0]
        cmax = float(sc_plasma.max())
        peak_i = int(np.argmax(sc_plasma))
        tmax = float(tt[peak_i])
        after = sc_plasma[peak_i:]
        t_after = tt[peak_i:]
        if len(after) > 1 and cmax > 0:
            t_half = float(t_after[int(np.argmin(np.abs(after - cmax / 2)))] - tmax)
        else:
            t_half = float("nan")
        return {
            "BA": float(BA),
            "AUC_SC": float(auc_sc),
            "AUC_IV": float(auc_iv),
            "drug": drug,
            "model": model,
            "pk_params": dict(k10=k10, k12=k12, k21=k21, k_input=k_input),
            "lymph_rate_params": dict(A=Al, B=Bl, C=Cl, D=Dl),
            "vessel_rate_params": dict(A=Av, B=Bv, C=Cv, D=Dv),
            "sc_curve": sc[:, 0],
            "iv_curve": iv[:, 0],
            "central_curve": sc[:, 0],
            "time_hr": tt,
            "Cmax": cmax,
            "Tmax_hr": tmax,
            "t_half_hr": t_half,
        }
    return float(BA)


def compute_bioavailability_ci(c_lymph, c_vessel, t_min, oof_resid,
                               dose=DEFAULT_DOSE, drug="IgG", model=3,
                               m5_params=None, n_mc=300, levels=(95.0, 99.0),
                               seed=42, lymph_idx=0, vessel_idx=1):
    """
    대리모델 불확실성을 전파한 BA 예측구간(신뢰구간).

    BA는 대리모델이 예측한 농도곡선(c_lymph, c_vessel)으로 계산되는데, 그 곡선에는
    honest OOF 잔차로 정량화된 불확실성(conformal 예측구간)이 있다. 이 함수는 그
    잔차를 곡선에 몬테카를로로 주입해 여러 개의 그럴듯한 곡선을 만들고, 각 곡선의
    BA를 계산해 BA 분포 → percentile 예측구간을 얻는다.

    인자:
      c_lymph, c_vessel : 누적 % mass 곡선 (대리모델 base 예측), 각 (n_t,)
      t_min             : 시간축 (분)
      oof_resid         : OOF 잔차 배열 (n_cases, n_t, n_out) = 실제 − 앙상블예측.
                          시점 간 상관을 보존하기 위해 잔차 '행'(케이스) 전체를
                          부트스트랩 샘플링한다.
      dose, drug, model, m5_params : compute_bioavailability와 동일
      n_mc              : 몬테카를로 표본 수
      levels            : 예측구간 수준(%) 튜플, 예 (95, 99)
      seed              : 난수 시드
      lymph_idx,vessel_idx : oof_resid의 구획 축에서 lymph·vessel 위치
    반환:
      dict: base, median, std, ci({level:(lo,hi)}), samples, n_valid,
            curve_bands({time_hr, sc_base, iv_curve, sc_lo95/hi95, sc_lo99/hi99})
            — curve_bands는 SC 혈장곡선의 시점별 95%/99% 예측 밴드(그래프용).
    """
    c_lymph = np.asarray(c_lymph, float)
    c_vessel = np.asarray(c_vessel, float)
    oof_resid = np.asarray(oof_resid, float)
    n_cases = oof_resid.shape[0]

    def _detail(cl, cv):
        return compute_bioavailability(cl, cv, t_min, dose=dose, drug=drug,
                                       model=model, m5_params=m5_params,
                                       return_detail=True)

    base_d = _detail(c_lymph, c_vessel)
    base = base_d["BA"]
    time_hr = np.asarray(base_d["time_hr"], float)     # SC/IV 곡선 시간축(hr)
    iv_curve = np.asarray(base_d["iv_curve"], float)
    sc_base = np.asarray(base_d["sc_curve"], float)

    rng = np.random.default_rng(seed)
    il = rng.integers(0, n_cases, n_mc)
    iv = rng.integers(0, n_cases, n_mc)
    samples = []
    sc_curves = []                                     # 표본별 SC 혈장곡선
    for i in range(n_mc):
        cl = np.clip(c_lymph + oof_resid[il[i], :, lymph_idx], 0.0, None)
        cv = np.clip(c_vessel + oof_resid[iv[i], :, vessel_idx], 0.0, None)
        try:
            d = _detail(cl, cv)
        except Exception:
            continue
        b = d["BA"]
        if np.isfinite(b):
            samples.append(b)
            sc_curves.append(np.asarray(d["sc_curve"], float))
    samples = np.asarray(samples)
    sc_curves = np.asarray(sc_curves)                  # (n_valid, n_tt)

    ci = {}
    for lv in levels:
        a = (100.0 - lv) / 2.0
        ci[lv] = (float(np.percentile(samples, a)),
                  float(np.percentile(samples, 100.0 - a)))

    # 시점별 SC 혈장곡선 밴드 (예측 그래프와 동일한 95%/99% 표현)
    curve_bands = {"time_hr": time_hr, "sc_base": sc_base, "iv_curve": iv_curve}
    if len(sc_curves):
        for lv in levels:
            a = (100.0 - lv) / 2.0
            curve_bands[f"sc_lo{int(lv)}"] = np.percentile(sc_curves, a, axis=0)
            curve_bands[f"sc_hi{int(lv)}"] = np.percentile(sc_curves, 100.0 - a, axis=0)

    return {
        "base": float(base),
        "median": float(np.median(samples)),
        "std": float(samples.std()),
        "ci": ci,
        "samples": samples,
        "n_valid": int(len(samples)),
        "curve_bands": curve_bands,
    }


def compute_pk_metrics_ci(c_lymph, c_vessel, t_min, conf_half,
                          dose=DEFAULT_DOSE, drug="IgG", model=3,
                          m5_params=None):
    """
    conformal 예측밴드를 그대로 전파한 약동학 지표 구간 (몬테카를로 없음, 빠름).

    농도 예측 그래프에 쓰는 것과 '동일한' 시점별 conformal 반폭(conf_half)을
    림프·혈관 곡선에 적용해 상한/하한 흡수곡선을 만들고, 각각으로 SC 혈장곡선을
    적분한다. 그 결과가 BA 그래프에 보이는 SC 밴드(sc_lo/hi)다. 지표는:
      - BA, Tmax : 상·하한 곡선의 detail에서 직접
      - Ka, t½   : 상·하한 SC '혈장곡선'을 흡수-소실 이중지수로 피팅해 추출
                   (SC 혈장곡선을 피팅하므로 흡수곡선 피팅과 달리 안정적)
    모든 구간은 base 값을 포함하도록 [min(base,lo,hi), max(base,lo,hi)]로 만든다
    (밴드가 시간에 따라 비대칭이라 상·하한 지표가 한쪽으로 쏠릴 수 있으므로).

    인자:
      c_lymph, c_vessel : 대리모델 base 예측 곡선 (n_t,)
      t_min             : 시간축 (분)
      conf_half         : {level(%): (q_lymph, q_vessel)} 시점별 반폭 배열 dict.
      dose, drug, model, m5_params : compute_bioavailability와 동일
    반환:
      dict:
        base    : base 곡선의 detail
        metrics : {"BA","Ka","Tmax_hr","t_half_hr"} 각각 base 점추정 값
        ci      : {metric: {level: (lo, hi)}}  — 4개 지표 모두 구간 부여
        curve_bands : {time_hr, sc_base, iv_curve, sc_lo/hi95, sc_lo/hi99}
    """
    c_lymph = np.asarray(c_lymph, float)
    c_vessel = np.asarray(c_vessel, float)

    def _detail(cl, cv):
        return compute_bioavailability(cl, cv, t_min, dose=dose, drug=drug,
                                       model=model, m5_params=m5_params,
                                       return_detail=True)

    def _sc_fit(t_hr_curve, sc):
        """SC 혈장곡선을 흡수-소실 이중지수 A*(e^{-ke t}-e^{-ka t})로 피팅.
        반환 (Ka=ka, t_half=ln2/ke). 실패 시 (nan, nan)."""
        def _biexp(t, A, ka, ke):
            return A * (np.exp(-ke * t) - np.exp(-ka * t))
        ip = int(np.argmax(sc)); cmax = float(sc[ip])
        if cmax <= 0:
            return float("nan"), float("nan")
        try:
            popt, _ = curve_fit(_biexp, t_hr_curve, sc,
                                p0=[cmax * 2, 0.1, 0.01], maxfev=10000,
                                bounds=([0, 1e-4, 1e-5], [np.inf, 50, 10]))
            _, ka, ke = popt
            ka, ke = max(ka, ke), min(ka, ke)          # ka(흡수) > ke(소실)
            return float(ka), float(np.log(2) / ke)
        except Exception:
            return float("nan"), float("nan")

    base_d = _detail(c_lymph, c_vessel)
    t_hr_curve = np.asarray(base_d["time_hr"], float)
    Ka0, th0 = _sc_fit(t_hr_curve, np.asarray(base_d["sc_curve"], float))
    base_metrics = {"BA": base_d["BA"], "Ka": Ka0,
                    "Tmax_hr": base_d["Tmax_hr"], "t_half_hr": th0}

    metric_names = ["BA", "Ka", "Tmax_hr", "t_half_hr"]
    ci = {m: {} for m in metric_names}
    curve_bands = {"time_hr": t_hr_curve, "sc_base": base_d["sc_curve"],
                   "iv_curve": base_d["iv_curve"]}

    for lv, (qL, qV) in conf_half.items():
        qL = np.asarray(qL, float); qV = np.asarray(qV, float)
        cl_lo = np.clip(c_lymph - qL, 0, None); cl_hi = c_lymph + qL
        cv_lo = np.clip(c_vessel - qV, 0, None); cv_hi = c_vessel + qV
        d_lo = _detail(cl_lo, cv_lo)      # 흡수곡선 하한 → SC 혈장곡선 하한
        d_hi = _detail(cl_hi, cv_hi)      # 흡수곡선 상한 → SC 혈장곡선 상한
        sc_lo = np.asarray(d_lo["sc_curve"], float)
        sc_hi = np.asarray(d_hi["sc_curve"], float)
        Ka_lo, th_lo = _sc_fit(t_hr_curve, sc_lo)   # SC 밴드 곡선 직접 피팅
        Ka_hi, th_hi = _sc_fit(t_hr_curve, sc_hi)
        vals = {
            "BA":        (d_lo["BA"], d_hi["BA"]),
            "Ka":        (Ka_lo, Ka_hi),
            "Tmax_hr":   (d_lo["Tmax_hr"], d_hi["Tmax_hr"]),
            "t_half_hr": (th_lo, th_hi),
        }
        for m in metric_names:
            a, b = vals[m]
            base_v = base_metrics[m]
            cand = [x for x in (a, b, base_v) if np.isfinite(x)]  # base 포함
            ci[m][lv] = (float(min(cand)), float(max(cand)))
        curve_bands[f"sc_lo{int(lv)}"] = sc_lo
        curve_bands[f"sc_hi{int(lv)}"] = sc_hi

    return {"base": base_d, "metrics": base_metrics,
            "ci": ci, "curve_bands": curve_bands}


# ── 자체 검증 (원본 SimBiology BA=96.8% 재현, model 3 / Fig 2B) ──
if __name__ == "__main__":
    dose = 40.0
    def kle(t): return (10.02*np.exp(-8.868*t) + 2.407*np.exp(-0.035*t)) * dose * 0.01
    def kve(t): return (0.919*np.exp(-0.035*t) + 1.061*np.exp(-1.576*t)) * dose * 0.01
    tt = np.linspace(0, 1000, 20000)
    sc = odeint(_sc_ode_m3, [0, 0, 0], tt,
                args=(kle, kve, K_10, K_12, K_21, K_INPUT))
    iv = odeint(_iv_ode_2pool, [dose, 0], tt, args=(K_10, K_12, K_21))
    ba = _trapz(sc[:, 0], tt) / _trapz(iv[:, 0], tt) * 100
    print(f"검증(model 3, Fig 2B): BA={ba:.1f}% (원본 SimBiology 96.8%, 오차 {abs(ba-96.8):.2f}%p)")