All articles

Nimotuzumabの実測PKから考える、曝露比較に必要なモデルの複雑さ

Daiki Kumakura · 2026-10-04 · Technical note · 1.0

On this page

研究・学習目的の二次解析。未査読。

公開された個人別濃度データから、線形1区画、線形2区画、非線形消失を含む2区画のPopPKモデルを推定した。問いは「最も複雑なモデルはどれか」ではなく、投与条件による曝露の違いを比較するために、どの構造が必要で、どこまで結果を信用できるかである。

今回、複雑なモデルは学習データへの適合を改善した。しかし、被験者を除外した内部検証での改善は小さく、再標本化推定にも収束上の問題が残った。投与条件の順位は一致しても、絶対曝露と用量比は一致しない。したがって、非線形モデルを採用する根拠も、線形モデルで十分と断定する根拠も揃わなかった。この結論に至る計算と、その先に必要な情報を示す。

問いと用途

Nimotuzumabの原研究では、進行乳がん患者のデータに対して半機序的な非線形PKモデルが検討されている[1]。本解析は原論文モデルの完全再現ではない。公開パッケージの配布データを使って、用途に対するモデルの複雑さを比較する二次解析である。

  • QOI:50–400 mgの反復投与について、モデル構造を変えると、初回・第10回の曝露と投与条件間の比較がどれだけ変わるか。
  • COU:公開データを用いた研究・学習上の曝露比較、構造の必要性と追加測定の検討。
  • 対象外:治療用量の選択、個人の安全性・有効性予測、EGFR占有率、臨床試験の用量決定。

「同じ順番に並ぶ」と「数値が用途に十分な精度で一致する」は別の判断である。解析計画では、複雑化を支持する条件として、被験者単位の予測RMSEの10%以上の改善に加え、較正と識別性を要求した。簡単な構造で十分とするには、妥当な代替モデルに対して主要曝露・比の差が10%以内であることと、不確実性の評価を要求した。この10%は本研究の事前基準であり、一般的な規制基準ではない。

データを先に監査する

固定commitのnlmixr2data 2.0.10にあるnimoDataを使った[2]。原本RDAのSHA-256は896d256aee8919061b76358c5e723a786571182f88e5069bf5aa79fc24b82352。RDAからCSVに書き出し、投与と観測を分けて確認した。

項目 配布実体と扱い
被験者 12人。50 / 100 / 200 / 400 mgに各3人
イベント 441行:投与120、濃度観測321
反復投与 全員10回。実TIMEを連続使用
DV 自然対数濃度。exp(DV)がµg/mL(mg/Lと同値)
時刻・投与量・速度 h、mg、mg/h
点滴時間 記録AMT/RATEでは0.37–2.91 h。原研究の公称0.5 hとの対応は未解決
TAD 143/321観測で直前投与からの時間と不整合。ODEには使用しない
BLQ 専用フラグなし。未収録BLQの有無は判断できない
独立な情報量 321行ではなく、まず12人。各用量の独立被験者は3人

対数尺度は公式モデル例のIPRED = log(conc)と整合する[3]。投与行のDV=0を濃度観測に混ぜない。投与間隔は記録上93.0–312.6 hであり、推定時に名目168 hへ丸めていない。

A01の原論文監査では443濃度観測と記録したが、配布実体は321であり、抽出・加工の対応表はない[1,2]。解析対象は配布された321観測である。元試験の全データを再現したとは書けない。記録RATEを優先して推定し、全点滴を0.5 hへ置き換える場合は別の感度解析として再推定した。

3つの候補を同じ枠組みで推定する

濃度を C = A1 / V1、点滴入力を R_in(t) とする。M1は線形1区画、M2は以下の線形2区画である。

\[ \frac{dA_1}{dt}=R_{\mathrm{in}}(t)-CL\,C-Q\frac{A_1}{V_1}+Q\frac{A_2}{V_2},\qquad \frac{dA_2}{dt}=Q\frac{A_1}{V_1}-Q\frac{A_2}{V_2}. \]

M3はM2の中心区画消失に、次のMichaelis–Menten項を追加する。

