Ensemble Modeling of Biomolecules by Small-angle X-ray Scattering and Coarse-grained Molecular Dynamics Simulation

Abstract

小角散乱プロファイルは液中に存在するあらゆる分子の構造情報を反映するため，分子運動や揺らぎの包括的理解に寄与する．本稿ではマルチドメインタンパク質を対象に，X線小角散乱測定と分子動力学計算による分子運動解析手法の現状と課題，さらに他の実験との統合的解析を目指した新戦略について解説する．

1.  はじめに

人類が明らかにした生体分子の立体構造は現在も加速度的に増加し続けている．こうした高解像度立体構造研究の進展は，同時に時間軸も加えた生体分子ダイナミクス研究への興味を引く原動力ともなっている．分子モーターを筆頭に生体分子の運動は機能と密接に関わっており，核磁気共鳴法（NMR），Fluorescence Resonance Energy Transfer（FRET），高速原子間力顕微鏡（HS-AFM），小角散乱法（SAS）等，様々な手法によって研究が進められている．動的なタンパク質の多くはマルチドメイン構造を有しており，ドメインの連動機構を可視化することが構造機能相関の理解に重要である．本記事で紹介するX線および中性子小角散乱法（SAXS, SANS）は，同一の分子種であっても分子構造の違いが散乱プロファイルに反映されるため，既知構造と異なるコンフォーメーションの存在を鋭敏に検知する．従って適切に利用すれば水溶液中における分子運動の解明に強力なツールとなる．一方，小角散乱プロファイルの解釈は容易でないが，その方法は多くの研究者によって模索されている．そこで，様々な実験手法の連携による高度な分子構造解析の発展を見据え，小角散乱プロファイル解釈の現状と，マルチドメインタンパク質の分子運動理解に向けた我々の試みを紹介する¹⁾．

2.  小角散乱プロファイルによるアンサンブル解析

小角散乱法はその名の通り，試料溶液にビーム（X線，中性子）を入射し，その散乱を記録する研究手法である．大きな利点は生体分子に何も束縛を与えずに水溶液中で様態観測が可能な点にある．得られるのは直接的な原子座標ではなく，散乱ベクトルQと呼ばれる量についての散乱強度関数である．そこには1 nmから100 nm程度の範囲の分子構造情報が含まれている．小角散乱法の原理の詳細は文献2を参照されたい．分子の原子配置から小角散乱プロファイルを簡便に理論計算することが可能であるため³⁾，例えば既知結晶構造の理論散乱曲線と実測の散乱データとの比較により，既知構造が溶液中の主要構造か否かを容易に判別可能である．さらには，小角散乱プロファイルから分子形状をモデリングするab initio法も確立している⁴⁾．

結晶構造は，あくまで結晶場によって選択された一つの準安定構造を解析によって可視化したもので，溶液中では本来分子は動的に揺らいでいるため，構造は一つではなく，多数の構造状態が混在する状態（構造アンサンブル）になっているはずである．しかしながら，SASプロファイルは以下の式（1）で表わされるように，溶液中に存在するすべての分子，つまり構造アンサンブルからの平均散乱強度になっている．

  

$$I\left(Q\right)={\sum }_j^{\text{全分子構造}}{w}_j{\times I}_j\left(Q\right)$$（1）

なお式（1）では，生体分子が溶液中でとりうるすべての構造に通し番号を振ったときの，j番目の構造の体積分率をw_j，散乱強度をI_j(Q)とした．本来の構造機能相関解析のためにはそのアンサンブルの実像を明らかにしたいのだが，式（1）から想像されるように，アンサンブルを構成する構造一つ一つのSASプロファイルと体積分率w_j × I_j(Q)を，実験で得られたSASプロファイルI(Q)だけから直接的に同定することは不可能であった．しかし近年，SASプロファイルを活用して溶液中の構造アンサンブルをモデリングする手法やソフトウェアがいくつか公開され，それらを活用したチャレンジングな解析が進められている．そのような手法はアプローチの仕方によって2つに分類することができ，それぞれmaximum parsimony（MP），maximum entropy（ME）と呼称される⁵⁾．前者のMP法では，可能な限り少数の分子構造で小角散乱プロファイルを説明する立場をとる．MP法では離散的な少数の分子構造が解として得られるため，構造変化を伴うタンパク質のうちでも，複数の準安定構造（例えばopen，closeの2状態）の存在がある程度予想されているものに対して特に有効である．なお日本生物物理学会年会等で比較的見かけるensemble optimization method（EOM）法⁶⁾もこの戦略に基づく．後者のME法は，分子動力学（MD）計算等の理論計算から求まる構造出現分布（自由エネルギー地形）に対し，可能な限り微小な変形によって実験データを再現させる考え方である．「可能な限り微小」の程度は，理論計算と実験データの乖離の程度に依存する．またME法は，自由エネルギー地形を事前分布とみなし，小角散乱プロファイルを観測として取り込むことにより事後分布を求めるとの解釈が可能であり，これはベイズ統計の考え方そのものである．そのためME法ではベイズの式を用いたアプローチや，ベイズ統計を活用する手法が提案されている．ME法では，①予め生体分子の構造を十分に探索しておき，後から構造分布を変形する方法，あるいは②実験データとの整合度をポテンシャル関数の形で力場に組み込み，それを用いて複数のMD計算を並列で実施することにより達成する．どちらを採用しても，ME法では連続的かつ多様に構造変化する天然変性領域の解析が可能であり，さらに離散的な自由エネルギー地形に対してMP法と類似の解を出力できる意味では，MP法を包含している．統計モデリング技術の発展と並行し，ME法が今後さらに発展するだろう．

