ベイズ的アプローチによる高速原子間力顕微鏡データと分子シミュレーションの融合

Abstract

データ同化では，ベイズ統計学に基づいて観測データとシミュレーションを融合させることで，観測対象の詳細を推定することができる．このデータ同化手法の1つである粒子フィルタを高速原子間力顕微鏡データと分子シミュレーションに適用するためのプロトコルを開発し，観測不可能な特性が推定できるかどうかを検証した．

1.  はじめに

近年，観測データのデータ不足や解像度不足を補うために高精度な数値シミュレーションが開発され，さらに両者をベイズ統計学に基づいて融合させることによって観測対象の詳細に迫ることが可能となってきている．このような観測データとシミュレーションの融合はデータ同化と呼ばれ，気象学や海洋学をはじめとする多くの分野で活発に適用されている¹⁾．私たちは，観測データとして，日本発の技術である高速原子間力顕微鏡（HS-AFM）で得られるデータに着目した．

HS-AFMでは，生体分子の構造とその動きを一分子レベルで観察することができる．これまでに数多くの生体分子の興味深い動的挙動が明らかとなり，機能発現機構に関する重要な知見が得られている²⁾．しかし，HS-AFM動画の分解能は時間に関して約100 ms，空間については約1 nmであることから，原子レベルの構造や機能発現の動的機構を明らかにすることは難しい．

一方，コンピュータによる全原子分子動力学（MD）シミュレーションは高い時間空間分解能（約1 fsと約1 Å）を持ち，近年の計算機の発展と力場の改良により，様々な実験結果を再現し，生体分子の構造ダイナミクスの詳細を解析できるようになってきている³⁾．しかし，標準的な計算機で到達可能な時間スケールは10 μs程度であり（MD専用計算機Anton 3を使えば昼食までの数時間で終わってしまうが⁴⁾），生体分子の動的挙動を実験結果と直接比較することは未だ困難である．

この時間スケールのギャップを埋めるためには，生体分子の構造を粗視化することが有効である．約10個の非水素原子を1つの粒子に粗視化し，溶媒分子を明示的に扱わない粗視化分子動力学（CG-MD）シミュレーションでは，溶媒分子を明示的に扱う全原子MDシミュレーションに比べて10⁵から10⁷倍程度高速化され，実験結果との直接比較が可能となる⁵⁾．しかし，粗視化モデルでも定量的な再現は容易ではない．

そこで，私たちはデータ同化手法の1つである粒子フィルタを用いたHS-AFM計測とCG-MDシミュレーションとのデータ同化プロトコル（粒子フィルタ・シミュレーション）を開発し，リンカーDNA付きヌクレオソームを対象系とした双子実験によりプロトコルの検証を行った⁶⁾．本稿では，その概要を紹介する．

2.  粒子フィルタによるデータ同化

粒子フィルタは，実験で得られる観測量の時系列y_1:t = {y₁, y₂, ..., y_t}とシミュレーションで得られる状態量の時系列x_1:t = {x₁, x₂, ..., x_t}とを融合するベイズ型のデータ同化手法である¹⁾．本稿で紹介する研究では，生体分子のHS-AFM動画（観測量の時系列）と粗視化モデルのトラジェクトリ（状態量の時系列）とを融合させた．このような融合によって，観測値の補間やシミュレーションの高精度化，観測不能な特性の推定などが可能となる．

時刻tにおける状態量x_tは，観測量の時系列y_1:tの条件付き確率分布p(x_t|y_1:t)として推定することができる．このp(x_t|y_1:t)は事後分布と呼ばれ，ベイズの定理を用いると，

  

$$p\left(x_t|y_{1:t}\right)~ p\left(y_t|x_{t,}{ y}_{1:t-1}\right)p\left(x_t|y_{1:t-1}\right)$$

と変形することができ，尤度p(y_t|x_t, y_1:t–1)と事前分布p(x_t|y_1:t–1)の積に比例することがわかる．この式による事後分布の導出を「フィルタリング」と呼ぶ．一方，事前分布p(x_t|y_1:t–1)は，時刻t – 1の事後分布p(x_t–1|y_1:t–1)のシミュレーションによる時間発展によって算出することができる；

  