\[ \frac{dA_1}{dt}=R_{\mathrm{in}}(t)-CL\,C-Q\frac{A_1}{V_1}+Q\frac{A_2}{V_2} -\frac{V_{\max}C}{K_m+C}. \]

M1ではQと末梢区画を除く。3候補共通でCLとV1に対数正規の個体間変動(IIV)を置き、共分散を0に固定した。

\[ CL_i=CL_{\mathrm{pop}}\exp(\eta_{CL,i}),\quad V_{1,i}=V_{1,\mathrm{pop}}\exp(\eta_{V,i}),\quad \eta_i\sim N(0,\Omega). \]

残差はlog濃度への加算誤差とした。

\[ y_{ij}=\log\{C_i(t_{ij})+10^{-12}\}+\epsilon_{ij},\qquad \epsilon_{ij}\sim N(0,\sigma^2). \]

M3のKmはEGFR結合親和性や受容体占有率ではない。公式例のTMDDモデルとも異なる。今回の非線形項はPK曲線を表す候補構造であり、生物学的な機序を同定したとは解釈しない。

推定にはR 4.5.1、nlmixr2est 7.1.0、rxode2 5.1.8を使用した。A02では3開始点ずつのSAEM、M2/M3では各3開始点のFOCEi精密化を行った。尤度の積分設定による違いを検出したため、同じデータ尺度で独立した適応Gauss–Hermite求積(25点)によりAICを比較した。15点との差も保存した。A03の再推定はFOCEi、上限outer 1000 / inner 300、ODE atol 1e-9、rtol 1e-6 とした[4]。

選択した候補解の推定値を次に示す。CL/QはL/h、V1/V2はL、Vmaxはmg/h、Kmはµg/mL。値は指数変換後の自然尺度である。sigmaはlog濃度のSD、omegaはlogパラメータの分散である。

Model CL V1 Q V2 Vmax Km sigma omega_CL omega_V1
M1 0.008629 1.775 — — — — 0.7459 0.2394 0.2786
M2 0.003797 1.47 0.004676 60.65 — — 0.7035 0.8028 0.2611
M3 0.001595 1.473 0.001861 8.279 2.372 396.7 0.6932 2.133 0.2803

共通尤度でのAICはM1 790.42、M2 751.25、M3 747.47だった。M3はM2より約3.78低い。しかし、A02のM3ではパラメータ相関の絶対値が最大約0.990であり、共分散とoptimizerの警告が残る。AICの順番だけで採用しない。

適合の良さと、新しい被験者への予測を分ける

モデルごとの適合と残差

図1. 学習データへのPREDとDV、IRESと時間・用量。対数尺度の図である。IRESは個人効果を推定した後の残差で、学習データへの当てはまりを示す。独立した予測検証ではない。

濃度・投与回数に対する残差

図2. IRESとIPRED、投与occasion。単純な平均残差だけでは、時間や反復投与に依存する偏りを見逃す可能性がある。図1と合わせて確認する。

個人効果が0へ縮む程度も確認した。SDで定義したeta shrinkageは次のとおり(%)。分散で定義した値もCSVに保存した。定義を混ぜない。低いshrinkageでも、モデルの識別性や予測妥当性が保証されるわけではない[5]。

model parameter SD_shrinkage
M1 eta.cl 0.3
M1 eta.v 4.792
M2 eta.cl 5.921
M2 eta.v 7.557
M3 eta.cl 13.75
M3 eta.v 6.2

被験者を1人ずつ除外し、残る11人から固定効果・IIV・残差を推定し直した(LOSO、12fold)。除外した被験者のDVから個人効果を推定せず、IIVを積分した周辺予測を評価した。各fold2048回のMonte Carlo近似を使用した。下表は、各被験者で計算したlog濃度RMSE・被覆率・点ごとの対数予測密度を、12人で同じ重みで平均した結果である。予測密度は反復観測全体のjoint densityではない。

model RMSE RMSE improvement vs M1 (%) coverage mean_lpd
M1 0.9458 0 0.9186 -1.386
M2 0.9061 4.195 0.9218 -1.347
M3 0.9147 3.288 0.9189 -1.332

被験者を除外した予測の比較

図3. 同じ被験者について、モデルごとのRMSEと90%予測区間の被覆率を比較。点の連結は対応関係を示す。観測数が多い被験者だけを重視しない。