アンサンブルモデリング技術の発展の一方で，依然として熟考すべき難点も存在する．小角散乱プロファイルが有する溶液構造情報は他の実験では得られないものであるが，その情報量は少ない点である．これは，アンサンブルモデリングの結果が理論モデリングのパラメータ（あるいは力場）に大きく依存することを意味しており，例えばME法を利用する場合には適切なMD力場選択が重要となる．この点は，膨大なデータを取り込むことで最適解に近づく，いわゆる「ビッグデータ解析」とは対照的な状況となっている．

3.  小角散乱プロファイルの情報量を考慮したアンサンブルモデリング

以上の背景を踏まえつつ，我々は新たな構造アンサンブルモデリングを考案した．上述した既存のME法は単一の構造アンサンブルを同定しているが，MD計算では特定の力場の下で行われた結果になっている．我々は，単一の構造アンサンブルではなく，複数の力場を使用して一定の論理的整合性のある構造アンサンブル候補を可能な限り多数提示するほうが有効ではないか？と考えた．原理的には①力場パラメータを少しずつ変えながら各力場で構造サンプリングを実施し，②その結果からそれぞれ少しずつ異なる自由エネルギー地形を得た後に，③各々に対しME法によるモデリングを実施すれば実現可能である．ただし膨大かつ非現実的な計算コストが予想されるため，それに代わる簡便な近似手法を検討した．我々は，「平衡状態を充分に探索し尽くしてはいないが，ナノスケールの分子構造変化が繰り返し観測される程度の長時間MD計算」から得られる分子構造分布は，平衡状態における自由エネルギー地形とは異なりつつも，その特徴を保持していることに着目した．そしてその分子構造分布を，「本来の力場からパラメータを少し変更した力場より得られる自由エネルギー地形」の一種と考えることにした．このアイデアに基づけば，異なる乱数を用いながら複数回の等温MD計算を実施して得られる構造分布群の各々は，それぞれ異なる力場でのMD計算から得られた自由エネルギー地形とみなせる．また，MD計算データに対するME法の適用は，初めから実験データに合致するMD計算の時系列を用意することで代用可能である．以上を踏まえ，最終的に我々は，①スケーリングの余地がある力場パラメータの値を様々に変更しつつ多数の等温MD計算を実施し，②各等温MD計算の時系列から，実測の小角散乱プロファイルを再現し，かつ充分長い領域を切り出し，目的の構造アンサンブルとする戦略に至った．実用上，MD計算による構造探索効率が高いほうが多くの構造アンサンブル候補を提示できることになるため，粗視化（CG）モデルによるMD計算（CGMD計算）を採用した．CGモデルは原子ではなく原子団，例えばアミノ酸残基，あるいはアミノ酸残基に含まれる官能基を一つの質点として扱うモデルであり，計算時間を短縮できる．

4.  マルチドメインタンパク質ER-60の分子運動解析

