용량-반응 곡선 피팅 — numpy로 EC50과 신뢰구간 자동 계산

이 토픽을 마치면

교과서에서 배운 numpy, matplotlib, 탐욕 최적화(greedy) 를 엮어서, 약물 스크리닝의 용량-반응 데이터에서 EC50과 신뢰구간을 자동 계산하는 도구를 직접 만들 수 있습니다. 96웰 플레이트 결과를 하나씩 손으로 계산하던 시간이 30초로 압축됩니다.

이 글은 교육용 일반 예제 입니다. 실전에서는 GraphPad Prism, R의 drc 패키지 등을 씁니다. 여기서는 그 도구들의 내부 원리를 정확히 다룹니다.

"이걸 다 손으로 하나?" — 반복 계산의 함정

여러분이 신약 스크리닝을 하고 있다고 합시다. 96웰 플레이트 하나에서 다음 데이터가 나옵니다.

각 웰에 서로 다른 농도의 후보 화합물을 넣음 (예: 0.001, 0.01, 0.1, 1, 10, 100 μM)
각 농도에서 세포 생존율(%)을 측정

여러분이 알고 싶은 것: EC50 (반응의 절반을 일으키는 농도) 와 최대 반응.

가장 순진한 접근: Excel에 데이터 붙여넣고, 산점도를 그리고, 눈으로 절반 지점을 찾음.

이 접근의 실전 문제:

문제 1: 정확성. 눈으로 그은 곡선은 사람마다 다릅니다. 재현성 X.

문제 2: 규모. 후보 화합물 100개 × 각각 8개 농도 = 800 데이터 포인트. Excel 클릭으로는 하루 종일.

문제 3: 신뢰구간. EC50이 예를 들어 2.3 μM 이라고 나왔을 때, 이 값이 "실제로 얼마나 정확한가"가 알 수 없음.

진짜 접근은 non-linear curve fitting입니다. Hill 방정식(시그모이드 함수)을 데이터에 최적 피팅하는 파라미터를 찾고, 그 파라미터의 통계적 신뢰구간까지 계산합니다.

블랙박스에서 부품으로

부품 1: Hill 방정식

용량-반응 곡선의 표준 수학적 모델이 Hill 방정식입니다.

text

response(dose) = bottom + (top - bottom) / (1 + (EC50 / dose)^hill)

파라미터 4개:

bottom: 최소 반응 (약물 없을 때)
top: 최대 반응 (충분히 많이 넣었을 때)
EC50: 반응이 (top+bottom)/2가 되는 농도
hill: 기울기 (곡선의 가파른 정도)

파이썬으로:

python

import numpy as np


def hill_equation(dose: np.ndarray, bottom: float, top: float, ec50: float, hill: float) -> np.ndarray:
    return bottom + (top - bottom) / (1 + (ec50 / dose) ** hill)

로그 스케일에서 편하게 시각화되므로, 일반적으로 dose를 log 변환 후 그립니다.

부품 2: 탐욕 최적화

파라미터 4개를 어떻게 찾는가? 데이터와 예측 사이 오차(residual)를 최소화하는 것. 최소 제곱법.

text

minimize sum_i (y_i - y_hat_i)^2

Non-linear 최소 제곱은 닫힌 해가 없어 반복 알고리즘으로 풉니다. Levenberg-Marquardt, trust-region 등. 이 알고리즘들은 greedy 접근으로 작동합니다 — 현재 위치에서 오차가 가장 줄어드는 방향으로 조금 이동, 다시 그 위치에서 반복.

scipy.optimize.curve_fit이 이 알고리즘의 구현입니다.

python

from scipy.optimize import curve_fit


def fit_hill(doses: np.ndarray, responses: np.ndarray) -> dict:
    # 초기값 추측
    p0 = [
        responses.min(),  # bottom
        responses.max(),  # top
        np.median(doses), # EC50
        1.0               # hill
    ]
    
    popt, pcov = curve_fit(hill_equation, doses, responses, p0=p0, maxfev=5000)
    perr = np.sqrt(np.diag(pcov))  # 표준오차
    
    return {
        "bottom": popt[0],
        "top": popt[1],
        "ec50": popt[2],
        "hill": popt[3],
        "bottom_se": perr[0],
        "top_se": perr[1],
        "ec50_se": perr[2],
        "hill_se": perr[3]
    }

부품 3: matplotlib 시각화

피팅 결과를 확인하는 표준 그래프.

python

import matplotlib.pyplot as plt


def plot_dose_response(doses: np.ndarray, responses: np.ndarray, fit: dict) -> None:
    fig, ax = plt.subplots(figsize=(7, 5))
    
    ax.scatter(doses, responses, color="steelblue", s=50, alpha=0.7, label="Data")
    
    dose_grid = np.logspace(np.log10(doses.min() * 0.5), np.log10(doses.max() * 2), 200)
    fit_curve = hill_equation(dose_grid, fit["bottom"], fit["top"], fit["ec50"], fit["hill"])
    ax.plot(dose_grid, fit_curve, color="darkred", linewidth=2, label=f"Fit (EC50={fit['ec50']:.3g})")
    
    ax.axvline(fit["ec50"], linestyle="--", color="gray", alpha=0.5)
    ax.axhline((fit["top"] + fit["bottom"]) / 2, linestyle="--", color="gray", alpha=0.5)
    
    ax.set_xscale("log")
    ax.set_xlabel("Dose (μM)")
    ax.set_ylabel("Response (%)")
    ax.legend()
    plt.tight_layout()
    plt.show()