$$p\left(x_t|y_{1:t-1}\right)=\int p\left(x_t|x_{t-1}\right)p\left(x_{t-1}|y_{1:t-1}\right)d{x}_{t-1}.$$

ここで，x_tはマルコフ性を持つことを仮定した．この事前分布の導出を「予測」と呼ぶ．この「予測」と「フィルタリング」のステップを1ラウンドとし，これを繰り返すことで，p(x_t|y_1:t)の時系列を得ることができる．粒子フィルタでは，多数のサンプル（「粒子」と呼ぶ）からなるアンサンブルによって確率分布を近似する．「予測」ステップでは各粒子について一定時間のシミュレーションを実行し，「フィルタリング」ステップでは観測量との比較によって推定される各粒子の尤度に基づいて粒子の再抽出を行う．

私たちはこの粒子フィルタを生体分子のHS-AFM計測データに適用するために，CG-MDシミュレーションを用いたプロトコル，粒子フィルタ・シミュレーション，を開発した（図1）⁶⁾．ここでは，実験の観測量は生体分子のAFM像であり，シミュレーションの状態量は生体分子の粗視化モデルとなる．予測ステップでは，粗視化生体分子シミュレータCafeMol⁷⁾を用いてCG-MDシミュレーションを実行した．一方，フィルタリングステップでは，粗視化モデルからafmize⁸⁾を用いて作成した疑似的なAFM像を観測量であるAFM像と比較して尤度を推定した．方法の詳細は文献6を参照してほしい．

図1

HS-AFM計測データとCG-MDシミュレーションを融合した粒子フィルタ・シミュレーションにおける1ラウンドのプロトコル．文献6より改変．

3.  粒子フィルタ・シミュレーションの検証

開発したプロトコルを検証するために，リンカーDNA付きヌクレオソームを対象系とし，10⁶ MDステップ間隔で疑似的に作成したAFM像の時系列を観測データとする双子実験を実施した．まず，粒子フィルタ・シミュレーションにおけるサンプルの数（粒子数）を512とし，1ラウンドを観測データと同じ10⁶ MDステップにして合計10ラウンドの計算を実行した．各ラウンドの予測ステップで得られたアンサンブルはリンカーDNAが多様な配置を取っており，その影響でフィルタリングステップにおける尤度は幅広い分布を示した（図2左）．特に，ラウンド1と10では他のラウンドに比べてより幅広い分布となっていたが，これは観測データにおけるリンカーDNAがより高い位置にあり，サンプリングが難しいことに起因していると考えられる．幅広い尤度分布の結果，再抽出される粒子が1つだけとなる退化の問題が発生したものの，最大尤度の粒子のAFM像は観測データからのズレが小さく，よく再現できていた（図2右）．

図2

粒子数512の粒子フィルタ・シミュレーションの結果．各ラウンドの予測ステップで得られた対数尤度分布（左；赤矢印は最大尤度の位置）と最大尤度の粒子と観測データのAFM像（右）．文献6より改変．

続いて，粒子フィルタ・シミュレーションの粒子数依存性，特に退化の問題が解決されるか，を調べるため，粒子数を8，32，128，512，2048，8192と変化させてシミュレーションを実行した．尤度最大のAFM像は，いずれの粒子数においても，観測データとよく似たものとなっていた．さらに，最大尤度の時系列について10ラウンドの累積尤度を比較したところ，粒子数の増加にしたがって着実に尤度が大きくなるという期待通りの結果が得られた（図3a）一方，各ラウンドの尤度を比較したところ，尤度と粒子数の関係が逆転するケースが確認された（図3b）．この結果は，異なる粒子数の粒子フィルタ・シミュレーションの結果を比較する場合，1ラウンドのデータだけよりも複数ラウンドのデータを用いる方がより適切であることを意味する．また，粒子数を8192まで増やしても退化の問題は解決されなかった．

図3

異なる粒子数の粒子フィルタ・シミュレーションの結果．最大尤度の時系列のa）10ラウンドの累積対数尤度とb）各ラウンドの寄与の内訳．文献6より改変．

4.  実験で観測不能な特性のベイズ推定

