affinityfit

滴定データ（濃度とシグナルの 2 列）から Kd などを求め、その推定値を信用してよいかを 診断する Python ライブラリ。

観測量が結合率と線形の関係にあれば、測定手段を問わず同じ数式（飽和曲線）が使えます。 NMR のピーク強度、SPR の定常状態応答、酵素反応の初速度がその例です。

signal = baseline + Bmax * [L] / (Kd + [L])

複数のデータセットを同時にフィッティングし、パラメータをデータセット間で共有したり 定数に固定したりできます。測定範囲が Kd を挟んでいないために単独では決定不能な曲線 でも、共有によって推定が成立します。

Usage

Python >= 3.12

uv add affinityfit

from affinityfit import fit, load_csv

conc, signal = load_csv("titration.csv")
res = fit(conc, signal, unit="nM")

print(res.report())
print(res.params["kd"], res.intervals["kd"].format("nM"))
for warning in res.warnings:
    print(warning)

入力 CSV は 1 列目が濃度、2 列目がシグナル。ヘッダ行・コメント行・空行は読み飛ばします。 同じ濃度の行が複数あればレプリケートとしてそのまま使われます。

欠測（空欄）、nan、inf は行番号つきのエラーになります。空欄の行を黙って捨てると 点数が変わり、自由度と診断が静かに変わってしまうためです。

concentration_nM,signal
0,0.0249
1,0.0963
3,0.2493
10,0.5090
30,0.7657
100,0.9249
300,0.9758
1000,1.0035

report() の出力:

model    : langmuir  (1:1 binding: signal = baseline + Bmax * L / (Kd + L))
interval : profile
Kd       = 10.1 +/- 0.8 nM
Bmax     = 0.995 +/- 0.017
baseline = 0.018 +/- 0.014
R^2      = 0.9998   (n = 8)
AICc     = -71.12   (AIC = -77.12)

診断チェック: 問題は検出されませんでした。

固定側（受容体・レクチンなど）の濃度がわかっている場合は渡してください。 リガンド枯渇のチェックが有効になります。

res = fit(conc, signal, receptor_conc=1.0, unit="nM")

モデルを変えるには model= を渡します。

from affinityfit import hill

res = fit(conc, signal, model=hill)  # 協同性を検定する
print(res.intervals["n"].contains(1.0))

複数データセットを同時にフィッティングし、パラメータを共有・固定する場合:

from affinityfit import Dataset, fit_global

res = fit_global(
    [Dataset("oxidized", conc, sig_ox), Dataset("reduced", conc, sig_red)],
    shared=["bmax"],  # 全データセットで 1 つの値を推定
    fixed={"baseline": 0.0},  # 定数に固定して推定しない
    unit="mM",
)
print(res.report())

# 1 データセット分を取り出す。そのデータセットの警告・注記も一緒に付いてくる。
sub = res.result_for("oxidized")
print(sub.report())
x, y = sub.curve()

warnings は問題、notes は「共有したおかげで推定できている」のような補足です。 GlobalFitResult.warnings / .notes は [データセット名] の接頭辞つきで全件を返し、 result_for() はそのデータセット分（およびフィッティング全体に関わるもの）を接頭辞 なしで返します。

Weighting

既定では全点を等価値に扱います。これは「誤差の大きさがどこでも同じ」という仮定で、 SPR の RU のような系では妥当ですが、蛍光・発光・吸光では誤差が信号に比例するため 成り立ちません。飽和曲線は信号が baseline から baseline+Bmax まで動くので、飽和側の 点の絶対誤差が低濃度側の数十倍になります。

res = fit(conc, signal, sigma=0.01 + 0.10 * signal)  # 点ごとの標準偏差

sigma は相対値として扱われ、全体のスケールは残差から推定されます。定数倍しても 結果は変わりません。グローバルフィットでは Dataset(..., sigma=...) でデータセット ごとに渡せ、データセット間の相対的な重み（片方が 10 倍ノイジー等）も反映されます。

sigma は全データセットに与えるか、どこにも与えないかのどちらかです。一部にだけ 与えると、sigma のないデータセットの重みが暗黙に 1 になり、両者の重み比が sigma の 絶対スケール（分数か % か ppm か）で決まってしまうためです。混在はエラーになります。

重みを渡さずに不均一分散のデータをフィッティングすると、精度が落ちるだけでなく 信頼区間が狭く出ます。Kd=10 に対し 30% の比例誤差、Kd が測定範囲上端にある設計 での実測です。

