容量-反応曲線フィッティング — NumPy を使用して EC50 と信頼区間を自動計算
このトピックを終えると
教科書で学んだ numpy、matplotlib、貪欲法 を組み合わせて、薬物スクリーニングの用量反応データから EC50 と信頼区間を自動計算するツールを自分で作成できます。 これまで 96 ウェルプレートの結果を一つずつ手で計算していた時間が、30 秒に短縮されます。
この記事は 教育目的の一般的な例 です。 実務では、GraphPad Prism、R の drc パッケージなどを使用します。 ここでは、これらのツールの内部原理を正確に説明します。
「これを全部手作業でやっているのか？」—繰り返し計算の落とし穴
あなたが新薬スクリーニングを行っているとしましょう。96ウェルプレートから、次のデータが得られます。
- 各ウェルに異なる濃度の候補化合物を入れる（例：0.001、0.01、0.1、1、10、100 μM）
- 各濃度で細胞生存率（%）を測定
あなたが知りたいこと：EC50（反応の半分を引き起こす濃度）と最大反応。
最も単純なアプローチ：Excelにデータを貼り付け、散布図を描き、目で見て半分になる点を特定します。
このアプローチの現実的な問題点：
問題1：正確性。目で描いた曲線は人によって異なります。再現性はありません。
問題2：規模。候補化合物100個 × それぞれ8つの濃度 = 800データポイント。Excelのクリック作業では、一日がかりです。
問題3：信頼区間。EC50が例えば2.3 μMと出た場合、この値が「実際にどれくらい正確なのか」がわかりません。
本当のアプローチは非線形曲線フィッティングです。ヒル方程式（シグモイド関数）をデータに最適にフィッティングするパラメータを見つけ、そのパラメータの統計的な信頼区間まで計算します。
ブラックボックスからコンポーネントへ
コンポーネント1：ヒル方程式
用量反応曲線の標準的な数学モデルは、ヒル方程式です。
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)通常、用量を対数変換してプロットすることで、対数スケールで視覚化しやすくなります。
コンポーネント2：貪欲法による最適化
4つのパラメータをどのように求めるか？ データと予測との誤差（残差）を最小化すること。最小二乗法。
minimize sum_i (y_i - y_hat_i)^2非線形最小二乗法は、閉じた解がないため、反復アルゴリズムで解きます。Levenberg-Marquardt法、trust-region法など。これらのアルゴリズムは、貪欲法のアプローチで動作します。現在の位置から誤差が最も減少する方向に少し移動し、その位置で繰り返します。
scipy.optimize.curve_fitは、このアルゴリズムの実装です。
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による視覚化
フィッティング結果を確認するための標準的なグラフ。
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()実践的な例
仮想データで確認します。
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)期待される出力：
EC50: 0.271 ± 0.0084
Hill: 0.98
Range: 1.2 → 99.5フェーディング — 埋めるべき3つの空白
空白1：95%信頼区間
標準誤差だけでは不十分です。95%信頼区間をt分布で計算します。
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）から、各化合物ごとに自動フィッティングを行います。
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名とともに追加します。 pass return pd.DataFrame(results)ヒント: doses = group["dose"].values; responses = group["response"].values; fit = fit_hill(doses, responses); results.append({"compound": compound, **fit}).
空白3：複数化合物の比較プロット
複数の化合物の曲線を、同じ軸で比較します。
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_fit は、すべてのデータを同等に扱います。 Huber損失やRANSACが、実用的なツールの代替手段です。
重み付きフィッティング：実用的なツールは、各データポイントの測定誤差が異なることを考慮します。例えば、低い反応の場合、相対誤差が大きくなります。sigma パラメータを使用して、各ポイントに重みを付けることができます。
モデル選択：ヒル方程式以外にも、さまざまなモデル（4パラメータロジスティック、2相モデル、シグモイドEmaxなど）があります。実用的なツールは、AIC/BICを使用して複数のモデルを比較し、最適なモデルを選択します。
GraphPad Prism：臨床薬理学分野における標準的なツールです。あなたのPythonによるアプローチは、自動化と再現性において優れていますが、GUIの使いやすさにおいてはPrismが優れています。
ブートストラップによる信頼区間：t分布に基づく信頼区間は、誤差が正規分布であるという仮定に基づいています。この仮定が成り立たない場合、ブートストラップリサンプリングによって信頼区間を得ることができます。実用的なツールでは、このアプローチも頻繁に使用されます。
拡張プロジェクト
1. Streamlitアプリ: ユーザーがCSVファイルをアップロードすると、自動的にフィッティング、可視化、PDFレポートを作成。
2. ブートストラップ信頼区間: データをリサンプリングしながら複数回フィッティングを行い、EC50の分布を得て信頼区間を計算。
3. 用量反応データベース: 複数の化合物の結果をSQLiteに蓄積し、経時的なSAR（構造活性相関）分析を行う。
4. 併用薬物相互作用: 2つの薬物の組み合わせにおいて、相加作用/相乗作用/拮抗作用を判別するBliss/Loeweモデルを実装。
この演習の部品マップ
- [F] numpy: 配列操作、対数スケール変換、標準誤差の計算。
- [F] matplotlib: 対数スケールの散布図 + フィッティング曲線 + 参照線。
- [F] 貪欲法による最適化:
scipy.optimize.curve_fit内部のLevenberg-Marquardtアルゴリズムを理解する。 - [W] ファイル入出力: CSVの解析（完成したスクリプトを提供する）。
[F] = 自分で実装 / [W] = 完成したコードとして提供。