M2/M3の平均改善は10%に届かない。別seedの8192点Sobol積分でも改善はそれぞれ4.10%、3.16%で、この判断は変わらなかった。ただし、この結果を「複雑化は不要」と結論しない。LOSOには収束警告が残っており、下表のとおり、全foldが収束スクリーニングを通過したわけではない。表の指標は全ての有限解による探索的な内部診断である。全体データで選んだ初期値を使うため、モデル選択を含むnested validationではなく、外部検証でもない。

Model Folds Finite Pass convergence screen
nlminb-M1 12 12 5
M2 12 12 0
M3 12 12 1

VPCは投与履歴を保って行う

実際の投与時刻・速度・採血時刻を保持し、モデルごとに1000データセットを生成した。同一被験者の反復観測では同じetaを使い、観測ごとに残差を加えた。初回と第10回を用量別に分け、直前投与からの時間で集計した。時間binは0–2、2–30、30–72、72–168、168–360、360–720 h。記録スケジュールの範囲を覆うように設定し、3観測未満のbinは表示しない。

用量別・初回と第10回のVPC

図4. 黒が観測中央値、色線が反復シミュレーションの中央値、帯が各binの中央値の90%シミュレーション区間。縦軸は対数濃度、横軸はbin上端であり採血点そのものではない。完全なVPC出力には10/50/90百分位も保存したが、図は中央値だけを示す。帯は個体間分布の90%区間でも、パラメータの信頼区間でもない。

各用量3人のため、binに複数観測があっても独立被験者は増えない。一部の用量・時期に中央値のずれがあり、広い帯への包含だけで適合を保証しない。今回はdose別の非補正VPCであり、prediction-corrected VPCではない。

bootstrapから信頼区間を無理に作らない

用量層を保ち、被験者を復元抽出してモデルごとに200回再推定した。観測行をばらばらに再抽出しない。3人しかいない層の経験分布であるため、少数集団の情報不足は残る。M1ではbobyqaに加えてnlminbへoptimizerを変え、さらに200回を実行した。追加試行を元の結果と置き換えず、両方残した。

group attempted finite qualified qualified_rate CI_defensible
M1 200 200 9 0.045 False
M2 200 200 9 0.045 False
M3 200 200 11 0.055 False
nlminb-M1 200 200 44 0.22 False

bootstrapの有限出力と収束判定

図5. 灰色は有限推定値が得られた数、青色は指定の収束スクリーニングを通過した数。not at minimum、最大反復、failure、convergence関連の記録を全件確認した。通過していてもgradient等の注意記録が残り得るため、通過数を「完全に警告なしの成功数」と解釈しない。

全候補が事前の80%基準に届かない。したがって、このbootstrapによる95%信頼区間は報告しない。成功した少数だけを選んで区間を作ると、失敗するデータ構成を除外した分布になる。有限出力、収束記録、推定値自体は監査用に保存した。

非線形消失パラメータの有限グリッド探索

図6. KmまたはVmaxを候補解の1/64〜64倍(7点)に固定して残るパラメータを再推定した。縦軸はA03基準FOCEi目的関数との差。これは有限グリッドの探索であり、AICに使った独立求積の曲線でも、正式なprofile信頼区間でもない。収束警告付きの点も隠さず白抜きで示す。曲線を滑らかな真のprofileとみなさず、区間の有界性を断定しない。

投与条件の順位が同じでも、数値は同じではない

推定時には実際の投与履歴を使った。一方、比較シミュレーションでは条件を揃え、50 / 100 / 200 / 400 mg、168 h間隔、0.5 h点滴、10回投与を設定した。固定効果を候補解に固定し、同じ標準正規乱数から各モデルのOmegaに従う10000人を生成した。各モデル内では同一仮想個体へ全条件を与えるpaired比較である。モデル間の同じ乱数は、同一実患者の同じ生物学的個体差を保証しない。

濃度から初回0–168 hと第10回1512–1680 hのAUC、Cmax、Ctrough、Cavgを計算した。第10回を定常状態と呼ばない。Cmaxは点滴終了を含む評価グリッド上の最大値で、AUCは数値積分である。IIVを示す曝露シミュレーションには残差を加えていない。

