2024 年 23 巻 4 号 p. 105-114
Recent efforts have focused on utilizing molecular interaction data from FMO calculations for phase separation simulations in materials design. Accurate prediction of phase separation, which is closely related to molecular affinity, has long been a challenge due to difficulties in calculating accurate interaction parameters. We have developed a framework for estimating effective interaction parameters between coarse-grained components using FMO calculations. In addition, a simulation scheme called FMO-DPD, which applies these parameters to dissipative particle dynamics (DPD) simulations, has demonstrated its effectiveness in various systems. These developments are discussed in this paper.
Recent efforts have focused on utilizing molecular interaction data from FMO calculations for phase separation simulations in materials design. Accurate prediction of phase separation, which is closely related to molecular affinity, has long been a challenge due to difficulties in calculating accurate interaction parameters. We have developed a framework for estimating effective interaction parameters between coarse-grained components using FMO calculations. In addition, a simulation scheme called FMO-DPD, which applies these parameters to dissipative particle dynamics (DPD) simulations, has demonstrated its effectiveness in various systems. These developments are discussed in this paper.
巨大分子系に対し非経験的な量子化学計算が可能なフラグメント分子軌道 (FMO) 法 [1,2,3]は,タンパク質・核酸・分子凝集系を主な対象に,分子間相互作用を中心としたナノレベルの詳細解析に用いられてきたが,近年,FMOから得られる分子間相互作用情報を,材料設計で重要となる相分離シミュレーションに活用する取り組みを推進している.材料機能の最適化において,メソスケールの分子集団の相分離構造の制御が重要な要素のひとつであり,相分離には分子間の親和性が密接に関わっている.分子計算による相分離の予測には,原子団を粗視化したシミュレーションが用いられるが,成分間の親和性を示す有効相互作用パラメータの高精度な算定は難しく,長い間大きな課題となっている.そこで私たちのグループでは,FMO計算に基づき粗視化成分間の有効相互作用パラメータを算定するフレームワークを確立した.更に,FMOから得られた相互作用パラメータで,相分離の予測に優れた散逸粒子動力学 (DPD) シミュレーション [4]を実施する一連のスキームFMO-DPDと称し,種々の系で先導的な実証応用計算を進めてきた [5,6,7,8,9,10,11,12,13,14,15,16,17,18].本稿では,これらの試みを概説するとともに,HPC資源を活用した産学の連携や,更なる高度化,機械学習を活用した高速化 [19]など,最近の取り組みについても言及する.
粗視化シミュレーションは,分子動力学 (MD) などの原子レベルのシミュレーションと比較して,大規模で長い時間スケールを処理することができるため,大規模構造の解析に広く用いられる.粗視化手法の中でも,一般的なMDの関数系を保ったまま原子団を粗視化する手法と比べ,系内の成分の濃度の時間発展を計算する動的平均場理論 [20, 21]や,流体の理論を取り込んだDPD [4, 22, 23]は,相分離構造の予測に優れている.中でもDPDは粒子の扱いを保った上で,流体効果を含んだ相挙動を低コストで追うことができるために,電解質膜 [24,25,26,27,28,29],界面活性剤 [30,31,32,33,34,35],液晶 [36],ドラッグデリバリーシステム (DDS) [37] など,化学構造に基づいた応用例が多数存在する.
平均場法やDPDシミュレーションでは,設定した粒子間の相互作用を記述するパラメータが重要となり,これらはFlory-Huggins理論の成分間の親和性を示すχパラメータに密接に関連している [38].Flory-Hugginsの格子理論において,2成分混合系の自由エネルギー変化ΔGは式(1)で示される.
(1) |
ここで,φは体積分率,Pi (i = 1,2) は各成分 (セグメント)の鎖長で,ポリマー分子量とセグメント分子量の商で表されるセグメントの連結数であり,Rは気体定数,Tは温度である.2成分間のχパラメータは次のように定義される.
(2) |
(3) |
化学構造に基づいたχの予測は古くからおこなわれているものの,信頼できるχの評価は困難であることがよく知られ,シミュレーションには実験値を元にした値や経験的な値が用いられることも多い.そのために,汎用性の問題が指摘されている.化学構造からのχの予測について代表的なものは2種類あり,1つ目は溶解度パラメータSP値からの予測である.SP値はHildebrand [39]によって考案され,原子団寄与モデル [40, 41],bicerano [42]の手法など様々な経験的な推算モデルが存在する.MDシミュレーションを用い,単一分子種のバルクの計算から凝集エネルギーを求める方法もある [43].2つ目は,異種分子間の相互作用からの算出である.主にFanらによるセグメント間の接触エネルギーからの予測 [44]と,混合系のMDにより求める方法 [45]がある.また,直接的に異種分子間の相互作用を考慮している訳ではないが,量子化学計算で得られた単一分子の表面電荷に基づいて算出するCOSMO法 [46]が挙げられる.FMO法との連携は,上記Fanらの手法を踏襲し構築している.この手法はFlory-Hugginsの定式化に近く,接触する粒子間の相互作用を直接分子論的に求めることができる利点がある.Fanらの手法と,FMO法との連携の概要を以下に簡単に示す.詳細は文献 [5, 6]をご参照されたい.
Fanらの手法は,セグメントにあたる小分子の全原子モデルを準備し,各々を剛体とした2分子対の網羅的な配置について相互作用計算を行うことで
(4) |
ここで,
構築した手法におけるFanらの手法からの改良として,主に2点が挙げられる.まず,Fanらの手法は接触エネルギーを固定電荷,あるいは電荷を考慮しない古典力場により評価しているために,電荷移動が相互作用の本質である系では適当な平均的固定電荷を決定する上で信頼性が低く,また電子の非局在化の寄与を考慮できない問題があるため,FMO法を用いて
2点目は,粗視化成分として扱うセグメントの異方性の考慮である.前述のように,網羅的な孤立2分子対の配置から,Metropolis法で採用した配置の平均エネルギーを
(5) |
(6) |
ここで,
(7) |
上記を経て,最終的に,以下の2通りでχが算出される.
(8) |
(9) |
上記ワークフローの検証計算として,典型的な系の相転移上部臨界温度の評価を行った.
(10) |
本研究の手法により得られた温度ごとの
開発した上記ワークフローは,任意の分子間のχパラメータを算定できるシステムとして整備し,FMO-based Chi parameter Evaluation Workflow System (FCEWS) という名称で2018年より公開している[6].プログラムはPython (3系)とFortran 90からなっており,ABINIT-MP [47, 48]が使用可能であることが前提となっている.ユーザーはパラメータを作成したいセグメントのxyzファイルを作成するのみで,構造の作成,ABINIT-MPによる計算,χパラメータの算出を一括で行うことができる.計算コストの観点から,Linuxでの運用を想定しているが,Windowsでの使用も可能である.各種ジョブ管理システム(PBS, Torque, LSF, Lava)のほか,FOCUSや「富岳」といったスーパーコンピュータへのジョブ投入にも対応している.FCEWSはフリーで利用可能であり,計算工学ナビにて公開されているほか,JSOLが開発するマルチスケールシミュレーションソフトウェアJ-OCTA [49]において,一部機能を除いてGUIを用いて使用することが可能である.2024年現在のFCEWS最新版はVer. 2.10で,動作の安定性や,各種計算機の対応が改良されている.
前述のχパラメータ算定システムFCEWSを基に,粗視化セグメント間の相互作用パラメータをFMO法で決定し,DPDを実施する一連のスキームをFMO-DPDと称し,種々の応用計算を進めている.本章では,応用例の概要を紹介する.紙面の都合で概要のみにとどめるが,詳細は各項目に記載する文献を参照されたし.セットアップや解析の多くはJ-OCTAソフトウェア
[49]で行い,DPDの実行はOCTAのCOGNACモジュール [50]で実施した.なおDPD内では粒子間の反発パラメータである
(11) |
高分子電解質膜は環境負荷の低いエネルギーとして期待される燃料電池のイオン交換膜として使われている.内部のプロトン伝導が電池性能に大きく影響するため,水和構造のつくる水ネットワークが重要とされているが,電子レベルの相互作用が構造形成の鍵となるため,パラメータ算定は大きな課題となっていた.本研究ではフッ素系骨格を有し燃料電池膜に一般的に使われるNafion ® と 芳香族炭化水素骨格を持つSPEEK膜に対して,Figure 1のように分子内を小部位に分けてパラメータを算定し(MP2/6-31G(d')レベル) ,複数の等価質量のモデルで,含水量を変えたFMO-DPDシミュレーションを行った [7].Nafionについて,水クラスタの直径とクラスタ間距離を水粒子の動径分布関数から算出したところ,実験値 (それぞれ約4 nm, 5 nm) と良い一致を示し,含水率により増減する傾向を再現した.更に,プロトン伝導に直結する指標として,水粒子の連結評価によるpercolation解析を実施した.実験だと,Nafionでは含水量が10 vol%ほどでプロトン伝導が生じるのに対し,SPEEKでは30 vol%程度まで生じない旨が報告されているが,DPDシミュレーション結果を用いた連結評価おいても,両分子種の差について近しい値で傾向を再現することができた.
Segment structures of Nafion and the results of DPD simulations. Adapted from Ref. [7] under the Creative Commons Attribution 3.0 Unported (CC BY 3.0) license.
リン脂質は生体内で二重膜構造をとり,膜に嵌入したタンパク質が分子を認識し透過することで複雑な機能を保っている.応用の一例として,生体膜構造を模したセンサーなどのナノバイオデバイスのほか,DDSのキャリアとしても研究がなされている.リン脂質電荷をもつ頭部と,尾部の疎水部の微細構造が種々の相分離を引き起こすため,FMOによるパラメータ算定と相性がよい.まず代表的なリン脂質であるPOPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) を対象に構造再現を試みたところ,エネルギー計算に溶媒効果を付与することで適切なχパラメータが得られ,脂質濃度に応じてベシクルや二重膜構造が得られた.構造検証として,二重膜状態の構造において,単一脂質分子が占める面積である膜面積と,垂直方向の膜厚を実験値と比較した.二重膜平面に対して,平行方向と垂直方向の圧力が釣り合った体積分率の分子数から膜面積を算出したところ69.4 Å2となり,実験値 62-68 Å2と良好な一致を示した.また,膜に垂直な方向の一次元濃度分布から膜厚を算出したところ2.8 nmであり,H-NMRの実測値2.58 nmと良い一致を示した [8].最近では多価不飽和脂質とコレステロール (Chol) 混合系よって誘起されるミクロドメイン形成をシミュレーションした.POPC, コレステロール,多価不飽和ホスファチジルエタノールアミン (PE) の混合系のDPDシミュレーションを実施したところ,片側の尾部に6個の不飽和部位を持つ0+6PE,両方の尾部に3個ずつ不飽和部位がある3+3PE,両方の尾部に6個の不飽和部位を持つ6+6PEの順にマイクロドメインサイズが大きくなる実験事実が再現された [9].
更にナノバイオデバイスへの応用端緒として,基盤となるシリカ表面への膜の吸着の再現を行った.界面の電子状態を考慮し,バルク中の水と異なる性質を持つ水粒子(界面水)を定義することで,シリカ基盤への膜の吸着に成功した [10].その他,電位差を印加することで膜に細孔形成させるエレクトロポレーションの検討も行っている.DDS 目線では,siRNAキャリアとしてのリポソームを意図した正電荷脂質,リン脂質の混合膜についての構造検証も実施した.DPPC, DOPC, TAP, DOTAPについてパラメータを算定し,混合・濃度の条件を変えながらDPDシミュレーションを行うことで,膜表面で分子がクラスタ化し,偏球状の構造となることを定性的に再現した [11].
3.3 脂質ナノ粒子 (LNP)LNPもリン脂質と同様,核酸医薬においてRNAのキャリアとして研究が盛んにおこなわれている.COVID-19のmRNAワクチンでも使用され注目を浴びているが,送達性能を最適化するため,構造と物性の解明が求められている.国際医療福祉大学の米持悦生先生,大阪大学の福澤薫先生らのグループを中心に応用が進められており,ssPalmと呼ばれるLNPを対象に,芳香環の付与による遺伝子導入効率の向上の解明のために,分子構造によるLNPの構造変化をFMO-DPDにより解析した [12].DPDの結果,ベンゼン環の有無によってLNPのコンフォメーションが大きく異なることが確認され,エンドソーム膜との相互作用影響する可能性が示唆された.コレステロールや核酸を考慮したより詳細な検討や,COVID-19のワクチン組成を模したDPD計算も実施され,実験の解釈に役立つ結果が得られつつある.
3.4 溶解性改善製剤製剤分野で,難水溶性の薬剤を高分子担体中に非晶質状態で分散させることで過飽和状態を作り,経口吸収を改善する非晶質固体分散体 (ASD) 技術が注目されており,シミュレーションによる設計の最適化が期待されている.米持・福澤先生らにより応用が進められ,pH依存で水溶性のEUD-E,脂溶性のEUD-NE担体を用いて,カルバマゼピン (CBZ) 原薬の分散性を評価した [13].FMO法により担体ポリマーと原薬間の相互作用エネルギーを計算し,原薬が担体内でどのように分散するかをDPDで解析したところ,EUD-E,CBZの混合系に対し,少量のEUD-NEを添加すると分散が促進され,NE含量が一定量を超えると相分離が生じるという実験結果と良好に一致することが確認された. 現在は界面活性剤を含めたASDの吸湿シミュレーションなど,さらなる応用を進めている.
3.5 タンパク質およびペプチドFMO-DPDの展開として,元来FMOで計算されていたタンパク質やペプチド系にも応用を行っている.タンパク質は古くから全原子MDで無数のシミュレーションが実施されているが,大規模系の効率の良い計算のために,粗視化による折り畳みも多く研究されている.我々のグループでは,土居らを中心にDPDに特化し,タンパク質に重要な骨格情報を容易に導入できるCAMUSプログラムが開発されており [16],構造再現の端緒として,10残基からなる最小タンパク質であるChignolinの折り畳みを実行した.粒子同士の平衡距離や角度,二面角に相当するポテンシャルを適切に付与することで,骨格がヘアピンカーブを描く構造の再現に成功し,水素結合やCH/πといった,特徴的な残基の相互作用からなる構造も再現した [14]. 最近ではぺプチド様分子であるぺプトイドのナノシート構造に関する応用計算を実施した.フェニル基や-NH3+, -COO-の側鎖を持つペプトイドについて,実験的には1分子の鎖長が16残基以上で安定したナノシートを作ることが確認されており,DPDシミュレーションにおいてもシート状初期構造から,種々の鎖長において構造維持の検証を行った.配向秩序パラメータや動径分布関数から構造の維持度合いを算出したところ,実験と同様の結果が再現された [15].
3.6 エマルションの安定性解析慶應義塾大学の荒井規允先生らのグループで,化粧品等で扱われるエマルションの安定性の解明にFMO-DPDが用いられている [17].水中油型 (O/W) エマルションはポリビニルアルコール (PVA) のような界面活性効果を持つポリマーが油-水界面に吸着することで安定性が強化されるが,PVAのけん化度がポリメタクリル酸メチル(PMMA)表面への吸着に与える影響を解析した.PMMA/水の2層分離系の界面にPVAを配置し,異なるけん化度条件でPVAの吸着を比較したところ,部分的に加水分解されたPVAは完全に加水分解されたPVAよりもPMMAに強く吸着し,エマルションの安定性を向上させることが確認された.更にPVAのアセチル基が,油-水界面での粒子の吸着を促進する役割を果たすことが示唆された.
2018年より公開したFCEWS,FMO-DPDは上記のように一定の成果を上げているが,精度向上や高度化,速度の問題を解決するために,更なる改良に取り組んでいる.
4.1 リバースマッピングDPDシミュレーションは系の大域的な配置情報を得るのには有効だが,そのままでは原子レベルの詳細解析を行うことは難しい.そこで,DPDの結果に対し原子を割り当てるリバースマッピングの手法を用い,FMOに繋げることでナノレベルでの詳細解析を行うDPD-based Structure Reverse Mapping System (DSRMS) システムを構築した [18].一般的に粗視化サイズが大きくなると,粗視化粒子内の構造の自由度が増すために全原子化は難しくなるが,官能基の性質を重視したDPDは比較的小さなレベルの粗視化なため,単純なスキームでのリバースマッピングが可能である.Figure 2にFCEWSとDSRMSを用いた解析スキームの概要を示したが,FMO法は10万原子程度のサイズを全電子計算することができ,この連携は,局所的な詳細解析と共に得られたパラメータの精度の確認にもなりうる.応用例で示したリン脂質,高分子電解質膜について小規模のDPDモデルを作成し,リバースマッピングののちのFMO計算を実行したところ,スルホン酸-水,リン酸-水間の界面に強い水素結合が確認された.
Overview of the FCEWS and DSRMS schemes. Adapted from Ref. [18] under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.
FCEWSのエネルギー計算は高精度なものの,前述のように剛体の孤立2分子対モデルから相互作用を得るため,接触時の構造緩和や会合状態の相互作用の考慮など,モデルとして不十分な点が複数残っている.そこで複数分子が直接配位した状態,かつコンフォメーションを最適化した上で相互作用エネルギーを評価するために,MDにより多体配置を生成し,FMO法で分子凝集体の相互作用を直接算定するシステムを整備した.Benzene-水, Cyclohexane–水の2系について検証し,得られたχの結果を用いたDPD計算で界面張力を計算したところ,実験値と良好な結果が得られている.従来の2体モデルとの相互作用を比較すると,2体では特に水分子の相互作用を過大評価する傾向が見受けられた.
4.3 機械学習による高速化pre_fcews上記のようにFCEWSシステムは汎用的な分子種の算定が可能なものの,小分子に対して大量のFMO計算を行うため,計算コストの削減が大きな課題である.「富岳」のような超並列スパコンを同時並行実行(capacity computing) 的に用いることで,高速な算定が可能だが,研究室レベルの小規模サーバを使う場合,多成分系では1週間以上の時間を要することもある.そこで土居らを中心に,機械学習を用い,精度を担保しつつ計算の総数を削減できるpre_fcewsが開発された [19].記述子にはセグメントのRESP電荷や距離に関連するパラメータ群を用い,セグメント対のエネルギーの一部を予測値で置き換える.系に依存するものの,現状総計算コストを1/3程度に抑えることができている.
本稿では,FMO法を基盤とした粗視化相互作用パラメータの算定システムFCEWSと,それらを用いたマルチスケールシミュレーションの展開について述べた.FMO-DPDの応用では,高分子電解質膜やリン脂質,LNP,タンパク質,エマルションなど,分子間用作用が重要な様々な分野での応用が進められている.高度化も継続して推進しており,リバースマップや多体効果,機械学習を活用した高速化についても,順次公開予定である.また,FMO-DPDの発展と活用のため,「富岳」の計算力を活かした産学連携のプロジェクトも進められている.2021年の秋から機動的課題 (一般) のhp210261として1年間行い,2023年度から産業課題{hp230016, hp240013}として継続している.JSOLを幹事に,各業界から5社に参加いただき,参画各社の課題への適用とともにシステムの改良を継続している.
粗視化シミュレーションへの連携において,現在はDPD法を中心としているが,DPD法が適用できる対象が,比較的短い鎖の溶液系に近いものに限られること,メソスケールの構造予測に限られることは注意が必要な点である.FMO法はχパラメータやDPDのポテンシャルに限らず,種々の粗視化手法のポテンシャル推算に寄与できる可能性を持っており,産業界の発展に寄与できるよう,今後もFCEWSプログラムの機能開発と種々の応用計算を進めていく予定である.
FMO関係の研究開発は,立教SFRからの支援を受けています.また,一連の計算はHPCI システム利用研究課題 (課題番号;hp210261, hp230016, hp230375, hp240013) を通じて,スーパーコンピュータ「富岳」の計算資源の提供を受け実施しました.また,東北大学の「MASAMUNE」も利用し,{20S0024, 2012SC0008} 課題にて実施しました.
また本研究の一連の実施において,FCEWSの理論面では慶應義塾大学の泰岡顕治先生,DPDの技術では慶應義塾大学の荒井規允先生や「富岳」FCEWS課題にご参加の方々に多数のご助言をいただきました.また,製剤分野では国際医療福祉大学の米持悦生先生,大阪大学の福澤薫先生,千葉大学の東顕二郎先生ならびに星薬科大学の学生の皆様に多大なご助言,ご協力をいただきましたことを深謝いたします.