""" 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)")