1. はじめに：プレプリント公開後に直面した「AI共創ならではの落とし穴」
先日、40年前の生物有機化学・ホストゲスト化学の知見と最新のAIエージェントを掛け合わせ、グラフ理論（ラプラシアン行列のFiedlerベクトル）と動的極性ロッキングによるタンパク質フォールディング・シミュレータを開発し、Zenodoにてプレプリント論文（Version 1）および GitHubリポジトリ（fiedler-water-folding-sim） を公開しました。
論文タイトル
Graph-Spectral Protein Folding: Simulating Peptide Self-Assembly via Laplacian Fiedler Vector Optimization
（グラフスペクトル・タンパク質フォールディング：ラプラシアンFiedlerベクトル最適化によるペプチド自己組織化のシミュレーション）
科学哲学者カール・ポパーが提唱した**「反証可能性（Falsifiability）」を道標に、AIと対話しながら仮説を立て、自ら反証実験を繰り返して組み上げた研究です。しかし、V1を公開した直後、AIを活用した研究特有の「データとコードのトレーサビリティ（追跡可能性）の危機」**に直面しました。
長期間にわたるAIとのセッションでは、対話ログ（JSONL）が数万行・数十MBへと肥大化し、コンテキストの圧縮（Compaction）や複数セッションへの分岐が発生します。その中で、以下のような問題が潜在していることに気づいたのです。
- プロトタイプ段階の仮データと最終計算データの混在：初期にグラフのレイアウト確認用として書いたスクリプトの数値や、感度分析のモックアップコードが、最終原稿の図表に一部残っていないか？
-
複数の実験条件・対象分子の取り違え：野生型Chignolin（
1UAO）と超安定変異体CLN025（5AWL）の配列や慣性半径（$R_g$）の定義が、本文中で混同されていないか？ - 開発環境のクォータ制限というタイミング：ちょうどメインで使っていたAI開発環境（Antigravity）が月次の利用制限に入り、28日まで復活しないという状況に直面しました。
「公開した論文に1ミリでも曖昧な数字や検証不足の箇所を残したくない」——そう考え、手元に残っていた全実験データ（PDB構造ファイル、トラジェクトリJSON、シミュレータ本体コード）を別の検証環境へ持ち込み、コード1行・数値1桁単位での徹底的なデータ監査（Data Audit）と再計算を実施しました。
本記事では、その監査で見つかった「5つの課題」と、それらをすべて実計算で置き換えて Zenodo Version 2（V2） へと昇華させたプロセス、そして厳密に再検証したからこそ見えてきた**新しい生物物理学的発見?**を包み隠さず共有します。
Zenodo : https://zenodo.org/records/23667251
2. 徹底監査で見つかった5つの課題（自己反証とコード監査）
自分自身の論文とスクリプトに対して「最も厳しい査読者」としてメスを入れた結果、V1には以下の5つの改善・修正すべき点があることが判明しました。
課題①：野生型 Chignolin（1UAO）と変異体 CLN025（5AWL）の配列・ID表記の混同
V1のSection 3では、テストケースとして PDB ID: 1UAO を挙げながら、アミノ酸配列として YYDPETGTWY を記載していました。
しかし厳密には、本田ら（2004）が報告したオリジナルの野生型Chignolin（1UAO）の配列は GYDPETGTWG（両末端がグリシン、重原子数77）であり、両末端をチロシンに置換して結晶化・超安定化させた変異体（本田ら, 2008）が CLN025（5AWL、配列：YYDPETGTWY、重原子数93） です。長いセッションの中で、両者の表記が1つのセクション内で混ざっていました。
課題②：慣性半径（$R_g$）における対象原子サブセット（$\text{C}\alpha$ / 重原子 / 全原子）の未明記
V1では実験値のネイティブ慣性半径を単一の数値（$5.17\text{ \AA}$ または $5.26\text{ \AA}$）として引用していましたが、それが「$\text{C}\alpha$炭素のみ」なのか「水素を除く重原子（Heavy atoms）」なのか「水素を含む全原子（All atoms）」なのかを明記していませんでした。10残基程度の短いペプチドでは、側鎖の芳香環や水素を含めるかどうかで $R_g$ が $0.6\text{ \AA}$ 以上変動するため、比較の基準を統一する必要がありました。
課題③：図表生成スクリプト内のハードコード配列と、全体RMSD指標の不足
V1の図生成スクリプト（generate_paper_figures.py）を点検したところ、トラジェクトリ曲線を描画する部分に、初期のプロトタイプ段階で作成した例示用の数値配列（ハードコードされた配列）が残っていました。また、構造の正確性を主に $R_g$ とターン部位（Asp3）の局所偏差（$1.67\text{ \AA}$）で議論しており、構造生物学の標準指標である Kabsch法（特異値分解・SVD）による全鎖 $\text{C}\alpha$ RMSD の明示が不足していました。
課題④：感度分析（Table 2）のプロトタイプコードに残っていた擬似乱数処理
V1でハイパーパラメータの頑健性を示すために作成した感度分析スクリプトを確認したところ、ベースラインの計算結果に対して np.random.uniform(0.98, 1.02) を掛け合わせるモックアップ処理が残っていました。論文としてパラメータの頑健性を主張する以上、乱数による近似ではなく、全パラメータ組み合わせに対して物理シミュレータを1ステップ目から決定論的に再実行した実測値を載せる必要があります。
課題⑤：「最適化ループの中でPDBの正解座標をカンニングしていないか？」のコード証明
これはデータの誤りではなく、独自のアルゴリズムを公表する上で最も重要な**「カンニングなし（Zero Target-Coordinate Leakage）の証明」**です。サーバスクリプト（server_overall_fiedler_polar_priority.py）の冒頭では、比較参照用および3Dビューア表示用に 1UAO.pdb を読み込んでいます。第三者から見て「最適化ループの中で、正解のPDB座標に近づくようなペナルティ項がこっそり入っているのではないか？」という疑念を完全に晴らすため、コードの行単位での監査が必要でした。
3. すべてを実データ・実再計算で置き換える（V1 → V2 修正の全貌）
上記の5つの課題を洗い出した後、一切の妥協を排してすべての数値を生データと決定論的再計算から再構築しました。
3.1 PDB実験構造の厳密な $R_g$ 再計算（1UAO 全18モデル ＆ 5AWL）
まず、RCSB PDBから取得した 1UAO.pdb（NMR全18モデル）および 5awl.pdb（X線結晶構造 Chain A、水分子等のHETATMを除外）から、原子選択定義ごとの非加重慣性半径（$R_g$）を直接計算しました。
| PDB構造・モデル | 配列 | $\text{C}\alpha$ $R_g$ ($\text{\AA}$) | 重原子（Heavy-Atom）$R_g$ ($\text{\AA}$) | 全原子（All-Atom）$R_g$ ($\text{\AA}$) |
|---|---|---|---|---|
| 1UAO Model 1（野生型代表構造） | GYDPETGTWG |
4.66 (4.6564) |
5.26 (5.2559) |
5.45 (5.4497) |
| 1UAO NMR Ensemble（全18モデル平均±SD） | GYDPETGTWG |
4.99 ± 0.19 (4.65–5.25) |
5.49 ± 0.16 (5.14–5.66) |
5.66 ± 0.15 (5.35–5.83) |
| 5AWL Chain A（変異体 CLN025） | YYDPETGTWY |
4.85 (4.8489) |
6.07 (6.0688) |
6.07（PDB内Hなし） |
| 1L2Y Model 1（Trp-cage 20残基参考） | NLYIQWLKDGGPSSGRPPPS |
6.88 (6.8832) |
7.41 (7.4050) |
7.59 (7.5931) |
この計算により、V1で引用していた「$5.26\text{ \AA}$」という数値は 1UAO Model 1の重原子 $R_g$（5.2559 Å）そのものであり、一方で $\text{C}\alpha$ のみの $R_g$ は $4.66\text{ \AA}$（全18モデル範囲 $4.65\text{--}5.25\text{ \AA}$）、またCLN025（5AWL）は両末端に嵩高いチロシン側鎖（Tyr1, Tyr10）を持つため重原子 $R_g$ が $6.07\text{ \AA}$（$\text{C}\alpha\ R_g = 4.85\text{ \AA}$）になることが明確に整理されました。
3.2 4つのラプラシアン定式化の決定論的再計算と Kabsch $\text{C}\alpha$ RMSD
ハードコードされていた図表データを完全に撤廃し、まったく同じアンフォールド状態の初期ランダムコイル構造（$\text{C}\alpha\ R_g = 8.44\text{ \AA}$、重原子 $R_g = 8.46\text{ \AA}$、$\text{RMSD} = 7.30\text{ \AA}$）から、4つのラプラシアン定式化（① Binary Cutoff、② Continuous $1/d^2$、③ Contact Locking、④ Polar-Priority Directional）を実コードで走らせ、さらに全トラジェクトリに対して Kabsch SVD アライメントによる $\text{C}\alpha$ RMSD を算出しました。
| 対象配列・シミュレーション条件 | 最終 $\text{C}\alpha$ $R_g$ ($\text{\AA}$) | 最終 重原子 $R_g$ ($\text{\AA}$) | 代数的連結度 $\lambda_2$ | 全鎖 $\text{C}\alpha$ RMSD ($\text{\AA}$) | ターン部位（残基3–7）局所 RMSD ($\text{\AA}$) |
|---|---|---|---|---|---|
| 1UAO 初期ランダムコイル（Cycle 0） | 8.44 | 8.46 | 0.029（Binary） | 7.30 | 1.98 |
| ① Binary Cutoff ($1/0$)（10サイクル後） | 8.44 | 8.46 | 0.029 | 7.30 | 1.98 |
| ② Continuous ($1/d^2$, 記憶なし) | 6.29 | 6.72 | 0.094 | 4.29 | 2.13 |
| ③ Contact Locking（10サイクル再計算） | 4.50 | 4.92 | 3.746 | 3.84 | 2.09 |
| ④ Polar-Priority Directional（10サイクル再計算） | 4.73 | 5.19 | 3.562 | 3.50（vs Model 11: 3.22） | 1.92 |
1UAO 保存トラジェクトリ（1uao_folded_directional.json） |
5.15 | 5.55 | 3.106 | 3.97（vs Model 14: 3.00） | 1.15（Asp3偏差: 1.67） |
1UAO 連続勾配降下（gydpetgtwg_differentiable） |
4.66 | 4.98 | 3.611 | 3.56（vs Model 11: 3.25） | 1.91 |
CLN025 (5AWL) Polar-Priority（cln025_folded） |
4.81 | 5.35 | 3.377 | 2.22（最小値: 2.03） | 1.25 |
この実データ再計算から、グラフ理論的フォールディングの本質が鮮やかに浮かび上がりました。
- Binary Cutoff（距離閾値による $1/0$ 隣接行列）が1ミリも動かない理由：閾値（$6.5\text{ \AA}$）より遠い原子間の重みが完全にゼロになるため、二面角 $\Delta\phi, \Delta\psi$ を微小変化させてもラプラシアン行列が変化せず、勾配消失（$\Delta\lambda_2 = 0$） に陥って初期構造（$\text{C}\alpha\ R_g = 8.44\text{ \AA}$）のまま完全にフリーズします。
- Continuous（記憶なし連続重み $1/d^2$）が途中で止まる理由：長距離の勾配は伝わるものの、一度形成された水素結合を保持するヒステリシス（記憶機構）がないため、$\text{C}\alpha\ R_g = 6.29\text{ \AA}$、$\text{RMSD} = 4.29\text{ \AA}$ の緩い局所解で停滞します。
- Contact Locking と Polar-Priority の劇的な凝縮：距離減衰勾配と極性ロッキングを組み合わせることで、代数的連結度 $\lambda_2$ が $0.029 \to 3.562$ へと2桁跳ね上がり、ネイティブ同等のコンパクトさ（$\text{C}\alpha\ R_g = 4.73\text{ \AA}$、重原子 $R_g = 5.19\text{ \AA}$）へと一気に自己組織化します。
3.3 乱数なしの9条件グリッド・決定論的感度分析（Table 2の完全刷新）
課題④であった感度分析についても、np.random.uniform を完全に削除し、共有結合重み $w_{\text{cov}} \in {8.0, 10.0, 12.0}$（$\pm 20%$）× 非共有結合スケール $s \in {0.9, 1.0, 1.1}$（$\pm 10%$）の 全9条件について、10サイクルの二面角最適化を最初から最後まで決定論的に再実行 しました。
| 共有結合重み $w_{\text{cov}}$ | 非共有結合スケール $s$ | 最終 $\text{C}\alpha$ $R_g$ ($\text{\AA}$) | 最終 重原子 $R_g$ ($\text{\AA}$) | 最終 $\lambda_2$ | 全鎖 $\text{C}\alpha$ RMSD vs 1UAO Model 1 ($\text{\AA}$) |
|---|---|---|---|---|---|
| 8.0 (-20%) | 0.9 (-10%) | 4.86 | 5.29 | 3.168 | 3.75 |
| 8.0 (-20%) | 1.0 (Baseline) | 5.27 | 5.66 | 3.249 | 3.92 |
| 8.0 (-20%) | 1.1 (+10%) | 5.01 | 5.39 | 3.635 | 2.94 |
| 10.0 (Baseline) | 0.9 (-10%) | 5.17 | 5.58 | 3.413 | 4.53 |
| 10.0 (Baseline) | 1.0 (Baseline) | 4.73 | 5.19 | 3.562 | 3.50 |
| 10.0 (Baseline) | 1.1 (+10%) | 4.87 | 5.22 | 3.920 | 3.85 |
| 12.0 (+20%) | 0.9 (-10%) | 4.54 | 5.13 | 4.027 | 2.97 |
| 12.0 (+20%) | 1.0 (Baseline) | 4.74 | 5.13 | 4.303 | 2.89 |
| 12.0 (+20%) | 1.1 (+10%) | 4.66 | 5.24 | 4.469 | 3.04 |
| 全9条件 平均 ± SD | — | 4.87 ± 0.25 | 5.31 ± 0.20 | 3.750 ± 0.414 | 3.49 ± 0.61 |
全9条件すべてにおいて、初期の $8.44\text{ \AA}$ から $\text{C}\alpha\ R_g = 4.54\text{--}5.27\text{ \AA}$（平均 $4.87 \pm 0.25\text{ \AA}$）のネイティブ領域へと確実に収束し、特に $w_{\text{cov}}=12.0, s=1.0$ や $w_{\text{cov}}=8.0, s=1.1$ では 全鎖 $\text{C}\alpha\ \text{RMSD} = 2.89\text{--}2.94\text{ \AA}$ に到達することが実証されました。
3.4 コード監査による「カンニングなし（Zero Target Leakage）」の証明
server_overall_fiedler_polar_priority.py を1行ずつ監査した結果、以下の事実がコードレベルで完全に裏付けられました。
-
1UAO.pdbの座標（native_coords等）は参照メトリクス計算およびフロントエンドの3Dビューア比較表示用にのみ読み込まれており、二面角最適化ループ（Lines 672–995）の内部には一切渡されていない。 - 最適化の目的関数（Line 842）は、現在の候補座標から構築した重み付きラプラシアン行列の第2最小固有値 $\lambda_2$、立体反発ペナルティ（$d < 3.2\text{ \AA}$）、および主鎖極性原子対（N–O）の方向性水素結合ボーナス（$< 4.0\text{ \AA}$）の3項のみで構成されている：
$$\text{Score} = \lambda_2 - P_{\text{clash}} + 0.05 \sum_{(i,j) \in \text{polar}} w_{\text{dir}}(i,j)$$
つまり、このシミュレータは実験構造の座標を一切知ることなく、トポロジー（グラフの連結度最大化）と局所的な立体・極性ルールだけでペプチドを折りたたんでいることが100%確認できました。
4. 厳密に検証したからこそ見えた新しい物理的発見
今回のV2改訂で最も面白かったのは、不都合かもしれないデータ（野生型 1UAO の全鎖 RMSD が $3.50\text{--}3.97\text{ \AA}$ であること）を隠さず正確に計算した結果、むしろ分子物理学的に非常に筋の通ったストーリーが浮かび上がったことです。
なぜ野生型 1UAO（GYDPETGTWG）はターンが $1.15\text{ \AA}$ で組めるのに、全鎖 RMSD が $3.5\text{ \AA}$ になるのか？
野生型 1UAO のトラジェクトリ（1uao_folded_directional.json）を詳しく解析すると、中央の5残基ターン領域（Asp3–Pro4–Glu5–Thr6–Gly7）の局所 $\text{C}\alpha\ \text{RMSD}$ は $1.15\text{ \AA}$（Asp3偏差 $1.67\text{ \AA}$）という極めて高い精度で核形成（Nucleation）していました。
それにもかかわらず全鎖 $\text{C}\alpha\ \text{RMSD}$ が $3.50\text{--}3.97\text{ \AA}$ に留まる理由は、両末端の Gly1 と Gly10（側鎖を持たないアキラルで柔軟なグリシン）が閉じきらず、末端間距離が $12.72\text{ \AA}$ 開いたままの「ターン核形成・部分フォールド中間体（Compact, Turn-Nucleated Partially Folded Intermediate）」で止まっているためでした。
変異体 CLN025（5AWL, YYDPETGTWY）が $\text{RMSD} = 2.22\text{ \AA}$ まで劇的に収束する理由
一方、両末端のグリシンを芳香族のチロシン（Tyr1, Tyr10）に置換した超安定変異体 CLN025（5AWL） を同じアルゴリズムで折りたたむと、中央ターンの形成（局所 RMSD $1.25\text{ \AA}$）に続いて、両末端の嵩高いチロシン側鎖同士が疎水性パッキングによって引き寄せられ、ヘアピン全体のジッパーが閉じて全鎖 $\text{C}\alpha\ \text{RMSD} = 2.22\text{ \AA}$（フレーム30で最小 $2.03\text{ \AA}$）のネイティブ $\beta$-ヘアピン構造へと収束したのです。
そもそも実験化学の世界でも、本田ら（2004, 2008）が野生型 Chignolin（1UAO、水溶液中でのフォールド率約60%）の両末端 Gly1/Gly10 を Tyr1/Tyr10 に置換して CLN025（5AWL、フォールド率90%以上・X線結晶構造解析に成功）を設計した理由は、「末端の芳香環パッキングによってヘアピンの末端フaying（ほどけ）を抑え、熱力学的に決定論的なネイティブ構造へロックする」 ためでした。
つまり、私たちのグラフスペクトル・シミュレータが「1UAO では中央ターンは正確にできるが末端のグリシンが揺らいで中間体（RMSD $\approx 3.5\text{ \AA}$）になり、CLN025 では末端チロシンのパッキングが効いてネイティブ構造（RMSD $= 2.22\text{ \AA}$）まで折りたたまれる」という挙動を示したことは、アルゴリズムの欠陥ではなく、実際のペプチド鎖が持つ末端配列依存の熱力学的安定性の違いをグラフラプラシアンが見事に捉えている証拠だったのです。
ごまかさずに全データを再計算し、ネガティブに見える数値とも真っ向から向き合ったからこそ、この発見に辿り着くことができました。
5. 再発防止のために策定した「AI研究データ紐付け鉄則（GEMINI.md）」
今回の教訓を踏まえ、今後のAIエージェントとの共同研究において「どのスクリプトが、どのパラメータで、どの図表の数値を出力したか」が二度と迷子にならないよう、プロジェクト設定ファイル（GEMINI.md）に以下の**データ・コード紐付け鉄則（Provenance Rules）**を組み込みました。
-
実験ごとの「3点セット」自動保存：シミュレーションや解析を実行する際は、①実行スクリプトのコピー、②使用パラメータと乱数シード（
params.json）、③出力生データ（JSON/CSV）を必ずタイムスタンプ付きの同一フォルダにセットで保存する。 -
論文・Markdown表への
[Provenance]コメント埋め込み：原稿内に数値を記載する際は、直前に<!-- [Provenance] Script: recompute_trajectories_and_sensitivity.py | Input: 1UAO.pdb -->のように生成元スクリプトと入力データをHTMLコメントとして強制記録する。 - ハードコード配列・擬似乱数モックアップの禁止：図表生成スクリプト内で数値を直接手打ちすることや、感度分析で乱数近似を使うことをルールとして禁止し、必ず生データファイル（JSON/PDB）を読み込んで描画・集計させる。
- 長大セッションの定期分割と凝縮引き継ぎ：AIとのチャットログが肥大化する前に、専用スクリプトでMarkdownへ抽出・分割し、検証済みの事実だけを凝縮した引き継ぎドキュメントで新セッションを開始する。
6. おわりに：Zenodo V2 と検証コード一式の公開
AIは驚異的なスピードでコードと原稿を生成してくれますが、「その数値が本当に物理計算から出たものか」「実験事実と1対1で対応しているか」を最終的に担保し、責任を持つのは人間の研究者自身です。
今回、開発環境の制限という思わぬインターバルをきっかけに、全データの徹底監査と決定論的再計算を行い、論文の信頼性と物理的洞察を大幅に強化した Version 2（V2） を公開できたことは、研究者として非常に大きな財産になりました。
改訂版論文（PDF）および、誰でも手元のPython環境（NumPy / Matplotlib）で数分以内に全トラジェクトリと図表・感度分析を100%再現できる検証コード・生データ一式（Supplementary Code & Data ZIP）は、ZenodoおよびGitHubにて公開しています。ぜひお手元でコードを動かし、グラフ理論によるタンパク質フォールディングの世界を追体験していただければ幸いです！
- GitHub リポジトリ: blackredfoxcrow7/fiedler-water-folding-sim
- Zenodo Preprint (Version 2): Graph-Spectral Protein Folding: Simulating Peptide Self-Assembly via Laplacian Fiedler Vector Optimization (Revised & Fully Verified Edition)