残差の大きさがフィッティング値とともに増えている場合は診断が警告します（誤報率 2%、 比例誤差の検出率 56〜71%）。

Confidence Intervals

ci= で 3 つの方法を選べます。

反復実験から誤差を出す慣行はこの分野で広く使われています。反復測定を渡せばその 流儀に揃えられます。

res = fit(conc, signal, ci="bootstrap", replicates=reps, n_boot=2000)

Dataset では反復測定だけを渡せます（signal は平均が使われます）。

Dataset("oxidized", conc, replicates=reps)

n_boot は百分位区間に必要な下限（100）以上でなければエラーになります。下回った値を 渡すと区間が作れないためです。

報告される桁数は不確かさに合わせて丸められます。相対誤差 17% の測定から Kd = 4.70e-8 と 3 桁で書くのは過剰なので、(4.7 +/- 0.8)e-08 になります。

Model

langmuir と hill の bmax は符号を制限していません。結合によって減少する観測量 （蛍光クエンチ、強度の減少）でも Kd の意味は変わらず、負の応答係数として同じ式で表せます。 michaelis の vmax は反応速度なので非負に留めているため、減少するデータには langmuir か hill を使ってください。噛み合わないモデルを指定した場合は、それらしい 数値を返さずに警告します。

Km = (k_off + k_cat)/k_on であり、k_cat が無視できるときだけ Kd に一致します。

Kd や Km は対数スケールで最適化されます。線形に扱うとフィッティング結果が濃度の単位に 依存し、M 単位のまま pM 親和性（Kd = 1e-12 M）を扱うと相対誤差が 6800% に達します。 現在は fM から 100 M まで 18 桁にわたって相対誤差 1e-6 以内で回収し、同じ実験を M / nM / pM のどの単位で書いても同じ答えになります。 モデルの追加は models.py に Model を 1 つ定義するだけで、フィッティングと診断の コードには手を入れません。

Model Selection

モデルやパラメータ共有の比較には aicc（小標本補正つき赤池情報量規準）を使ってください。 aic も参照用に持っていますが、補正なしの AIC は漸近的にしか妥当でなく、滴定実験の 規模（点数 6〜15、パラメータ 3〜6 で n/k が 2〜5）ではパラメータの多い側に偏ります。

実測: 振幅が真に共通なデータ（12 点、共有で 3 パラメータ / 非共有で 4）40 通りで、 補正なし AIC が共有を選んだのは 32 回、AICc は 38 回でした。10 点・5 対 6 パラメータの 構成では 12 通り中 7 通りで結論が食い違い、いずれも AIC が過剰パラメータ側を選びました。

shared.aicc < free.aicc
fit_global(ds, model=hill).aicc < fit_global(ds, model=langmuir).aicc

標本が小さすぎて補正が定義できない場合（n − k − 1 ≤ 0）、aicc は無限大になります。 report() は AICc を先に、参照用の AIC を括弧内に表示します。

Diagnostics

R² が 0.99 を超えていても Kd が決まっていないことは珍しくありません。以下を自動で警告します。 Kd / Km のようにモデルごとに呼び名が変わる量は、そのモデルの表記で表示されます。

診断は2種類あります。フィッティング自体が健全かと、測定点の配置が適切かです。

フィッティング自体が壊れている場合（R² < 0.5、振幅が潰れた）、「共有したおかげで 推定できている」という NOTE は出しません。矛盾した助言になるためです。 | 自由度が 0（点数 ≤ パラメータ数） | 曲線は定義上すべての点を通る。信頼区間は算出できない | | ヤコビ行列がランク落ち | パラメータが区別できず値が一意に決まらない。濃度の種類を増やす | | 信頼区間の片側が決定できない | 点推定値を報告してはいけない。共有か測定範囲の拡張が必要 |

グローバルフィットでは、共有・固定によって前提が消えた指摘は抑制されます。たとえば bmax を共有していれば「Kd と Bmax を分離できない」は当てはまらないので、代わりに 「共有しているため推定できている」という NOTE が出ます。

examples/titration_unsaturated.csv は R² = 0.9936 と当てはまりが良く見えても、最高濃度が Kd の 1/3 しかないため上限が決まらず Kd > 76 nM としか言えない例です。

信頼区間が算出できない場合は数値を作りません。自由度が 0 のとき、ヤコビ行列が ランク落ちしているとき、プロファイル探索の内側の再フィットが収束しなかったとき、 ブートストラップの成功回数が足りないときは、区間を「決定不能」として返し、理由を 警告で説明します。点推定値そのものは返します。