考案法をマルチドメインタンパク質に適用した場合に，その分子運動についてどれだけの情報が得られるかを検証した．我々は構造アンサンブル解析のモデル系として，マルチドメインタンパク質であるER-60を対象とした．ER-60はprotein disulfide isomerase familyに属し，ジスルフィド結合の架け換えによってタンパク質の折れたたみを助ける機能を持つ⁷⁾．構造の観点からは，4つの球状ドメイン（a, b, bʹ, aʹ）から構成される（図1A）．ER-60のSAXSデータは我々のグループが過去に測定したものを用いた⁸⁾．本研究ではMartini力場⁹⁾によるCGMD計算を実施し，続き考案法によるアンサンブルモデリングを行った．CGMD力場においては，アミノ酸－水分子間の相互作用とアミノ酸－アミノ酸間の相互作用のバランスに検討の余地があるため，このバランスに関するパラメータを走査しながら多数の等温CGMD計算を実施した．各MD計算データからの時系列切り出しでは，実測のSAXSプロファイルを再現し，かつ1 μs以上であることを判断基準とした．なお計算から得られる原子座標に基づく理論小角散乱プロファイルと，実測のプロファイルとの合致度は式（2）で定義されるχ²値を用いて評価した．

図1

考案法によるER-60のアンサンブルモデリング．A．マルチドメインタンパク質ER-60のドメイン構成．ER-60はa（青），b（緑），bʹ（黄），aʹ（赤）の4ドメインから構成される．B．CGMD計算中におけるER-60の構造変化に伴うχ²値の変化の一例．考案法でSAXSプロファイルに合致する領域として切り出した区間の背景を水色で着色した．C．切り出した構造アンサンブルの理論SAXSプロファイル（水色実線）と実測プロファイル（○）．エラーバーは標準偏差を示す．

  

$${\chi }^2=\frac{1}{N-1}\sum _{j=1}^N{\left(\frac{I_{exp}\left(Q_j\right)-I_{mdl}\left(Q_j\right)}{{\sigma }_j}\right)}^2$$（2）

ここでNはデータ点の数，I_exp(Q_j)，I_mdl(Q_j)はそれぞれQ_jにおける実測散乱強度あるいはモデルからの理論散乱強度，σ_jは実測値の標準偏差である．両プロファイルの合致度が高い場合にχ²値は小さく，経験的に約3以下であれば合致しているとみなせるため，本研究でも基準に用いた．CGMD計算の時系列データに対して，各時点におけるスナップショット構造からの理論散乱プロファイルと実測プロファイルとの一致度を比較したχ²値を図1Bにプロットしている．χ²値は3.0から303.0という広い範囲に分布し，χ² < 5.0となることは稀であった．しかしながら，背景を水色で示した領域におけるすべての理論プロファイルを平均すると，その平均プロファイルは実測プロファイルと一致した（図1C）．このように，考案法でSAXSプロファイルを再現する構造アンサンブルが得られた．併せて，各々では実測のSAXSプロファイルを全く再現しない多数の分子構造によって構造アンサンブルが構成されうることが改めて示された．SAXSプロファイルは全分子の平均散乱強度であるため，このような結果であっても問題はない．

我々はこのアプローチを多数のMD計算データに適用することにより，SAXSプロファイルを再現する複数の構造アンサンブルを得たが，ここではそのうち3つについて概観する．それぞれについてaドメインとaʹドメインの幾何重心間距離の分布を図2に示した．構造アンサンブル#1では，この距離が40 Åと50 Åに存在する確率が高い傾向があるものの，全体としてなだらかな分布となっている．一方#2では，40 Åと60 Åにそれぞれピークがあり，aドメインとaʹドメインの距離に関して，2タイプの構造の存在を示唆する．さらに#3では50 Å付近で出現頻度が最大だが，#1と比べて70 Å以上となる確率は非常に小さい．なお，構造アンサンブル中の全構造の理論散乱プロファイルを計算し，その平均散乱強度を実験値と比較した場合のχ²値はそれぞれ2.19，2.32，2.03となっており，いずれも実験データを良く再現する．つまり，「SAXSプロファイルとCGMD計算に基づき，論理的整合性のある構造アンサンブル候補を複数提示する」という目的は達成した．逆に，考案法ではこれ以上の絞り込みの手がかりをSAXSプロファイルとCGMD力場に求めない．すべての構造アンサンブルの共通の特徴としては，aおよびaʹドメインは大きく揺らぐことが示唆された．

図2

考案法により得られたアンサンブルの内3つについて，aドメインとaʹドメインの幾何重心間距離の分布を示した．