データ同化の目的の1つとして，実験では観測できない隠れた特性の推定がある．HS-AFM計測の場合，試料周辺の局所的なイオン濃度や試料ステージ表面・探針との相互作用などが生体分子の挙動に大きな影響を与えるが，その影響を観測データから評価することは一般に困難である．しかし，観測データに基づいて確率的な推定を行うベイズ推定の枠組みを用いると，観測不能な特性を推定することが可能である．そこで，粒子フィルタ・シミュレーションによって，このような推定が可能であるかどうかの調査・検証を行った．

まず，イオン濃度に着目し，CG-MDシミュレーションのイオン濃度を疑似的に作成した実験データの場合（0.2 M）の2倍（0.4 M）と1/2倍（0.1 M）にして，粒子フィルタ・シミュレーションを実行した．解析対象としたリンカーDNA付きヌクレオソームは電荷に富んでおり，その挙動はイオン濃度の違いによって大きく影響を受けると予想される．シミュレーションで得られた最大尤度の時系列について10ラウンドの累積尤度を比較したところ，異なるイオン濃度において明瞭な差が観察された（図4a）．疑似的実験データと同じ場合（0.2 M）が最大となり，イオン濃度が正しく推定できることが確認された．

図4

粒子フィルタ・シミュレーションによる観測不能な特性の推定．a）異なるイオン濃度とb）異なる時間スケールのシミュレーションで得られた最大尤度の時系列における10ラウンドの累積対数尤度．文献6より改変．

続いて，実験と粗視化シミュレーションの時間スケールの対応付けが可能かどうかを検証するために，1ラウンドのシミュレーションのMDステップ数を10分の1の10⁵とした粒子フィルタ・シミュレーションを実行した．尤度最大の時系列について10ラウンドの累積尤度の比較を行ったところ，正しい時間スケールの対応付けの場合が有意に大きくなり（図4b），時間スケールの対応付けが可能であることが示唆された．

5.  おわりに

粒子フィルタにおける退化は一般には望ましいことではないが，AFM像が大次元データであり，解析対象の挙動が本質的に確率的であることを踏まえると，やむを得ないことと考えられる．一方で，退化があっても隠れた特性のベイズ推定は可能であることが示され，HS-AFM動画と完全に整合する高い時間空間分解能の構造データを得ることには成功していることから，退化の問題は深刻な問題ではないと考えている．

本稿で紹介した研究では1つの構造モデルから作成した疑似的なAFM像を用いたが，実際のHS-AFM計測で得られるAFM像ではピクセルごとに観測時間が異なる（非斉時性）．そこでこの問題に対処すべく，非斉時性を考慮した粒子フィルタ・シミュレーションへの拡張を行い，双子実験による検証でその有効性を確認済みである．今後は，実際のHS-AFM計測データへの適用を進めていきたいと考えている．

文献

• 1)  樋口知之ら (2011) データ同化入門 ―次世代のシミュレーション技術―，朝倉書店，東京．
• 2)   Ando,  T. (2019) Curr. Opin. Chem. Biol. 51, 105-112. DOI: 10.1016/j.cbpa.2019.05.010.
• 3)   Bottaro,  S.,  Lindorff-Larsen,  K. (2018) Science 361, 355-360. DOI: 10.1126/science.aat4010.
• 4)   Shaw,  D. E. et al. (2021) Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC '21), 1, pp. 1-11. DOI: 10.1145/3458817.3487397.
• 5)   高田 彰二（2021）生物物理 61, 144-151. DOI: 10.2142/biophys. 61.144.
• 6)   Fuchigami,  S. et al. (2020) J. Chem. Theory Comput. 16, 6609-6619. DOI: 10.1021/acs.jctc.0c00234.
• 7)   Kenzaki,  H. et al. (2011) J. Chem. Theory Comput. 7, 1979-1989. DOI: 10.1021/ct2001045.
• 8)   Niina,  T. et al. (2020) J. Chem. Theory Comput. 16, 1349-1358. DOI: 10.1021/acs.jctc.9b00991.

Biographies

渕上壮太郎（ふちがみ　そうたろう）

京都大学大学院理学研究科特定准教授