""" clamp_core.py — GlucoClampTracer 순수 계산 모듈 ================================================= 하이퍼인슐린혈성-정상혈당 클램프(hyperinsulinemic-euglycemic clamp) 실험의 in vivo 인슐린 감수성 정량을 위한 순수 Python 계산 모듈. Streamlit에 의존하지 않으므로 단독 import / 단위 테스트가 가능하다. 의존성: numpy, scipy, pandas (statsmodels는 선택 — 혼합효과모형에만 사용). 참고용 · 연구용(for reference / research use only). 임상 진단 목적 사용 금지. """ from __future__ import annotations import numpy as np import pandas as pd from scipy import stats # statsmodels는 선택적 의존성 — 혼합효과모형(mixed-effects)에만 필요 try: # pragma: no cover - 환경 의존 import statsmodels.formula.api as smf # noqa: F401 _HAS_STATSMODELS = True except Exception: # pragma: no cover _HAS_STATSMODELS = False # ===================================================================== # 1. 정상상태(steady-state) 윈도우 자동 선택 # ===================================================================== def detect_steady_state( time_min, glucose, target_glucose: float = 150.0, target_tol: float = 20.0, max_cv_pct: float = 10.0, min_window_min: float = 30.0, ): """혈당 시계열에서 정상상태 윈도우를 자동 탐지한다. 기준(사용자 정의): - 윈도우 내 혈당 변동계수(CV%) <= max_cv_pct - 윈도우 평균 혈당이 target_glucose +/- target_tol 범위 내 - 윈도우 길이 >= min_window_min (분) 조건을 만족하는 가장 긴 윈도우를 반환한다. 동일 길이면 CV가 가장 낮은 것. Returns ------- dict : start_idx, end_idx, start_min, end_min, mean_glucose, cv_pct, n_points, found(bool) """ time_min = np.asarray(time_min, dtype=float) glucose = np.asarray(glucose, dtype=float) n = len(time_min) best = { "start_idx": None, "end_idx": None, "start_min": None, "end_min": None, "mean_glucose": None, "cv_pct": None, "n_points": 0, "found": False, } best_score = (-1.0, np.inf) # (윈도우 길이, -CV) — 길수록 / CV 낮을수록 좋음 for i in range(n): for j in range(i + 1, n + 1): seg_t = time_min[i:j] seg_g = glucose[i:j] if len(seg_g) < 2: continue win_len = seg_t[-1] - seg_t[0] if win_len < min_window_min: continue mean_g = float(np.mean(seg_g)) sd_g = float(np.std(seg_g, ddof=1)) cv = (sd_g / mean_g * 100.0) if mean_g else np.inf if cv > max_cv_pct: continue if abs(mean_g - target_glucose) > target_tol: continue score = (win_len, -cv) if score > best_score: best_score = score best = { "start_idx": i, "end_idx": j - 1, "start_min": float(seg_t[0]), "end_min": float(seg_t[-1]), "mean_glucose": mean_g, "cv_pct": cv, "n_points": len(seg_g), "found": True, } return best def mean_gir_in_window(time_min, gir, start_min, end_min): """정상상태 윈도우 내 평균 GIR을 계산한다(시간 가중 평균).""" time_min = np.asarray(time_min, dtype=float) gir = np.asarray(gir, dtype=float) mask = (time_min >= start_min) & (time_min <= end_min) if mask.sum() == 0: return float("nan") return float(np.mean(gir[mask])) # ===================================================================== # 2. 단위 변환 / 정규화 # ===================================================================== def gir_to_per_kg(gir_value, unit: str, body_weight_g: float): """GIR 절대값을 mg/kg/min 으로 변환한다. Parameters ---------- gir_value : float — 펌프 로그 GIR unit : str — 'mg/min', 'uL/min', 'mL/h' 중 하나 body_weight_g : float — 체중(g) 주입액은 표준 20% 덱스트로스(0.2 g glucose / mL)로 가정. """ bw_kg = body_weight_g / 1000.0 if bw_kg <= 0: return float("nan") glucose_conc_mg_per_uL = 0.2 # 20% dextrose = 200 mg/mL = 0.2 mg/uL if unit == "mg/min": mg_per_min = gir_value elif unit == "uL/min": mg_per_min = gir_value * glucose_conc_mg_per_uL elif unit == "mL/h": uL_per_min = gir_value * 1000.0 / 60.0 mg_per_min = uL_per_min * glucose_conc_mg_per_uL else: raise ValueError(f"지원하지 않는 GIR 단위: {unit}") return mg_per_min / bw_kg def normalize_by_mass(rate_per_kg_bw, body_weight_g, lean_mass_g=None): """체중 기준 rate를 lean mass 기준으로 환산한다(lean mass 제공 시).""" if lean_mass_g is None or lean_mass_g <= 0: return rate_per_kg_bw return rate_per_kg_bw * (body_weight_g / lean_mass_g) # ===================================================================== # 3. 트레이서 동역학 — Steele 방정식 # ===================================================================== def steele_steady_state(tracer_infusion_rate, plasma_enrichment): """정상상태 Steele 방정식으로 포도당 출현율 Ra를 계산한다. 정상상태에서는 Ra = Rd 이며: Ra = F / SA (방사성 트레이서, SA = specific activity) Ra = F / (E/100) (안정동위원소, E = enrichment %) Parameters ---------- tracer_infusion_rate : float — 트레이서 주입율(dpm/min 또는 mole/min) plasma_enrichment : float — 혈장 비방사능(dpm/mg) 또는 enrichment 분율 비율 Returns ------- float : Ra (= Rd), 단위는 입력에 의존(보통 mg/kg/min) """ if plasma_enrichment <= 0: return float("nan") return tracer_infusion_rate / plasma_enrichment def steele_non_steady_state( tracer_infusion_rate, enrichment_t1, enrichment_t2, glucose_t1, glucose_t2, dt_min, pool_fraction: float = 0.65, distribution_volume_mL_per_kg: float = 200.0, ): """비정상상태 Steele 1조구획(one-compartment) 방정식으로 Ra를 계산한다. Steele(1959) 비정상상태 식: Ra = [ F - pV * ( (G1+G2)/2 ) * (E2-E1)/(dt) ] / ( (E1+E2)/2 ) 여기서 E는 비방사능 또는 enrichment, G는 혈당, pV는 유효 분포용적이다. Parameters ---------- tracer_infusion_rate : float — 트레이서 주입율 enrichment_t1, enrichment_t2 : float — 두 시점의 비방사능/enrichment glucose_t1, glucose_t2 : float — 두 시점의 혈당(mg/dL) dt_min : float — 두 시점 간 간격(분) pool_fraction : float — pool correction factor p (보통 0.65) distribution_volume_mL_per_kg : float — 포도당 분포용적(mL/kg) Returns ------- float : Ra (mg/kg/min) """ if dt_min <= 0: return float("nan") mean_E = (enrichment_t1 + enrichment_t2) / 2.0 if mean_E <= 0: return float("nan") mean_G = (glucose_t1 + glucose_t2) / 2.0 # pV: 유효 분포용적 (mL/kg) -> dL/kg 환산 (혈당이 mg/dL 이므로) pV_dL_per_kg = pool_fraction * distribution_volume_mL_per_kg / 100.0 dE_dt = (enrichment_t2 - enrichment_t1) / dt_min ra = (tracer_infusion_rate - pV_dL_per_kg * mean_G * dE_dt) / mean_E return float(ra) def hot_ginf_correction(tracer_infusion_rate, hot_ginf_rate): """Hot-GINF 보정: 주입 포도당에 트레이서가 섞인 경우 유효 주입율 보정. 클램프 중 주입하는 비표지(cold) 포도당에 트레이서가 함께 들어가면 내인성 출현율을 과대평가하게 된다. 유효 트레이서 주입율에서 hot-GINF 기여분을 차감한다. """ return max(tracer_infusion_rate - hot_ginf_rate, 0.0) def compute_kinetics( tracer_infusion_rate, plasma_enrichment_basal, plasma_enrichment_clamp, gir_per_kg, mode: str = "steady", hot_ginf_rate: float = 0.0, nss_params: dict | None = None, ): """기저(basal) 및 클램프(clamp) 상태의 포도당 동역학을 한번에 계산한다. 핵심 정의 --------- Ra (rate of appearance) : 포도당 출현율 Rd (rate of disappearance): 포도당 소실율 = 전신 포도당 처리율 HGP (hepatic glucose production) : 간 포도당 생성 - 기저: HGP_basal = Ra_basal (외인성 포도당 없음) - 클램프: HGP_clamp = Ra_clamp - GIR HGP 억제율(%) = (HGP_basal - HGP_clamp) / HGP_basal * 100 클램프 Rd = Ra_clamp = HGP_clamp + GIR Returns ------- dict : ra_basal, ra_clamp, hgp_basal, hgp_clamp, rd_clamp, hgp_suppression_pct, gir """ eff_tir = hot_ginf_correction(tracer_infusion_rate, hot_ginf_rate) if mode == "steady": ra_basal = steele_steady_state(eff_tir, plasma_enrichment_basal) ra_clamp = steele_steady_state(eff_tir, plasma_enrichment_clamp) elif mode == "non-steady": p = nss_params or {} ra_basal = steele_steady_state(eff_tir, plasma_enrichment_basal) ra_clamp = steele_non_steady_state( eff_tir, p.get("enrichment_t1", plasma_enrichment_clamp), p.get("enrichment_t2", plasma_enrichment_clamp), p.get("glucose_t1", 150.0), p.get("glucose_t2", 150.0), p.get("dt_min", 10.0), p.get("pool_fraction", 0.65), p.get("distribution_volume_mL_per_kg", 200.0), ) else: raise ValueError("mode 는 'steady' 또는 'non-steady' 여야 합니다.") hgp_basal = ra_basal # 기저상태: 외인성 포도당 없음 hgp_clamp = ra_clamp - gir_per_kg rd_clamp = ra_clamp # 정상상태 클램프: Rd = Ra if hgp_basal and not np.isnan(hgp_basal) and hgp_basal != 0: hgp_supp = (hgp_basal - hgp_clamp) / hgp_basal * 100.0 else: hgp_supp = float("nan") return { "ra_basal": float(ra_basal), "ra_clamp": float(ra_clamp), "hgp_basal": float(hgp_basal), "hgp_clamp": float(hgp_clamp), "rd_clamp": float(rd_clamp), "hgp_suppression_pct": float(hgp_supp), "gir": float(gir_per_kg), } # ===================================================================== # 4. 조직 특이적 2-deoxyglucose 흡수 (Rg') # ===================================================================== def plasma_2dg_auc(time_min, plasma_2dg_dpm): """혈장 2-DG 소실곡선의 AUC를 사다리꼴 적분으로 계산한다. Rg' 계산 시 분모로 사용되는 적분 비방사능(integrated specific activity). Returns ------- float : AUC (dpm·min / mL) """ time_min = np.asarray(time_min, dtype=float) plasma_2dg_dpm = np.asarray(plasma_2dg_dpm, dtype=float) if len(time_min) < 2: return float("nan") return float(np.trapz(plasma_2dg_dpm, time_min)) def tissue_rg( tissue_2dg6p_dpm_per_g, plasma_2dg_auc_dpm_min_per_mL, mean_plasma_glucose_mg_per_mL, ): """조직 특이적 포도당 대사지수 Rg' 를 계산한다. Rg' = [조직 내 2-DG-6-P 축적량] / [혈장 2-DG AUC] * [평균 혈장 포도당] Parameters ---------- tissue_2dg6p_dpm_per_g : float — 조직 인산화 2-DG-6-P (dpm/g tissue) plasma_2dg_auc_dpm_min_per_mL : float — 혈장 2-DG AUC (dpm·min/mL) mean_plasma_glucose_mg_per_mL : float — 클램프 평균 혈장 포도당 (mg/mL) Returns ------- float : Rg' (umol/100g/min 스케일 — 입력 단위에 따름; 보통 nmol/g/min) """ if plasma_2dg_auc_dpm_min_per_mL <= 0: return float("nan") return ( tissue_2dg6p_dpm_per_g / plasma_2dg_auc_dpm_min_per_mL * mean_plasma_glucose_mg_per_mL ) def compute_tissue_panel(tissue_df, plasma_auc, mean_plasma_glucose_mg_per_mL): """여러 조직의 Rg' 를 한번에 계산한다. Parameters ---------- tissue_df : DataFrame — 'tissue', 'tissue_2dg6p_dpm_per_g' 컬럼 필요 plasma_auc : float — 혈장 2-DG AUC mean_plasma_glucose_mg_per_mL : float Returns ------- DataFrame : tissue, rg_prime 컬럼 """ out = tissue_df.copy() out["rg_prime"] = out["tissue_2dg6p_dpm_per_g"].apply( lambda x: tissue_rg(x, plasma_auc, mean_plasma_glucose_mg_per_mL) ) return out # ===================================================================== # 5. 모델 참조 범위 + QC 플래그 # ===================================================================== # 마우스/랫 모델별 전형적 클램프 지표 참조 범위(문헌 기반 근사치). # 단위: GIR mg/kg/min, HGP 억제율 %, 클램프 Rd mg/kg/min MODEL_REFERENCE = { "C57BL/6": { # 정상 인슐린 감수성 "gir": (25.0, 60.0), "hgp_suppression_pct": (70.0, 100.0), "rd_clamp": (35.0, 75.0), "note": "정상 야생형 — 높은 GIR, 강한 HGP 억제", }, "DIO": { # diet-induced obese — 인슐린 저항성 "gir": (8.0, 25.0), "hgp_suppression_pct": (30.0, 65.0), "rd_clamp": (15.0, 40.0), "note": "고지방식이 비만 — 중등도 인슐린 저항성", }, "db/db": { # leptin receptor 결손 — 심한 저항성 "gir": (0.0, 12.0), "hgp_suppression_pct": (0.0, 40.0), "rd_clamp": (5.0, 25.0), "note": "렙틴수용체 결손 — 심한 인슐린 저항성/고혈당", }, "ob/ob": { # leptin 결손 "gir": (0.0, 15.0), "hgp_suppression_pct": (10.0, 45.0), "rd_clamp": (8.0, 28.0), "note": "렙틴 결손 — 심한 비만/인슐린 저항성", }, "ZDF": { # Zucker diabetic fatty rat "gir": (2.0, 14.0), "hgp_suppression_pct": (5.0, 45.0), "rd_clamp": (6.0, 26.0), "note": "Zucker 당뇨비만 랫 — 베타세포 부전 동반", }, "GK": { # Goto-Kakizaki rat — 비비만 2형 당뇨 "gir": (10.0, 28.0), "hgp_suppression_pct": (25.0, 60.0), "rd_clamp": (14.0, 36.0), "note": "Goto-Kakizaki — 비비만 2형 당뇨 모델", }, } def reference_range(model: str, metric: str): """모델/지표에 대한 참조 범위 튜플 (low, high) 를 반환한다. 없으면 None.""" m = MODEL_REFERENCE.get(model) if not m: return None return m.get(metric) def qc_flags( glucose_cv_pct, steady_state_found: bool, tracer_enrichment_cv_pct, hot_ginf_applied: bool, model: str | None = None, gir: float | None = None, hgp_suppression_pct: float | None = None, rd_clamp: float | None = None, max_glucose_cv: float = 10.0, max_tracer_cv: float = 15.0, ): """클램프 실험 품질관리(QC) 경고 목록을 생성한다. Returns ------- list[dict] : 각 항목 {level: 'error'/'warning'/'info', message: str} """ flags = [] if not steady_state_found: flags.append({ "level": "error", "message": "정상상태 윈도우를 찾지 못했습니다 — 혈당이 수렴하지 않았습니다.", }) if glucose_cv_pct is not None and glucose_cv_pct > max_glucose_cv: flags.append({ "level": "warning", "message": f"혈당 변동계수 {glucose_cv_pct:.1f}% > 기준 {max_glucose_cv:.0f}% — 불안정한 정상혈당.", }) if tracer_enrichment_cv_pct is not None and tracer_enrichment_cv_pct > max_tracer_cv: flags.append({ "level": "warning", "message": f"트레이서 enrichment CV {tracer_enrichment_cv_pct:.1f}% > 기준 {max_tracer_cv:.0f}% — 트레이서 정상상태 불안정.", }) if not hot_ginf_applied: flags.append({ "level": "info", "message": "Hot-GINF 보정이 적용되지 않았습니다 — 표지 포도당 주입 시 HGP 과대평가 가능.", }) # 모델 참조 범위 대비 이탈 검사 if model: for metric, value in ( ("gir", gir), ("hgp_suppression_pct", hgp_suppression_pct), ("rd_clamp", rd_clamp), ): if value is None or (isinstance(value, float) and np.isnan(value)): continue rng = reference_range(model, metric) if rng and not (rng[0] <= value <= rng[1]): flags.append({ "level": "warning", "message": ( f"{metric} = {value:.1f} 이 {model} 모델 참조범위 " f"{rng[0]:.0f}-{rng[1]:.0f} 를 벗어났습니다." ), }) if not flags: flags.append({"level": "info", "message": "QC 경고 없음 — 모든 점검 통과."}) return flags # ===================================================================== # 6. 코호트 통계 # ===================================================================== def one_way_anova(df, group_col: str, value_col: str): """그룹 간 일원배치 분산분석(one-way ANOVA)을 수행한다. Returns ------- dict : f_stat, p_value, groups, group_means, group_n, k(그룹수) """ df = df[[group_col, value_col]].dropna() groups = [] samples = [] means = {} counts = {} for name, sub in df.groupby(group_col): vals = sub[value_col].to_numpy(dtype=float) if len(vals) < 1: continue groups.append(str(name)) samples.append(vals) means[str(name)] = float(np.mean(vals)) counts[str(name)] = int(len(vals)) result = { "f_stat": float("nan"), "p_value": float("nan"), "groups": groups, "group_means": means, "group_n": counts, "k": len(groups), } if len(samples) >= 2 and all(len(s) >= 2 for s in samples): f_stat, p_value = stats.f_oneway(*samples) result["f_stat"] = float(f_stat) result["p_value"] = float(p_value) return result def tukey_pairwise(df, group_col: str, value_col: str): """그룹 쌍별 비교 — 독립표본 t검정 + Bonferroni 보정. (statsmodels 미설치 환경에서도 동작하도록 scipy 기반으로 구현) Returns ------- DataFrame : group1, group2, mean_diff, t_stat, p_raw, p_bonferroni """ df = df[[group_col, value_col]].dropna() grp = {str(k): v[value_col].to_numpy(dtype=float) for k, v in df.groupby(group_col)} names = list(grp.keys()) rows = [] pairs = [(names[i], names[j]) for i in range(len(names)) for j in range(i + 1, len(names))] m = max(len(pairs), 1) for g1, g2 in pairs: a, b = grp[g1], grp[g2] if len(a) >= 2 and len(b) >= 2: t_stat, p_raw = stats.ttest_ind(a, b, equal_var=False) else: t_stat, p_raw = float("nan"), float("nan") rows.append({ "group1": g1, "group2": g2, "mean_diff": float(np.mean(a) - np.mean(b)), "t_stat": float(t_stat), "p_raw": float(p_raw), "p_bonferroni": float(min(p_raw * m, 1.0)) if not np.isnan(p_raw) else float("nan"), }) return pd.DataFrame(rows) def mixed_effects_available() -> bool: """혼합효과모형 사용 가능 여부(statsmodels 설치 여부).""" return _HAS_STATSMODELS def mixed_effects_model(df, value_col: str, group_col: str, subject_col: str): """선형 혼합효과모형 — 개체를 random effect로 둔 그룹 효과 추정. statsmodels가 설치된 경우에만 동작. 미설치 시 None 반환. Returns ------- dict | None : params, pvalues, summary_text """ if not _HAS_STATSMODELS: return None import statsmodels.formula.api as smf d = df[[value_col, group_col, subject_col]].dropna().copy() d = d.rename(columns={value_col: "y", group_col: "grp", subject_col: "subj"}) model = smf.mixedlm("y ~ C(grp)", d, groups=d["subj"]) fit = model.fit(reml=True) return { "params": fit.params.to_dict(), "pvalues": fit.pvalues.to_dict(), "summary_text": str(fit.summary()), } def group_summary(df, group_col: str, value_cols): """그룹별 mean / sd / sem / n 요약표를 생성한다.""" if isinstance(value_cols, str): value_cols = [value_cols] rows = [] for name, sub in df.groupby(group_col): row = {"group": str(name), "n": int(len(sub))} for col in value_cols: vals = sub[col].to_numpy(dtype=float) vals = vals[~np.isnan(vals)] if len(vals) == 0: row[f"{col}_mean"] = float("nan") row[f"{col}_sd"] = float("nan") row[f"{col}_sem"] = float("nan") continue row[f"{col}_mean"] = float(np.mean(vals)) row[f"{col}_sd"] = float(np.std(vals, ddof=1)) if len(vals) > 1 else 0.0 row[f"{col}_sem"] = ( float(np.std(vals, ddof=1) / np.sqrt(len(vals))) if len(vals) > 1 else 0.0 ) rows.append(row) return pd.DataFrame(rows) # ===================================================================== # 자체 점검(모듈 단독 실행 시) # ===================================================================== if __name__ == "__main__": # 간단한 자기검증 — 의존성/수식 sanity check t = np.arange(0, 121, 10) g = np.array([300, 250, 200, 165, 152, 149, 151, 150, 148, 150, 151, 149, 150], dtype=float) ss = detect_steady_state(t, g, target_glucose=150, target_tol=20, max_cv_pct=10, min_window_min=30) print("steady-state:", ss["found"], ss["start_min"], "-", ss["end_min"], "min, CV=%.2f%%" % ss["cv_pct"]) kin = compute_kinetics( tracer_infusion_rate=1.0e6, # dpm/min 예시 plasma_enrichment_basal=4.0e4, plasma_enrichment_clamp=2.5e4, gir_per_kg=30.0, mode="steady", ) print("kinetics:", {k: round(v, 2) for k, v in kin.items()}) rg = tissue_rg(1.2e5, 3.5e6, 1.5) print("tissue Rg':", round(rg, 4)) print("statsmodels available:", mixed_effects_available()) print("clamp_core self-check OK")