最後に，単一の小角散乱プロファイルのみではER-60の構造アンサンブルを同定できない原因を検討した．図3Aには構造アンサンブル#3に含まれる10構造の理論SAXSプロファイルをそれぞれ示しているが，興味深いことにQ = 0.12 Å^–1近傍の2か所に等散乱点が存在していた．分光学においては，「等吸収点が存在する場合，吸光度と線形関係の量が存在する」¹⁰⁾という近似法則が適用できるのだが，SAXSにおける等散乱点も同様に「ER-60の分子構造に散乱強度I(Q)と線形関係となる構造量が存在するのではないか」と推測された．実際，aドメインの幾何重心とaʹドメインの幾何重心の距離とI(Q)が近似的な線形関係にあることが見出された．例として図3BにはQ = 0.07 Å^–1におけるI(Q)とaドメイン－aʹドメイン間距離の関係を示す．こうした関係が成立した場合には話は非常に単純になり，分子構造毎の散乱プロファイルの違いは，aドメイン－aʹドメイン間距離の違いと解釈できる．Q = 0.07 Å^–1における散乱強度の実測値は6.4 × 10^–3 cm² mg^–1であり，図3Bによればそれには約50 Åのドメイン間距離が対応する．つまり，アンサンブル解析結果の要約は，「溶液中ではaドメインとaʹドメインの重心間距離の平均はおよそ50 Åであり，かつその妥当なドメイン間距離分布のいくつかをCGMD計算が提示した．具体的には，40-50 Åあたりにピークを持つ単峰の分布や，40 Åおよび60 Åにそれぞれピークを持つ二峰の分布等である．」となる．ただしMD計算を組み合わせているため，ER-60の構造アンサンブルは単なるドメイン間距離の情報のみならず，それを実現している各ドメインの配向や相互作用面の情報をも有している．単なるドメイン間距離情報以上の情報を持つ点に関しては，線形近似が成立しないQ > 0.2 Å^–1領域においても構造アンサンブルが実験結果を再現することからも支持される．

図3

構造アンサンブルに含まれる各構造とI(Q)との関係．A．構造アンサンブル#3に含まれる分子構造のうち10構造の散乱プロファイルを重ね書きした．交点が集まる箇所を矢印で示した．B．aドメインとaʹドメインの幾何重心間距離とQ = 0.07 Å^–1におけるI(Q)との関係を，出現確率のヒートマップにより示した．

5.  考案法の今後の展望

今回考案した手法は，MD計算力場のパラメータが必ずしも精確ではない状況における，実験データと理論計算の統合解析の選択肢の一つである．「実際の分子運動がどのようであるか？」を特定するのではなく，「次にどのような実験を行えばさらに理解が進むか？」への手がかりを得ることを主目的としている．ER-60に関しても，得られた多数のアンサンブル候補の中でどれが最も正解に近いかを調べるためのSANS実験，FRET実験等を，本記事で紹介したデータに基づいて設計可能である．念のため，考案法で得られた構造アンサンブルのあらゆる線形結合も実験データを再現するために分子運動の候補となることを補足しておく．研究の過程で思いもかけず現れる新発見の前では，次にすべき実験がすぐには思いつかないことも多い．大抵の場合に複数のグループによる共同研究を行う現代においては，様々な手法を連携させるための戦略の確立も重要である．本記事の考案手法がその一助となり，タンパク質ダイナミクスの解明に役立てば幸いである．

文献

• 1)   Shimizu,  M. et al. (2022) Sci. Rep. 12, 9970. DOI: 10.1038/s41598-022-13982-9.
• 2)   Bernado,  P. et al. (2018) Biochim. Biophys. Acta Gen. Subj. 1862, 253-274. DOI: 10.1016/j.bbagen.2017.10.015.
• 3)   Grudinin,  S. et al. (2017) Acta Crystallogr. D Struct. Biol. 73, 449-464. DOI: 10.1107/S2059798317005745.
• 4)   Svergun,  D. I. (1999) Biophys. J. 76, 2879-2886. DOI: 10.1016/S0006-3495(99)77443-6.
• 5)   Bonomi,  M. et al. (2017) Curr. Opin. Struct. Biol. 42, 106-116. DOI: 10.1016/j.sbi.2016.12.004.
• 6)   Bernadó,  P. et al. (2007) J. Am. Chem. Soc. 129, 5656-5664. DOI: 10.1021/ja069124n.
• 7)   Kozlov,  G. et al. (2010) FEBS J. 277, 3924-3936. DOI: 10.1111/j.1742-4658.2010.07793.x.
• 8)   Okuda,  A. et al. (2021) Sci. Rep. 11, 5655. DOI: 10.1038/s41598-021-85219-0.
• 9)   Souza,  P. C. T. et al. (2021) Nat. Methods 18, 382-388. DOI: 10.1038/s41592-021-01098-3.
• 10)   Greger,  M. et al. (2013) Phys. Rev. B 87, 195140. DOI: 10.1103/PhysRevB.87.195140.

Biographies

清水将裕（しみず　まさひろ）

計算モデリング担当

杉山正明（すぎやま　まさあき）

京都大学複合原子力科学研究所粒子線物性学分野教授