실전 예시

가상 데이터로 확인:

python

doses = np.array([0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 100])
responses = np.array([2, 3, 5, 12, 28, 55, 78, 90, 96, 98, 99])

fit = fit_hill(doses, responses)
print(f"EC50: {fit['ec50']:.3g} ± {fit['ec50_se']:.3g}")
print(f"Hill: {fit['hill']:.2f}")
print(f"Range: {fit['bottom']:.1f} → {fit['top']:.1f}")

plot_dose_response(doses, responses, fit)

기대 출력:

text

EC50: 0.271 ± 0.0084
Hill: 0.98
Range: 1.2 → 99.5

페이딩 — 여러분이 채워야 할 세 빈칸

빈칸 1: 95% 신뢰구간

표준 오차만으로는 부족합니다. 95% 신뢰구간을 t-분포로 계산합니다.

python

from scipy import stats


def compute_ci(fit: dict, n_data: int, alpha: float = 0.05) -> dict:
    """
    95% 신뢰구간 계산.
    """
    dof = n_data - 4  # 4 파라미터
    t_val = stats.t.ppf(1 - alpha / 2, dof)
    
    # TODO: 각 파라미터에 대해 popt ± t_val * se
    # 반환: {"ec50_ci": (low, high), "hill_ci": (low, high), ...}
    pass

힌트: ec50_low = fit["ec50"] - t_val * fit["ec50_se"], ec50_high = fit["ec50"] + t_val * fit["ec50_se"].

빈칸 2: 96웰 배치 처리

플레이트 파일(CSV)에서 각 화합물별로 자동 피팅.

python

def batch_fit_plate(csv_path: str) -> "pd.DataFrame":
    """
    CSV: columns = [compound, dose, response]
    각 compound별로 fit_hill 호출, 결과 DataFrame으로 반환
    """
    import pandas as pd
    df = pd.read_csv(csv_path)
    
    results = []
    for compound, group in df.groupby("compound"):
        # TODO: doses, responses 추출, fit_hill 호출
        # 결과에 compound 이름과 함께 append
        pass
    
    return pd.DataFrame(results)

힌트: doses = group["dose"].values; responses = group["response"].values; fit = fit_hill(doses, responses); results.append({"compound": compound, **fit}).

빈칸 3: 다중 화합물 비교 플롯

여러 화합물의 곡선을 한 축에서 비교.

python

def plot_multi_compounds(fits: dict, doses_dict: dict) -> None:
    """
    fits: {compound_name: fit_result_dict}
    doses_dict: {compound_name: (doses, responses)}
    각 화합물의 피팅 곡선을 다른 색으로 겹쳐 그림
    """
    fig, ax = plt.subplots(figsize=(9, 6))
    
    # TODO: 각 화합물마다 다른 색으로 산점도 + 피팅 곡선
    # 범례에 각 화합물의 EC50 표시
    pass

힌트: colors = plt.cm.tab10.colors; for i, (name, fit) in enumerate(fits.items()): ax.scatter(..., color=colors[i]); ax.plot(..., color=colors[i], label=f"{name} EC50={fit['ec50']:.3g}").

성찰 — 실전 curve fitting 도구와의 차이

Robust regression: 실전은 outlier에 강한 fitting을 씁니다. 여러분의 curve_fit은 모든 데이터를 동등 취급. Huber loss나 RANSAC이 실전 대안.

Weighted fitting: 실전은 각 데이터 포인트의 측정 오차가 다릅니다. 예를 들어 낮은 반응에서 상대 오차가 큼. sigma 파라미터로 각 포인트에 가중치를 줄 수 있습니다.

모델 선택: Hill 방정식 외에도 여러 모델(4-parameter logistic, biphasic, sigmoid emax 등)이 있습니다. 실전은 AIC/BIC로 여러 모델을 비교해 최적을 선택.

GraphPad Prism: 임상 약리학 분야의 표준 도구. 여러분의 파이썬 접근은 자동화·재현성에서 강하지만, GUI 접근성에서는 Prism이 우위.

bootstrap 신뢰구간: t-분포 기반 신뢰구간은 오차가 정규분포라는 가정. 이 가정이 깨질 때 bootstrap resampling 으로 신뢰구간을 얻을 수 있습니다. 실전은 이 접근도 자주 씁니다.

확장 프로젝트

1. Streamlit 앱: 사용자가 CSV 업로드하면 자동 피팅 + 시각화 + PDF 리포트.

2. bootstrap 신뢰구간: 데이터를 리샘플링하며 여러 번 피팅해 EC50 분포를 얻어 신뢰구간 계산.

3. Dose-response DB: 여러 화합물의 결과를 SQLite에 축적하고 시간별 SAR(structure-activity relationship) 분석.

4. 병용 약물 상호작용: 두 약물의 조합에서 addition/synergy/antagonism을 판별하는 Bliss/Loewe 모델 구현.

이 편의 부품 지도

[F] numpy: 배열 조작, log 스케일 변환, 표준오차 계산.
[F] matplotlib: 로그 스케일 산점도 + 피팅 곡선 + 참조선.
[F] greedy 최적화: scipy.optimize.curve_fit 내부의 Levenberg-Marquardt 알고리즘 이해.
[W] 파일 I/O: CSV 파싱 (완성 스크립트 제공).

[F] = 여러분이 직접 구현 / [W] = 완성 코드로 제공.