alopecia/js/twin-engine.js

/* ============================================================
   twin-engine.js — 모낭 디지털 트윈 ODE 엔진 (브라우저 포팅)
   digital_twin/follicle_model.py 의 정확한 포팅. RK4 적분.
   라이브 'what-if' 모드용. 정밀 시나리오는 twin_scenarios.json 사용.
   ============================================================ */
(function (global) {
  'use strict';

  const STATE_NAMES = ["Wnt", "BMP", "SHH", "DP", "HFSC", "APO", "Hair"];

  // Python Params 와 1:1 동일
  const P = {
    kWp: 0.45, kWd: 0.32, kBp: 0.30, kBd: 0.38, kSp: 0.50, kSd: 0.40,
    kDp: 0.50, kDd: 0.30, kHp: 0.58, kHd: 0.24, kAd: 0.50,
    kHairP: 0.80, kHairD: 0.16, Wbasal: 1.0, Bbasal: 0.4, Dbasal: 1.0,
    kDKK: 0.65, Kd: 1.0, Kb: 1.2, kBMPand: 0.30, Kw: 0.7, Ksw: 0.5, Kb2: 0.9,
    kIGFand: 0.55, kMiniAND: 0.40, hpow: 0.6, Hcap: 3.5,
    kDinf: 0.45, kDapo: 0.5, Kh: 0.35, Ks: 0.35, Kb3: 1.5,
    kHinf: 0.55, kHapo: 0.6, kApAND: 0.30, kApINF: 0.85, kSurv: 1.6,
    Ksh: 0.4, kHairApo: 0.55,
  };

  const DISEASE_PRESETS = {
    "Healthy": { AND: 0.10, INF: 0.05, chemo_amp: 0.0 },
    "Androgenetic Alopecia": { AND: 0.62, INF: 0.08, chemo_amp: 0.0 },
    "Alopecia Areata": { AND: 0.12, INF: 0.82, chemo_amp: 0.0 },
    "Chemotherapy-induced Alopecia": { AND: 0.10, INF: 0.06, chemo_amp: 1.7 },
  };
  const CHEMO_TIMES = [10, 31, 52, 73];

  // 개입 → drive 수정자
  const INTERVENTIONS = {
    none: d => d,
    finasteride: d => (d.AND *= 0.40, d),
    dutasteride: d => (d.AND *= 0.25, d),
    AR_antagonist: d => (d.AND *= 0.45, d),
    minoxidil: d => (d.uDP += 0.80, d.uWnt += 0.15, d),
    anti_DKK1: d => (d.uWnt += 0.45, d),
    wnt_agonist: d => (d.uWnt += 0.50, d),
    JAK_inhibitor: d => (d.INF *= 0.18, d),
    corticosteroid: d => (d.INF *= 0.55, d),
    CDK46_inhibitor: d => (d.chemo_protect = 0.70, d),
    scalp_cooling: d => (d.chemo_protect = Math.max(d.chemo_protect, 0.45), d),
    PTH_CBD: d => (d.uDP += 0.35, d),
    exosome_MSC: d => (d.uDP += 0.35, d.uWnt += 0.25, d),
  };

  function buildDrive(disease, interventions, overrides) {
    const pre = DISEASE_PRESETS[disease] || DISEASE_PRESETS["Healthy"];
    let d = { AND: pre.AND, INF: pre.INF, chemo_amp: pre.chemo_amp || 0,
              uWnt: 0, uDP: 0, uNog: 0, chemo_protect: 0 };
    (interventions || []).forEach(iv => { if (INTERVENTIONS[iv]) d = INTERVENTIONS[iv](d); });
    if (overrides) Object.assign(d, overrides);   // 라이브 슬라이더
    d.AND = Math.max(0, d.AND); d.INF = Math.max(0, d.INF);
    return d;
  }

  function chemo(d, t) {
    if (!d.chemo_amp) return 0;
    let s = 0;
    for (const tc of CHEMO_TIMES) s += Math.exp(-((t - tc) ** 2) / (2 * 1.5 * 1.5));
    return d.chemo_amp * s * (1 - (d.chemo_protect || 0));
  }

  function rhs(t, y, d) {
    let [Wnt, BMP, SHH, DP, HFSC, APO, Hair] = y.map(v => Math.max(0, v));
    const AND = d.AND, INF = d.INF, ch = chemo(d, t);
    const DKK = P.kDKK * AND;
    const dWnt = P.kWp * (P.Wbasal + d.uWnt) / (1 + DKK / P.Kd) / (1 + BMP / P.Kb) - P.kWd * Wnt;
    const dBMP = P.kBp * (P.Bbasal + P.kBMPand * AND) / (1 + d.uNog) / (1 + Wnt / P.Kw) - P.kBd * BMP;
    const dSHH = P.kSp * (Wnt / (P.Ksw + Wnt)) / (1 + BMP / P.Kb2) - P.kSd * SHH;
    const IGF = 1 / (1 + P.kIGFand * AND);
    const dDP = P.kDp * (P.Dbasal + d.uDP) * IGF - P.kDd * DP - DP * (P.kDinf * INF + P.kDapo * APO);
    const dHFSC = P.kHp * (Wnt / (P.Kh + Wnt)) * (SHH / (P.Ks + SHH)) / (1 + BMP / P.Kb3)
                  - P.kHd * HFSC - HFSC * (P.kHinf * INF + P.kHapo * APO);
    const dAPO = (P.kApAND * AND + P.kApINF * INF + ch) / (1 + P.kSurv * DP) - P.kAd * APO;
    // 분수 지수(hpow=0.6)의 음수 밑 → Math.pow(neg,0.6)=NaN 방지: 밑을 명시적으로 0 하한.
    // (위 y.map(Math.max(0,·)) 클램프와 중복이나, 적분기 단계 값이 음수를 흘려도 NaN 전파 차단)
    const HFSCp = Math.pow(Math.max(0, HFSC), P.hpow);
    const DPp = Math.pow(Math.max(0, DP), P.hpow);
    const dHair = P.kHairP * HFSCp * DPp * (SHH / (P.Ksh + SHH))
                  / (1 + P.kMiniAND * AND) * Math.max(0, 1 - Hair / P.Hcap)
                  - P.kHairD * Hair - P.kHairApo * APO * Hair;
    return [dWnt, dBMP, dSHH, dDP, dHFSC, dAPO, dHair];
  }

  function rk4(y0, d, days, dt) {
    dt = dt || 0.5;
    const steps = Math.round(days / dt);
    const out = { t: [], y: y0.map(() => []) };
    let y = y0.slice();
    for (let i = 0; i <= steps; i++) {
      const t = i * dt;
      if (Math.abs(t - Math.round(t)) < 1e-9) {       // 일 단위 샘플만 저장
        out.t.push(Math.round(t));
        y.forEach((v, j) => out.y[j].push(v));
      }
      const k1 = rhs(t, y, d);
      const k2 = rhs(t + dt / 2, y.map((v, j) => v + dt / 2 * k1[j]), d);
      const k3 = rhs(t + dt / 2, y.map((v, j) => v + dt / 2 * k2[j]), d);
      const k4 = rhs(t + dt, y.map((v, j) => v + dt * k3[j]), d);
      y = y.map((v, j) => Math.max(0, v + dt / 6 * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j])));
    }
    return out;
  }

  let _healthy = null;
  function healthyState() {
    if (!_healthy) {
      const d = buildDrive("Healthy", []);
      const sol = rk4([1.0, 0.45, 0.7, 1.0, 0.8, 0.1, 1.0], d, 400);
      _healthy = sol.y.map(arr => arr[arr.length - 1]);
    }
    return _healthy;
  }

  function proteinReadouts(states, d) {
    const n = states.Wnt.length, ones = new Array(n);
    const AND = d.AND, INF = d.INF;
    return {
      "AR/DHT (AR)": ones.fill(AND).slice(),
      "DKK1": new Array(n).fill(P.kDKK * AND),
      "β-catenin (CTNNB1)": states.Wnt.slice(),
      "BMP4": states.BMP.slice(),
      "SHH": states.SHH.slice(),
      "IGF1": new Array(n).fill(1 / (1 + P.kIGFand * AND)),
      "VEGFA (DP)": states.DP.slice(),
      "JAK-STAT (STAT1)": new Array(n).fill(INF),
      "p53/apoptosis (TP53)": states.APO.slice(),
    };
  }

  // 질환의 '무처치 평형상태' = 환자가 이미 탈모를 가진 출발점(치료 전).
  // 건강 상태에서 질환 구동을 충분히 오래 적분해 정착시킨 종단 상태.
  const _diseq = {};
  function diseaseEquilibrium(disease) {
    if (!_diseq[disease]) {
      const d = buildDrive(disease, []);
      const sol = rk4(healthyState(), d, 600);
      _diseq[disease] = sol.y.map(arr => arr[arr.length - 1]);
    }
    return _diseq[disease].slice();
  }

  // 메인 API: 시뮬레이션 실행. opts.y0 = 커스텀 초기상태(미지정 시 건강 평형).
  function run(disease, interventions, opts) {
    opts = opts || {};
    const d = buildDrive(disease, interventions, opts.overrides);
    const y0 = opts.y0 || healthyState();
    const sol = rk4(y0, d, opts.days || 240);
    const states = {};
    STATE_NAMES.forEach((nm, i) => states[nm] = sol.y[i]);
    const hss = healthyState()[6];
    const rel = states.Hair.map(h => Math.min(100, 100 * h / hss));
    states.HairDensity = rel;
    const proteins = proteinReadouts(states, d);
    const metrics = {
      final_hair_density_pct: +rel[rel.length - 1].toFixed(1),
      min_hair_density_pct: +Math.min(...rel).toFixed(1),
      anagen_fraction: +(rel.filter(v => v > 70).length / rel.length).toFixed(3),
      AND_load: +d.AND.toFixed(3), INF_load: +d.INF.toFixed(3),
    };
    let tracked = [];
    (interventions || []).forEach(iv => {
      // genes 는 scenarios JSON 에서 보강; 엔진은 비움
    });
    return { disease, interventions, drive: d, t: sol.t, states, proteins, metrics };
  }

  global.TwinEngine = { run, diseaseEquilibrium, STATE_NAMES, DISEASE_PRESETS, INTERVENTIONS, buildDrive };
  if (typeof module !== 'undefined' && module.exports) module.exports = global.TwinEngine;
})(typeof window !== 'undefined' ? window : globalThis);