ブートストラップは収束しなかった再抽出の回数を数えて報告します。収束する再抽出は フィッティングしやすいものに偏るため、黙って捨てると区間が系統的に狭くなります。

残差が厳密に 0 のフィッティング（合成データなど）では区間幅も 0 になりますが、裸の 数値ではなく 10.000 +/- 0 (no residual scatter) と理由つきで表示します。

残差の検定にはラン検定（Wald–Wolfowitz）と隣接残差の自己相関の両方を使い、どちらかが 有意なら警告します。15 点の滴定を 300 通りのノイズで試した実測で、誤ったモデル （協同性データを 1:1 でフィッティングする、R² は 0.94〜0.98）の検出率は 100%、 正しいモデルでの誤検出率は 0.3〜1.7% でした。

統計量

警告の判定は固定の有意水準（例: 残差検定は片側 2.5%、不均一分散は片側 1%）で行われて います。これは 1 件のフィッティングには適切ですが、複数のデータセットや複数の論文図に わたる検定を積み重ねる場合には、その有意水準は検定の数だけ増えた偽陽性率を勝手に 受け入れていることになります。何本の検定になるか、それをどう補正すべきかはこちら側では 決められないので、判定の元になった統計量と p 値をそのまま返します。

res = fit(conc, signal)
for s in res.statistics:
    print(s.name, s.statistic, s.p_value)

複数データセットの場合は GlobalFitResult.statistics_per["データセット名"] で データセットごとに取得できます。集めた p 値に対して、scipy.stats.false_discovery_control や statsmodels.stats.multitest.multipletests で Bonferroni や Benjamini–Hochberg 補正をかけられます。

from scipy.stats import false_discovery_control

p_values = [
    s.p_value
    for name, stats in res.statistics_per.items()
    for s in stats
    if s.name == "heteroscedasticity" and s.p_value is not None
]
q_values = false_discovery_control(p_values, method="bh")

residual_autocorrelation は固定閾値（0.3）による判定で、p 値を持ちません （p_value は None）。

Plotting

図を返す API はありません。論文用の図はフォント・色・パネル構成を自分で制御する必要が あるので、affinityfit は描画に必要な数値だけを返します。

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(*res.curve(), "-")
ax.plot(conc, signal, "o")
ax.set_xscale("log")

濃度軸は対数にしてください。線形軸だと低濃度側が潰れ、Kd 近傍の曲率と「測定点が 飽和領域に偏っている」といった設計上の問題が見えなくなります。

examples/plot_fit.py に、測定点・フィッティング曲線・残差パネルを 1 枚に描くサンプル があります。配布物には含まれません。

uv run python examples/plot_fit.py examples/titration_good.csv --unit nM

examples/titration_good.csv にこのサンプルを実行した結果:

Development

uv run pytest -q
uv run ruff check .
uv run ruff format --check .
uv run ty check

テストは既知のパラメータ（Kd は 1e-16 〜 1e2 の 18 桁、Hill 係数は 0.4 〜 3）を回収 できるか、各警告が意図した条件でのみ発火するか、共有パラメータが識別不能なデータセット を救うかを検証しています。信頼区間については、良条件のデータで 3 手法が一致すること、 プロファイル区間の端が F 検定の閾値に厳密に一致すること、識別不能なデータで片側限界に なることを確認しています。examples/ のサンプルもテストから実際に呼んでいるので、 API を変えれば落ちます。

Roadmap

未実装のもの。上の各節に書かれている機能はすべて動きます。

モデル

• 時間経過モデル A(1 - exp(-k·t))。プラトーに達していないと振幅と速度定数が 分離できないため、その判定つきで
• 競合阻害モデル（阻害剤存在下の見かけの Kd から真の Ki を出す）
• tight-binding（二次式）モデル。受容体濃度が Kd と同程度の場合に必要
• 山型のプロファイル（金属濃度やイオン強度の最適値、pH 最適）。単調でないため 飽和曲線とは別の族になる

誤差モデル

• sigma を自前で用意できない場合の分散モデル（proportional / poisson）。 当てはめ値から重みを作るため反復再重み付けが必要で、当てはめ値が 0 の点で 重みが発散する。σ の下限を含めた定式化が前提
• 反復測定から sigma を推定する。反復数が少ないと標準偏差自体がばらつくので、 そのばらつきを織り込む必要がある

入出力

• 複数データセットの結果表を CSV / JSON で出力するヘルパ
• 装置固有のエクスポート形式（SPR のセンサーグラムなど）のインポータ
• LOD / LOQ の算出（IUPAC / CLSI 準拠）