下表は第10回AUC(µg·h/mL)の10/50/90百分位と、各仮想個体の200 mg条件に対する比の中央値である。「中央値の比」と「個体ごとの比の中央値」を区別する。

model amount_mg p10 median p90 Paired median ratio vs 200 mg
M1 50 3102 5781 1.046e+04 0.25
M1 100 6205 1.156e+04 2.092e+04 0.5
M1 200 1.241e+04 2.312e+04 4.184e+04 1
M1 400 2.482e+04 4.625e+04 8.369e+04 2
M2 50 3107 6213 9010 0.25
M2 100 6213 1.243e+04 1.802e+04 0.5
M2 200 1.243e+04 2.485e+04 3.604e+04 1
M2 400 2.485e+04 4.97e+04 7.208e+04 2
M3 50 2891 5905 7086 0.2141
M3 100 5879 1.248e+04 1.534e+04 0.4524
M3 200 1.212e+04 2.758e+04 3.563e+04 1
M3 400 2.541e+04 6.497e+04 9.158e+04 2.351

第10回の絶対曝露とpaired用量比

図7. 左はIIVによるAUCの中央値と10–90百分位、右は同一仮想個体内のAUC比。帯はパラメータ推定の信頼区間ではない。

M1/M2では線形性により400 mg対200 mgのAUC比は厳密に2となる。M3のpaired中央値は約2.35となった。この比の差は事前の10%基準を超える。一方、M3は非線形項の識別と再推定の安定性に問題を残す。したがって、「2.35が正しい用量比」とも、「線形の2で十分」とも判断できない。モデル依存性を確認できたことが、ここで支持できる成果である。

200 mgを144 / 168 / 192 h間隔で10回投与する条件も比較した。総投与量は同じだが、総期間が異なるため、AUCだけでなく第10回区間のCavgを比較する。短い間隔で濃度が上がることを、安全性・有効性上の利点と読み替えない。

投与間隔と第10回平均濃度

図8. 各候補解のIIVによる中央値と10–90百分位。10回目の区間長自体が条件によって異なる。固定効果の推定不確実性は帯に含まない。

同じパラメータのまま0.5 h点滴を1 hにすると、初回Cmaxのpaired中央値の変化は3モデルとも約−0.1〜−0.2%だった。これは「同じ薬物動態で点滴だけを変更した」比較であり、次の入力履歴感度とは異なる。

入力履歴と体重の仮定を変更するとどうなるか

被験者を除いた再推定解による曝露も調べた。次は第10回の典型個体AUC200 mgと、400/200 mg比の12除外解における範囲である。個体間分布や信頼区間ではない。

model endpoint minimum maximum screen_pass_folds
M1 typical_AUC200 2.126e+04 2.513e+04 5
M1 typical_ratio400_200 2 2 5
M2 typical_AUC200 2.342e+04 2.723e+04 0
M2 typical_ratio400_200 2 2 0
M3 typical_AUC200 2.558e+04 3.113e+04 1
M3 typical_ratio400_200 2.187 2.466 1

M3の用量比は1人ずつ除外しても約2.19–2.47で、線形の2と異なった。ただし、収束スクリーニングを通過したM3は1/12解だけである。単一被験者だけで変化するかを調べた感度として残し、頑健な集団予測とはしない。Km/Vmaxの除外解も別のCSVに保存した。

全ての記録点滴を0.5 hへ置き換えた再推定、初回体重を使う固定allometry(CL指数0.75、V1指数1、基準70 kg)の再推定を、元の候補と別に保存した。指数は推定していない。次は各再推定解の200 mg週1回・10回・0.5 h点滴における典型個体の第10回AUCである。allometryは70 kgで比較する。

group allometry-fixed baseline infusion-0.5h infusion-0.5h / baseline allometry-fixed / baseline
M1 2.164e+04 2.299e+04 2.293e+04 0.9973 0.9412
M2 2.406e+04 2.483e+04 2.499e+04 1.007 0.9693
M3 2.759e+04 2.751e+04 2.797e+04 1.017 1.003
nlminb-M1 2.159e+04 2.317e+04 2.28e+04 0.9839 0.9319

各感度再推定にも収束警告が残るため、差を頑健な共変量効果や正しい点滴時間の証明として扱わない。さらに、記録された初回体重を復元抽出した10000人について、体重とIIVを加えた曝露も計算した(allometry-exposure.csv)。これは観測された12人の体重分布に条件付けた感度であり、一般患者集団の分布ではない。未知の患者背景のjoint distributionを合成したとは主張しない。

数値検証で確かめたこと

独立したPython実装(線形モデルは解析的遷移、M3はDOP853)とRの予測を、基準・点滴感度・allometryの9推定解で比較した。最大log予測差は3.3e-06だった。積分グリッド337→673点の変更によるAUCの最大相対差は3.6e-07(32人の同一乱数コホート)だった。M1/M2のpaired用量比が線形理論値と一致すること、LOSOが36被験者fold・963観測評価、曝露表が168行であることも確認した。

これらは実装と数値積分の確認であり、データの由来、推定の収束、モデルの生物学的妥当性を解決する検証ではない。

この解析から何を判断できるか

判断 今回支持できる範囲 支持できない範囲
実測データによるPopPK 個体間変動を含む3候補の推定・比較を実行した 元試験全体の完全再現
構造の追加 学習適合は改善するが、予測改善は小さく不安定性が残る 非線形モデルの採用確定
線形モデルの十分性 単純な基準としての比較に使える 用量比・絶対曝露が十分正確との保証
投与条件比較 固定した候補モデル内の曝露・比の違いを数値化できる 治療用量の最適化、臨床リスク予測
不確実性 IIV、構造差、再推定失敗を区別できる bootstrap由来の信頼できる95%CI

次に必要なのは、単に別の複雑な式を追加することではない。まず、点滴開始・終了と採血時刻の原記録、配布対象の選択・加工履歴を確認する。そのうえで、初期分布相と低濃度の末端相を区別できる測定、用量ごとの独立被験者、選択に使っていない検証データを増やす必要がある。

既存データへの追加測定を設計する場合には、低濃度域と高濃度域の候補モデル予測差を、測定誤差・採血負担と並べて評価するのが次の問いになる。今回の解析だけでは最適採血時刻や必要人数を求めていないため、具体的な設計値は提示しない。

MIDDの観点で得た結論は、モデルを複雑にするかどうかは、適合の改善だけでなく、その構造が比較したい数値を変えるか、そしてデータから安定して推定できるかで決めるということである。今回、数値への影響は確認できたが、安定した推定と独立検証は揃わなかった。採用保留と追加情報の特定までが、このデータで支えられる判断である。

再実行資料と出典

解析記録は実行済みa03-diagnostics.ipynbに保存した。配布用には、患者別データなしで動くレビューNotebook、再実行手順、モデル定義・R再推定スクリプト、Pythonの独立予測・診断・シミュレーション、集計CSV、図を準備した。レビューNotebookは保存結果から比較と検算を再実行する。長時間のR推定そのものは別コマンドである。配布用資料に患者別データやRDSは含めず、固定版の取得元と原本hashを記載した。コードの利用条件と、元データの利用権限・原論文の著作権は分けて扱う。

  1. Rodríguez-Vera L, et al. Semimechanistic model to characterize nonlinear pharmacokinetics of nimotuzumab in patients with advanced breast cancer. J Clin Pharmacol. 2015;55:888–898. DOI:10.1002/jcph.496. A01で監査した著者公開本文.
  2. nlmixr2data:解析に用いた固定commit. パッケージ2.0.10、原本hashは本文参照。
  3. nlmixr2公式Nimotuzumab例. データ尺度と原研究モデルへの参照。今回のM3をこのTMDDモデルの再現と呼ばない。
  4. nlmixr2est公式foceiControlリファレンス. optimizerと推定制御。
  5. Savic RM, Karlsson MO. Importance of shrinkage in empirical Bayes estimates for diagnostics: problems and solutions. AAPS J. 2009;11:558–569. DOI:10.1208/s12248-009-9133-0、公開全文.

本記事の表・図は配布データを用いた本解析から作成した。引用論文の図を転載していない。

解析資料一式(コード・集計結果・実行済みNotebook)